Jack J. Dongarra：他用标准数值库和性能基准，为高性能计算建立了统一语言，让超级计算机从竞速工具变成科学基础设施

图灵奖第五十六届（2021）| Jack J. Dongarra：他用标准数值库和性能基准，为高性能计算建立了统一语言，让超级计算机从竞速工具变成科学基础设施

一句话概括：当科学家们用各自的方式编写数学计算程序、超级计算机以性能数字互相攀比时，Dongarra创建了LINPACK、BLAS、LAPACK等标准数值库，设计了TOP500基准，让全世界的高性能计算有了共同的算法基础和评价标准，从天气预报到基因测序，从核武器模拟到深度学习，让”算得快且算得准”成为科学发现的引擎。

🏆 获奖简介

Jack Joseph Dongarra（杰克·约瑟夫·唐加拉，1950-）是高性能数值计算的架构师，线性代数软件库的奠基人，超级计算机性能基准的设计者。

个人信息：

出生时间：1950年7月18日
出生地点：美国伊利诺伊州芝加哥
主要成就：
- LINPACK（线性方程组求解库）开发者
- BLAS（基本线性代数子程序）标准化推动者
- LAPACK（线性代数包）主要设计者
- ScaLAPACK（可扩展LAPACK）创建者
- MAGMA（GPU加速线性代数）先驱
- TOP500超级计算机排行榜联合创始人
- Netlib数值软件库创始人

教育背景：

1972年：芝加哥州立大学数学学士
1973年：伊利诺伊理工学院计算机科学硕士
1980年：新墨西哥大学应用数学博士

职业经历：

1980-1989年：阿贡国家实验室研究员
1989-至今：田纳西大学杰出教授
1989-至今：橡树岭国家实验室杰出科学家
获奖年份：2021年
获奖原因：对数值算法和库的开创性贡献，使这些方法成为高性能计算的基石

为什么他是第五十六位？ 1970年代，科学计算领域混乱不堪：每个实验室用自己的算法求解线性方程组，每台超级计算机用不同方式报告性能，程序从Cray移植到IBM需要重写。Dongarra改变了这一切——LINPACK提供了标准求解器和性能测试，BLAS定义了矩阵运算的通用接口，LAPACK继承LINPACK但针对现代内存层次优化，ScaLAPACK扩展到分布式集群，MAGMA适应GPU时代。更重要的是，TOP500排行榜用LINPACK基准统一了性能评价标准，Netlib网站让这些库免费传播全球。从量子化学到气候模拟，从石油勘探到机器学习，今天90%以上的科学计算都依赖他创建的数值库。他让”高性能计算”从各自为战的艺术变成标准化的工程，让超级计算机从实验室的炫耀品变成推动人类认知边界的基础设施。

🚀 主要贡献

1. LINPACK：线性代数的第一个标准库

背景（1970年代末）：

问题：求解线性方程组 Ax = b 是科学计算的核心
- 结构力学：有限元网格生成数百万方程
- 流体力学：Navier-Stokes方程离散化
- 量子力学：薛定谔方程特征值问题
现状混乱：
- 每个实验室自己实现高斯消元
- 代码质量参差不齐（数值稳定性问题）
- 无法在不同机器间移植

LINPACK的诞生（1979）：

设计目标：

正确性：数值稳定的算法（列主元LU分解）
高效性：充分利用向量处理器（Cray-1时代）
可移植性：Fortran代码，最小机器依赖
文档化：每个函数有详细说明

核心算法：

LU分解（求解 Ax = b）：

1. 分解：A = PLU
   - P：行交换矩阵（列主元策略）
   - L：下三角矩阵
   - U：上三角矩阵

2. 前向替换：解 Ly = Pb
3. 后向替换：解 Ux = y

为什么需要列主元：

避免除以接近零的数（数值不稳定）

例子：

0.0001x + y = 1
x + y = 2

不选主元：第一步消元会放大误差
选主元：交换两行，避免小除数

函数命名规范（沿用至今）：

