简单来说，计算机由计算单元、储存单元和二者之间的通信构成。

• 下载 源码

向量是 PETSc 最基本的结构，用于保存计算结果等

PETSc 支持多种矩阵类型，包括稠密矩阵、稀疏矩阵等

使用线性求解器 KSP 解线性方程 \(Ax=b\)

1. KSPCreate(MPI_Comm comm,KSP *ksp);
   1. 初始化 KSP
2. KSPSetOperators(KSP ksp,Mat Amat,Mat Pmat);
   1. 设置矩阵
3. KSPSetFromOptions(KSP ksp);
   1. 设置求解器选项
4. KSPSolve(KSP ksp,Vec b,Vec x);
   1. 运行求解器
5. KSPDestroy(KSP ksp);
   1. 运算结束销毁求解器

只能通过编译选项 --with-scalar-type=complex 将 PetscScalar 设成复数

参考

debug 模式：

   ./configure PETSC_ARCH=arch-complex-debug --with-scalar-type=complex

release 模式:

   ./configure PETSC_ARCH=arch-complex-release --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --with-debugging=0 COPTFLAGS='-O3 -march=native -mtune=native' CXXOPTFLAGS='-O3 -march=native -mtune=native' FOPTFLAGS='-O3 -march=native -mtune=native' --with-scalar-type=complex --download-mpich

64 位支持：

   --with-64-bit-indices

python 支持：

   –with-petsc4py=1

PETSc 是纯 MPI 框架，不提供混合线程的函数，同时所有函数也不保证线程安全。

在调用 PetscInitialize() 之前设置

   PETSC_COMM_WORLD = MPI_COMM_WORLD; // 或其它 comm

依赖越复杂，莫名其妙的问题就越多，只需要 c 的话就简单不少。

./configure PETSC_ARCH=arch-complex-release --with-cc=mpicc --with-cxx=0 --with-fc=0 --with-debugging=0 COPTFLAGS='-O3 -march=native -mtune=native'  --with-scalar-type=complex --with-64-bit-indices  --with-make-dir=$HOME/app/make --with-blas-lapack-dir=/opt/intel/composer_xe_2015.1.133/mkl/lib/intel64

• COO: Coordinate or triplet
• CSR: Compressed Sparse Row
• CSC: Compressed Sparse Column
• DIA: Diagonal
• BSR: Block Sparse Row

1. 也叫做 triplet 格式
2. 最简单和基本的格式，通常用做向其它格式转换的入门格式
3. 顺序是任意的
4. 保存三列数据：值、行指标、列指标

A:
   [[ 1,  0,  0,  2,  0],
    [ 3,  4,  0,  5,  0],
    [ 6,  0,  7,  8,  9],
    [ 0,  0, 10, 11,  0],
    [ 0,  0,  0,  0, 12]]

data:   [1, 3, 6, 4, 7, 10, 2, 5, 8, 11, 9, 12]
row(i): [0, 0, 1, 1, 1,  2, 2, 2, 2,  3, 3,  4]
col(j): [0, 3, 0, 1, 3,  0, 2, 3, 4,  2, 3,  4]

A[ i[k], j[k] ] = data[k]

1. 按行压缩，即coo格式中的数据值、列指标不变，但行指标改成指向这一行开始的列的指标
2. 一般来说，Fortran 里使用这种格式

A:
   [[ 1,  0,  0,  2,  0],
    [ 3,  4,  0,  5,  0],
    [ 6,  0,  7,  8,  9],
    [ 0,  0, 10, 11,  0],
    [ 0,  0,  0,  0, 12]]

data:   [1, 3, 6, 4, 7, 10, 2, 5, 8, 11, 9, 12]
row(i): [0, 2, 5, 9, 11, 12]
col(j): [0, 3, 0, 1, 3,  0, 2, 3, 4,  2, 3,  4]

1. 按列压缩，与 CSR 类似，不过行列换一下
2. c 中常用这种格式

A:
   [[ 1,  0,  0,  2,  0],
    [ 3,  4,  0,  5,  0],
    [ 6,  0,  7,  8,  9],
    [ 0,  0, 10, 11,  0],
    [ 0,  0,  0,  0, 12]]

data:   [1, 3, 6, 4, 7, 10, 2, 5, 8, 11, 9, 12]
row(i): [0, 1, 2, 1, 2,  3, 0, 1, 2,  3, 2,  4]
col(j): [0, 3, 4, 6, 10, 12]

1. 按对角线压缩，保存偏移值，负的偏移忽略头元素，正的偏移忽略尾元素

