日安，
我正在开发一个C库（针对我自己，代码：https://github.com/BattlestarSC/matrixLibrary.git)以处理矩阵函数。这主要是一种学习/实践活动。我的一个挑战是有效地计算矩阵的行列式。由于我目前的尝试失败了，我想采取不同的方法。我从麻省理工学院的文件中读到了这个方法：http://web.mit.edu/18.06/www/Spring17/Determinants.pdf这很有道理。我面临的问题是如何达到上述目的。由于高斯消去法适用于多变量方程组，所以我的矩阵不是由方程建立的，因此不是系统的一部分。如中所示，每个方程都没有设定结果，也不符合本文的形式：https://math.oregonstate.edu/home/programs/undergrad/CalculusQuestStudyGuides/vcalc/gauss/gauss.html
从这一点来说，我不知道如何继续这个方法。
从麻省理工学院论文中描述的每组方程中获取枢轴点是很有意义的，但我应该如何设置矩阵以使所述结果有效？
当您执行高斯消去法时，您交换行并重复地从一行中减去另一行的倍数，以产生上三角形式。
当您在方程组或“增广矩阵”上执行此操作时，您不使用结果列中的任何信息。无论结果列中的数字是什么，关于交换哪些行以及用哪个乘数减去哪些行的决定都是完全相同的。
因为不使用结果列，所以您可以对普通方阵执行相同的过程。由于这些操作不会改变行列式（如果您在交换时否定一行），因此您最终会得到一个与原始值相同的上三角矩阵。
麻省理工学院的作者在开始附近的示例中调用了一个函数lu来实现这一点。这对矩阵进行L-U分解，返回“U”部分的高斯消去矩阵：https://en.wikipedia.org/wiki/LU_decomposition.
L-U分解非常酷。这就像做高斯消去一次求解所有具有相同“矩阵部分”的系统，这同样可以做到，因为该过程根本不需要看到结果列。
从矩阵M开始，你得到L和U，使得LU=M。这意味着，如果你想解决：
Mx=y
...其中（x an y 是列向量），您有：
勒克斯 = y
解决Lv=y，这很容易（只是替换），因为L是下三角形的。然后你有：
Ux = v
...这很容易解决，因为 U 是上三角形的。
GEM 对计算机来说不是很好，因为它需要对行重新排序，因此算法会导致一个有效的结果，从而增加相对较大的开销和潜在的不稳定性（如果排序不当）。GEM 更适合人类和纸/铅笔，因为我们本能地重新排序/选择行......
所以你应该按照你最初的想法，使用(次)行列式方法。更快更安全。我知道从纸上学习有点棘手。如果有帮助的话，这是我的古老的matrix.h class(但是用C语言写的)我还是个新手的时候写的(所以可能有一些我不知道的隐藏的bug，我已经很久没用这个了):
//--- matrix ver: 2.1 -------------------------------------------------------
#ifndef _matrix_h
#define _matrix_h
//---------------------------------------------------------------------------
double fabs(double x)
{
if (x<0) x=-x;
return x;
}
//---------------------------------------------------------------------------
class matrix
{
private:double **p;
int xs,ys;
double zeroacc;
public: matrix() { p=NULL; xs=0; ys=0; resize(1,1); zeroacc=1e-10; }
~matrix() { free(); }
void free();
int resize(int _xs,int _ys);
matrix& operator=(const matrix &b);
matrix& operator+();
matrix& operator-();
matrix& operator+(matrix &b);
matrix& operator-(matrix &b);
matrix& operator*(matrix &b);
matrix& operator+=(matrix &b);
matrix& operator-=(matrix &b);
matrix& operator*=(matrix &b);
matrix& operator!();
double& operator()(int y,int x);
double* operator[](int y) { return p[y]; }
void one();
int get_xs() { return xs; }
int get_ys() { return ys; }
double get_zeroacc() { return zeroacc; }
void set_zeroacc(double _zeroacc) { zeroacc=_zeroacc; if (zeroacc<0) zeroacc=-zeroacc; }
void ld(int y,double x0=0.0,double x1=0.0,double x2=0.0,double x3=0.0,double x4=0.0,double x5=0.0,double x6=0.0,double x7=0.0,double x8=0.0,double x9=0.0);
void prn(TCanvas *scr,int x0,int y0);
void lxch(int y1,int y2);
void lcom(int y1,int y2,double k);
void lmul(int y,double k);
void ldiv(int y,double k);
int gaus(matrix &b);
matrix& matrix::submatrix(int _x,int _y);
double determinant();
double subdeterminant();
matrix& inv_det();
matrix& inv_gaus();
};
//---------------------------------------------------------------------------
void matrix::free()
{
int y;
if (p!=NULL)
for (y=0;y<ys;y++)
delete[] p[y];
delete[] p;
p=NULL;
xs=0;
ys=0;
}
//---------------------------------------------------------------------------
int matrix::resize(int _xs,int _ys)
{
int y;
free();
if (_xs<1) _xs=1;
if (_ys<1) _ys=1;
xs=_xs;
ys=_ys;
p=new double*[ys];
if (p==NULL)
{
xs=0;
ys=0;
return 0;
}
for (y=0;y<ys;y++)
{
p[y]=new double[xs];
if (p[y]==NULL)
{
if (y>0)
for (y--;y>=0;y--)
delete p[y];
delete p;
p=NULL;
xs=0;
ys=0;
return 0;
}
}
return 1;
}
//---------------------------------------------------------------------------
matrix& matrix::operator=(const matrix &b)
{
int x,y;
if (!resize(b.get_xs(),b.get_ys())) return *this;
if (b.p)
for (y=0;y<ys;y++)
for (x=0;x<xs;x++)
p[y][x]=b.p[y][x];
return *this;
}
//---------------------------------------------------------------------------
matrix& matrix::operator+()
{
static matrix c;
int x,y;
c.resize(xs,ys);
for (y=0;y<ys;y++)
for (x=0;x<xs;x++)
c.p[y][x]= p[y][x];
return c;
}
//---------------------------------------------------------------------------
matrix& matrix::operator-()
{
static matrix c;
int x,y;
c.resize(xs,ys);
for (y=0;y<ys;y++)
for (x=0;x<xs;x++)
c.p[y][x]=-p[y][x];
return c;
}
//---------------------------------------------------------------------------
matrix& matrix::operator+(matrix &b)
{
static matrix c;
int x,y;
c.free();
if (xs!=b.get_xs()) return c;
if (ys!=b.get_ys()) return c;
c.resize(xs,ys);
for (y=0;y<ys;y++)
for (x=0;x<xs;x++)
c.p[y][x]=p[y][x]+b.p[y][x];
return c;
}
//---------------------------------------------------------------------------
matrix& matrix::operator-(matrix &b)
{
static matrix c;
int x,y;
c.free();
if (xs!=b.get_xs()) return c;
if (ys!=b.get_ys()) return c;
c.resize(xs,ys);
for (y=0;y<ys;y++)
for (x=0;x<xs;x++)
c.p[y][x]=p[y][x]-b.p[y][x];
return c;
}
//---------------------------------------------------------------------------
matrix& matrix::operator*(matrix &b)
{
static matrix c;
int i,j,k,ii,jj,kk;
c.free();
ii=ys;
jj=b.get_xs();
kk=b.get_ys();
if (kk!=xs) return c;
if (!c.resize(jj,ii)) return c;
for (i=0;i<ii;i++)
for (j=0;j<jj;j++)
c.p[i][j]=0.0;
for (i=0;i<ii;i++)
for (j=0;j<jj;j++)
for (k=0;k<kk;k++)
c.p[i][j]+=p[i][k]*b.p[k][j];
return c;
}
//---------------------------------------------------------------------------
matrix& matrix::operator+=(matrix &b)
{
int x,y;
if (xs!=b.get_xs()) { free(); return *this; }
if (ys!=b.get_ys()) { free(); return *this; }
for (y=0;y<ys;y++)
for (x=0;x<xs;x++)
p[y][x]+=b.p[y][x];
return *this;
}
//---------------------------------------------------------------------------
matrix& matrix::operator-=(matrix &b)
{
int x,y;
if (xs!=b.get_xs()) { free(); return *this; }
if (ys!=b.get_ys()) { free(); return *this; }
for (y=0;y<ys;y++)
for (x=0;x<xs;x++)
p[y][x]-=b.p[y][x];
return *this;
}
//---------------------------------------------------------------------------
matrix& matrix::operator*=(matrix &b)
{
matrix c;
int i,j,k,ii,jj,kk;
c.free();
ii=ys;
jj=b.get_xs();
kk=b.get_ys();
if (kk!=xs) { *this=c; return *this; }
if (!c.resize(jj,ii)) { *this=c; return *this; }
for (i=0;i<ii;i++)
for (j=0;j<jj;j++)
c.p[i][j]=0.0;
for (i=0;i<ii;i++)
for (j=0;j<jj;j++)
for (k=0;k<kk;k++)
c.p[i][j]+=p[i][k]*b.p[k][j];
*this=c; return *this;
}
//---------------------------------------------------------------------------
matrix& matrix::operator!()
{
// return inv_det();
return inv_gaus();
}
//---------------------------------------------------------------------------
double& matrix::operator()(int y,int x)
{
static double _null;
if (x<0) return _null;
if (y<0) return _null;
if (x>=xs) return _null;
if (y>=ys) return _null;
return p[y][x];
}
//---------------------------------------------------------------------------
void matrix::one()
{
int x,y;
for (y=0;y<ys;y++)
for (x=0;x<xs;x++)
if (x!=y) p[y][x]=0.0;
else p[y][x]=1.0;
}
//---------------------------------------------------------------------------
void matrix::ld(int y,double x0,double x1,double x2,double x3,double x4,double x5,double x6,double x7,double x8,double x9)
{
int x;
if (y<0) return;
if (y>=ys) return;
x=0;
if (x<xs) p[y][x]=x0; x++;
if (x<xs) p[y][x]=x1; x++;
if (x<xs) p[y][x]=x2; x++;
if (x<xs) p[y][x]=x3; x++;
if (x<xs) p[y][x]=x4; x++;
if (x<xs) p[y][x]=x5; x++;
if (x<xs) p[y][x]=x6; x++;
if (x<xs) p[y][x]=x7; x++;
if (x<xs) p[y][x]=x8; x++;
if (x<xs) p[y][x]=x9; x++;
}
//---------------------------------------------------------------------------
void matrix::prn(TCanvas *scr,int x0,int y0)
{
int x,y,xx,yy,dx,dy;
dx=50;
dy=13;
yy=y0;
for (y=0;y<ys;y++)
{
xx=x0;
for (x=0;x<xs;x++)
{
scr->TextOutA(xx,yy,AnsiString().sprintf("%.4lf",p[y][x]));
xx+=dx;
}
yy+=dy;
}
}
//---------------------------------------------------------------------------
void matrix::lxch(int y1,int y2)
{
int x;
double a;
if (y1<0) return;
if (y2<0) return;
if (y1>=ys) return;
if (y2>=ys) return;
for (x=0;x<xs;x++) { a=p[y1][x]; p[y1][x]=p[y2][x]; p[y2][x]=a; }
}
//---------------------------------------------------------------------------
void matrix::lcom(int y1,int y2,double k)
{
int x;
if (y1<0) return;
if (y2<0) return;
if (y1>=ys) return;
if (y2>=ys) return;
for (x=0;x<xs;x++) p[y1][x]+=p[y2][x]*k;
}
//---------------------------------------------------------------------------
void matrix::lmul(int y,double k)
{
int x;
if (y<0) return;
if (y>=ys) return;
for (x=0;x<xs;x++) p[y][x]*=k;
}
//---------------------------------------------------------------------------
void matrix::ldiv(int y,double k)
{
int x;
if (y<0) return;
if (y>=ys) return;
if ((k> zeroacc)||(k<-zeroacc)) k=1.0/k; else k=0.0;
for (x=0;x<xs;x++) p[y][x]*=k;
}
//---------------------------------------------------------------------------
int matrix::gaus(matrix &b)
{
int x,y;
double a;
if (xs!=ys) return 0;
if (ys!=b.ys) return 0;
for (x=0;x<xs;x++)
{
a=p[x][x]; // je aktualny prvok (x,x) na diagonale = 0 ?
if (a<0) a=-a;
if (a<=zeroacc)
for (y=0;y<ys;y++) // ak hej najdi nejaky nenulovy riadok v aktualnom stlpci (x)
if (x!=y)
{
a=p[y][x];
if (a<0) a=-a;
if (a>=zeroacc) // ak sa nasiel tak ho pripocitaj k aktualnemu riadku co zrusi tu nulu
{
b.lcom(x,y,1.0);
lcom(x,y,1.0);
break;
}
}
a=p[x][x]; // este raz otestuj ci na diagonale neni nula
if (a<0) a=-a;
if (a<=zeroacc) return 0; // ak je tak koniec
b.ldiv(x,p[x][x]); // sprav na diagonale 1-tku
ldiv(x,p[x][x]);
for (y=0;y<ys;y++) // a vynuluj zvysne riadky v stlpci(x)
if (y!=x)
{
b.lcom(y,x,-p[y][x]);
lcom(y,x,-p[y][x]);
}
}
return 1;
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
matrix& matrix::submatrix(int _x,int _y)
{
static matrix c;
int x,y,xx,yy;
c.resize(xs-1,ys-1);
yy=0; for (y=0;y<ys;y++)
if (y!=_y)
{
xx=0; for (x=0;x<xs;x++)
if (x!=_x)
{
c.p[yy][xx]=p[y][x];
xx++;
}
yy++;
}
return c;
}
//---------------------------------------------------------------------------
double matrix::determinant()
{
double D;
matrix a;
int x,y,s;
D=0;
if (xs!=ys) return D;
if (xs==1) { D=p[0][0]; return D; }
y=0;
s=y&1;
for (x=0;x<xs;x++)
{
a=submatrix(x,y);
if (s) D-=a.determinant()*p[y][x];
else D+=a.determinant()*p[y][x];
s=!s;
}
return D;
}
//---------------------------------------------------------------------------
double matrix::subdeterminant()
{
double D;
matrix a,b;
int x,y,s;
D=0;
if (xs!=ys) return D;
if (xs==1) { D=p[0][0]; return D; }
b=this[0];
for (y=0;y<ys;y++)
for (x=0;x<xs;x++)
{
a=b.submatrix(x,y);
p[y][x]=a.determinant();
}
y=0;
s=y&1;
for (x=0;x<xs;x++)
{
if (s) D-=p[y][x]*b.p[y][x];
else D+=p[y][x]*b.p[y][x];
s=!s;
}
return D;
}
//---------------------------------------------------------------------------
matrix& matrix::inv_det()
{
int x,y,s;
double D;
static matrix a,b;
a=this[0];
b=this[0];
D=b.subdeterminant();
if (fabs(D)>zeroacc) D=1.0/D;
for (y=0;y<ys;y++)
for (x=0;x<xs;x++)
{
s=(x+y)&1;
if (s) a.p[y][x]=-b.p[x][y]*D;
else a.p[y][x]= b.p[x][y]*D;
}
return a;
}
//---------------------------------------------------------------------------
matrix& matrix::inv_gaus()
{
static matrix a,b;
a=*this;
b.resize(xs,ys);
b.one();
a.gaus(b);
return b;
}
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
GEM＜code＞inv_gaus＜/code＞和（子）行列式＜code＞inv_det＜/code>方法都存在，因此只需从中提取/比较所需内容即可。
顺便说一下，最近我需要一些N维空间的数学知识，一旦我做到了这一点，我还编写了一个方阵作为模板，其中（子）行列式方法作为递归模板nd_math.h完成：
//--- N-Dimensional math ver: 1.002 -----------------------------------------
#ifndef _ND_math_h
#define _ND_math_h
//---------------------------------------------------------------------------
#include <math.h>
//---------------------------------------------------------------------------
#ifndef _rep4d_h
double divide(double a,double b) { if (fabs(b)<1e-30) return 0.0; return a/b; }
#endif
//---------------------------------------------------------------------------
template <const DWORD N> class vector
{
public:
double a[N];
vector() {}
vector(vector& a) { *this=a; }
~vector() {}
vector* operator = (const vector<N> *a) { *this=*a; return this; }
//vector* operator = (vector<N> &a) { ...copy... return this; }
double& operator [](const int i) { return a[i]; }
vector<N> operator + () { return *this; } // =+v0
vector<N> operator - () { int i; vector<N> q; for ( i=0;i<N;i++) q.a[i]= -a[i]; return q; } // =-v0
vector<N> operator + (vector<N> &v) { int i; vector<N> q; for ( i=0;i<N;i++) q.a[i]=a[i]+v.a[i]; return q; } // =v0+v1
vector<N> operator - (vector<N> &v) { int i; vector<N> q; for ( i=0;i<N;i++) q.a[i]=a[i]-v.a[i]; return q; } // =v0-v1
double operator * (vector<N> &v) { int i; double q; for (q=0.0,i=0;i<N;i++) q +=a[i]*v.a[i]; return q; } // =(v0.v1) dot product
vector<N> operator + (const double &c) { int i; vector<N> q; for ( i=0;i<N;i++) q.a[i]=a[i]+c; return q; } // =v0+(c,c,c,c,...)
vector<N> operator - (const double &c) { int i; vector<N> q; for ( i=0;i<N;i++) q.a[i]=a[i]-c; return q; } // =v0-(c,c,c,c,...)
vector<N> operator * (const double &c) { int i; vector<N> q; for ( i=0;i<N;i++) q.a[i]=a[i]*c; return q; } // =v0*c
vector<N> operator / ( double c) { int i; vector<N> q; c=divide(1.0,c); for ( i=0;i<N;i++) q.a[i]=a[i]*c; return q; } // =v0/c
vector<N> operator +=(vector<N> &v) { this[0]=this[0]+v; return *this; }; // v0+=v1
vector<N> operator -=(vector<N> &v) { this[0]=this[0]-v; return *this; }; // v0-=v1
vector<N> operator +=(const double &c) { this[0]=this[0]+c; return *this; }; // v0+=(c,c,c,c,...)
vector<N> operator -=(const double &c) { this[0]=this[0]-c; return *this; }; // v0-=(c,c,c,c,...)
vector<N> operator *=(const double &c) { this[0]=this[0]*c; return *this; }; // v0*=c
vector<N> operator /=(const double &c) { this[0]=this[0]/c; return *this; }; // v0/=c
AnsiString str() { int i; AnsiString q; for (q="( ",i=0;i<N;i++) q+=AnsiString().sprintf("%6.3lf ",a[i]); q+=")"; return q; }
double len() { int i; double l; for (l=0.0,i=0;i<N;i++) l+=a[i]*a[i]; return sqrt(l); } // get size
double len2() { int i; double l; for (l=0.0,i=0;i<N;i++) l+=a[i]*a[i]; return l; } // get size^2
void len(double l) { int i; l=divide(l,len()); for (i=0;i<N;i++) a[i]*=l; } // set size
void unit() { len(1.0); } // set unit size
void zero() { int i; for (i=0;i<N;i++) a[i]=0.0; } // set zero vector
void rnd() { int i; for (i=0;i<N;i++) a[i]=(2.0*Random())-1.0; } // set random unit vector
void set(double c) { int i; for (i=0;i<N;i++) a[i]=c; } // (c,c,c,c,...)
// i x j = k | | i j k |
// j x k = i | a x b = det | a0 a1 a2 | = + i*det | a1 a2 | - j*det | a0 a2 | + k*det | a0 a1 |
// k x i = j | | b0 b1 b2 | | b1 b2 | | b0 b2 | | b0 b1 |
void cross(const vector<N> *v)
{
int i,j;
matrix<N> m0;
matrix<N-1> m;
for (i=1;i<N;i++)
for (j=0;j<N;j++)
m0.a[i][j]=v[i-1].a[j];
for (j=0;j<N;j++)
{
m=m0.submatrix(0,j);
if (int(j&1)==0) a[j]=+m.det();
else a[j]=-m.det();
}
}
void cross(vector<N> **v)
{
int i,j;
matrix<N> m0;
matrix<N-1> m;
for (i=1;i<N;i++)
for (j=0;j<N;j++)
m0.a[i][j]=v[i-1]->a[j];
for (j=0;j<N;j++)
{
m=m0.submatrix(0,j);
if (int(j&1)==0) a[j]=+m.det();
else a[j]=-m.det();
}
}
void cross(vector<N> &v0) { vector<N> *v[ 1]={&v0}; cross(v); }
void cross(vector<N> &v0,vector<N> &v1) { vector<N> *v[ 2]={&v0,&v1}; cross(v); }
void cross(vector<N> &v0,vector<N> &v1,vector<N> &v2) { vector<N> *v[ 3]={&v0,&v1,&v2}; cross(v); }
void cross(vector<N> &v0,vector<N> &v1,vector<N> &v2,vector<N> &v3) { vector<N> *v[ 4]={&v0,&v1,&v2,&v3}; cross(v); }
void cross(vector<N> &v0,vector<N> &v1,vector<N> &v2,vector<N> &v3,vector<N> &v4) { vector<N> *v[ 5]={&v0,&v1,&v2,&v3,&v4}; cross(v); }
void cross(vector<N> &v0,vector<N> &v1,vector<N> &v2,vector<N> &v3,vector<N> &v4,vector<N> &v5) { vector<N> *v[ 6]={&v0,&v1,&v2,&v3,&v4,&v5}; cross(v); }
void cross(vector<N> &v0,vector<N> &v1,vector<N> &v2,vector<N> &v3,vector<N> &v4,vector<N> &v5,vector<N> &v6) { vector<N> *v[ 7]={&v0,&v1,&v2,&v3,&v4,&v5,&v6}; cross(v); }
void cross(vector<N> &v0,vector<N> &v1,vector<N> &v2,vector<N> &v3,vector<N> &v4,vector<N> &v5,vector<N> &v6,vector<N> &v7) { vector<N> *v[ 8]={&v0,&v1,&v2,&v3,&v4,&v5,&v6,&v7}; cross(v); }
void cross(vector<N> &v0,vector<N> &v1,vector<N> &v2,vector<N> &v3,vector<N> &v4,vector<N> &v5,vector<N> &v6,vector<N> &v7,vector<N> &v8) { vector<N> *v[ 9]={&v0,&v1,&v2,&v3,&v4,&v5,&v6,&v7,v8}; cross(v); }
void cross(vector<N> &v0,vector<N> &v1,vector<N> &v2,vector<N> &v3,vector<N> &v4,vector<N> &v5,vector<N> &v6,vector<N> &v7,vector<N> &v8,vector<N> &v9) { vector<N> *v[10]={&v0,&v1,&v2,&v3,&v4,&v5,&v6,&v7,v8,v9}; cross(v); }
void ld(const double &a0) { a[0]=a0; }
void ld(const double &a0,const double &a1) { a[0]=a0; a[1]=a1; }
void ld(const double &a0,const double &a1,const double &a2) { a[0]=a0; a[1]=a1; a[2]=a2; }
void ld(const double &a0,const double &a1,const double &a2,const double &a3) { a[0]=a0; a[1]=a1; a[2]=a2; a[3]=a3; }
void ld(const double &a0,const double &a1,const double &a2,const double &a3,const double &a4) { a[0]=a0; a[1]=a1; a[2]=a2; a[3]=a3; a[4]=a4; }
void ld(const double &a0,const double &a1,const double &a2,const double &a3,const double &a4,const double &a5) { a[0]=a0; a[1]=a1; a[2]=a2; a[3]=a3; a[4]=a4; a[5]=a5; }
void ld(const double &a0,const double &a1,const double &a2,const double &a3,const double &a4,const double &a5,const double &a6) { a[0]=a0; a[1]=a1; a[2]=a2; a[3]=a3; a[4]=a4; a[5]=a5; a[6]=a6; }
void ld(const double &a0,const double &a1,const double &a2,const double &a3,const double &a4,const double &a5,const double &a6,const double &a7) { a[0]=a0; a[1]=a1; a[2]=a2; a[3]=a3; a[4]=a4; a[5]=a5; a[6]=a6; a[7]=a7; }
void ld(const double &a0,const double &a1,const double &a2,const double &a3,const double &a4,const double &a5,const double &a6,const double &a7,const double &a8) { a[0]=a0; a[1]=a1; a[2]=a2; a[3]=a3; a[4]=a4; a[5]=a5; a[6]=a6; a[7]=a7; a[8]=a8; }
void ld(const double &a0,const double &a1,const double &a2,const double &a3,const double &a4,const double &a5,const double &a6,const double &a7,const double &a8,const double &a9) { a[0]=a0; a[1]=a1; a[2]=a2; a[3]=a3; a[4]=a4; a[5]=a5; a[6]=a6; a[7]=a7; a[8]=a8; a[9]=a9; }
};
//---------------------------------------------------------------------------
template <DWORD N> class matrix // square matrix
{
public:
vector<N> a[N];
matrix() {}
matrix(matrix& a) { *this=a; }
~matrix() {}
matrix* operator = (const matrix<N> *a) { *this=*a; return this; }
//matrix* operator = (matrix<N> &a) { ...copy... return this; }
vector<N>& operator [](const int i) { return a[i]; }
matrix<N> operator + () { return *this; }
matrix<N> operator - () { matrix<N> q; int i,j; for (i=0;i<M;i++) for (j=0;j<N;j++) q[i][j]=-a[i][j]; return q; } // = -m0
matrix<N> operator * (const matrix &m)
{
matrix<N> q;
int i,j,k;
for (i=0;i<N;i++)
for (j=0;j<N;j++)
for (q.a[i][j]=0.0,k=0;k<N;k++)
q.a[i].a[j]+=a[i].a[k]*m.a[k].a[j];
return q;
}
vector<N> operator * (vector<N> &v)
{
vector<N> q;
int i,j;
for (i=0;i<N;i++)
for (q.a[i]=0.0,j=0;j<N;j++)
q.a[i]+=a[i][j]*v.a[j];
return q;
}
matrix<N> operator * (const double &c)
{
matrix<N> q;
int i,j;
for (i=0;i<N;i++)
for (j=0;j<N;j++)
q.a[i].a[j]=a[i].a[j]*c;
return q;
}
matrix<N> operator / (const double &c)
{
return this[0]*divide(1.0,c);
}
matrix<N> operator *=(matrix<N> &m) { this[0]=this[0]*m; return *this; };
vector<N> operator *=(vector<N> &v) { this[0]=this[0]*v; return *this; };
matrix<N> operator *=(const double &c) { this[0]=this[0]*c; return *this; };
matrix<N> operator /=(const double &c) { this[0]=this[0]/c; return *this; };
AnsiString str() { int i,j; AnsiString q; for (q="",i=0;i<N;i++,q+="\r\n") { for (q+="( ",j=0;j<N;j++) q+=AnsiString().sprintf("%6.3lf ",a[i][j]); q+=")"; } return q; }
void unit() { int i,j; for (i=0;i<N;a[i][i]=1.0,i++) for (j=0;j<N;j++) a[i][j]=0.0; } // set unit matrix
void zero() { int i,j; for (i=0;i<N;i++) for (j=0;j<N;j++) a[i][j]=0.0; } // set zero matrix
void rnd() { int i,j; for (i=0;i<N;i++) for (j=0;j<N;j++) a[i][j]=(2.0*Random())-1.0; } // set random <-1,+1> matrix
void set(double c) { int i,j; for (i=0;i<N;i++) for (j=0;j<N;j++) a[i][j]=c; } // (c,c,c,c,...)
void orthonormal() // convert to orthonormal matrix
{
int i,j;
vector<N> *pV[N],*pp;
for (i=0;i<N;i++) { a[i].unit(); pV[i]=a+i; }
for (i=1;i<N;i++)
{
pV[0]->cross(pV+1);
pp=pV[0]; for (j=1;j<N;j++) pV[j-1]=pV[j]; pV[N-1]=pp;
}
}
matrix<N> transpose()
{
int i,j;
matrix<N> M;
for (i=0;i<N;i++)
for (j=0;j<N;j++)
M[i][j]=a[j][i];
return M;
}
matrix<N> inverse()
{
return adjugate()/det();
}
matrix<N> adjugate()
{
matrix<N> C;
double s;
int i,j;
for (i=0;i<N;i++)
for ((i&1)?s=-1.0:s=+1.0,j=0;j<N;j++,s=-s)
C[j][i]=minor(i,j)*s;
return C;
}
matrix<N> cofactor()
{
matrix<N> C;
double s;
int i,j;
for (i=0;i<N;i++)
for ((i&1)?s=+1.0:s=-1.0,j=0;j<N;j++,s=-s)
C[i][j]=minor(i,j)*s;
return C;
}
double minor(int i,int j)
{
return submatrix(i,j).det();
}
matrix<N-1> submatrix(int i,int j)
{
matrix<N-1> m;
int i0,i1,j0,j1;
for (i0=0,i1=0;i1<N;i1++)
if (i1!=i){ for (j0=0,j1=0;j1<N;j1++)
if (j1!=j){ m.a[i0][j0]=a[i1][j1]; j0++; } i0++; }
return m;
}
double det();
};
//---------------------------------------------------------------------------
double matrix<1>::det() { return a[0][0]; }
double matrix<2>::det() { return (a[0][0]*a[1][1])-(a[0][1]*a[1][0]); }
template <DWORD N> double matrix<N>::det()
{
double d=0.0; int j;
matrix<N-1> m;
for (j=0;j<N;j++)
{
m=submatrix(0,j);
if (int(j&1)==0) d+=a[0][j]*m.det();
else d-=a[0][j]*m.det();
}
return d;
}
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
但是正如你所看到的，代码有点复杂，因为我现在处于不同的编码级别（查找inverse）…
如果您还需要结果，则将其计算为矩阵方程：
A*X = Y
X = inv(A)*Y
其中，X是未知数（矢量），Y为未知数（向量），A则是矩阵。