SGESV：单精度通用矩阵求解器
  S - Single precision（单精度）
  GE - GEneral matrix（通用矩阵）
  SV - SolVe（求解）

DGETRF：双精度通用矩阵LU分解
  D - Double precision
  GE - GEneral
  TRF - TRiangular Factorization

软件结构：

驱动程序（Driver routines）：高层接口，如SGESV
计算程序（Computational routines）：核心算法，如SGETRF（LU分解）
辅助程序（Auxiliary routines）：工具函数

LINPACK基准测试：

意外副产品：测试LINPACK性能时，Dongarra记录了各种机器的速度

1979年报告：列出了10多台计算机的性能
意义：
- 首次用统一方法比较不同架构
- 问题规模：1000×1000矩阵（当时很大）
- 测量指标：MFLOPS（百万次浮点运算/秒）

演变：这一基准后来成为TOP500排行榜的基础（见下文）。

影响：

标准化：LINPACK成为求解线性系统的事实标准
教育：数值分析课程的实验库
商业软件：MATLAB早期内核就是LINPACK

2. BLAS：数值计算的通用语言

问题（1970年代）：

LINPACK虽然可移植，但在不同机器上性能差异巨大
原因：
- Cray-1：向量处理器，适合向量操作
- IBM 3090：标量处理器
- 编译器无法自动优化到硬件特性

解决方案：定义标准接口，让硬件厂商优化实现

BLAS层次结构：

Level 1 BLAS（1973年，Charles Lawson等提出，Dongarra推广）：

向量操作：O(n)计算量

例子：

SAXPY：y = a*x + y（向量加法）
SDOT：计算向量点积 x·y
SNRM2：计算向量2-范数 ||x||

Level 2 BLAS（1988年，Dongarra等设计）：

矩阵-向量操作：O(n²)计算量

例子：

SGEMV：y = α*A*x + β*y（矩阵向量乘）
STRSV：解三角方程 A*x = b

动机：数据复用，减少内存访问

Level 3 BLAS（1990年，Dongarra等设计）：

矩阵-矩阵操作：O(n³)计算量

例子：

SGEMM：C = α*A*B + β*C（矩阵乘法）
STRSM：解多右端三角方程 A*X = B

关键：最高数据复用率，最适合优化

为什么BLAS重要？

性能差异（示例：1000×1000矩阵乘法）：

朴素三重循环：1 GFLOPS
优化BLAS（利用缓存、向量化）：100 GFLOPS
差距：100倍！

可移植性：

科学家调用SGEMM
Intel提供MKL（Math Kernel Library）优化实现
AMD提供BLIS
NVIDIA提供cuBLAS（GPU版本）
程序不改一行代码，在不同平台都能高效运行

实例（BLAS-3的威力）：矩阵乘法 C = A * B（n×n矩阵）：

朴素算法：
  for i in 1..n
    for j in 1..n
      for k in 1..n
        C[i,j] += A[i,k] * B[k,j]

问题：
  - 内循环每次从内存读取B[k,j]（缓存不友好）
  - 无向量化

BLAS-3优化：
  - 分块算法（Block Algorithm）：
    - 将矩阵分成小块（如64×64）
    - 每块放入L1缓存
    - 对块内做密集计算
  - 数据复用：每个矩阵元素被使用O(n)次，而内存访问只有O(1)
  - 向量化：利用SIMD指令

结果：现代科学计算90%的时间花在BLAS-3调用上。

3. LAPACK：继承LINPACK的现代化升级

背景（1980年代末）：

硬件变化：
- 1970年代：Cray向量机，长向量操作高效
- 1980年代：RISC工作站，缓存层次结构关键
LINPACK局限：
- 基于Level 1 BLAS（向量操作）
- 在缓存机器上性能差

LAPACK项目（1987-1992，Dongarra领导）：

设计原则：

基于Level 3 BLAS：最大化缓存利用
完全替代LINPACK：兼容接口
扩展功能：更多矩阵类型、特征值问题

核心功能：

线性系统求解：