A:
   [[ 1.,  0.,  2.,  0.,  0.],
    [ 3.,  4.,  0.,  5.,  0.],
    [ 0.,  6.,  7.,  0.,  8.],
    [ 0.,  0.,  9., 10.,  0.],
    [ 0.,  0.,  0., 11., 12.]]

data:
   [[ 3.,  6.,  9., 11., nan],
    [ 1.,  4.,  7., 10., 12.],
    [nan, nan,  2.,  5.,  8.]]

offsets: [-1,  0,  2]

1. 类似 CSR，只不过每个行列指的是一个小的矩阵，这些小矩阵大小都相等

1. 很容易使用，参考 scipy - sparse

1. 稀疏矩阵指的是那些有很多元素为 0 的矩阵
2. 一般而言 \(m\times n\) 的矩阵， 其非零元素的个数在 \(O(\mathrm{min}(m,n))\) 就可以作为稀疏矩阵来处理，一般来说就是每行或每列的元素个数是常数
3. 很多非零元是 \(O(\mathrm{log}(n))\) 的矩阵不能作为稀疏矩阵

1. 稀疏矩阵 \(A\) 的逆 \(A^{ -1}\) 一般来说是稠密的
2. 稀疏矩阵 \(A\) 的 LU 分解 \(L\) 和 \(U\) 矩阵可以也是稀疏的

矩阵满足其特征方程。

例：考虑矩阵 \(A = \begin{pmatrix} 1&2 \\ 3&4 \end{pmatrix}\) ，其特征方程为

那么有

\[ A^{2} - 5 A - 2 = 0 \]

稀疏矩阵可以用图来表示，进而可以用图论的方法分析。

图的定义 图 \(G\) 定义为一个集合 \(G = (V,E)\), \(E \subset V \times V\) ，其中 \(V\) 是顶点 (vertex) 的集合, \(E\) 是边 (edge) 的集合。

1. 图是无向图，当且仅当矩阵是对称的。
2. 矩阵的图表示，就是将行和列作为顶点的标号，非零矩阵元表示顶点间的连线。

典型的偏微分方程数值模拟的过程是：

1. 物理问题
2. 非线性偏微分方程
3. 离散化
4. 线性化近似
5. 变成稀疏矩阵线性代数问题

可以想象成将自变量作为矩阵指标，离散化的偏微分方程每一步的作用结果就是利用周围的几个点的数值得到下一步当前点的值，这个过程可以写成一个稀疏矩阵。

1. 由一个下划线和一个大写字母或两个下划线开头的标识符是标准库保留的
2. 由一个下划线开头的是文件作用域保留的
3. 标准库中的外部链接标识符都是保留的

流(stream) 表示任意输入/输出的源或目的地，一般程序都是从一个流获得输入，再通过另一个流输出。

当发生错误时，大多数函数会将一个错误码存储到 errno 变量中。此外，如果函数返回值大于double类型最大值会返回 HUGE_VAL 值。

1. 定义域错误。当函数的参数超出定义域，会将 EDOM 存储到 errno 中
2. 取值范围错误。当函数返回值超出double范围时，会将 ERANGE 存储到 errno 中

C99中新增加了融合乘加 (fused multiply-add) 函数，即

   a = b * c + d;
   a = fma(b, c, d);

这种合并可能会速度更快一点，编译器是否自动进行紧缩可以由 #pragma STDC FP_CONTRACT ON/OFF/DEFAULT 来控制。

提供了5种函数

1. 复制: memcpy, memmove, strcpy, strncpy
2. 拼接: strcat, strncat
3. 比较: memcmp, strcmp, strcoll, strncmp, strxfrm
4. 搜索: memchr, strchr, strcpn, strpbrk, strrchr, strspn, strstr, strtok
5. 其它（初始化、长度）: memset, strlen

assert 是一个宏。当参数值为 0 时，assert 会向 stderr 写消息，并调用 abort 函数中止程序。

   void assert(scalar expression);

错误代码存储在 errno 变量中，每次使用都要把它置零

errno 主要作用是说明错误类型，而不是发生错误

   /* 设置跳转位置 */
   int setjmp(jmp_buf env);

   /* 跳转到 val 位置 */
   void longjmp(jmp_buf env, int val);

可以实现跨函数的跳转

加载数据(load data)，需要完成如下几个工作

1. 读取和准备数据
2. 对数据归一化
3. 划分训练数据和测试数据

数据应该是有 \(x_{1},x_{2},\cdots,x_{n}\) 和一组 \(y\)

