# Matrix inverse by elimination

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 ?

• "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). – WooWapDaBug Jan 12 '18 at 8:39
• There is also the Eigen library, check this answer for computing the inverse with it – WooWapDaBug Jan 12 '18 at 8:43