通用矩阵：GESV（LU分解）
对称正定：POSV（Cholesky分解）
对称不定：SYSV（Bunch-Kaufman分解）
带状矩阵、三对角矩阵…

特征值与奇异值：

GEEV：通用矩阵特征值
- 应用：振动模式分析、量子力学
SYEV：对称矩阵特征值（Jacobi、分治法、MRRR）
GESVD：奇异值分解（SVD）
- 应用：主成分分析、推荐系统

最小二乘：

GELS：求解超定/欠定系统（QR或LQ分解）
- 应用：数据拟合、参数估计

算法创新：

分块算法示例（Cholesky分解）：

传统LINPACK（列算法）：
  for j = 1 to n
    A[j,j] = sqrt(A[j,j])  // 单个元素
    for i = j+1 to n
      A[i,j] /= A[j,j]     // 向量操作（BLAS-1）
    更新剩余矩阵            // BLAS-2

LAPACK分块算法：
  for j = 1 to n step NB（块大小）
    分解对角块（NB×NB）     // 递归调用
    更新当前列块            // BLAS-2
    更新剩余矩阵            // BLAS-3（GEMM）

性能提升：

LINPACK：约 n³/3 FLOPS，性能 < 10% 峰值
LAPACK：同样计算量，性能 > 80% 峰值（通过BLAS-3）

软件工程：

7000+行Fortran（后来增加C接口）
模块化设计：驱动、计算、辅助分离
严格测试：数值稳定性、边界条件

影响：

MATLAB：从LINPACK切换到LAPACK内核
NumPy/SciPy：调用LAPACK
R统计软件：线性代数基于LAPACK
Julia：性能依赖优化BLAS/LAPACK

4. ScaLAPACK：从单机到集群的跨越

挑战（1990年代）：

问题规模爆炸：
- 气候模拟：矩阵 10⁶ × 10⁶（无法放入单机内存）
- 分子动力学：原子数百万
并行机出现：
- 分布式内存：每个节点独立内存
- 消息传递：MPI（Message Passing Interface）

ScaLAPACK（Scalable LAPACK，1995）：

设计挑战：

数据分布：如何将矩阵切分到多个节点？
负载平衡：避免某些节点空闲
通信开销：减少节点间数据传输

2D块循环分布（2D Block-Cyclic Distribution）：

例：8×8矩阵，4个处理器（2×2网格），块大小2×2

原矩阵：
  A = [a11 a12 | a13 a14 | a15 a16 | a17 a18]
      [a21 a22 | a23 a24 | a25 a26 | a27 a28]
      [--------+--------+--------+--------]
      [a31 a32 | a33 a34 | a35 a36 | a37 a38]
      [a41 a42 | a43 a44 | a45 a46 | a47 a48]
      [--------+--------+--------+--------]
      [a51 a52 | a53 a54 | a55 a56 | a57 a58]
      [a61 a62 | a63 a64 | a65 a66 | a67 a68]
      [--------+--------+--------+--------]
      [a71 a72 | a73 a74 | a75 a76 | a77 a78]
      [a81 a82 | a83 a84 | a85 a86 | a87 a88]

分配：
  P0: [a11,a12,a21,a22] [a15,a16,a25,a26] [a51,a52,a61,a62]...
  P1: [a13,a14,a23,a24] [a17,a18,a27,a28]...
  P2: [a31,a32,a41,a42]...
  P3: [a33,a34,a43,a44]...

优势：

负载平衡：每个处理器获得均匀数据
通信局部性：相邻块通信频繁，分配到邻近处理器

并行算法示例（LU分解）：

for k = 1 to n（迭代消元列）
  1. 主列广播：拥有第k列的处理器发送给同一行的其他处理器
  2. 主行广播：拥有第k行的处理器发送给同一列的其他处理器
  3. 局部更新：每个处理器更新自己的子矩阵（GEMM）

性能：

可扩展性：从4个处理器到数千个
效率：80-90%（理想情况下100%）

应用：

核武器模拟：洛斯阿拉莫斯国家实验室
气候预报：全球模型需要百万级矩阵
天体物理：暗物质模拟

