SLISC 库简介

                     

贡献者: addis

预备知识 C++ 基础

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

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

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

   一个例子(编译方法见 “SLISC 的编译和测试”):

代码 1:intro.cpp
#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 的一些特性,我们以后会详细介绍。

   在 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] 的方式,定义自己的类型名称。

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

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


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

                     

© 小时科技 保留一切权利