硬件和科学计算的演变关系：

• 几十年来的硬件进步：计算机硬件不断快速发展，从提升单核速度，到多核并行。
• 科学计算的驱动力：科学计算需求推动硬件创新，比如需要更多计算能力、更高性能。
• 当前的解决方案是并行架构：为了继续提升性能，硬件设计转向多核、多线程、分布式计算等并行体系结构。
• 机器变得越来越复杂：随着并行和异构计算的发展，硬件架构愈加复杂，软件也要适应这种复杂性。

现代硬件复杂性的一个具体例子：

• 硬件背景：过去几十年硬件快速发展，现在科学计算成为推动硬件创新的主要力量。
• 并行架构：为满足计算需求，硬件采用了并行设计。
• 机器复杂性增加：即使是“简单”的笔记本，也包含多种并行技术和硬件加速器。

以这台笔记本为例：

• CPU：Intel Core i5-2410M，主频2.3 GHz，有4个逻辑核心（通常是2物理核+超线程）。
• SIMD扩展：支持SSE2到SSE4.2，AVX（高级向量扩展），能一次处理多个数据，提升单核性能。
• GPU：NVIDIA GeForce GT 520M，有48个CUDA核心，支持GPU并行计算。

总结

这说明现代计算设备包含多层次的并行：

• 多核CPU并行（多线程、多核）
• SIMD指令并行（单指令多数据，向量计算）
• GPU并行（大量小核心，适合大规模并行任务）

计算机硬件和软件的发展阶段，以及随之而来的性能提升和表达能力（expressiveness）变化：

1. 单核时代 (Single Core Era)

• 性能重点：通过提升单核频率和指令集优化来提升性能。
• 代表语言：C、Fortran、C++，多以顺序执行为主。
• 特点：顺序执行，编程相对简单，但性能提升受限于时钟频率。

2. 多核/SIMD时代 (Multi-Core/SIMD Era)

• 性能重点：
  • 利用多核并行（线程）
  • 利用SIMD指令集（单指令多数据）提高单核效率
• 代表语言：Java（支持多线程和并行模型）
• 特点：
  • 需要写并发代码，提高表达能力（expressiveness）以处理线程和数据并行。
  • 性能依赖于有效的线程管理和SIMD使用。

3. 异构计算时代 (Heterogeneous Era)

• 性能重点：
  • 利用多种并行架构，如GPU、SIMD、Intel Phi（加速器）
  • 线程级并行和分布式计算
• 特点：
  • 复杂度更高，软件必须支持多种硬件资源
  • 分布式系统和异构硬件成为主流，性能提升巨大，但编程复杂度也显著增加。

额外说明：

• 表达能力（Expressiveness）：随着硬件复杂度增加，软件需要更丰富的抽象和语言特性来表达并行、异构计算，保证代码可维护和高效执行。
• 你提到的“Sequential”（顺序执行）是在单核时代的标志，后来被多线程和异构计算取代或补充。

科学计算设计工具时面临的挑战和解决思路：

挑战（Challenges）

1. 非破坏性（Be non-disruptive）
   工具设计要尽量不打断用户已有的工作流程或代码架构，容易集成。
2. 领域驱动优化（Domain driven optimizations）
   针对具体科学计算领域做专门优化，而不是通用但低效。
3. 直观的用户API（Provide intuitive API for the user）
   让用户使用简单，易于理解和掌握。
4. 支持多样化硬件架构（Support a wide architectural landscape）
   适应各种并行硬件和计算平台，比如CPU、GPU、分布式系统等。
5. 高效性能（Be efficient）
   性能不能妥协，工具必须产生高性能代码。

我们的方案（Our Approach）

1. 设计成C++库（Design tools as C++ libraries）
   利用C++强大的抽象和性能能力，方便集成且高效。
2. 用领域特定嵌入式语言(DSEL)设计库（Design these libraries as Domain Specific Embedded Language）
   让库本身像一种专门的语言，提供符合领域思维的表达方式，满足挑战2和3。
3. 使用并行骨架（Parallel Skeletons）作为并行组件（Use Parallel Skeletons as parallel components）
   抽象并行模式，方便复用和跨平台，实现挑战4。
4. 利用生成式编程（Generative Programming）交付性能（Use Generative Programming to deliver performance）
   编译期代码生成和模板元编程，最大化运行时效率，满足挑战5。
   简单总结：他们选择用C++写一套类似“领域语言”的库，封装并行模式，借助编译期代码生成实现高性能且跨架构的科学计算工具，同时保持易用和兼容性。

