2
\$\begingroup\$

Because I had trouble in finding a reliable source on the internet for the matrix inverse algorithm using Gaussian elimination, I wrote one myself.

Matrix.h

#ifndef MAT4X4_H
#define MAT4X4_H

class Mat4x4{
    public:
        float mat4x4[4][4];

        Mat4x4();
        Mat4x4(const Mat4x4& m);
        Mat4x4(const float m[4][4]);
        Mat4x4(float t00, float t01, float t02, float t03 ,
            float t10, float t11, float t12, float t13 ,
            float t20, float t21, float t22, float t23 ,
            float t30, float t31, float t32, float t33 );

        static Mat4x4 transpose(const Mat4x4 &m);
        static Mat4x4 inverse(Mat4x4 m);
        bool operator==(const Mat4x4 &m);
        bool operator!=(const Mat4x4 &m);
        Mat4x4 operator*(const Mat4x4 &m);
        void print();
};

#endif

Matrix.cpp

#include "Mat4x4.h"
#include <cstring>
#include <cstdio>
#include <iostream>

Mat4x4::Mat4x4(){
    memset(mat4x4, 0, sizeof(mat4x4));
}


Mat4x4::Mat4x4(const Mat4x4& m){
    memset(mat4x4, 0, sizeof(mat4x4));
    for(int i = 0; i < 4; i++){
        for(int j = 0; j < 4; j++){
            mat4x4[i][j] = m.mat4x4[i][j];
        }
    }

}

Mat4x4::Mat4x4(const float m[4][4]) {
    for(int i = 0; i < 4; i++){
        for(int j = 0; j < 4; j++){
            mat4x4[i][j] = m[i][j];        
        }    
    }
}

Mat4x4::Mat4x4(float t00, float t01, float t02, float t03 ,
            float t10, float t11, float t12, float t13 ,
            float t20, float t21, float t22, float t23 ,
            float t30, float t31, float t32, float t33 )

    {

    mat4x4[0][0] = t00; mat4x4[0][1] = t01; mat4x4[0][2] = t02; mat4x4[0][3] = t03;
    mat4x4[1][0] = t10; mat4x4[1][1] = t11; mat4x4[1][2] = t12; mat4x4[1][3] = t13;
    mat4x4[2][0] = t20; mat4x4[2][1] = t21; mat4x4[2][2] = t22; mat4x4[2][3] = t23;
    mat4x4[3][0] = t30; mat4x4[3][1] = t31; mat4x4[3][2] = t32; mat4x4[3][3] = t33;
}

Mat4x4 Mat4x4::transpose(const Mat4x4 &m) {
    float temp = 0.0f;
    Mat4x4 t(m.mat4x4);
    for(int i = 0; i < 4; i++){
        for(int j = 0; j < i; j++){

            temp = t.mat4x4[i][j];
            t.mat4x4[i][j] = t.mat4x4[j][i];
            t.mat4x4[j][i] = temp;
        }    
    }

    return t;
}

Mat4x4 Mat4x4::inverse(Mat4x4 m){
    Mat4x4 r(1,0,0,0, 
            0,1,0,0,
            0,0,1,0,
            0,0,0,1);

    int flag = 0;

    for(int i = 0; i < 4; i++){
        //Row swap
        if(m.mat4x4[i][i] == 0){
            for(int j = i+1; (j < 4); j++){
                if(m.mat4x4[j][i] != 0){
                    flag = 1;

                    for(int k = 0; k < 4; k++){
                        float temp = m.mat4x4[i][k];
                        m.mat4x4[i][k] = m.mat4x4[j][k];
                        m.mat4x4[j][k] = temp;

                        temp = r.mat4x4[i][k];
                        r.mat4x4[i][k] = r.mat4x4[j][k];
                        r.mat4x4[j][k] = temp;
                    }

                    break;
                }
            }

            if(flag == 0) {
                //return a zero matrix. Failure.
                return Mat4x4();
            }
        }
        //Normalising the pivots    
        if(m.mat4x4[i][i] != 1){ 
            float f = m.mat4x4[i][i];
            for(int j = 0; j < 4; j++){
                m.mat4x4[i][j] /= f; 
                r.mat4x4[i][j] /= f;
            }
        }
        //Row reduction. 
        for(int j = 0; (j < 4); j++){
            if(j == i) continue;
            if(m.mat4x4[j][i] == 0) continue;

            float f = m.mat4x4[j][i];
            m.mat4x4[j][i] = 0;
            for(int k = i+1; k < 4; k++){
                m.mat4x4[j][k] -= f*m.mat4x4[i][k]; 
            }

            for(int l = 0; l < 4; l++){
                r.mat4x4[j][l] -= f*r.mat4x4[i][l];
            }
        }

        flag = 0;
    }

    return r;
}

bool Mat4x4::operator==(const Mat4x4 &m){
    for(int i = 0; i < 4; i++){
        for(int j = 0; j < 4; j++){
            if(mat4x4[i][j] != m.mat4x4[i][j]){
                return false;
            }
        }    
    }

    return true;
}

bool Mat4x4::operator!=(const Mat4x4 &m){
    for(int i = 0; i < 4; i++){
        for(int j = 0; j < 4; j++){
            if(mat4x4[i][j] != m.mat4x4[i][j]){
                return true;
            }
        }    
    }

    return false;
}

Mat4x4 Mat4x4::operator*(const Mat4x4 &m){
    Mat4x4 r;
    for(int i = 0; i < 4; i++){
        for(int j = 0; j < 4; j++){

            for(int k = 0; k < 4; k++ ){
                r.mat4x4[i][j] += mat4x4[i][k]*m.mat4x4[k][j];
            }

        }    
    }

    return r;
}

void Mat4x4::print(){
    printf("t00 = %f, t01 = %f, t02 = %f, t03 = %f\n"
        "t10 = %f, t11 = %f, t12 = %f, t13 = %f\n"
        "t20 = %f, t21 = %f, t22 = %f, t23 = %f\n"
        "t30 = %f, t31 = %f, t32 = %f, t33 = %f\n", 

            mat4x4[0][0], mat4x4[0][1], mat4x4[0][2], mat4x4[0][3],
            mat4x4[1][0], mat4x4[1][1], mat4x4[1][2], mat4x4[1][3],
            mat4x4[2][0], mat4x4[2][1], mat4x4[2][2], mat4x4[2][3],
            mat4x4[3][0], mat4x4[3][1], mat4x4[3][2], mat4x4[3][3]);
}

Please don't bother with method other than Mat4x4 Mat4x4::inverse(Mat4x4 m); I provided them so that it would be easier to run my program.

How can I make my function more efficient ?

\$\endgroup\$
  • \$\begingroup\$ "I had trouble in finding a reliable source on the internet for the matrix inverse algorithm using Gaussian elimination" Why? After a quick search I found this : GSL LU Descomposition The algorithm used in the decomposition is Gaussian Elimination with partial pivoting (Golub & Van Loan, Matrix Computations, Algorithm 3.4.1). \$\endgroup\$ – WooWapDaBug Jan 12 '18 at 8:39
  • \$\begingroup\$ There is also the Eigen library, check this answer for computing the inverse with it \$\endgroup\$ – WooWapDaBug Jan 12 '18 at 8:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.