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.

致读者: 小时百科一直以来坚持所有内容免费无广告,这导致我们处于严重的亏损状态。 长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。 因此,我们请求广大读者热心打赏 ,使网站得以健康发展。 如果看到这条信息的每位读者能慷慨打赏 20 元,我们一周就能脱离亏损, 并在接下来的一年里向所有读者继续免费提供优质内容。 但遗憾的是只有不到 1% 的读者愿意捐款, 他们的付出帮助了 99% 的读者免费获取知识, 我们在此表示感谢。

                     

友情链接: 超理论坛 | ©小时科技 保留一切权利