**并行骨架（Parallel Skeletons）**的理由，分为软件抽象和硬件抽象两个方面：

为什么用并行骨架（Why using Parallel Skeletons）

1. 软件抽象（Software Abstraction）

• 写代码时不用管并行细节
  程序员只需关注算法逻辑，底层的线程、同步、负载均衡等复杂细节被骨架封装。
• 代码可扩展且易维护
  代码结构清晰，易于理解和扩展，减少并行代码常见的bug。
• 易调试、可证明、可认证（Debuggable, Provable, Certifiable）
  抽象良好，行为可预测，更容易进行形式验证和质量保证。

2. 硬件抽象（Hardware Abstraction）

• 语义固定，具体实现可自由切换
  并行骨架定义了高层语义接口，底层可以针对不同硬件（CPU多核、GPU、分布式）选择最佳实现。
• 可组合性 ⇒ 层次化架构（Composability ⇒ Hierarchical architecture）
  不同骨架可以组合形成复杂并行模式，支持分层设计，方便构建大型并行系统。
  简单来说，并行骨架帮你把复杂的并行编程拆解成易管理、易重用的模块，同时屏蔽硬件细节，写出高效、可维护的并行代码。
  

这张图片展示了一个由三个主要部分组成的系统架构：领域特定应用程序描述、生成组件和具体应用程序。

1. 领域特定应用程序描述：以树形结构表示应用程序的结构，展示特定领域的需求或配置。
2. 生成组件：包括“Translator”和“Parametric Sub-components”。Translator负责将领域描述转换为参数化子组件（各种形状），而Parametric Sub-components表示系统中使用的基本构建模块（例如：矩形、圆形、三角形等）。
3. 具体应用程序：通过Translator转换后的子组件组合，生成具体的应用程序或结构（例如：复杂的几何形状组合）。

生成式编程（Generative Programming），以及它作为工具的相关技术和定义：

生成式编程作为工具（Generative Programming as a Tool）

可用技术（Available techniques）

• 专用编译器（Dedicated compilers）
  设计专门的编译器来支持生成式编程。
• 外部预处理工具（External pre-processing tools）
  在编译前对源码做自动生成和转换。
• 支持元编程的语言（Languages supporting meta-programming）
  利用语言本身支持的元编程特性来实现代码生成。

元编程定义（Definition of Meta-programming）

• 元编程是编写程序，这些程序可以分析、转换并生成其他程序或自身，以程序作为数据进行操作。

C++元编程（C++ meta-programming）

• 利用C++的**模板子语言（template sub-language）**实现。
• 处理**类型（types）和整数常量（integral constants）**的编译时计算。
• 已被证明是图灵完备的（Turing-complete），也就是说，模板元编程理论上可以完成任何可计算任务。
  总结：生成式编程通过元编程技术，实现了在编译期对代码的分析和生成，使得C++代码可以写得更加灵活高效，同时提升性能和表达力。

领域特定嵌入式语言（Domain Specific Embedded Languages, DSEL），主要内容如下：

什么是 DSEL？

• DSL（Domain Specific Language，领域特定语言）
  是一种针对特定领域设计的声明式语言，简洁且易用，专门为该领域问题建模。
• DSEL（Domain Specific Embedded Language）
  是将DSL嵌入到通用编程语言中的一种技术，比如C++里的DSEL。

C++ 中的 DSEL（EDSL）

• 依赖于运算符重载的“滥用”技术，也就是表达式模板（Expression Templates）。
• 可以在代码片段中携带语义信息，不仅仅是普通的计算，而是让代码本身“懂”它想表达的领域含义。
• 这样泛型实现（Generic implementation）能够自我感知优化机会。
• 利用静态抽象语法树（AST），在编译期对代码进行分析和生成：
  • 表达式层面（expression level）：进行代码生成
  • 函数层面（function level）：实现跨函数的优化（inter-procedural optimization）
    总结：DSEL 通过嵌入式技术，让开发者用符合领域习惯的方式写代码，同时编译器可以根据这些表达的语义进行高级优化，从而达到高效又易用的目的。

**嵌入式领域特定语言（EDSL）**在C++中的应用，重点和优势如下：

EDSL 在 C++ 中的关键点

• 依赖运算符重载的“滥用”
  例如通过 Boost.Proto 库，实现表达式模板的机制，把代码片段当成带语义的信息载体。
