SLISC 库简介

                     

贡献者： addis

　　 相对于 Fortran 或者 Matlab 等，用 C++ 做数值计算的一个缺陷就是语言本身没有矩阵以及高维矩阵类。但我们可以用第三方库或者自己写一个。我们选择后者，原因是为了教学需要我们需要保持代码的简单易读，避免复杂的 C++ 语法。

　　 本部分中我们将大量使用自编的 SLISC（Scientific Library in Simple C++） 矩阵库。代码可以从 GitHub 仓库下载，另外提供英文说明。该库的特点是尽量不使用 C++ 的复杂语法（如模板）和复杂的类结构，使代码便于阅读学习和修改，同时又保持相对较高的性能。SLISC 使用兼容性较高的 C++11 标准。

　　 SLISC 库在矩阵类的基础上还实现了一些科学计算常用的功能，如线性代数运算、插值、微分方程、计时、文件和目录操作、矩阵文件读写、字符串处理、特殊函数、任意精度计算、量子力学计算等。当然，一些功能基于其他项目如 Blas、Lapack、Intel MKL、Arb、GSL、Eigen 和 Arpack。在编译时可以通过选项来决定是否使用这些依赖。

　　 一个例子（编译方法见 “SLISC 的编译和测试”）：

#include "SLISC/arithmetic.h"
#include "SLISC/disp.h"
#include "SLISC/cut.h"
#include "SLISC/matb.h"
int main()
{
	using namespace slisc;

	// === 矢量和矩阵 ===
	VecDoub u(3), v(3); // 三个元素的双精度矢量
	linspace(u, 0, 2); // u = [0,1,2]
	cout << "u = \n"; disp(u); // 显示矩阵
	copy(v, 3.14); // 所有元素 = 3.14
	copy(u, v); // 复制
	u += v; // 矢量 + 矢量
	v += 2; // 矢量 + 标量
	CmatDoub a, b(2, 3); // 双精度矩阵（行主序）
	a.resize(2, 3); // 改变大小（重新分配内存）
	a(0, 0) = 1.1; // 双指标
	a[3] = 9.9; // 单指标
	a.end() = 5.5; a.end(2) = 4.4; // 最后和倒数第二个矩阵元
	cout << "a 有 " << a.n0() << " 行和 " << a.n1()
	<< " 列, 共计 " << a.size() << " 个元素。" << endl;
	disp(a); // display

	// === 矩阵切割 ===
	SvecDoub bs = cut0(b, 0); // 切割矩阵的第一列
	bs[0] += 1; // 切割类型的使用方法和矩阵类似
	Doub s2 = sum(bs); // 求和
	DvecDoub bd = cut1(b, 1); // 切割第二行
	copy(cut1(b, 0), cut1(b, 1)); // 复制第二行到第一行

	// === 文件读写 ===
	save_matb(b, "b", "data0.matb"); // 把矩阵存为二进制文件
	// 新建二进制文件储存多个变量
	Matb matb("data.matb", "w");
	save(u, "u", matb); save(v, "v", matb);
	save(a, "a", matb); save(b, "b", matb);
	matb.close();
}

1. 主要特点

　　 以下列出 SLISC 的一些特性，我们以后会详细介绍。

• 全部定义使用 namespace slisc，所有宏以 SLS_ 开头。
• 尽量不使用模板（尤其是高级的模板技巧），用 Matlab/Octave 生成代码。初学用户不需要了解代码如何生成，可以直接阅读或使用生成后的代码。
• 实现了密矩阵和一些稀疏矩阵。以对列主序的支持为主。
• 实现了密矩阵的切割，并用于函数接口。这会有少量的额外运算，但使用起来却比指针方便得多。用户可以自行选择使用哪种接口。
• Debug 模式实现了详尽的尺寸和指标越界检查，可以及时发现指标超出长度等错误。
• 函数尽量首先使用和 Lapack 类似的指针接口实现，然后再封装一层更友好的接口。
• 为了保证性能，函数内部尽量不改变矩阵尺寸，这是因为内存的动态分配往往耗时较多。

