稀疏矩阵线性解析库SPOOLES的简单应用

稀疏矩阵线性解析库SPOOLES的简单应用
烤鱼片(@eii.dlmu)
cleverysm@163.com

SPOOLES的全称是SParse Object Oriented Linear Equations Solver，中文名大概就是面向对象的稀疏线性等式解析器。顾名思义，就是可以用来解稀疏矩阵为参数的线性方程组的数学函数库。所谓面向对象是指的应用了面向对象的封装思想，但实际上SPOOLES是用非面向对象的C语言来写的。
SPOOLES的主页是http://www.netlib.org/linalg/spooles/spooles.2.2.html ，最新版是2.2，支持单线程，多线程和MPI三种计算模式。主页上的文档也比较齐全，对SPOOLES的使用介绍的也很清晰，仔细看看用不了很长时间就能解自己的方程组。如果说只是拿来用的话，Install.ps.gz ，ReferenceManual.ps.gz ，LinSol.ps.gz ， AllInOne.ps.gz 这四个文档就足够用了，Install.ps.gz 不用多说是安装指南，ps格式，down下来转换成pdf看就行，其他几个也一样。ReferenceManual.ps.gz 是SPOOLES的功能参考手册，AllInOne.ps.gz 是比较详细的介绍了在单线程，多线程和MPI环境下应用SPOOLES解线性方程AX=Y的步骤，包括LU和QR两种因式分解的方法。不过AllInOne.ps.gz 里的解决方法步骤比较多，调用的函数也多了一点，所以就有人作了一个LinSol的库，也就是LinSol.ps.gz 中介绍的库。我在本文中介绍的也是应用LinSol的解决方法。

Spooles的安装

首先是安装SPOOLES。从网站的down下来，解压到某文件夹下，比如我们就解压缩到/home/mpi/spooles下。编译代码之前要修改一下Make.inc文件的某些内容，Make.inc这个文件会在makefile里include进去。需要改的地方有设定编译器，比如设为CC = gcc，MPI_INSTALL_DIR为mpi安装的路径，如在我的机器上是MPI_INSTALL_DIR = /home/mpi/mpich2，mpi的库位置MPI_LIB_PATH，如MPI_LIB_PATH = -L$(MPI_INSTALL_DIR)/lib，mpi的库MPI_LIBS可以按照自己需要修改，mpi的头文件地址，MPI_INCLUDE_DIR = -I$(MPI_INSTALL_DIR)/include。然后就可以编译程序了。
如果需要使用单线程的spooles库，只需在存放spooles的根目录下，make lib，这个命令将编译所需要的代码在spooles的根目录下生成一个spooles.a的静态库。不过make lib需要使用一个perl脚本，所以机器上需要安装了perl。如果机器没有perl的话，可以使用make global。不过在我的机器上make global没有成功，而且官方的安装文档也说make global不太保险。如果需要多线程版的spooles的话，就需要到MT/src目录下，make spoolesMT.a，就可以在这个目录下生成库文件spoolesMT.a，如果想把多线程库编译到根目录的spooles.a中的话，可以用，make makeLib。但是，官方安装文档中推荐使用第一种方式。MPI版的情况与多线程的类似，存放目录是在MPI/src中。
不过，直接使用spooles的api来计算的话相对来说有点烦琐，需要调用的函数和各种参数比较繁杂一点，所以在spooles的发行版中还附带了一个封装好了的辅助库，代码放在在LinSol目录下，叫做Wrapper Objects for Spooles，或者也称作是一种spooles的bridge，我就是用的这个工具。同样的也分为单线程，多线程和MPI版，分别在LinSol目录下的srcST，srcMT和srcMPI目录下。我们只需要在LinSol目录下make lib就可以分别编译生成这三个版本对应的静态库文件，以单线程版为例，就是Bridge.a。

解方程
简单的讲，用spooles的LinSol解一个AX=B方程可分为四个基本的步骤，初始化，因式分解，解出结果，释放资源。

初始化

初始化的过程首先要构造矩阵A，B矩阵和X。系数矩阵A是按照稀疏矩阵的模式来存储的，也就是仅仅保存非零元素既可，而B和X则是稠密矩阵。

a) 构造矩阵A