5. TOP500与性能基准

动机（1993）：

问题：
- 超级计算机厂商各自声称”最快”
- 无统一标准比较
- 峰值速度 vs 实际性能差距大
需求：客观、可重复的基准测试

TOP500创立（1993年6月，Dongarra与Hans Meuer、Erich Strohmaier、Horst Simon）：

LINPACK基准详解：

HPL（High Performance LINPACK）：

问题：求解 Ax = b，A是n×n稠密矩阵
算法：LU分解 + 三角求解
自由度：
- 矩阵大小n（通常选择尽可能大，接近内存上限）
- 数据分布、块大小、算法参数
测量：
- 执行时间 T
- 计算量：(2/3)n³ + O(n²) FLOPS
- 性能 R_max = FLOPS / T

榜单规则：

每年6月、11月发布
按R_max排序
报告：
- R_max（实际性能）
- R_peak（理论峰值）
- 效率 = R_max / R_peak
- 功耗（近年新增Green500）

历史趋势：

1993年（首届）：

第1名：Thinking Machines CM-5/1024（59.7 GFLOPS）
第500名：0.4 GFLOPS

2008年（Petaflop时代）：

IBM Roadrunner：1.026 PFLOPS（10¹⁵）

2022年（Exaflop时代）：

Frontier（橡树岭）：1.102 EFLOPS（10¹⁸）
采用AMD EPYC + Radeon GPU

增长率：

约每13.5个月翻倍（Moore定律的体现）

批评与辩护：

批评：

“HPL不代表真实应用”（稠密矩阵 vs 稀疏矩阵）
“过度优化基准”（厂商针对HPL调优）

Dongarra的回应：

“HPL是最小共识”：所有系统都能跑
“其他基准欢迎补充”：HPCG（稀疏矩阵，2013）、AI基准（MLPerf）
“透明性”：算法公开，可验证

实际影响：

采购决策：政府/企业参考TOP500
技术驱动：厂商为排名投入研发
科研指标：国家科技实力象征

6. 后续工作：GPU时代的MAGMA

新挑战（2000年代末）：

GPU崛起：NVIDIA CUDA让GPU做通用计算
峰值性能：单GPU可达数TFLOPS，远超CPU
问题：
- LAPACK只支持CPU
- CUDA编程复杂

MAGMA项目（Matrix Algebra on GPU and Multicore Architectures，2008）：

设计思路：

异构计算：CPU + GPU协同
- CPU：控制流、小矩阵
- GPU：大规模并行计算（BLAS-3）
自动调度：运行时决定任务分配

示例（LU分解）：

传统CPU：
  for k = 1 to n
    列消元（CPU）
    矩阵更新（CPU）

MAGMA混合：
  for k = 1 to n
    列消元（CPU，依赖性强）
    矩阵更新（GPU，可并行，GEMM）

性能：

矩阵乘法：GPU版本比CPU快10-100倍
整体加速：5-20倍（取决于问题）

影响：

深度学习：TensorFlow、PyTorch底层调用cuBLAS/cuDNN（受MAGMA启发）
科学计算：地震成像、气候模拟用GPU加速
通用性：MAGMA支持AMD ROCm、Intel oneAPI

7. Netlib：开源数值软件的先驱

创建（1984年，Dongarra与Eric Grosse）：

目的：免费分发高质量数值软件
方式：
- FTP服务器（当时互联网早期）
- 电子邮件请求（发邮件自动回复代码）

内容：

LINPACK、LAPACK、BLAS
特殊函数库（Bessel函数、误差函数）
测试矩阵（Matrix Market）

意义：

民主化：任何人可免费获取
标准化：避免重复发明轮子
教育：学生可学习实际工业级代码

现代遗产：

今天GitHub、CRAN、PyPI继承精神
Netlib仍在运行（netlib.org）

🌍 对世界的深远影响

科学发现的加速器

案例1：天气预报

问题：求解Navier-Stokes方程（流体力学），矩阵 10⁸ × 10⁸
使用：ScaLAPACK在超级计算机上并行求解
结果：5天预报准确率从60%提升到85%（1990 vs 2020）