　　 在 SLISC 中，我们希望把矩阵进行一定程度的封装，但又几乎不损失性能。许多人使用流行的矩阵库例如 Eigen，但是其代码复杂，错误信息不容易懂，代码修改起来非常困难。例如 Eigen 的 MatrixXd 矩阵类有 6 次继承，一个 2x2 的矩阵在 gdb 里面调试的时候显示出来的变量信息是这样的：

<Eigen::PlainObjectBase<Eigen::Matrix<double,-1,-1,0,-1,-1>>> = 
{<Eigen::MatrixBase<Eigen::Matrix<double,-1,-1,0,-1,-1>>> = 
{<Eigen::DenseBase<Eigen::Matrix<double,-1,-1,0,-1,-1>>> = 
{<Eigen::DenseCoeffsBase<Eigen::Matrix<double,-1,-1,0,-1,-1>, 3>> = 
{<Eigen::DenseCoeffsBase<Eigen::Matrix<double,-1,-1,0,-1,-1>, 1>> = 
{<Eigen::DenseCoeffsBase<Eigen::Matrix<double,-1,-1,0,-1,-1>, 0>> = 
{<Eigen::EigenBase<Eigen::Matrix<double,-1,-1,0,-1,-1>>> =
{<No data fields>}, <No data fields>}, 
<No data fields>}, <No data fields>}, <No data fields>}, <No data fields>},
m_storage = {m_data = 0x855ceb0, m_rows = 2, m_cols = 2}}, <No data fields>

实际上这里面的重点只有最后一行，显示了矩阵在内存中的地址 m_data，行数 m_rows 以及列数 m_cols。这样的库只适合直接拿来用，不适合阅读、学习和修改，尤其是对非计算机专业的同学来说。

重新定义类型名称

　　 在 C++ 标准中，一些基本类型在内存中的长度在不同计算机上可能会不同。例如 long 有时候和 int 一样是 4 字节，而另一些时候则和 long long 一样是 8 字节。因此，SLISC 仿照 Numerical Recipe [1] 的方式，定义自己的类型名称。

• Char 是 char，即 1 字节字符或整数
• Uchar 是 unsigned char，即 1 字节无符号整数
• Short 是 2 字节整数
• Int 是 4 字节整数
• Long 是 Int 或者 Llong，取决于宏 SLS_USE_INT_AS_LONG 是否定义。SLISC 库中，Long 被用来作为矩阵的索引类型。
• Llong 是 8 字节整数
• Float 是 float（4 字节浮点数）
• Doub 是 double（8 字节浮点数）
• Ldoub 是 long double（16 字节浮点数，但并不一定是四精度，取决于机器和编译器）
• Qdoub 是 16 字节浮点数，标准的四精度（需要在编译时打开 quadmath 选项）。
• Fcomp 是 std::complex<float>（8 字节复数）
• Comp 是 std::complex<double>（16 字节复数）
• Lcomp 是 std::complex<long double>（32 字节复数）
• Qcomp 是 std::complex<Qdoub> 32 字节复数，标准的四精度。

　　 另外我们还定义一些附加类型，用于声明函数的参数。例如 Doub_I 是 const Doub，表示函数的输入参数。又如 Doub_O 是 Doub &，表示函数的输出参数（原来的值不会被使用）。Doub_IO 也是 Doub &，但表示既作为输入也作为输出（原来的值会被使用且被覆盖）。以上的标量类型都同理，在后面加 _I，_O，_IO 分别表示 in，out，in-out。

　　 以后还会看到矢量类和矩阵类，但和标量不同，为了避免不必要的复制，它们都必须 pass by reference 而不是 pass by value。例如矩阵类 CmatDoub 的 CmatDoub_I 是 const CmatDoub &，CmatDoub_O 和 CmatDoub_IO 都是 CmatDoub &。

---

[1] ^ W. H. Press, et al. Numerical Recipes 3rd edition.