矩阵A使用InpMtx对象来存储。代码如下。
InpMtx *mtxA;
mtxA = InpMtx_new() ;
InpMtx_init(mtxA, coordType, inputMode, 0, 0) ;
参数coordType是int型变量，指的是数据在对象中的存储模式，包括三种，INPMTX_BY_ CHEVRONS，INPMTX_BY_COLUMNS，INPMTX_BY_ROWS，一般来说用INPMTX_BY_ROWS就可以。
参数inputMode也是int型变量，指的是在对象中的存储的数据类型，包括三种，仅存储索引INPMTX_INDICES_ONLY，实数SPOOLES_REAL，复数SPOOLES_COMPLEX，在我的程序里使用的是一般的实数，所以就是用的SPOOLES_REAL。
然后，我们需要将我们的数据填充到mtxA对象中。针对三种不同的inputMode也就对应有三种不同的方法来给mtxA赋值，分别可以从方法的名称体现出来，如void InpMtx_inputEntry，void InpMtx_inputRealEntry，void InpMtx_inputComplexEntry。每一种模式下又有四种不同的操作方法，以实数为例，分别是
void InpMtx_inputRealEntry ( InpMtx *mtxA, int row, int col, double value ) ;
每次赋值一个元素，row和col分别指的行列坐标，value则是元素值
void InpMtx_inputRealRow ( InpMtx *mtxA, int row, int rowsize,int rowind[], double rowent[] ) ;
每次赋值一行元素，row指的行坐标，rowsize是非零元素个数，rowind是存储非零元素的对应索引的数组（从0开始），rowent是对应于rowind的元素值数组。
void InpMtx_inputRealColumn ( InpMtx *mtxA, int col, int colsize,int colind[], double colent[] ) ;
每次赋值一列元素，参数含义类似于InpMtx_inputRealRow
void InpMtx_inputRealMatrix ( InpMtx *mtxA, int nrow, int ncol, int rowstride,int colstride, int rowind[], colind[], double mtxent[] ) ;
没用过，大家自己研究吧：－）
以InpMtx_inputRealColumn为例，我们来插入下面这个矩阵的第一列
1 2 4
0 1 0
3 0 5
int col=0;
int colsize=2;
int colind[2]={0,2};
double colent[2]={1,3};
InpMtx_inputRealColumn (mtxA, col, colsize, colind, colent ) ;

b) 构造矩阵B和X

矩阵B和X使用DenseMtx对象来存储。代码如下。
int type, rowid, colid, nrow, ncol, inc1, inc2 ;
DenseMtx *mtx = DenseMtx_new() ;
DenseMtx_init(mtx, type, rowid, colid, nrow, ncol, inc1, inc2) ;
Type指的是数据类型，分实数SPOOLES_REAL，复数SPOOLES_COMPLEX两种。
rowid, colid是为了标示所创建的这个矩阵是某个矩阵的子矩阵，如果不是子矩阵，则全为0。
nrow, ncol表示在这个矩阵中的行数和列数
inc1, inc2的赋值则取决于矩阵是行主序存储还是列主序存储。行主序时inc1 = ncol 且 inc2 = 1，列主序时inc1 = 1 且 inc2 = nrow。
例如，我们可以如下来初始化一个列主序稠密矩阵。
DenseMtx_init(mtx, SPOOLES_REAL, 0, 0, 100, 5, 1, 100) ;
对DenseMtx对象中数据的读写，我们可以简单的使用如下几个方法。
double value, real, imag ;
int irow, jcol ;
DenseMtx_realEntry(mtx, irow, jcol, &value) ;
DenseMtx_complexEntry(mtx, irow, jcol, &real, &imag) ;
DenseMtx_setRealEntry(mtx, irow, jcol, value) ;
DenseMtx_setComplexEntry(mtx, irow, jcol, real , imag ) ;以实数情况为例，irow和jcol分别是行列坐标（从0开始），value就是具体的元素值了。

比如像下面这一段代码
mtxY = DenseMtx_new();
DenseMtx_init(mtxY, SPOOLES_REAL, 0, 0, n1*n2, 1, 1, n1*n2) ;
DenseMtx_zero(mtxY) ;
//注意这里，是给稠密矩阵mtxY所有元素赋值为0
//最好在每次对mtxY重新赋值之前都这么做
ii = n1 - 1 ;
for ( jj = 0 ; jj < n2 ; jj++ ) {
ij = ii + jj*n1 ;
DenseMtx_setRealEntry(mtxY, ij, 1, 1.0) ;
}
c) bridge对象的初始化
bridge的初始化比较简单，只需要下面几个方法。
bridge = Bridge_new() ;
Bridge_setMatrixParams(bridge, neqns, type,) ;
//neqns指的是等式的数量
//type等同于上面提到的inputMode
//symmetryflag是矩阵的对称性
//可以是SPOOLES_SYMMETRIC，SPOOLES_HERMITIAN或者SPOOLES_NONSYMMETRIC
//默认是SPOOLES_SYMMETRIC
Bridge_setMessageInfo(bridge, msglvl, msgFile) ;
//msglvl是指的spooles对外输出的辅助信息类型，一般1就可以
//msgFile是一个FILE指针，spooles会根据msglvl的类型，将一些辅助信息输出到msgFile中
//msgFile可以简单的设为stdout，就是打印到标准输出屏幕上
rc = Bridge_setup(bridge, mtxA) ;
因式分解
因式分解的过程很简单，只需要
int permuteflag(1), error;
int rc = Bridge_factor(bridge, mtxA, permuteflag, &error) ;
// mtxA系数矩阵A
//permuteflag是有关矩阵内部排序的参数，一般就是1，具体细节可以不用管它
//error，错误输出的信息
//如果正常完成因式分解过程，则rc为1
解得结果
也很简单。
rc = Bridge_solve(bridge, permuteflag, mtxX, mtxY) ;
// mtxY，B矩阵
// mtxX，X结果矩阵
//permuteflag含义同Bridge_factor中的
然后，方程的结果就会被存放在mtxX中，按照上面介绍的操作稠密矩阵的方法，就可以把答案取出来了。

