matlab稀疏矩阵

(zZ)

先说一下matlab如何存储稀疏矩阵。matlab使用三个数组表示一个稀疏矩阵，Ir，Jc，Pr，其中Ir的长度和稀疏矩阵非零元素个数一样，依次记录了所有元素的行号。Pr也是依次记录所有非零元素的取值，和Ir对应。一般情况下，我们会想到建立一个同样大小记录列号的数组，然而matlab没有这么做。matlab使用了一个更小的数组Jc，这个数组长度等于稀疏矩阵的列数加1. Jc中依次记录了Ir中对应每一列结束的位置。

举个例子，如果一个5*5的稀疏矩阵，Ir为[1 1 2 4]，表示四个非零元素的行号分别为1,1,2,4.  如果Jc等于 [0 1 3 4 4 4]. 就说明 Ir中第一个元素在第一列，元素2-3对应第二列，元素4对应第三列。再比如Jc等于[0 1 4 4 4 4]的话，就说明第一个元素位置为(1,1)， 其他几个元素位置为(2,1),(2,2),(2,4)。

在mex文件中，有个数据类型，mwIndex，稀疏矩阵存储中的Ir，Jc都是这个类型的数组。在32位系统中，这个类型为4个字节（32位），而在64位系统中，为了能存储更大的稀疏矩阵，这个类型升级到了8字节（64位）。然而，在64位windows中，考虑到用户习惯，int类型仍然是32位的，long类型为64位。

在一些老代码中，很多人习惯用int代替mwIndex，因为他们曾经是等价的。然而这些代码用到64位windows时，这种做法就会存在错误。而且这种非语法错误在编译的时候时不提示的，运行时候就会出错。

因此，如果要用64位windows编译涉及稀疏矩阵的mex文件，除了加上-largeArrayDims以外，记得检查所有Ir和Jc的类型，确定是mwIndex。避免像我这样在莫名奇妙的错误中浪费好几个小时

=========================

例子：

a=[0 1 0 0;
   1 0 0 0;
   0 0 0 0
   0 0 0 1];
a=sparse(a);
aaa(a);

输出：

>> test
1	0	3	
0	1	2	2	3

　　　

其中aaa.cpp为：

#include "mex.h"
 
void mexFunction( int nlhs, mxArray *plhs[],
		int nrhs, const mxArray *prhs[] )
{
    mwIndex *ir, *jc;
    double *samples;
   	samples = mxGetPr(prhs[0]);
	ir = mxGetIr(prhs[0]);
	jc = mxGetJc(prhs[0]);
    for(int i=0;i<3;i++)
    {printf("%d\t",ir[i]);
    }printf("\n");
    for(int i=0;i<5;i++)
    {printf("%d\t",jc[i]);
    }
}

　　matlab是按列存储！