• 携带语义信息
  代码不仅仅是执行动作，还在类型和表达式结构中携带了领域相关的语义。
• 泛型实现具备“自我感知”的优化能力
  因为表达的语义在类型系统中，所以编译期可以根据语义做特定优化。

优势（Advantages）

• 能引入 DSL 而不破坏开发流程
  不需要额外的工具链，只用现有的C++编译器和语言特性就能实现领域特定语言。
• 语义定义为类型信息，意味着可以在编译期解决
  编译器可以在编译阶段进行静态推导和优化，提升性能和安全性。
• 还可以利用丰富的运行时绑定机制
  兼顾灵活性和性能，结合编译期和运行时的优势。
  总结：通过 C++ 的 EDSL，开发者可以无缝地在已有C++环境里构建专门领域的语言抽象，既享受DSL的表达简洁性，又保留C++的高性能和工具链兼容性。

表达式模板（Expression Templates） 在矩阵计算中的应用，重点如下：

关键点

• 代码示例：

  matrix x(h,w), a(h,w), b(h,w);
  x = cos(a) + (b * a);
  

• 编译器不会直接执行普通的逐元素操作，而是构建了一个元表达式树（meta-AST），表示这个计算过程：

  expr<
    assign,
    expr<matrix&>,
    expr<
      plus,
      expr<cos, expr<matrix&>>,
      expr<multiplies, expr<matrix&>, expr<matrix&>>
    >
  
  >(x,a,b);
  
  

• 这个树可以被编译器/库用来做任意变换，比如：
  • 合并计算，减少临时变量
  • 调度并行执行（比如用 OpenMP ）
  • 针对特定硬件做优化

运行时展开示例（伪代码）：

#pragma omp parallel for
for(int j=0; j < h; ++j) {
  for(int i=0; i < w; ++i) {
    x(j,i) = cos(a(j,i)) + (b(j,i) * a(j,i));
  }
}

总结

表达式模板通过静态构建表达式树，让编译器能够：

• 避免创建临时矩阵，减少内存开销
• 实现自动向量化、并行化等优化
• 在编译期对代码做重写，提高性能
  这是一种非常强大且高效的元编程技术，广泛应用于线性代数库如 Eigen、Blaze 等。

“并行领域专用嵌入式语言（Parallel DSEL）”在实际中的应用，核心内容总结如下：

目标 (Objectives)

• 利用DSEL生成技术，针对不同硬件进行代码生成
• 展示抽象层带来的低性能开销（低成本抽象）
• 演示并行骨架（skeletons）的可用性和效果

贡献 (Our contribution)

• BSP++
  通用C++ BSP（Bulk Synchronous Parallel）框架，支持共享内存和分布式内存架构
• Quaff
  针对并行骨架编程的领域专用语言
• Boost.SIMD
  便携式SIMD编程的DSEL，方便使用SIMD指令集实现向量化
• NT2
  类MATLAB的科学计算DSEL，简洁高效，支持向量化和并行化

总结

这些项目结合了DSEL和并行编程技术，解决了如何用高层抽象同时保证效率的问题。它们通过模板元编程和代码生成技术，支持多种硬件，方便科学计算领域的开发者写出高性能且可维护的代码。

NT2，这是一个科学计算库，核心特点总结如下：

NT2 概述

• 科学计算库，目标是简化数值计算，提供类似 MATLAB 的简洁接口
• 高性能：提供高效的计算基本单元和操作原语，保证计算性能
• 易扩展：设计灵活，方便用户和开发者扩展功能

组成部分

• Boost.SIMD
  利用Boost.SIMD实现底层向量化优化，充分发挥CPU SIMD指令集优势
• 递归并行骨架
  利用并行骨架（parallel skeletons），递归分解任务，实现高效的并行执行
• 架构与运行时无关
  代码设计上与硬件架构和具体运行时解耦，提高移植性和灵活性

总结

NT2通过现代C++模板元编程和DSEL技术，将用户友好的接口和底层高性能结合起来，非常适合科研人员和工程师进行大规模科学计算。

The Numerical Template Toolbox (NT2) 的核心设计理念和使用流程：

核心原则

• table<T,S> 类型
  一个简单的多维数组类模板，行为和 MATLAB 的数组一模一样
• 丰富的函数库
  超过500个函数，既可以作用于 table（数组），也可以作用于标量，接口和 MATLAB 非常接近

使用流程

1. 拿一个 MATLAB 的 .m 文件
2. 复制改写成 .cpp 文件
3. 在代码开头添加：

   #include <nt2/nt2.hpp>
   

   并做必要的格式（语法）调整
