记得第一个跟我讲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 3
4 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 0
0 1 0
（和MATLAB里的eye相同）

3，矩阵的基本运算
由于这里的矩阵元素类型不一定为double型，所以不一定是用来做数值运算的，或者你也可以用其它你定义的类型，但有时是需要有+,-,*运算定义以及整数0，1的构造函数构造出0元和单位元的构造函数。有时你也可以把它想成一种容器（它现在的底层实现还不是很好），一种矩阵的容器。

3.1 选择运算
选择运算是最重要的，拿到一个矩阵要很容易选定它的元素或部分。
a,operator[] 选行（列）运算
矩阵的行和列是对等的，我这里定义的是行标为正数列标为负数，如下：
a -1 -2 -3 -4
1  1
2     1
3        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 2
3 4
转置运算：
~ma;
其结果是一个引用：
r1,r2
r3,r4

r1->1,r2->3,r3->2,r4->4

所以代表了矩阵：
1 3
2 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 5
3 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 2
3 4
5 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; // OK
1+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

评论