c++矩阵转置_<读书笔记4> 稀疏矩阵基础算法

本篇为稀疏矩阵求解算法经典论著<Direct Methods for Sparse Linear System>的<读书笔记 4>

Chapter 2 Basic algorithms

上一篇主要介绍了稀疏矩阵的数据结构，那么在接下来本章主要包括了其几个基础的算法。

1. Matrix-vector multiplication

对于矩阵-向量乘法操作

,其中

和

是稠密向量，

为稀疏矩阵。

其伪代码为

这个伪代码是基于CSC格式的，第一层for循环实际上是以x向量展开的。以CSR格式进行计算更容易理解。

C代码为：

int 

2.3 Utilities

此部分主要是一些用于内存空间处理的应用函数。不展开介绍，比较简单。这里的函数仅仅是创建了内存空间，但是没有存储空间中没有实际内容。

函数	用途
void *cs_malloc (int n, size_t size)	内存分配
void *cs_calloc (int n, size_t size)	内存分配并清零
void *cs_free (void *p)	清楚分配空间
void *cs_realloc (void *p, int n, size_t size, int *ok)	改变一块存储空间p的size
cs *cs_spalloc (int m, int n, int nzmax, int values, int triplet)	创建m*n稀疏矩阵，用于保存nzmax个entries,参数triplet用于设定是COO还是CSC格式
int cs_spfree (cs *A, int nzmax)	释放创建的稀疏矩阵空间
int cs_sprealloc (cs *A, int nzmax)	依据nzmax,重新对稀疏矩阵A分配空间

2.4 Triplet form

一般情况下，需要从文件中或许原始稀疏矩阵数据，常用的数据文件，例如mtx文件等，采用的是COO（或者在本书中称 Triplet form）。程序需要读入此类文件，获取数据，然后再得到这个数据的CSC，或者CSR格式。

如下面的代码，首先将COO格式的数据，存放再Ti，Tj，Tx中，然后再转成CSC格式的数据，得到对应的p，以及i。如果存在具有相同行和列索引的多个entries，则对应的数值是所有这些重复entries的和。

cs 

函数	用途
int cs_entry (cs *T, int i, int j, double x)	将aij=x的数据存入到T中
cs *cs_compress (const cs *T)	将COO转变为CSC格式
cs *cs_done (cs *C, void *w, void *x, int ok)	释放w和x的内存空间
double cs_cumsum (int *p, int *c, int n)	计算累加和

2.5 Transpose

稀疏矩阵的转置，与稠密矩阵不同，更像是将矩阵的从CSC格式转变为CSR格式的转换。

函数	用途
cs *cs_transpose (const cs *A, int values)	稀疏矩阵A的转置

2.6 Summing up duplicate entries

有限元素方法生成的矩阵是元素（elements）的集合或密集的子矩阵。完备矩阵是元素(elements)的总和。如果两个元素构成同一个条目，则应该对它们的值进行求和。

函数	用途
int cs_dupl (cs *A)	对重复elements求和

2.7 Removing entries from a matrix

CSparse并不要求它的稀疏矩阵不包含数值上的零项，但它的MATLAB接口要做到这一点。cs_fkeep函数使用了一个函数指针fkeep (i ,j ,aij , other)作为参数，fkeep 用于评估A中的每一个entry。如果fkeep是true，则此entry保持。

函数	用途
int cs_fkeep (cs *A, int (*fkeep) (int, int, double, void *), void *other)	移除零entry
static int cs_nonzero (int i, int j, double aij, void *other)	判断entry非零
int cs_dropzeros (cs *A)	移除零

2.8 Matrix multiplication

在CSparse中，采用CSC格式，矩阵乘

，其中

为

，A为

，B为

，取A和B的列，每次计算出C的列。

,

,

分别代表矩阵

的列向量，则有

。将A分解为k列，将

分解为k个独立的entries。

矩阵C的非零特征由以下的定理得到：

定理2.1：

的非零特征是

的非零特征的并集，其中

需要满足

为非零。如果

表示,

的非零entries的行索引集合，则有

矩阵乘法均需要计算

以及

。

函数	用途
int cs_scatter (const cs *A, int j, double beta, int *w, double *x, int mark, cs *C, int nz)	对于给定的j，实现上述两个公式
cs *cs_multiply (const cs *A, const cs *B)	计算稀疏矩阵乘法C=AB

稀疏矩阵乘法的计算量为

，其中

为浮点运算性能。

2.9 Matrix addition

矩阵加

等价于

。所以函数

cs_add实际上看起来很像cs_multiply。

函数	用途
cs *cs_add (const cs *A, const cs *B, double alpha, double beta)	计算C=alpha A+beta B

2.10 Vector permutation

置换矩阵P实际上也是一个稀疏矩阵，它的每一行和列均有只有一个非零元素并且为1。所以置换矩阵P可以用一个稀疏矩阵P来表示。或者用一个长度为n的整型向量p来表示，p即称为置换向量，其中p[k]=i，意味着

。矩阵P是正交的，所以

。同样，置换矩阵的逆被定义为

pinv[i]=k，如果

。

函数	用途
int cs_pvec (const int *p, const double *b, double *x, int n)	计算x=Pb
int cs_ipvec (const int *p, const double *b, double *x, int n)	计算x=PTb
int *cs_pinv (int const *p, int n)	计算p的逆

2.11 Matrix permutation

函数cs_permute用于计算矩阵的置换，

(在MATLAB中为C=A(p,q))。其中，长度为n的向量q为列转置向量，长度为m的

pinv(不是p)为行转置向量，其中A为m*n。如果 pinv[i]=k，则A的第i行变为C的第k行；同理，A的第j列通过q进行置换。

函数cs_symperm用于对称矩阵的置换。如果一个对称矩阵，仅保存了上三角矩阵或者下三角元素，则利用cs_cs_permute则会使得到C不是一个对称矩阵，所以需要cs_symperm实现此功能，并且返回值C是与A同样的三角矩阵形式。

函数	用途
cs *cs_permute (const cs *A, const int *pinv, const int *q, int values)	计算一般矩阵的置换矩阵
cs *cs_symperm (const cs *A, const int *pinv, int values)	计算对称矩阵的置换矩阵

2.12 Matrix norm

矩阵范数：

• 1-范数：
  ，是最容易计算的，利用
  cs_norm函数得到。
• 2-范数：A的最大奇异值，稀疏矩阵的2-范数并不重要，此处没有给计算函数。
• ∞-范数：最大的行的和。

函数	用途
double cs_norm (const cs *A)	计算稀疏矩阵的1-范数