Fortran 稀疏矩阵

Fortran 处理稀疏矩阵

稀疏矩阵Ax=b

在Fortran里面使用稀疏矩阵最基础的是用BLAS（Basic Linear Algebra Subprograms），但是在后来的MKL库中有集成BLAS。如果是解Ax=b的线性方程组，建议使用Pardiso，同样在MKL库中有集成，可以去官网查找资料。
这里给出两篇博文，介绍了如何在Fortran里面解线性方程组以及BLAS库函数的示例
https://blog.csdn.net/chd_lkl/article/details/83011186
https://blog.csdn.net/chd_lkl/article/details/107433757

目前查到的除上述表格中的MKL库函数之外，还有SPARSEKIT库函数，但是该库函数是在GNU下运行的，运行之前还需要安装该函数包，可以参考SPARSEKIT网站。

以及SOL（Systems Optimization Laboratory）求解器，参考网站http://stanford.edu/group/SOL/download.html给出一个解线性方程组的包，里面还有LU分解函数。

这里再附一个老师提供的稀疏矩阵求解器GMERS
参考网站，主要功能就是两个：

• 一个是对COO格式稀疏矩阵直接迭代求解Ax=b 使用函数 mgmres_st
• 另一个是对CSR格式稀疏矩阵通过ILU方法求解Ax=b，使用函数pmgmres_ilu_cr
  测试结果来看，明显后者速度要快不少，如果只有COO格式的稀疏矩阵，可以通过MKL的库函数 mkl_dcsrcoo转换，具体的使用规则参考MKL库函数，还是之前给的链接。这个求解器里面还有一个求稀疏矩阵乘向量的小函数比较好用 ax_st ax_cr注意这里的向量是稠密的向量，st cr分别表示稀疏矩阵的存储格式是COO还是CSR格式。

稀疏矩阵存储方式及运算

稀疏矩阵的存储方式一般有两种：

• COO 存三元组（行标、列标、值） 这种方法就很直观
• CSR 存三元组（值，列标，行指标）
  其中CSR的rowindex表示每行第一个非0元素是第几个非0元素，rowindex最后一个表示最后一个非0元素是第几个。
  于是rowindex的长度为行数+1，rowindex最后一个元素为非0元素的个数。

需要用到的mkl处理稀疏矩阵的函数主要参考文档：
里面给出了各种函数的用法，主要包括：

功能	函数
稀疏矩阵乘一般向量	mkl_?csrgemv
解线性方程组Ax=b	mkl_?csrsv
稀疏矩阵乘一般矩阵	mkl_?csrmm

需要注意的是，稀疏矩阵之间的加法、乘法这里没有给出具体的函数处理，需要自己去补充写。。。

这里补充一个讲稀疏矩阵乘法的算法博客
该算法就是两个嵌套的for循环分别遍历两个矩阵的非0元素

关于稀疏矩阵的乘法目前在网上能查到的基本都是做两次遍历，外循环遍历矩阵A的所有非0元，内循环遍历B的非0元，这样时间复杂度就是$O(nu{m}^2)$但是对于我们的问题（二维非局域热传导方程DG方法），$num=C{N}_x{N}_y{k}^3$，当网格加密一倍，即$N_x,N_y$变为原来的两倍，则时间复杂度变为原来的16倍，这对于大规模矩阵运算是难以接受的。但实际上我们的稀疏矩阵是接近于三对角阵甚至对角阵，实际计算起来应该不太难，但是由于存储稀疏矩阵用COO格式或者CSR格式，并不能及时找到对应行列的元素，所以一般的稀疏矩阵乘法对于我们这的特殊矩阵并没有明显的优势。

补充

再补充几个编程中遇到的问题：

• Fortran开辟动态数组的时候，需要给一个指定的长度，然后传参都是传地址，所以对于长度不固定的情况，目前的处理方法就是最简单无脑的给一个足够长的空间，由于传参传的是地址，如果有函数计算A-B并赋给C
  Minus(A, B, C)
  但此时如果需要计算A-B并赋给A，写成Minus(A, B, A)就会出问题，最好还是别这么写
• Fortran里面声明完一个变量之后如果用的话尽量给它初始化！！！因为Fortran默认会给一些很奇怪的数，不管是integer变量还是real变量！！！