4. 编译这个 .cpp 文件并链接到 libnt2.a 静态库
5. 完成，程序即可高效运行，且几乎和 MATLAB 代码一样简单

总结

NT2 就像是 MATLAB 的 C++ 版本，让用户可以用接近 MATLAB 的风格写高性能 C++ 数值计算程序，兼顾易用性和速度。

这段代码展示了如何用 MATLAB 代码表达数值计算，具体是：

MATLAB代码含义：

A1 = 1:1000;                           % 生成向量 A1，包含1到1000的整数
A2 = A1 + randn(size(A1));            % A2 = A1 + 正态分布随机噪声，长度与A1相同
X = lu(A1*A1');                       % 计算矩阵 A1*A1' 的 LU 分解
rms = sqrt(sum(sqr(A1(:)-A2(:))) / numel(A1));  % 计算 A1 和 A2 差的均方根误差

这里每步操作：

• 1:1000 生成 1 到 1000 的行向量
• randn(size(A1)) 生成与 A1 同尺寸的正态分布随机值
• A1*A1' 是向量与其转置的乘积，生成一个 1000×1000 矩阵
• lu(...) 是矩阵的 LU 分解
• A1(:) 和 A2(:) 都是把矩阵转换成列向量
• sqr(...) 是逐元素平方
• sum(...) 是求和
• numel(A1) 是元素个数，1000
• sqrt(...) 是平方根

NT2 里写法（伪代码风格）

#include <nt2/nt2.hpp>
nt2::table<double> A1 = nt2::_(1,1000);  // 创建1到1000的数组
nt2::table<double> A2 = A1 + nt2::randn(nt2::size(A1)); // 加噪声
auto X = nt2::lu(A1 * nt2::trans(A1));   // LU分解
auto diff = A1 - A2;
double rms = nt2::sqrt(nt2::sum(nt2::sqr(nt2::reshape(diff, nt2::numel(A1), 1))) / nt2::numel(A1));

C++ 代码用 NT2 库，实现了与之前 MATLAB 代码几乎相同的功能。详细说明如下：

table<double> A1 = _(1., 1000.);                   // 创建一个从1到1000的数组（向量）
table<double> A2 = A1 + randn(size(A1));          // 给A1加上正态随机噪声，生成A2
table<double> X = lu(mtimes(A1, trans(A1)));      // 计算A1乘以转置A1的矩阵乘积，并对结果做LU分解
double rms = sqrt(sum(sqr(A1(_) - A2(_))) / numel(A1)); // 计算A1与A2差值的均方根误差

各部分解释：

• _(1., 1000.)：生成从1到1000的数组，类似 MATLAB 的 1:1000。
• randn(size(A1))：生成和 A1 同样大小的正态分布随机数组。
• mtimes(A1, trans(A1))：矩阵乘法，trans(A1) 是转置。
• lu(...)：LU分解函数。
• A1(_)、A2(_)：NT2 里“全部元素”的语法，类似 A1(:)。
• sqr：逐元素平方。
• sum：求和。
• numel(A1)：元素个数。

整体来说：

NT2通过重载操作符和提供类似 MATLAB 的API，让C++也能像MATLAB一样写矩阵和向量运算，且支持高效编译。

Motion Detection

table<char> sigma_delta(table<char>& background, table<char> const& frame, table<char>& variance) {
    // Estimate Raw Movement
    background = selinc(background < frame, seldec(background > frame, background));
    table<char> diff = dist(background, frame);
    // Compute Local Variance
    table<char> sig3 = muls(diff, 3);
    var = if_else(diff != 0, selinc(variance < sig3, seldec(var > sig3, variance)), variance);
    // Generate Movement Label
    return if_zero_else_one(diff < variance);
}

这段代码描述的是一个 Sigma-Delta运动检测算法 的 NT2 实现，它是一种用于视频分析的轻量级背景建模方法，适合在低功耗、嵌入式设备中使用（例如机器人视觉、安防监控等）。下面是代码分段讲解：

背景

• 该算法用于检测运动物体。
• 不使用浮点数，仅使用小范围整数（如 char），以便提升性能并减少资源使用。
• 目标是用 SIMD 向量化库（NT2）进行高效实现。

代码功能分解：

background = selinc( background < frame,
                     seldec(background > frame, background) );

1. 背景模型更新

• selinc / seldec 是 条件自增/自减。
• background 向当前帧 frame 靠拢：
  • 如果当前帧更亮，则背景变亮；
  • 如果当前帧更暗，则背景变暗；
• 这是 Sigma-Delta 的核心建模步骤。

table<char> diff = dist(background, frame);

2. 差值计算

• dist 计算每个像素点当前帧与背景之间的差值（绝对值）。

table<char> sig3 = muls(diff,3);

3. 扩大差值

• 将差值乘3，用于调整容差窗口（threshold）；
• muls 是饱和乘法，防止溢出。

var = if_else( diff != 0,
               selinc(variance < sig3,
                      seldec(var > sig3, variance)),
               variance);

4. 更新局部方差（variance）

• 如果有差异（diff != 0），则：
  • 如果当前方差过小，则递增；
  • 如果过大，则递减；
• 否则，保持不变。

return if_zero_else_one( diff < variance );

5. 输出运动掩码

• 比较差值是否小于方差（在容差范围内 → 认为是背景）；
• if_zero_else_one 产生 0/1 掩码（1 表示运动，0 表示无变化）。

总结

步骤	含义
背景建模	background 向 frame 靠拢
差值计算	估计当前帧与背景差异
方差调整	根据信号变化动态调整容差
运动检测	若变化超出容差，则为“运动”
该算法结构简单但实用，且易于并行和 SIMD 化，适用于资源受限场景。

NT2（Numerical Template Toolbox）库编写的 Black-Scholes 期权定价公式 的两种版本：

1. 普通版（直接表达式计算）
2. Loop Fusion 优化版（缓存友好、更高性能）

背景：Black-Scholes 模型公式

欧式看涨期权定价公式如下：
$$C=S\cdot N(d_1)-X\cdot e^{-rT}\cdot N(d_2)$$
其中：
$$d_1=\frac{\ln (S/X)+(r+\frac{1}{2}{v}^2)T}{v\sqrt{T}},d_2=d_1-v\sqrt{T}$$

• $S$: 当前资产价格 (Sa)
• $X$: 执行价格 (Xa)
• $T$: 到期时间 (Ta)
• $r$: 无风险利率 (ra)
• $v$: 波动率 (va)
• $N(\cdot )$: 标准正态累积分布函数 (normcdf)

普通版代码：直接链式表达式

table<float> blackscholes(table<float> const& Sa,
                          table<float> const& Xa,
                          table<float> const& Ta,
                          table<float> const& ra,
                          table<float> const& va) {
    table<float> da = sqrt(Ta);
    table<float> d1 = log(Sa / Xa) + (sqr(va) * 0.5f + ra) * Ta / (va * da);
    table<float> d2 = d1 - va * da;
    return Sa * normcdf(d1) - Xa * exp(-ra * Ta) * normcdf(d2);
}

优点：简洁、可读性强
缺点：每一步可能创建中间临时变量，性能可能受影响，特别是 多次遍历内存，对 CPU cache 不友好。

Loop Fusion 优化版

table<float> blackscholes(table<float> const& Sa, table<float> const& Xa,
                          table<float> const& Ta, table<float> const& ra,
                          table<float> const& va) {
    // 预分配表（避免反复申请内存）
    table<float> da(extent(Ta)), d1(extent(Ta)), d2(extent(Ta)), R(extent(Ta));
    // tie + expression templates 实现 loop fusion
    tie(da, d1, d2, R) = tie(
        sqrt(Ta),
        log(Sa / Xa) + (sqr(va) * 0.5f + ra) * Ta / (va * sqrt(Ta)),  // 注意此处da是并行绑定的
        d1 - va * da,
        Sa * normcdf(d1) - Xa * exp(-ra * Ta) * normcdf(d2)
    );
    return R;
}

优化点详解

• tie(...) = tie(...)：
  • 这是 NT2 的一个 loop fusion 技巧，可以让多个数组的元素计算在一个循环内完成。
• 提高缓存局部性（cache locality）和流水线效率
• 减少临时变量和内存分配开销
• 显式管理中间变量 (da, d1, d2) 可以在后续重用或进一步优化

总结对比

版本	优点	缺点
普通表达式	简洁、易读	中间变量多、可能缓存不友好
Loop Fusion	高性能、单遍多算	可读性差一点、维护难度略高
如果你想深入实践这个例子，我还可以给你：

• OpenMP 或 CUDA 后端扩展方式（基于 NT2）
• 对比测试结果模板（性能评估）
• 更详细解释 tie / extent / normcdf 的底层实现