Гауссово исключение без результата для ускорения

Question

Добрый день,

Я работаю над библиотекой C (для себя код: https://github.com/BattlestarSC/matrixLibrary.git) для работы с матричными функциями. В основном это учебная / практическая деятельность. Одна из моих задач - эффективно использовать определитель матрицы. Поскольку мои текущие попытки не увенчались успехом, я хотел выбрать другой подход. Я читал этот метод из документов MIT: http://web.mit.edu/18.06/www/Spring17/Determinants.pdf, и он имел большой смысл. Проблема, которую я имею, состоит в том, как получить к указанному пункту. Поскольку метод исключения Гаусса хорош для систем уравнений с несколькими переменными, мои матрицы не построены из уравнений и поэтому не являются частью системы. Как и в случае, каждое уравнение не имеет заданного результата и не вписывается в Форма из этой статьи здесь: https://math.oregonstate.edu/home/programs/undergrad/CalculusQuestStudyGuides/vcalc/gauss/gauss.html

С этого момента я не знаю, как поступить с этим методом.

Имеет смысл брать точку опоры из каждого набора уравнений, как описано в документе MIT, но как мне настроить свои матрицы, чтобы сделать указанный результат действительным?

Answer 1

Когда вы выполняете исключение по Гауссу, вы меняете строки и многократно вычитаете кратные из одной строки из другой, чтобы получить верхнюю треугольную форму.

Когда вы делаете это в системе уравнений или «расширенной матрице», вы не используете любую информацию из столбца результата. Решения о том, какие строки поменять местами, а какие вычесть с каким множителем, будут одинаковыми, независимо от того, какие числа находятся в столбце результата.

Поскольку «столбец результата» не используется, вы можете выполнить ту же процедуру для обычной квадратной матрицы. Поскольку операции не изменяют определитель (если вы отменяете одну строку при каждом обмене), в результате получается верхняя треугольная матрица с той же самой точкой det, что и в оригинале.

Автор MIT вызывает функцию lu, чтобы сделать это в примере рядом с началом. Это выполняет L-U-разложение на матрицу, которая возвращает исключенную по Гауссу матрицу в части 'U': https://en.wikipedia.org/wiki/LU_decomposition.

L-U разложение довольно круто. Это все равно что выполнять исключение по Гауссу для одновременного решения всех систем с одной и той же «матричной частью», что опять же можно сделать, потому что процессу вообще не нужно видеть столбец результата.

Начиная с матрицы M , вы получаете L и U , такие что LU = M . Это означает, что если вы хотите решить:

Mx = y

... где (x a y - векторы столбцов), у вас есть:

LUx = y

Решить Lv = y , что легко (просто подстановка), потому что L имеет нижнюю треугольную форму. Тогда у вас есть:

Ux = v

... который легко решить, потому что U имеет верхнюю треугольную форму.

Answer 2

GEM не очень хорош для компьютеров, так как он требует переупорядочения строк, поэтому алгоритм приводит к достоверному результату, который добавляет относительно большие накладные расходы и потенциальную нестабильность (если заказан плохо). GEM гораздо лучше подходит для людей и бумаги / карандаша, поскольку мы инстинктивно переупорядочиваем / выбираем строки ...

Итак, вы должны использовать подход (sub) Determinant , как вы и хотели. Это быстрее и безопаснее. Я знаю, что это немного сложно узнать из бумаги. Если это поможет, это моя древняя matrix.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 и (суб) определитель 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 - матрица.