Gaussian elimination without result for acceleration

Question

I'm working on a C library (for myself, code: https://github.com/BattlestarSC/matrixLibrary.git) to handle matrix functions. This is mostly a learning/practice activity.

One of my challenges is to take the determinant of a matrix efficiently. As my current attempts have failed, I wanted to take a different approach. I was reading though this method from MIT documentation: Determinants and it made a lot of sense.

The issue I'm having is how to get to said point. As the Gaussian elimination method is good for multi-variable systems of equations, my matrices are not built from equations, and therefore are not part of a system. As in, each equation doesn't have any set result and does not fit into the form from this paper.

From this point, I'm at a loss as far as how to proceed with this method.

It makes a lot of sense to take the pivot point from each set of equations as described in the MIT paper, but how should I set up my matrices to make said result valid?

Beauharnais

If you only have problems with the math, there's math.stackexchange.com for that. If you have problems with C code, you will have to include more info, like a minimal reproducible example. Also note that links to external resources are bad, because they can go out of scope easily. BTW: Use test-driven development (TDD). When applied, it makes sure you catch errors immediately. – Ladonnalady

The math.oregonstate.edu link is broken: "The requested page could not be found." – Daphie

Answer 1

When you perform a Gaussian elimination, you swap rows and repeatedly subtract a multiple of one row from another to produce an upper triangular form.

When you do this on a system of equations or an "augmented matrix", you do not use any information from the result column. The decisions about which rows to swap and which to subtract with what multiplier would be exactly the same no matter what numbers are in the result column.

Because the "result column" is not used, you can perform the same procedure on a normal square matrix. Since the operations don't change the determinant (if you negate one row whenever you swap), you end up with an upper triangular matrix with the same det as the original.

The MIT author calls a function lu to do this in the example near the start. This does L-U decomposition on the matrix, which returns the Gaussian-eliminated matrix in the 'U' part: https://en.wikipedia.org/wiki/LU_decomposition.

L-U decomposition is pretty cool. It's like doing Gaussian elimination to solve all systems with the same "matrix part" all at once, which again you can do because the process doesn't need to see the result column at all.

Starting with a matrix M, you get L and U such that LU = M. That means, if you want to solve:

Mx = y

... where (x an y are column vectors), you have:

LUx = y

Solve Lv=y, which is easy (just substitution) because L is lower-triangular. Then you have:

Ux = v

... which is easy to solve because U is upper-triangular.

Applicative

Answer 2

GEM needs to reorder the rows so the algorithm leads to a valid result which adds relatively big overhead and potential instability (if ordered badly). The GEM is much better suited for humans and paper/pencil as we instinctively reorder/chose rows ... but for bigger matrices is faster.

For small/medium sized matrices you should go with the (sub)determinant approach as you wanted in the first place. It is faster and safer. I know it’s a bit tricky to learn it from papers. If it helps, this is my ancient matrix.h class (but in C++) I wrote when I was still a rookie (so there might be some hidden bugs I do not know of; I haven't used this for ages):

//--- 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
//---------------------------------------------------------------------------

Both GEM inv_gaus and (sub)determinant inv_det approaches are present so just extract/compare from it what you need.

BTW, lately I needed some math stuff for N-dimensional space and once I was at it, I also coded a square matrix as a template where the (sub)determinant approach is done as recursive template 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
//---------------------------------------------------------------------------

But as you can see, that code is a bit more complicated to follow as I am in a different coding level now (look for inverse)...

If you also need results, then compute it as a matrix equation:

A*X = Y
  X = inv(A)*Y

Where X are unknowns (vector) , Y are knowns (vector) and A is the matrix.

Algar

@UlrichEckhardt yes I know I do not code in C but my intention was to extract the info OP author needs from code instead of papers which is in some cases like this faster (as in OP (s)he states (s)he got lost ... ) especially for programmers and C is very similar to C++ (appart classes and templates of coarse). The inverse function in C will be very similar to first example ... second would be possible only using macros I think ... but the second example is cleaner to understand the recursive math behind ... that is why I posted it too (also it contains the terminology) – Algar

While Laplace's expansion might be faster for very small matrices (eg 6x6 or smaller) it is of order n! and so will be catastrophically slow for moderate (eg 100x100) matrices and utterly infeasible for larger ones. – Appreciable

@Appreciable good point I do not use such large matrices (at least not for matrix math / linear algebra) – Algar