案例2：药物设计

问题：分子对接（Docking），计算蛋白质-药物结合能
使用：LAPACK求解特征值（分子轨道）
影响：疫苗开发时间从10年缩短到1年（COVID-19疫苗）

案例3：引力波探测

LIGO实验：需要实时处理海量数据
使用：BLAS/LAPACK做信号处理（FFT、滤波）
成果：2016年首次探测到引力波，验证广义相对论

产业应用的基础设施

金融工程：

期权定价：Black-Scholes方程数值解（有限差分法 → 线性系统）
风险管理：Monte Carlo模拟（依赖随机数矩阵运算）

石油勘探：

地震成像：反演问题，求解大规模线性系统
使用：商业软件（如Schlumberger）内核是LAPACK

机器学习：

训练神经网络：反向传播 = 大量矩阵乘法（GEMM）
TensorFlow/PyTorch：底层调用cuBLAS（GPU版BLAS）
事实：80%训练时间在GEMM上

教育与人才培养

大学课程：

数值分析：教材配套LAPACK实验
高性能计算：学生用MPI + ScaLAPACK完成项目

培训项目：

Dongarra在田纳西大学开设HPC暑期课程
培养数千名科学计算人才

开源文化：

代码公开→学生可阅读”真实工业代码”
影响一代程序员的软件工程观

🏆 获奖理由（通俗版）

ACM官方表彰：“对数值算法和库的开创性贡献，这些算法和库实现了高性能计算软件的长期发展。”

更通俗的理解：

Jack Dongarra证明：

标准化的力量：BLAS统一接口让软件在任何硬件上都能高效运行
开源的价值：Netlib让全世界共享高质量代码，避免浪费
基准的必要：TOP500用客观数据衡量进步，推动竞争和创新

他让”高性能计算”从各自为战的混乱变成协同发展的生态系统，让超级计算机从实验室玩具变成推动科学、工程、商业的引擎。没有他的工作，今天的天气预报不会准确，药物研发不会快速，深度学习不会爆发。

👤 个人生平与传奇

成长经历

童年（1950-1960年代）：

芝加哥工人家庭：父母非学术背景
数学天赋：高中时期表现出色
大学选择：芝加哥州立大学（非顶尖，但负担得起）

早期职业（1970年代）：

1972年毕业：数学学士
第一份工作：阿贡国家实验室程序员
- 机缘巧合：被分配到数值软件组
- 导师：Jim Pool，EISPACK项目（特征值软件包）

学术之路

研究生（1970年代末）：

边工作边读博：新墨西哥大学应用数学
论文（1980）：矩阵计算的并行算法
动力：不为学位，纯粹对问题感兴趣

阿贡实验室（1980-1989）：

LINPACK项目：首次主导大型软件开发
学到教训：
- 软件工程的重要性（文档、测试）
- 社区合作（多机构协作）

田纳西大学（1989-至今）：

职位：杰出教授（Distinguished Professor）
双重身份：
- 大学教授：教学、指导研究生
- 国家实验室科学家（橡树岭）：大型项目
优势：学术自由 + 实际问题接触

性格与风格

实用主义：

“理论必须能解决真实问题”
“性能数字说话，不吹嘘”

合作精神：

LAPACK有40+合作者
主张”站在巨人肩上”，继承EISPACK/LINPACK

持久力：

LINPACK(1979) → LAPACK(1992) → ScaLAPACK(1995) → MAGMA(2008)
40年持续改进同一领域

幽默感：

回忆TOP500早期：“厂商请我们喝酒，希望排名高一点”
自嘲：“我只是写了些Fortran代码，没想到会得图灵奖”

荣誉与认可

主要奖项：

图灵奖（2021）
IEEE计算机先驱奖（2020）
SIAM/ACM计算科学与工程奖（2019）
美国国家工程院院士（2006）

其他：

Google Scholar引用：20万+（高被引科学家）
h-index：100+（极高影响力）

趣闻：

