Determinates And Inverse Matrix

[steve_bank]

The last piece of my math library.

The last piece of my math library, determinants and inverse matrix.

This one I did not write, I cleaned it up and made useful from code on the the net.

Inverse Matrix in C - Sanfoundry

This C program sorts a given array of integer numbers using Bubble Sort technique. The algorithm gets its name from the way smaller elements “bubble” to the top of the list. Because it only uses comparisons to operate on elements, it is a comparison sort. Time Complexity of this algorithm is...

www.sanfoundry.com

The are different ways to do it.

11.4: Determinants and Cramer's Rule for n x n Matrices

math.libretexts.org

I tested it against Scilab an established math tool.

class MATINV{
    public:
    double det(int n,double *x);
    void minv(int n,double *x,double *xinv);
    private:
    double Det(int k,double a[25][25]);
    void cofactor(int f,double num[25][25],double fac[25][25]);
};

double MATINV::det(int n,double *x){
    int r,c,k=0;
    double a[25][25];
    for(r=0;r<n;r++)for(c=0;c<n;c++) a[r][c] = x[k++];
    return Det(n,a);
}//det()

double MATINV::Det(int k,double a[25][25]){
    double s = 1, determ = 0;
    double b[25][25];
    int i, j, m, n, c;
    if (k == 1)return (a[0][0]);
    determ = 0;
    for (c = 0; c < k; c++){
        m = 0;
        n = 0;
        for (i = 0;i < k; i++){
            for (j = 0 ;j < k; j++){
                b[i][j] = 0;
                if (i != 0 && j != c){
                    b[m][n] = a[i][j];
                    if (n < (k - 2)){n++;}
                       else{n = 0;m++;}
                }
              }
           }
            determ = determ + s * (a[0][c] * Det(k-1,b));
            s = -1. * s;
        }
    return determ;
} //Det()

void MATINV::cofactor(int f,double num[25][25],double fac[25][25]){
    double b[25][25];
    int p, q, m, n, i, j;
    for (q = 0;q < f; q++){
        for (p = 0;p < f; p++){
            m = 0;
            n = 0;
             for (i = 0;i < f; i++){
                for (j = 0;j < f; j++){
                    if (i != q && j != p){
                        b[m][n] = num[i][j];
                        if (n < (f - 2))
                            n++;
                         else{
                            n = 0;
                            m++;
                            }
                    }//if
              }//for j
              fac[q][p] = pow(-1, q + p) * Det(f - 1,b);
            }//for i
        }//for p
    }//for q
}//cofactor()

void MATINV::minv(int n,double *x,double *xinv){
    int i, j,k=0;
    double a[25][25],inverse[25][25],fac[25][25], d;
    for(i=0;i<n;i++)for(j=0;j<n;j++) a[i][j] = x[k++];
    cofactor(n,a,fac);
    d = Det(n,a);
    k = 0;
    for (i = 0;i < n; i++)for (j = 0;j < n; j++){
        inverse[i][j] = fac[i][j] / d;
        xinv[k++] = inverse[i][j];
        }
}//minv()


int main(){
    int i;
    MATINV mi;
   
    double de,xi[3][3],x[3][3] = {{3,5,2},{1,5,8},{3,9,2}};
    de = mi.det(n,&x[0][0]);
    cout<<de<<endl;
    mi.minv(3,&x[0][0],&xi[0][0]);
    for(int i = 0;i<3;i++){
        for(int j=0;j<3;j++)printf("%.10f\t",xi[i][j]);
        printf("\n");
        }
    return 0;
}

 

[lpetrich]

lkpetrich/Semisimple-Lie-Algebras: Does "semisimple" Lie algebras, important for describing elementary-particle and space-time symmetries. Sets up algebras and representation basis sets in them. Finds products and powers of representations ("plethysms"), also subalgebras and representations in them ("branching rules"). Has several examples, like Grand Unified Theories.

I wrote a solver of systems of linear equations because I wanted to invert some matrices, and I wanted to do so with exact rational numbers, not floating-point ones. To invert a matrix, run one's solver, but with columns taken from an identity matrix: {1, 0, ..., 0}, {0, 1, ..., 0}, ..., {0, 0, ..., 1}. Each column's result is a column in the inverse matrix.

To keep it simple, I used Gaussian elimination with pivoting, to avoid dividing by zero. I took the original matrix and extended each row with the appropriate row of an identity matrix, and then solved. The extra rows then contained the inverse. I could have used LU decomposition, but I didn't bother. LU has an advantage over Gauss in that one can construct a quick solver with it, one that runs O(n^2) instead of O(n^3) for each solution vector. But both algorithms take O(n^3) for doing inverses.

My code runs in Mathematica, Python, and C++. Mathematica has rational numbers as a built-in type, Python has rational numbers in its standard library, though they are much slower than Python integers, and I found a C++ rational-number object class online.