稀疏矩阵算法
1. 稀疏矩阵的运算
import java.util.*;
//稀疏矩阵算法。
//稀疏态昌矩阵算法是为了在大型矩阵中非零元素少时，减少存贮空间，并提高矩阵运算速度的。
//但本例中的矩阵只是为了演示算法，都比较小，时间告备和空间效率提升可以忽袜闭毁略。
public class SparseMatrix{
public static void main(String[] args){
TripleSMatrix tsm=new TripleSMatrix(7,4);
//tsm.printTriple();
tsm.printMatrix();
TripleSMatrix tsm2=new TripleSMatrix(7,4);
System.out.println("矩阵a:");
tsm.printMatrix();
System.out.println("矩阵b:");
tsm2.printMatrix();
int[][] matrixSum=addSMatrix(tsm,tsm2);
System.out.println("矩阵a+矩阵b:");
for(int i=0;i <matrixSum.length;i++){
for(int j=0;j <matrixSum[i].length;j++){
System.out.print(" "+matrixSum[i][j]);
}
System.out.println("");
}
}
public static int[][] addSMatrix(TripleSMatrix t1,TripleSMatrix t2){ //计算两个三元组表示的矩阵之和，返回结束数组
if(t1.rows!=t2.rows||t1.columns!=t2.columns){
System.out.println("这两个矩阵不能相加");
return null;
}
int[][] c=new int[t1.rows][t2.columns];
int i=1,j=1;
while(i <=t1.nonzeroElements||j <=t2.nonzeroElements){
if(t1.triple[i][0] <t2.triple[j][0]&&i <=t1.nonzeroElements){
c[t1.triple[i][0]-1][t1.triple[i][1]-1]=t1.triple[i][2];
i++;
}else if(t2.triple[j][0] <t1.triple[i][0]&&j <=t2.nonzeroElements){
c[t2.triple[j][0]-1][t2.triple[j][1]-1]=t2.triple[j][2];
j++;
}else{
if(t1.triple[i][1] <t2.triple[j][1]&&i <=t1.nonzeroElements){
c[t1.triple[i][0]-1][t1.triple[i][1]-1]=t1.triple[i][2];
i++;
}else if(t1.triple[i][1]>t2.triple[j][1]&&j <=t2.nonzeroElements){
c[t2.triple[j][0]-1][t2.triple[j][1]-1]=t2.triple[j][2];
j++;
}else{
c[t1.triple[i][0]-1][t1.triple[i][1]-1]=t1.triple[i][2]+t2.triple[j][2];
i++;j++;
}
}
}
return c;
}
}
//下面的类定义不一定是最好的，比如其中的属性大多是包访问权限，可以改进。
class TripleSMatrix{ //定义了一个三元组的类。
int[][] triple=new int[2001][3]; //三元组数组，假设稀疏矩阵的值都是整数。最多可以有2000个非零元素。第零行没有用。
int rows,columns,nonzeroElements; //稀疏矩阵的行列数和非零元素个数。
TripleSMatrix(int rows,int columns){ //构造方法，rows是稀疏矩阵的行数，columns是稀疏矩阵的列数。
Scanner input=new Scanner(System.in);
System.out.println("请输入稀疏矩阵三元组");
System.out.println("以行 列 值的形式输入，如：1 2 4表示第1行第2列元素的值为4，当输入的行为999时结束：");
int count=1;
int i=0,j,v; //i行j列,值v
while(i!=999&&input.hasNext()){
i=input.nextInt();
j=input.nextInt();
v=input.nextInt();
if(i>rows||i <1||j>columns||j <1){
System.out.println("刚才的行，列值错，将被忽略");
continue;
}
triple[count][0]=i;
triple[count][1]=j;
triple[count][2]=v;
count++;
}
this.rows=rows;
this.columns=columns;
this.nonzeroElements=count-1;
sortTriple(triple,1,count); //对输入的三元组排序。
}
static void sortTriple(int[][] triple,int first,int end){ //对三元组排序方法，按行排，行一样按列排。
Arrays.sort(triple,first,end,new Comparator <int[]>(){
public int compare(int[] t1,int[] t2){
if(t1[0]>t2[0]) return 1;
if(t1[0] <t2[0]) return -1;
if(t1[0]==t2[0]) return t1[1]-t2[1];
return 0; //没有用的一个语句，但没有它编译通不过。
}
});
}
public void printMatrix(){ //打印出当前三元组表示的稀疏矩阵。
int row=1,column=1; //row当前要打印的行，column当前要打印的列。
for(int t=1;t <=nonzeroElements;t++){
while(triple[t][0]>row){ //三元组中的行比当前行大
if(column!=1){ //前面打印的行没有打印完，继续打印完
for(;column <=columns;column++) System.out.print(" "+0);
column=1; //新的一行列从1开始。
}else{ //当前行全为0
for(int i=1;i <=columns;i++){
System.out.print(" "+0);
}
}
System.out.println(""); //换行
row++; //下一行
}
for(;column <triple[t][1];column++){ //当前打印的列小于三元组中的列，前面要补零。
System.out.print(" "+0);
}
System.out.print(" ".substring(0,6-(String.valueOf(triple[t][2])).length())+triple[t][2]); //打印三元组对应的元素。
column++;
}
if(column!=1){ //前面打印的行没有打印完，继续打印完
for(;column <=columns;column++) System.out.print(" "+0);
System.out.println("");
column=1;
row++ ;
}
for(;row <=rows;row++){ //三元组中没有对应的值了，矩阵后面的元素全为0
for(column=1;column <=columns;column++){
System.out.print(" "+0);
}
System.out.println("");
}
}
public void printTriple(){ //打印三元组
for(int i=1;i <=nonzeroElements;i++){
for(int j=0;j <3;j++){
System.out.print(triple[i][j]+" ");
}
System.out.println("");
}
}
}
为了节省存储空间并且加快处理速度，需要对这类矩阵进行压缩存储，压缩存储的原则是：不重复存储相同元素；不存储零值元素。稀疏矩阵，有三元组表示法、带辅助行向量的二元组表示法（也即行逻辑链表的顺序表），十字链表表示法等。算法基本思想：num[col]:第col列的非零元素个数；cpot[col]:第col列第一个非零元在b.data中的恰当位置；在转置过程中，指示该列下一个非零元在b.data中的位置。
3. 稀疏矩阵的运算
M AT L A B中对满矩阵的运算和函数同样可用在稀疏矩阵中.结果是稀疏矩阵还是满矩阵,
这取决于运算符或者函数及下列的操作数:当函数用一个矩阵作为输入参数,输出参数为一个标量或者一个给定大小的向量时,输出参数的格式总是返回一个满阵形式,如命令s i z e.
当函数用一个标量或者一个向量作为输入参数,输出参数为一个矩阵时,输出参数的格式也总是返回一个满矩阵,如命令e y e.还有一些特殊的命令可以得到稀疏矩阵,如命令s p e y e.
对于单参数的其他函数来说,通常返回的结果和参数的形式是一样的,如d i a g.
对于双参数的运算或者函数来说,如果两个参数的形式一样,那么也返回同样形式的结果.在两个参数形式不一样的情况下,除非运算的需要,均以满矩阵的形式给出结果.
两个矩阵的组和[A B],如果A或B中至少有一个是满矩阵,则得到的结果就是满矩阵.
表达式右边的冒号是要求一个参数的运算符,遵守这些运算规则.
表达式左边的冒号不改变矩阵的形式. 假设有:
这是一个5×5的单位满矩阵和相应的稀疏矩阵.
(a) C = 5*B,结果为:
这是一个稀疏矩阵.
(b) D = A + B,给出的结果为:
这是一个满矩阵.
(c) x = B h,结果为:
这是一个满向量.
有许多命令可以对非零元素进行操作.
命令集8 9矩阵的非零元素
n n z ( A )求矩阵A中非零元素的个数.它既可求满矩阵也可求稀疏矩阵.
s p y ( A )画出稀疏矩阵A中非零元素的分布.也可用在满矩阵中,在
这种情况下,只给出非零元素的分布.
s p y ( A , c s t r , s i z e )用指定的颜色c s t r(见表1 3 - 1 )和在s i z e规定的范围内画出稀疏
矩阵A中非零元素的分布.
n o n z e r o s ( A )按照列的顺序找出矩阵A中非零的元素.
s p o n e s ( A )把矩阵A中的非零元素全换为1.
s p a l l o c ( m , n ,产生一个m×n阶只有n z m a x个非零元素的稀疏矩阵.这样可以
n z m a x )有效地减少存储空间和提高运算速度.
n z m a x ( A )给出为矩阵A中非零元素分配的内存数.不一定和n n z ( A )得
到的数相同;参见s p a r s e或者s p a l l o c.
i s s p a r s e ( A )如果矩阵A是稀疏矩阵,则返回1;否则返回0.
s p f u n ( f c n , A )用A中所有非零元素对函数f c n求值,如果函数不是对稀疏矩
阵定义的,同样也可以求值.
s p r a n k( A )求稀疏矩阵A的结构秩.对于所有的矩阵来说,都有
s p r a n k ( A）≥r a n k ( A ). 用下面的命令定义稀疏矩阵:
创建一个大矩阵:
Big=kron(A, A)
这个矩阵B i g是什么样子呢?
K r o n e c k e r张量积给出一个大矩阵,它的元素是矩阵A的元素之间可能的乘积.因为参量都是稀疏矩阵,所以得到的矩阵也是一个稀疏矩阵.可以用命令 w h o s和i s s p a r s e来确认一下.
查看矩阵B i g的结构图,可输入s p y ( B i g ),结构图如右图所示. 从图中可以看出B i g是一个块双对角矩阵. MATLAB中有四个基本稀疏矩阵,它们是单位矩阵,随机矩阵,对称随机矩阵和对角矩阵.
命令集9 0单位稀疏矩阵
s p e y e ( n )生成n×n的单位稀疏矩阵.
s p e y e ( m , n )生成m×n的单位稀疏矩阵.
命令speye(A) 得到的结果和s p a r s e ( e y e ( A ) )是一样的,但是没有涉及到满阵的存储.
命令集9 1随机稀疏矩阵
s p r a n d ( A )生成与A有相同结构的随机稀疏矩阵,且元素服从均匀分布.
s p r a n d ( m , n , d e n s )生成一个m×n的服从均匀分布的随机稀疏矩阵,有d e n s×m×
n个非零元素,0≤d e n s≤1.参数d e n s是非零元素的分布密度.
s p r a n d ( m , n , d e n s ,生成一个近似的条件数为1 /rc,大小为m×n的随机稀疏矩
r c )阵.如果rc=rc是一个长度为l≤l ( m i n (m, n) )的向量,那么
矩阵将rci作为它l个奇异值的第一个,其他的奇异值为0.
s p r a n d n ( A )生成与A有相同结构的随机稀疏矩阵,且元素服从正态分布.
s p r a n d n ( m , n , d e n s ,生成一个m×n的服从正态分布的随机稀疏矩阵,和sprand
r c )一样.
s p r a n d s y m ( S )生成一个随机对称稀疏矩阵.它的下三角及主对角线部分与S的结构相同,矩阵元素服从正态分布.
s p r a n d s y m ( n , d e n s )生成一个m×n的随机对称稀疏矩阵.矩阵元素服从正态分布,分布密度为d e n s.
s p r a n d s y m ( n , d e n s ,生成一个近似条件数为1 /rc的随机对称稀疏矩阵.元素以0r c )对称分布,但不是正态分布.如果rc=rc是一个向量,则矩阵有特征值rci.也就是说,如果rc是一个正向量,则矩阵是正定矩阵.
s p r a n d s y m ( n , d e n s ,生成一个正定矩阵.如果k= 1,则矩阵是由一正定对称矩r c , k )阵经随机J a c o b i旋转得到的,其条件数正好等于1 /rc;如果k= 2,则矩阵为外积的换位和,其条件数近似等于1 /rc.
s p r a n d s y m ( S , d e n s ,生成一个与矩阵S结构相同的稀疏矩阵,近似条件数为1 /rc.
r c , 3 )参数d e n s被忽略,但是这个参数在这个位置以便函数能确认最后两个参数的正确与否. (a) 假设有矩阵A:
输入R a n d o m = s p r a n d n ( A ),可得到随机稀疏矩阵:
矩阵中随机数的位置和矩阵A中非零元素的位置相同.
(b) 对于( a )中的矩阵A,输入:
B = s p r a n d s y m ( A )
结果为:
这是一个用矩阵A的下三角及主对角线部分创建的对称矩阵,在非零元素的位置用随机数作为元素值.
用命令s p d i a g s可以取出对角线元素,并创建带状对角矩阵.假设矩阵A的大小为m×n,
在p个对角线上有非零元素.B的大小为m i n (m×n)×p,它的列是矩阵A的对角线.向量d的长度为p,其整型分量给定了A的对角元:
di0 主对角线上的对角线
命令集9 2对角稀疏矩阵
[ B , d ] = s p d i a g s ( A )求出A中所有的对角元,对角元保存在矩阵B中,它们的下标保存在向量d中.
s p d i a g s ( A , d )生成一个矩阵,这个矩阵包含有矩阵A中向量d规定的对角元.
s p d i a g s ( B , d , A )生成矩阵A,用矩阵B中的列替换d定义的对角元.
A = s p d i a g s ( B , d , m , n )用保存在由d定义的B中的对角元创建稀疏矩阵A.
例11 . 4给出了如何使用s p d i a g s命令来解普通微分方程组. 在许多实际应用中要保留稀疏矩阵的结构,但是在计算过程中的中间结果会减弱它的稀疏性,如L U分解.这就会导致增加浮点运算次数和存储空间.为了避免这种情况发生,在第稀疏矩阵1 2 9
M AT L A B中用命令对矩阵进行重新安排.这些命令都列在下面的命令集9 3中.通过h e l p命令
可以得到每个命令更多的帮助信息,也可见h e l p d e s k.
命令集9 3矩阵变换
c o l m m d ( A )返回一个变换向量,使得矩阵A列的秩为最小.
s y m m m d ( A )返回使对称矩阵秩为最小的变换.
s y m r c m ( A )矩阵A的C u t h i l l - M c K e e逆变换.矩阵A的非零元素在主对角线附近.
c o l p e r m ( A )返回一个矩阵A的列变换的向量.列按非零元素升序排列.有时这是L U因式分解前有用的变换:lu(A(:, j)).如果A是一个对称矩阵,对行和列进行排序,这有利于C h o l e s k y分解:chol(A(j, j)).
r a n d p e r m ( n )给出正数1,2,. . .,n的随机排列,可以用来创建随机变换矩阵.
d m p e r m ( A )对矩阵A进行D u l m a g e - M e n d e l s o h n分解,输入help dmperm可得更多信息. 创建一个秩为4的变换矩阵,可输入:
一旦运行p e r m = r a n d p e r m ( 4 ),就会得到:
给出的变换矩阵为:
如果矩阵A为:
输入命令:
运行结果为:
有两个不完全因式分解命令,它们是用来在解大线性方程组前进行预处理的.用h e l p d e s k命令可得更多信息.命令集9 4不完全因式分解c h o l i n c ( A , o p t )进行不完全C h o l e s k y分解,变量o p t取下列值之一:
d r o p t o l指定不完全分解的舍入误差,0给出完全分解.
m i c h o l如果m i c h o l = 1,则从对角线上抽取出被去掉的元素.
r d i a g用s q r t ( d r o p t o l*n o r m ( X ( : , j ) ) )代替上三角分
解因子中的零元素,j为零元素所在的列.
[ L , U , P ]=返回矩阵X的不完全分解得到的三个矩阵L,U和P,变量o p t取
l u i n c ( X , o p t )下列值之一:
d r o p t o l指定分解的舍入误差.
m i l u改变分解以便从上三角角分解因子中抽取被去掉的列元素.
u d i a g用d r o p t o l值代替上三角角分解因子中的对角线上的零元素.
t h r e s h中心极限.
解稀疏线性方程组既可用左除运算符解,也可用一些特殊命令来解.
命令集9 5稀疏矩阵和线性方程组
s p p a r m s ( k e y s t r , o p )设置稀疏矩阵算法的参数,用help spparms可得详细信息.
s p a u g m e n t ( A , c )根据[ c*l A; A' 0 ]创建稀疏矩阵,这是二次线性方程组的最
小二乘问题.参见7 . 7节.
s y m b f a c t ( A )给出稀疏矩阵的C h o l e s k y和L U因式分解的符号分解因子.
用help symbfact可得详细信息.
稀疏矩阵的范数计算和普通满矩阵的范数计算有一个重要的区别.稀疏矩阵的欧几里德范数不能直接求得.如果稀疏矩阵是一个小矩阵,则用n o r m ( f u l l ( A ) )来计算它的范数;但是对于大矩阵来说,这样计算是不可能的.然而M AT L A B可以计算出欧几里德范数的近似值,在计算条件数时也是一样.
命令集9 6稀疏矩阵的近似欧几里德范数和条件数
n o r m e s t ( A )计算A的近似欧几里德范数,相对误差为1 0-6.
n o r m e s t ( A , t o l )计算A的近似欧几里德范数,设置相对误差t o l,而不用缺省时的1 0-6.
[ n r m , n i t ] =计算近似n r m范数,还给出计算范数迭代的次数n i t.
n o r m e s t ( A )
c o n d e s t ( A )求矩阵A条件数的1 -范数中的下界估计值.
[ c , v ]=求矩阵A的1 -范数中条件数的下界估计值c和向量v,使得
c o n d e s t ( A , t r )| |Av| | = ( | |A| | | |v| | ) / c.如果给定t r,则给出计算的过程.t r= 1,
给出每步过程;t r=-1,给出商c / r c o n d ( A ). 假设给出:
用n o r m A p p r o x = n o r m e s t ( S p r s )计算出:
用t h e N o r m = n o r m ( f u l l ( S p r s ) )得:
为了找到它们之间的差别,计算d i f f e r e n c e = t h e N o r m - n o r m A p p r o x,得:
在许多应用中,n o r m e s t计算得到的近似值是一个很好的近似欧几里德范数,它的计算步数要比n o r m要少得多;可参见7 . 6节.
用e t r e e命令来找到稀疏对称矩阵的消元树,用向量f来描述消元树,还可用e t r e e p l o t命令画出来.元素fi是矩阵的上三角C h o l e s k y分解因子中i行上第1非零元素的列下标.如果有非零元素,则fi= 0.消元树可以这样来建立:
节点i是fi的孩子,或者如果fi= 0,则节点i是树的根节点.
命令集9 7矩阵的消元树
e t r e e ( A )求A的消元树向量f,这个命令有可选参数;输入h e l p
e t r e e获取帮助.
e t r e e p l o t ( A )画出向量f定义的消元树图形.
t r e e p l o t ( p , c , d )画出指针向量p的树图形,参数c和d分别指定节点的颜色和分支数.e t r e e p l o t可以调用这个命令.
t r e e l a y o u t显示树的结构,t r e e p l o t可以调用这个命令. 假设有对称稀疏矩阵B:
运行命令b t r e e = e t r e e ( B ),结果为:
开始的数字2不难理解,它是矩阵的第1列上第1个非零元素的行数,它决定了在C h o l e s k y分解因子的第1行第2列处有一个非零元素.当缩减第1列的元素时就得到第2列的数字5.B在缩减后,在( 5 , 2 )位置的元素是非零的,这样消元树向量中第2个元素的值为5.
s p y ( c h o l ( B ) )给出了C h o l e s k y分解因子的结构图,如图9 - 2所示:
图9-2 Cholesky分解结构图
图9-3 矩阵B的消元树
这个向量消元树可以这样来建立:上三角中只有一行有非零元素,节点8,因此这就是树
唯一的根.节点1是节点2的孩子,节点2和3是节点5的孩子,而节点5是节点6的孩子.节点4和6是节点7的孩子,而节点7又是节点8的孩子,即根的孩子.
命令e t r e e p l o t ( B )给出了树的结构图,如图9 - 3所示.
消元树的形状取决于列和行序,它可以用来分析消元过程.
用g p l o t命令可以画出坐标和矩阵元素间的联系图形.必须在n×2的矩阵中给出n个坐标,
矩阵的每一行作为一个点.这样就创建出点点之间连接的n×n矩阵,如果点4连接到点8,则(4, 8)的值为1.由于是一个大矩阵,而且非零元素较少,所以它应该被建成稀疏矩阵.
这个图可以说明网络问题,如传递问题.它还包含有线性方程组中未知量之间的相关信息.
命令集9 8网络图形
g p l o t ( A , K )如果矩阵A的a(i, j)不为0,则将点ki连接到点kj.K是一个n×
2的坐标矩阵,A是一个n×n的关联矩阵.
g p l o t ( A , K , s t r )用字符串s t r给定的颜色和线型画出的同上图形.字符串s t r的
取值参见表1 3 - 1.
[ X , A ] = u n m e s h ( E )求网格边界矩阵E的L a p l a c e矩阵A和网格点的坐标矩阵X.
例七
假设有下面的坐标矩阵K和关联矩阵A:
矩阵A在稀疏化后,用命令g p l o t ( A , K )画出图9 - 4,给出了点(0, 1)和点(4, 1)之间所有可能的路径.
4. c++ 十字链表实现稀疏矩阵的转置
mplate <class Type>SparseMatrix<Type>
SparseMatrix<Type>::SparseMatrix::FastTranspose()
{
int *RowSize = new int[Cols]; // 统计各列非零元素个数
int *RowStart = new int[Cols]; // 预计转置后各行存放位置
SparseMatrix b; // 存放转置结果
b.Rows = cols; b.cols = Rows; b.Terms = Terms;
if (Terms>0) //nonzero matrix
{
for (int i=0;i<Cols;i++) RowSize[i] = 0; //initialize
for (i=0;i<terms;i++) RowSize[smArray[i].col]++;
//RowStart[i]=starting position of row i in b
RowStart[0]=0;
for (i=1;i<Cols;i++) RowStart[i] = RowStart[i-1]+
RowSize[i-1];
for (i=0;i<Terms;i++) //move from a to b
{
int j=Rowstart[smArray[i].col];
b.smArray[j].row=smArray[i].col;
b.smArray[j].col= smArray[i].row;
b.smArray[j].value= smArray[i].value;
RowStart[smArray[i].col]++;
} //end of for
}
delete [ ]Rowsize;
delete [ ]Rowstart;
return b;
}
另外,站长团上有产品团购,便宜有保证
5. 稀疏矩阵的乘法的算法思想
/*Multiplicate part*/
//C = A * B
/*算法分析：
首先，由于楼主没有给出输入函数，也没有对三元组的稀疏矩阵的数据结构做完整的说明，
所以我只能猜测这个稀疏矩阵是以行为主序存储的。（后面的乘法函数佐证了我的猜测，
但是如果我不幸猜错了，还请楼主告知）
另一个猜测是楼主的程序中设定矩阵和数组的下标都是从1开始，而非从0开始。
接下来，我们说一下普通矩阵的乘法，这个在线性代数里面有定义，无需赘言。
我要说的是稀疏矩阵的乘法也是用这个公式来计算，但却有一个问题——效率。
我们举一个例子来说明：
假设我们已知A:m*n和B:n*l,要计算C:m*l，那么C(i,j)的计算公式就是：
C(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + …… + A(i,n)*B(n,j) 公式1
如果A、B都是普通矩阵（并且以二维数组存储），那么直接用循环语句就可以完成。
如果A、B都是稀疏矩阵（并且乱序存储，即是说没有以行序或列序存储），那么计算就会很麻烦。
我们需要首先遍历整个A矩阵去查找是否存在A(i,1)（若有则取其值，若无则其值为0），然后再去
遍历B矩阵查找B(1,j)，并将两者相乘；接着又是A(i,2)和B(2,j)，以此类推。
这样做当然是可行的，但是显然效率太低了，一种改进的方法（就是楼主的程序中所用的方法）如下：
首先我们要优化稀疏矩阵的存储，不能乱序存储，而是以行序或列序为主序来存储，
比如这里我们以行序为主序，以列序为次序。
具体来说，就是将稀疏矩阵的第一行中的非零元素排列在前面，然后才是第二行、第三行……
在各行的非零元素中又以列序来排列，这样存储的稀疏矩阵就是有序的。
接下来，在真正用公式来计算之前，我们要做一个预处理。
这个预处理是对矩阵B所做的，其实就是要计算得到矩阵B所对应的pos数组。
那么这个pos数组代表什么意思呢？
我们从代码中不难看出pos数组的长度是矩阵B的行数（也就是B->m）加1。
pos[i]（1 <= i<= B->m）表示矩阵B的第i行的第一个非零元素在三元数组B的位置。
pos[i]（i == B->m + 1）是为了方便后面的计算而增加一个标记，类似于监视哨。
有了这个pos数组之后，我们再来计算公式1就可以提高效率了。
效率的提高表现在两个方面：
1、我们看到代码中（如下）虽然有两层while循环，但是其实总共只执行了A->t次。
p=1;
while(p<=A->t)
{
……
while (p<=A->t&&A->data[p].row==crow)
{
……
p=p+1;
}
……
}
2、在第二层while循环中的for循环，由于pos数组的作用使得循环次数减少很多。
for(q=pos[brow];q<=pos[brow+1]-1;q++)
{
ccol=B->data[q].col;
ctemp[ccol]=ctemp[ccol] + A->data[p].v * B->data[q].v;
}
好了，以上就是这个稀疏矩阵乘法的大致过程，如有错误请指正。
*/
void MultS(TriTable *A,TriTable *B,TriTable *C)
{
int k,p,q,brow,crow,ccol;
int num[MAXSIZE],pos[MAXSIZE],ctemp[MAXSIZE];
if (A->n==B->m)
{
for(k=1;k<=B->m;k++)
num[k]=0;
for(k=1;k<=B->t;k++)
num[B->data[k].row]++;
pos[1]=1;
for(k=2;k<=B->m;k++) //这里将k<=B->t改为k<=B->m
pos[k]=pos[k-1]+num[k-1];
pos[1+B->m]=pos[B->m]+1; //这里将k<=B->t改为k<=B->m
C->m=A->m;
C->n=B->n;
C->t=0;
p=1;
while(p<=A->t)
{
crow=A->data[p].row;
for(k=1;k<=C->n;k++)
ctemp[k]=0;
while (p<=A->t&&A->data[p].row==crow)
{
brow=A->data[p].col;
for(q=pos[brow];q<=pos[brow+1]-1;q++)
{
ccol=B->data[q].col;
ctemp[ccol]=ctemp[ccol] + A->data[p].v * B->data[q].v;
}
p=p+1;
}
for(ccol=1;ccol<=B->n;ccol++)
if(ctemp[ccol]!=0)
{
C->t=C->t+1;
C->data[C->t].row=crow;
C->data[C->t].col=ccol;
C->data[C->t].v=ctemp[ccol];
}
}
}else
printf("these two arrat can't Multiplicate");
return;
}
6. 对稀疏矩阵进行压缩存储的目的是什么
对稀疏矩阵进行压缩存储目的是节省存储空间。
存储矩阵的一般方法是采用二维数组，其优点是可以随机地访问每一个元素，因而能够较容易地实现矩阵的各种运算。
但对于稀疏矩阵而言，若用二维数组来表示，会重复存储了很多个0了，浪费空间，而且要花费时间来进行零元素的无效计算。所以必须考虑对稀疏矩阵进行压缩存储。
(6)稀疏矩阵算法扩展阅读
优点
稀疏矩阵的计算速度更快，因为MATLAB只对非零元素进行操作，这是稀疏矩阵的一个突出的优点。假设矩阵A,B中的矩阵一样，计算2*A需要一百万次的浮点运算，而计算2*B只需要2000次浮点运算。
因为MATLAB不能自动创建稀疏矩阵，所以要用特殊的命令来得到稀疏矩阵。算术和逻辑运算都适用于稀疏矩阵。对于一个用二维数组存储的稀疏矩阵Amn，如果假设存储每个数组元素需要L个字节，那么存储整个矩阵需要m*n*L个字节。
7. 稀疏矩阵
// algo5-1.cpp 实现算法5.2的程序
// c1.h (程序名)
#include<string.h>
#include<ctype.h>数念
#include<malloc.h> // malloc()等
#include<limits.h> // INT_MAX等
#include<stdio.h> // EOF(=^Z或F6),NULL
#include<stdlib.h> // atoi()
#include<io.h> // eof()
#include<math.h> // floor(),ceil(),abs()
#include<process.h> // exit()
#include<iostream.h> // cout,cin
// 函数结果状态代码
#define TRUE 1
#define FALSE 0
#define OK 1
#define ERROR 0
#define INFEASIBLE -1
// #define OVERFLOW -2 因为在math.h中已定义OVERFLOW的值为3,故去掉此行
typedef int Status; // Status是函数的类型,其值是函数结果状态代码，如OK等
typedef int Boolean; // Boolean是布尔类型,其值是TRUE或FALSE
typedef int ElemType;
// c5-2.h 稀疏矩阵的三元组顺序表存储表示
#define MAXSIZE 100 // 非零元个数的最大值
struct Triple
{
int i,j; // 行下标,列下标
ElemType e; // 非零元素值
};
struct TSMatrix
{
Triple data[MAXSIZE+1]; // 非零元三元组表,data[0]未用
int mu,nu,tu; // 矩阵的行数、列数和非零元个数核哪
};
// bo5-2.cpp 三元组稀疏矩阵的基本操作,包括算法5.1(9个)
Status CreateSMatrix(TSMatrix &M)
{ // 创建稀疏矩阵M
int i,m,n;
ElemType e;
Status k;
printf("请输入矩阵的行数,列数,非零元素数：");
scanf("%d,%d,%d",&M.mu,&M.nu,&M.tu);
M.data[0].i=0; // 为以下比较顺序做准备
for(i=1;i<=M.tu;i++)
{
do
{
printf("请按行序顺序输入第%d个非零元素所在的行(1～%d),列(1～%d),元素值：",i,M.mu,M.nu);
scanf("%d,%d,%d"改毕码,&m,&n,&e);
k=0;
if(m<1||m>M.mu||n<1||n>M.nu) // 行或列超出范围
k=1;
if(m<M.data[i-1].i||m==M.data[i-1].i&&n<=M.data[i-1].j) // 行或列的顺序有错
k=1;
}while(k);
M.data[i].i=m;
M.data[i].j=n;
M.data[i].e=e;
}
return OK;
}
void DestroySMatrix(TSMatrix &M)
{ // 销毁稀疏矩阵M
M.mu=0;
M.nu=0;
M.tu=0;
}
void PrintSMatrix(TSMatrix M)
{ // 输出稀疏矩阵M
int i;
printf("%d行%d列%d个非零元素。\n",M.mu,M.nu,M.tu);
printf("行 列 元素值\n");
for(i=1;i<=M.tu;i++)
printf("%2d%4d%8d\n",M.data[i].i,M.data[i].j,M.data[i].e);
}
Status CopySMatrix(TSMatrix M,TSMatrix &T)
{ // 由稀疏矩阵M复制得到T
T=M;
return OK;
}
int comp(int c1,int c2) // 另加
{ // AddSMatrix函数要用到
int i;
if(c1<c2)
i=1;
else if(c1==c2)
i=0;
else
i=-1;
return i;
}
Status AddSMatrix(TSMatrix M,TSMatrix N,TSMatrix &Q)
{ // 求稀疏矩阵的和Q=M+N
Triple *Mp,*Me,*Np,*Ne,*Qh,*Qe;
if(M.mu!=N.mu)
return ERROR;
if(M.nu!=N.nu)
return ERROR;
Q.mu=M.mu;
Q.nu=M.nu;
Mp=&M.data[1]; // Mp的初值指向矩阵M的非零元素首地址
Np=&N.data[1]; // Np的初值指向矩阵N的非零元素首地址
Me=&M.data[M.tu]; // Me指向矩阵M的非零元素尾地址
Ne=&N.data[N.tu]; // Ne指向矩阵N的非零元素尾地址
Qh=Qe=Q.data; // Qh、Qe的初值指向矩阵Q的非零元素首地址的前一地址
while(Mp<=Me&&Np<=Ne)
{
Qe++;
switch(comp(Mp->i,Np->i))
{
case 1: *Qe=*Mp;
Mp++;
break;
case 0: switch(comp(Mp->j,Np->j)) // M、N矩阵当前非零元素的行相等,继续比较列
{
case 1: *Qe=*Mp;
Mp++;
break;
case 0: *Qe=*Mp;
Qe->e+=Np->e;
if(!Qe->e) // 元素值为0，不存入压缩矩阵
Qe--;
Mp++;
Np++;
break;
case -1: *Qe=*Np;
Np++;
}
break;
case -1: *Qe=*Np;
Np++;
}
}
if(Mp>Me) // 矩阵M的元素全部处理完毕
while(Np<=Ne)
{
Qe++;
*Qe=*Np;
Np++;
}
if(Np>Ne) // 矩阵N的元素全部处理完毕
while(Mp<=Me)
{
Qe++;
*Qe=*Mp;
Mp++;
}
Q.tu=Qe-Qh; // 矩阵Q的非零元素个数
return OK;
}
Status SubtSMatrix(TSMatrix M,TSMatrix N,TSMatrix &Q)
{ // 求稀疏矩阵的差Q=M-N
int i;
for(i=1;i<=N.tu;i++)
N.data[i].e*=-1;
AddSMatrix(M,N,Q);
return OK;
}
Status MultSMatrix(TSMatrix M,TSMatrix N,TSMatrix &Q)
{ // 求稀疏矩阵的乘积Q=M*N
int i,j,h=M.mu,l=N.nu,Qn=0;
// h,l分别为矩阵Q的行、列值,Qn为矩阵Q的非零元素个数，初值为0
ElemType *Qe;
if(M.nu!=N.mu)
return ERROR;
Q.mu=M.mu;
Q.nu=N.nu;
Qe=(ElemType *)malloc(h*l*sizeof(ElemType)); // Qe为矩阵Q的临时数组
// 矩阵Q的第i行j列的元素值存于*(Qe+(i-1)*l+j-1)中，初值为0
for(i=0;i<h*l;i++)
*(Qe+i)=0; // 赋初值0
for(i=1;i<=M.tu;i++) // 矩阵元素相乘，结果累加到Qe
for(j=1;j<=N.tu;j++)
if(M.data[i].j==N.data[j].i)
*(Qe+(M.data[i].i-1)*l+N.data[j].j-1)+=M.data[i].e*N.data[j].e;
for(i=1;i<=M.mu;i++)
for(j=1;j<=N.nu;j++)
if(*(Qe+(i-1)*l+j-1)!=0)
{
Qn++;
Q.data[Qn].e=*(Qe+(i-1)*l+j-1);
Q.data[Qn].i=i;
Q.data[Qn].j=j;
}
free(Qe);
Q.tu=Qn;
return OK;
}
Status TransposeSMatrix(TSMatrix M,TSMatrix &T)
{ // 求稀疏矩阵M的转置矩阵T。算法5.1
int p,q,col;
T.mu=M.nu;
T.nu=M.mu;
T.tu=M.tu;
if(T.tu)
{
q=1;
for(col=1;col<=M.nu;++col)
for(p=1;p<=M.tu;++p)
if(M.data[p].j==col)
{
T.data[q].i=M.data[p].j;
T.data[q].j=M.data[p].i;
T.data[q].e=M.data[p].e;
++q;
}
}
return OK;
}
Status FastTransposeSMatrix(TSMatrix M,TSMatrix &T)
{ // 快速求稀疏矩阵M的转置矩阵T。算法5.2
int p,q,t,col,*num,*cpot;
num=(int *)malloc((M.nu+1)*sizeof(int)); // 生成数组（[0]不用）
cpot=(int *)malloc((M.nu+1)*sizeof(int)); // 生成数组（[0]不用）
T.mu=M.nu;
T.nu=M.mu;
T.tu=M.tu;
if(T.tu)
{
for(col=1;col<=M.nu;++col)
num[col]=0; // 设初值
for(t=1;t<=M.tu;++t) // 求M中每一列含非零元素个数
++num[M.data[t].j];
cpot[1]=1;
for(col=2;col<=M.nu;++col) // 求第col列中第一个非零元在T.data中的序号
cpot[col]=cpot[col-1]+num[col-1];
for(p=1;p<=M.tu;++p)
{
col=M.data[p].j;
q=cpot[col];
T.data[q].i=M.data[p].j;
T.data[q].j=M.data[p].i;
T.data[q].e=M.data[p].e;
++cpot[col];
}
}
free(num);
free(cpot);
return OK;
}
void main()
{
TSMatrix A,B;
printf("创建矩阵A: ");
CreateSMatrix(A);
PrintSMatrix(A);
FastTransposeSMatrix(A,B);
printf("矩阵B(A的快速转置): ");
PrintSMatrix(B);
DestroySMatrix(A);
DestroySMatrix(B);
}
转置的，你参考下
8. 稀疏矩阵的压缩存储只需要存储什么
非零元素。
对于一个用二维数组存储的稀疏矩阵Amn，如果假设存储每个数组元素需要L个字节，那么存储整个矩阵需要m*n*L个字节。但是，这些存储空间的大部分存放的是0元素，从而造成大量的空间浪费。为了节省存储空间，可以只存储其中的非0元素。
(8)稀疏矩阵算法扩展阅读
稀疏矩阵算法的最大特点是通过只存储和处理非零元素从而大幅度降低存储空间需求以及计算复杂度，代价则是必须使用专门的稀疏矩阵压缩存储数据结构。稀疏矩阵算法是典型的不规则算法，计算访存比很低，并且计算过程中的访存轨迹与稀疏矩阵的稀疏结构相关。