记得第一个跟我讲MATLAB的老师是数字信息处理的老师，他大概的意思是说MATLAB如何好，举个例子在MATLAB里做个转置就只要a'就行了，说用C语言还要循环套循环，麻烦~我当时就想，MATLAB好用还不是底层把事情都做好了，你只调用当然方便啦，我还是觉得C/C++才是最好的语言。 我想模仿MATLAB做事，有了我写的库，很多事会变得同样简单。 主要功能，用C++模板将已知类型数据组织成矩阵的形式，并提供基本的矩阵运算。 1，矩阵定义：matrix<类型名> <变量名>; 比如定义一个double双精度的矩阵：matrix<double> a; 2，矩阵的初始化：从数组初始化：double a[]={1,2,3,4,5,6};matrix<double> ma(a,2,3);构造出一个2 x 3的矩阵ma，其元素依次是：1 2 34 5 6 从单个元素初始化：matrix<double> a=0.0;构造出一个1 x 1的矩阵，其元素是0 用size(h,w)成员函数指定存储大小：matrix<double> a;a.size(2,3);构造出一个2 x 3的矩阵，元素值未知。可以接着用a=0.0;将所有元素赋成0或者其它值 用提供函数identify(T(),h,w)构造单位矩阵：matrix<double> a=identify(double(),2,3);第一个参数随便指定一个相同类型的值就行了，结果构造一个2 x 3的单位矩阵：1 0 00 1 0（和MATLAB里的eye相同） 3，矩阵的基本运算由于这里的矩阵元素类型不一定为double型，所以不一定是用来做数值运算的，或者你也可以用其它你定义的类型，但有时是需要有+,-,*运算定义以及整数0，1的构造函数构造出0元和单位元的构造函数。有时你也可以把它想成一种容器（它现在的底层实现还不是很好），一种矩阵的容器。 3.1 选择运算选择运算是最重要的，拿到一个矩阵要很容易选定它的元素或部分。a,operator[] 选行（列）运算矩阵的行和列是对等的，我这里定义的是行标为正数列标为负数，如下：a -1 -2 -3 -41  12     13        1如果要选第2行，a[2]如果要选第2列，a[-2]选择的结果是得到一个向量（同样也是一个矩阵） 那如果它本身就是一个向量呢？那它将选到第i个元素，但返回类型还是矩阵，即1 x 1的矩阵。比如 a[2][2]就是一个1 x 1的矩阵，值为1，这和C语言二维数组很像，但是代价不同，一般只用它来选行（列）。 b,select(h,w) 多功能选择和operator[]类似，如果h=0，则选第w列，如果w=0，则选第h行，如果都不为0则选到第h行第w列的单元素矩阵。 c,elem(h,w) 选元素a.elem(2,2)选择第2行第2列的值，结果不是矩阵，而是标量数值。速度比[][]和select()要快。且结果是引用返回，可以用来赋值：a.elem(2,2)=0; 3.2 矩阵重组运算矩阵的重组是指将己知阵的数值不变，而位置改变成另外的形式构成的新的矩阵，而原矩阵的实际存储形式未改变。类似于C语言的指针，这里的重组运算的结果只是原矩阵的一个引用，而不是一个有实际存储的矩阵，比如：int a[]={1,2,3,4};matrix<int> ma(a,2,2);内部存储：1 23 4转置运算：~ma;其结果是一个引用：r1,r2r3,r4 r1->1,r2->3,r3->2,r4->4 所以代表了矩阵：1 32 4即转置矩阵。 但是原a并未做任何改变，还是原存储形式。 这样的矩阵定义为引用矩阵，对引用矩阵的改变会改变它所引用的对象。比如： (~a)=identify(int(0),2,2); 那a的值会发生变化，因为它的引用被赋值了。 a,转置运算operator~~a b,左右并operator,int a[]={1,2,3,4};int b[]={5,6};matrix<int> ma(a,2,2),mb(b,2,1);(a,b)=foo23();(a,b)将得到一个2 x 3的引用矩阵，引用原来的a和b，结果为：1 2 53 4 6建议在左右并的时候加上一对括号，为了和C++中的,区别开。 c,上下并operator()int a[]={1,2,3,4};int b[]={5,6};matrix<int> ma(a,2,2),mb(b,1,2);(a)(b)=foo32();(a)(b)的结果将a和b上下合并成3 x 2的矩阵，合并后结果为：1 23 45 6 d,自定义重组运算函数为了使自定义函数能返回重组引用矩阵，需要一个类似C语言中指针的运算，称引用运算：operator>>a>>b;将使a成为b的引用矩阵，如果a原来是实矩阵，则会被释放掉，重组构成和b有同样结构的引用矩阵，如果b也是一个引用矩阵，则得到和b同样的矩阵。这个时候和a=b是不同的，赋值运算不同于引用，它会要求二者的大小是相同的，而引用不需要。 在自定义函数中，使用引用矩阵重组原矩阵，并将引用矩阵返回。 引用矩阵的确定化：operator!它将一个引用矩阵确定化为实阵，其值为引用阵对应的值。 例如： matrix<int> foo(){  matrix<int> a=1,b=2;  return a,b;} a,b得到的是一个引用矩阵，其值指向两个局部变量a,b，它们在函数返回时会被析构，所以返回的引用矩阵将变得无效，相当于成了野指针，所以需要先将他们确定化再返回。 matrix<int> foo(){  matrix<int> a=1,b=2;  return !(a,b);} 3.3 基本代数运算简单的说就是+，-和*，其运算结果依赖于元素类型定义的对应运算。a,矩阵和矩阵要求矩阵的大小是一致的，*运算要求前者的列数和后者的行数相一致。返回是一个临时实矩阵，在未赋值时的下一条语句前将被析构。 b,矩阵和数称为标量运算，标量的+，-和*，为对应元素的值相应运算。注意只能是矩阵和数不能数和矩阵，比如 a+1; // OK1+a; // error 如果想实现：1-a这样的运算，需要先a*(-1)+1，或者先构造一个全为1的矩阵再-a。 c,方阵的乘幂运算operator^(int n) a^100; 计算方阵a*a*...*a，100个a的连乘。其运算复杂度为O(m^3*logn)，m为a的阶，m^3是一次乘法的代价。 3.4 赋值运算 operator=a=b;如果a是实矩阵，则赋值给a，如果a是引用矩阵，则赋值给a所引用的矩阵；如果b是实矩阵，直接拷贝，如果b是引用矩阵，拷贝b所引用的矩阵的值。 其它几个：operator+=operator-=operator*=+=和-=比a=a+b;和a=a-b;的效率要高，而*=和a=a*b;的效率是相同的。 另外还有对应的标准赋值运算，其中operator=为将矩阵的每一个元素都赋值为该标量，比如： matrix<int> m;m.size(3,3);m=0; 将构造一个3 x 3的0矩阵。 3.5 其它运算 a,交换矩阵的值 operator| a|b; 交换a和b的值，要求a,b的大小一致。 b,自定义函数运算只简单写了几个数值运算的函数，有rref（和MATLAB里一样），inv（用rref实现的）,det（也是消元实现的，列主元消元），代码如下： template<class T>void rref(matrix<T> &ab){ matrix<T> m; if(ab.height>ab.width)  m>>~ab; else  m>>ab; int i,j; for(i=1;i<=m.height;++i) {  //从a[i,i]到a[n,i]找出最大元素所在行  int max=i;//max指向最大列主元素所在行  for(j=i+1;j<=m.height;++j)  {   if(fabs(m.elem(j,i))>fabs(m.elem(max,i)))    max=j;  }  //ab.swap(i,max);//交换行  if(i!=max)   m[i] | m[max];  if(m.elem(i,i)==0.0)//det=0，计算停止   return;  //将a[i,i]化成1.0  //ab.mul(i,1.0/ab.elem(i,i));  m[i]*=1.0/m.elem(i,i);  //消元  for(j=i+1;j<=m.height;++j)  {   m[j]-=m[i]*m.elem(j,i);   //ab.pt(i,j,-(ab.elem(j,i)/ab.elem(i,i)));  } } for(i=m.height;i>=1;--i) {  //消第i列上三角元素  for(j=1;j<i;++j)  {   m[j]-=m[i]*m.elem(j,i);   //ab.pt(i,j,-ab.elem(j,i));  } }} template<class T>matrix<T> inv(matrix<T> m){ matrix<T> result; if(m.height<=0 || m.width<=0 || m.height!=m.width) {  return result; } result=identity(T(0),m.height,m.width); rref((m,result)); return result;} template<class T>T det(matrix<T> m){ if(m.height<=0 || m.width<=0 || m.height!=m.width) {  return T(0); } int i,j; T s(1),sign(1); for(i=1;i<=m.height;++i) {  //从a[i,i]到a[n,i]找出最大元素所在行  int max=i;//max指向最大列主元素所在行  for(j=i+1;j<=m.height;++j)  {   if(fabs(m.elem(j,i))>fabs(m.elem(max,i)))    max=j;  }  //ab.swap(i,max);//交换行  if(i!=max)  {   m[i] | m[max];   sign=-sign;  }  if(m.elem(i,i)==T(0))//det=0，计算停止   return T(0);  //将a[i,i]化成1.0  //ab.mul(i,1.0/ab.elem(i,i));  //m[i]*=1.0/m.elem(i,i);  s*=m.elem(i,i);  //消元  for(j=i+1;j<=m.height;++j)  {   m[j]-=m[i]*(m.elem(j,i)/m.elem(i,i));   //ab.pt(i,j,-(ab.elem(j,i)/ab.elem(i,i)));  } } return s*sign;} 好了，主要的接口已经介绍完了，没有MATLAB里的函数多，那是当然的，听说MATLAB卖2万多USD，我这么点怎么可能和他比，如果你懂一些数值计算的原理，自己算或许更好，我试过用这个inv计算Herimit严重的变态矩阵，效果比MATLAB要好一些。你也可以自己试；） 另外一个在设计初没有料到的事，可以用它实现多值返回，比如前面提到的： mint foo(){  mint a=1,b=2;  return !(a,b);} 调用时，mint a=0,b=0;(a,b)=foo();之后，a就是1，b就是2，这样的多值返回就像MATLAB里一样！当然这不是它本来的用意，他是用来做数值处理的，但我没有定义死它非要是double型，你可以先定义一个比如分数类，再重载了+，-和*，再装到矩阵里也可以动作了。再比如，原来做的巨数类hugeint，将它封到矩阵里一样工作。 注1 ：目前版本2，在VC6上实现，如果有更好的C++编译器，将重载2维数组直接构造到矩阵的构造函数，使构造更简单。 注2 ：暂不提供实现源码，因为还没有严格测试，如果有兴趣可以来email交流。 注3 ：现在的实现只适合小型矩阵，一些维度上千的大型稀疏矩阵没办法处理，以后有机会改进。 rickone 2007/7/14

评论