没有加速结果的高斯消元

Question

早安，

我正在开发一个 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

从这一点来看，我不知道如何继续使用这种方法。

从麻省理工学院论文中描述的每组方程中获取枢轴点很有意义，但是我应该如何设置矩阵以使所述结果有效？

Answer 1

当您执行高斯消元法时，您交换行并反复从另一行减去一行的倍数以生成上三角形式。

当您在方程组或“增广矩阵”上执行此操作时，您不会使用结果列中的任何信息。无论结果列中的数字是什么，关于交换哪些行以及用什么乘数减去哪些行的决定都是完全相同的。

由于不使用“结果列”，您可以对普通方阵执行相同的过程。由于这些操作不会改变行列式（如果你每次交换时都取反），你最终会得到一个与原始矩阵具有相同 det 的上三角矩阵。

麻省理工学院的作者在开头附近的示例中调用了一个函数lu 来执行此操作。这会对矩阵进行 L-U 分解，在“U”部分返回高斯消除矩阵：https://en.wikipedia.org/wiki/LU_decomposition。

L-U 分解非常酷。这就像使用高斯消元法一次性解决所有具有相同“矩阵部分”的系统，您也可以这样做，因为该过程根本不需要查看结果列。

从矩阵 M 开始，您会得到 L 和 U，使得 LU = M。这意味着，如果你想解决：

Mx = y

...其中（x 和 y 是列向量），您有：

LUx = y

求解 Lv=y，这很容易（只是替换），因为 L 是下三角形。那么你有：

Ux = v

...这很容易解决，因为 U 是上三角形。

Answer 2

GEM 对计算机来说不是很好，因为它需要对行重新排序，因此算法会产生有效结果，这会增加相对较大的开销和潜在的不稳定性（如果排序不当）。 GEM 更适合人类和纸/铅笔，因为我们本能地重新排序/选择行...

因此，您首先应该按照自己的意愿使用 (sub)Determinant 方法。更快更安全。我知道从论文中学习它有点棘手。如果它有帮助，这是我的古老 ma​​trix.h class（但在 C++ 中）我还是个菜鸟时写的（所以我可能有一些隐藏的错误不知道好久没用这个了）：

//--- 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 inv_gaus 和 (sub)determinant inv_det 方法都存在，因此只需从中提取/比较您需要的内容。

顺便说一句，最近我需要一些用于 N 维空间的数学知识，一旦我使用它，我还将一个方阵编码为模板，其中 (sub)Determinant 方法作为递归模板完成 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 是矩阵。