Dongarra的Erdős数是3（数学家协作网络的中心度指标）
他每年要飞行10万英里参加会议（疫情前）

💭 为什么他值得纪念？

1. 他建立了科学计算的”通用语言”

之前：

每个实验室用自己的矩阵求解器
代码无法共享、无法比较

之后：

BLAS定义接口 → 全世界程序员用同一套API
LAPACK成为事实标准 → 科学论文可复现

类比：就像USB接口统一了硬件连接，BLAS/LAPACK统一了数值计算。

2. 他证明了”标准化”不等于”停滞”

批评： “标准会扼杀创新”

Dongarra的实践：

LINPACK → LAPACK：重写以适应新硬件
LAPACK → ScaLAPACK：扩展到并行
ScaLAPACK → MAGMA：拥抱GPU

启示：好的标准是”演进的框架”，而非”僵化的教条”。

3. 他让”开源”成为科学计算的规范

Netlib精神（1984）：

代码免费
文档齐全
欢迎贡献

影响：

今天GitHub、arXiv延续这一传统
开放科学（Open Science）运动的先驱

4. 他的工作”无处不在但不可见”

悖论：

99%的程序员用过他的代码（通过NumPy、MATLAB、R）
但只有1%知道LAPACK是什么

价值：最好的基础设施是”用户无需关心的基础设施”。

🔍 技术深度：为什么BLAS-3这么快？

缓存局部性分析

现代CPU架构：

寄存器：~100 bytes，1周期
L1缓存：~32KB，4周期
L2缓存：~256KB，12周期
L3缓存：~8MB，40周期
主内存：~8GB，200周期

矩阵乘法 C = A * B（n×n）：

朴素算法内存访问：

for (i = 0; i < n; i++)
  for (j = 0; j < n; j++)
    for (k = 0; k < n; k++)
      C[i][j] += A[i][k] * B[k][j];

内循环每次读 B[k][j]（不同j → 缓存失效）
总内存访问：O(n³)（与计算量相同）

BLAS-3分块算法：

// 假设 n = 1000，块大小 NB = 100
for (i = 0; i < n; i += NB)
  for (j = 0; j < n; j += NB)
    for (k = 0; k < n; k += NB)
      // 对 NB×NB 子矩阵做密集计算
      for (ii = i; ii < i+NB; ii++)
        for (jj = j; jj < j+NB; jj++)
          for (kk = k; kk < k+NB; kk++)
            C[ii][jj] += A[ii][kk] * B[kk][jj];

关键优化：

子矩阵放入缓存：
- 3个 NB×NB 矩阵 = 3 × 100² × 8字节 = 240KB（放入L2）
- 在缓存中完成 NB³ 次乘法
- 数据复用：每个元素被访问 NB 次，而非1次
向量化：
- 内循环用SIMD指令（如AVX2：一次8个双精度）
指令级并行：
- 现代CPU可同时执行多条无关指令
- 展开内循环增加并行度

性能差距：

朴素：5 GFLOPS（5%峰值）
BLAS-3：90 GFLOPS（90%峰值）

ScaLAPACK的通信复杂度

LU分解通信量分析（n×n矩阵，P个处理器，P = p×q网格）：

算法：

for k = 1 to n（迭代消元列）
  1. 列广播：第k列发送给同行处理器（q-1个目标）
  2. 行广播：第k行发送给同列处理器（p-1个目标）
  3. 局部更新：每个处理器计算自己的子矩阵

通信量：

每步：
- 列广播：(n/p) × (q-1) 个元素
- 行广播：(n/q) × (p-1) 个元素
- 总计：O(n/√P) 每步
总共 n 步：O(n²/√P)

计算量：

每个处理器：O(n³/P)

通信/计算比：

O((n²/√P) / (n³/P)) = O(√P / n)
结论：n越大，通信占比越小（可扩展性好）

实际性能：

P = 64，n = 10000：效率 > 90%
P = 1024，n = 100000：效率 > 85%

🧪 实践意义：如何使用BLAS/LAPACK

Python示例（NumPy/SciPy）