记住要对数据归一化, 可以选择基本的移动归一化，即减去最小值再除以最大最小之差 \[ x' = \frac{x - \mathrm{min}(x)}{\mathrm{max}(x) - \mathrm{min}(x)} \]

或者减平均值除方差

\[ x' = \frac{x - \mathrm{mean}(x)}{\mathrm{std}(x)} \]

前向传播(forward propagation) 就是从输入参数计算得到输出参数。

线性回归公式是 \[ y = \sum_{j=1}^{M} x_{j}w_{j} + b \]

所以前向传播很容易写

   def forward(x, w, b):
       result = np.dot(x, w) + b
       return result

损失函数(loss function) 是衡量前向计算结果与数据目标结果之间偏差的函数。

简单的损失函数是均方误差 \[ \mathrm{Loss} = \frac{1}{N} \sum_{j=1}^{N} (y_{i} - z_{i})^{2} \] 其中 \(y_{i}\) 是目标结果， \(z_{i}\) 是计算结果。

   def loss(z, y):
       result = np.mean((z - y) ** 2)
       return result

梯度(gradient) 是损失函数随参数变化的改变。

对于这里的线性回归模型和均方误差损失函数，梯度定义如下 \[ \nabla \mathrm{Loss} = \left( \frac{\partial \mathrm{Loss}}{\partial w_{i}}, \dots \right) \]

对于我们要计算的损失函数 \[ L = \frac{1}{2N} \sum_{i} (y_{i} - z_{i})^{2} \] 预测值为 \[ z_{i} = \sum_{j} x_{i}^{j} w^{j} + b \] 其中 \(i\) 是样本的指标，\(j\) 是参数的指标。 那么损失函数 \(L\) 对参数 \(w_{j}\) 的微分为 \[ \frac{\partial L}{\partial w_{j}} = \frac{ 1 }{ N } \sum_{i} (z_{i}- y_{i}) \frac{\partial z_{i }}{\partial w_{j}} =\frac{ 1 }{ N } \sum_{i} (z_{i}- y_{i}) x_{i}^{j} \]

   def gradient(x, y, z, w, b):
       gradient_w = np.mean((z - y) * x, axis=1)
       gradient_b = np.mean((z - y))
       return gradient_w, gradient_b

梯度下降(gradient descent) 就是根据梯度来更新参数，最终到达收敛的过程。

在这部分要实现两个函数

1. 更新(update)，即从这一步梯度和学习率(learning rate)计算下一步的梯度
2. 训练(train)，即迭代执行前向计算-梯度-更新，并根据损失函数判断是否收敛

梯度下降的更新函数就是用上一步的参数减去学习率与梯度之积 \[ w_{i+1}^{j} = w_{i}^{j} - \eta \frac{\partial \mathrm{Loss}}{\partial w_{i}^{j}} \]

   def update(w, b, gradient_w, gradient_b, eta):
       w_new = w - eta * gradient_w
       b_new = b - eta * gradient_b
       return w_new, b_new

   def train(x, y, w, b, eta, N):
       L_list = []
       for i in range(N):
           z = forward(x, w, b)
           L = loss(z, y)
           gradient_w, gradient_b = gradient(x, y, z, w, b)
           w, b = update(w, b, gradient_w, gradient_b, eta)
           L_list.append(L)
       return w, b, L_list

随机梯度下降(Stochastic Gradient Descent, SGD) 是指每步迭代时从总训练集抽取一小部分来计算梯度，通过这种方式能够加快训练。 一些概念：

• mini batch: 每次迭代时用到的那一小部分数据集
• batch size: 一个 mini batch 中的样本数
• epoch: 迭代时不断抽取 mini batch，当取过一遍整个数据集时就叫做一轮 epoch

具体的操作：

1. 将总的数据集随机打乱
2. 将打乱后的数据集划分成若干个 mini batch
3. 用每个 mini batch 进行一次训练
4. 用所有 mini batch 训练过一遍后，返回第 1 步，开启下一轮 epoch

在一个 git 项目目录中运行

   dvc init
   git commit -m "Initialize DVC"

就自动创建好了 dvc 需要东西。

dvc 是寄生在git中的，其本身不提供版本管理能力，完全依靠 git 。

如果有一个大的数据文件 data/data.h5 ，那么就首先将它添加给 dvc， 然后把dvc创建的记录 *.dvc 和自动生成的 .gitignore 一起加到 git 中，之后再用 dvc 上传数据即可

   dvc add data/data.h5
   git add data/data.h5.dvc data/.gitignore
   git commit -m "Add raw data"

dvc 支持多种远程库，甚至是本地的远程库

   mkdir -p /tmp/dvcstore
   dvc remote add -d myremote /tmp/dvcstore
   git commit .dvc/config -m "Configure local remote"

如果已经配置好了远程库，那么上传和下载数据操作就类似于 git

   dvc pull

   dvc push