释放资源
就是把通过*_new()创建的对象都释放调，如
Bridge_free(bridge) ;
InpMtx_free(mtxA) ;
DenseMtx_free(mtxX) ;
DenseMtx_free(mtxY) ;
编译链接
编译没什么特殊的，链接的时候记得把需要的静态库加进去，再就是如果程序是用C++写的，在引入spooles的头文件的时候需要加上extern “C”的标示，因为spooles是用C语言写的，用C++链接的时候就必须生命按照C外部函数来链接。
如：
extern "C"
{
#include <LinSol/Bridge.h>
#include <LinSol/BridgeMPI.h>
}

例程，结方程AX=B
A矩阵
10 -1 -2
-1 10 -2
-1 -1 5
B矩阵
72
83
42
X矩阵
11
12
13

int AXB()
{
InpMtx *mtxA;
Bridge *bridge;
DenseMtx *mtxX,*mtxY;
int error,rc,count(3);
int msglvl(1);
int inputMode(SPOOLES_REAL),coordType(INPMTX_BY_COLUMNS);
int permuteflag(1) ;
int indexs[3];
double values[3];

mtxA = InpMtx_new() ;
InpMtx_init(mtxA, coordType, inputMode, 0, 0) ;
mtxY = DenseMtx_new() ;
DenseMtx_init(mtxY, inputMode, 0, 0, count, 1, 1, count) ;
mtxX = DenseMtx_new() ;
DenseMtx_init(mtxX, inputMode, 0, 0, count, 1, 1, count) ;
//给矩阵A赋值
indexs[0]=0;
indexs[1]=1;
indexs[2]=2;
values[0]=10;
values[1]=-1;
values[2]=-1;
InpMtx_inputRealColumn ( mtxA, 0, count,indexs, values) ;
indexs[0]=0;
indexs[1]=1;
indexs[2]=2;
values[0]=-1;
values[1]=10;
values[2]=-1;
InpMtx_inputRealColumn ( mtxA, 1, count,indexs, values) ;
indexs[0]=0;
indexs[1]=1;
indexs[2]=2;
values[0]=-2;
values[1]=-2;
values[2]=5;
InpMtx_inputRealColumn ( mtxA, 2, count,indexs, values) ;
//将A矩阵数据打印到标准输出
InpMtx_writeForHumanEye(mtxA, stdout) ;
//给矩阵B赋值
values[0]=10;
values[1]=-1;
values[2]=-1;
DenseMtx_zero(mtxY) ;
for(long i=0;i<count;i++)
{
DenseMtx_setRealEntry(mtxY, i, 0, values[i]) ;
}
//将Y矩阵数据打印到标准输出
DenseMtx_writeForHumanEye(mtxY, stdout) ;
//初始化bridge
bridge = Bridge_new() ;
Bridge_setMatrixParams(bridge, count, inputMode, SPOOLES_NONSYMMETRIC) ;
Bridge_setMessageInfo(bridge, msglvl, stdout) ;
rc = Bridge_setup(bridge, mtxA) ;
//因式分解
rc = Bridge_factor(bridge, mtxA, permuteflag, &error) ;
//解方程
DenseMtx_zero(mtxX) ;
rc = Bridge_solve(bridge, permuteflag, mtxX, mtxY) ;
//将X矩阵数据打印到标准输出
DenseMtx_writeForHumanEye(mtxX, stdout) ;
//释放资源
Bridge_free(bridge) ;
InpMtx_free(mtxA) ;
DenseMtx_free(mtxX) ;
DenseMtx_free(mtxY) ;

return 0;
}