矩阵乘法（调用BLAS-3）：

import numpy as np

A = np.random.rand(1000, 1000)
B = np.random.rand(1000, 1000)

# 底层调用 DGEMM（双精度通用矩阵乘法）
C = A @ B

# 等价于：
# C = scipy.linalg.blas.dgemm(1.0, A, B)

线性系统求解（调用LAPACK）：

from scipy.linalg import solve

A = np.random.rand(1000, 1000)
b = np.random.rand(1000)

# 底层调用 DGESV（LU分解+求解）
x = solve(A, b)

# 验证：
print(np.allclose(A @ x, b))  # True

特征值：

from scipy.linalg import eigh

A = np.random.rand(1000, 1000)
A = (A + A.T) / 2  # 对称化

# 底层调用 DSYEVD（对称矩阵特征值，分治法）
eigenvalues, eigenvectors = eigh(A)

性能调优技巧

选择正确的BLAS实现：

# 检查当前BLAS
python -c "import numpy; numpy.show_config()"

# Linux：安装优化版本
# Intel MKL（最快，商业）
pip install mkl

# OpenBLAS（开源，接近MKL）
sudo apt install libopenblas-dev

# BLIS（AMD优化）

矩阵布局：

# C-contiguous（行主序，NumPy默认）
A = np.array([[1,2,3], [4,5,6]], order='C')

# Fortran-contiguous（列主序，LAPACK偏好）
A = np.array([[1,2,3], [4,5,6]], order='F')

# 建议：预先转换避免运行时拷贝
A = np.asfortranarray(A)

避免不必要的拷贝：

# 坏：创建副本
C = A + B
C = C @ D

# 好：就地操作
C = A + B
C @= D  # 避免分配新数组

C/C++直接调用BLAS

示例（矩阵乘法）：

#include <cblas.h>

void matrix_multiply(double *A, double *B, double *C, int n) {
    // C = A * B
    // cblas_dgemm(布局, A转置?, B转置?, M, N, K,
    //             alpha, A, LDA, B, LDB, beta, C, LDC)
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                n, n, n,        // M, N, K
                1.0,            // alpha
                A, n,           // A, LDA
                B, n,           // B, LDB
                0.0,            // beta
                C, n);          // C, LDC
}

编译：

gcc -o program program.c -lopenblas -lm

📚 延伸阅读

书籍

Dongarra等：《Templates for the Solution of Linear Systems》（1994）
- 迭代法综述（CG、GMRES等）
Golub & Van Loan：《Matrix Computations》（经典教材，第4版2013）
- 矩阵算法的数学基础
Petersen & Arbenz：《Introduction to Parallel Computing》（并行计算入门）

软件文档

LAPACK Users’ Guide（netlib.org）
- 每个函数的详细说明
BLAS Technical Forum Standard
- BLAS接口的官方规范
ScaLAPACK Users’ Guide

在线资源

Netlib（netlib.org）：
- 源代码、文档、测试矩阵
TOP500（top500.org）：
- 历史榜单、趋势分析
Jack Dongarra主页（icl.utk.edu/~dongarra）：
- 论文、报告、演讲视频

历史文献

LINPACK Users’ Guide（1979）：
- Dongarra的早期工作
“The LINPACK Benchmark: Past, Present, and Future”（2003）：
- Dongarra回顾TOP500

🌟 精神遗产

”软件是科学的基础设施”

理念：

科学论文依赖软件，软件应该像实验设备一样可靠
数值不稳定的代码 = 损坏的望远镜

实践：

LAPACK经过几十年完善，极少bug
测试矩阵库（Matrix Market）：已知答案的标准问题

对今天的启示：

开源科学软件应该得到资助和认可
“可重复性危机”的部分原因是软件质量差

”标准让每个人都能站在巨人肩上”

Dongarra名言：

“如果每个人都重新发明高斯消元，我们永远到不了月球。”

BLAS的哲学：

硬件厂商优化底层（Intel MKL、NVIDIA cuBLAS）
应用科学家调用标准接口
分工合作，而非各自为战