Below is the code to calculate matrix inverse of a matrix of arbitrary size (order) by using analytic solution. This method is known to be slow for very large matrix because of the recursion. However, I used this mainly for calculating inverse of 4×4 matrices and it worked just fine. You can also use CalcDeterminant and GetMinor to calculate determinant of a matrix.

**MatrixInversion: **the main function to calculate the inverse matrix.**
CalcDeterminant**: calculate the determinant of a matrix.

**Get minor matrix.**

GetMinor:

GetMinor:

// matrix inversioon // the result is put in Y void MatrixInversion(float **A, int order, float **Y) { // get the determinant of a double det = 1.0/CalcDeterminant(A,order); // memory allocation float *temp = new float[(order-1)*(order-1)]; float **minor = new float*[order-1]; for(int i=0;i<order-1;i++) minor[i] = temp+(i*(order-1)); for(int j=0;j<order;j++) { for(int i=0;i<order;i++) { // get the co-factor (matrix) of A(j,i) GetMinor(A,minor,j,i,order); Y[i][j] = det*CalcDeterminant(minor,order-1); if( (i+j)%2 == 1) Y[i][j] = -Y[i][j]; } } // release memory //delete [] minor[0]; delete [] temp; delete [] minor; } // calculate the cofactor of element (row,col) int GetMinor(float **src, float **dest, int row, int col, int order) { // indicate which col and row is being copied to dest int colCount=0,rowCount=0; for(int i = 0; i < order; i++ ) { if( i != row ) { colCount = 0; for(int j = 0; j < order; j++ ) { // when j is not the element if( j != col ) { dest[rowCount][colCount] = src[i][j]; colCount++; } } rowCount++; } } return 1; } // Calculate the determinant recursively. double CalcDeterminant( float **mat, int order) { // order must be >= 0 // stop the recursion when matrix is a single element if( order == 1 ) return mat[0][0]; // the determinant value float det = 0; // allocate the cofactor matrix float **minor; minor = new float*[order-1]; for(int i=0;i<order-1;i++) minor[i] = new float[order-1]; for(int i = 0; i < order; i++ ) { // get minor of element (0,i) GetMinor( mat, minor, 0, i , order); // the recusion is here! det += (i%2==1?-1.0:1.0) * mat[0][i] * CalcDeterminant(minor,order-1); //det += pow( -1.0, i ) * mat[0][i] * CalcDeterminant( minor,order-1 ); } // release memory for(int i=0;i<order-1;i++) delete [] minor[i]; delete [] minor; return det; }

If you need to calculate matrix inverse frequently in your codes, you might want to optimize the code by declaring the matrix such as **minor **and **temp** as global variables and allocate memories for them only once.

implementation

matrix inversion in C++

Hi,

Found your webpage in google.

I think there is a little mistake in your source code l.47 :

dest[a][b] = src[i][j]; I don’t see any a or b but maybe is it i and j ?

Hi, Fab,

Thanks a lot for finding the bug. It should be dest[rowCount][colCount] instead of dest[a][b].

Hi,

I found your code for matrix inversion in google, too. But unfortunately I get a runtime error while trying to execute the code. I suppose the problem is conserned to dynamic memory allocation. Would you take a look at the main function, pls?

I use Dev C++ 4.9.9.1. Here is the code of the main function:

int main()

{ float cov[6][6]={ {228.585, 120.904, 147.073, 147.073, 228.585, 356.428},

{164.066, 85.4136, 104.528, 104.528, 164.066, 257.445},

{130.128, 66.7447, 82.1484, 82.1484, 130.128, 205.378},

{130.128, 66.7447, 82.1484, 82.1484, 130.128, 205.378},

{228.585, 120.904, 147.073, 147.073, 228.585, 356.428},

{286.483, 152.753, 185.253, 185.253, 286.483, 445.253}};

float **cov1 = new float *[6];

for(int p=0;p<6;p++)cov1[p]=new float [6];

for (int i=0; i<6; i++)

for (int j=0; j<6; j++) cov1[i][j]=cov[i][j];

for (int i=0; i<6; i++)

{for (int j=0; j<6; j++) cout<<cov1[i][j]<<” “;

cout<<endl;

}

float **res= new float *[6];

for(int p=0;p<6;p++)res[p]=new float [6];

MatrixInversion(cov1,6,res);

for (int i=0; i<6; i++)

{for (int j=0; j<6; j++) cout<<res[i][j]<<” “;

cout<<endl;

}

for(int i=0;i<6;i++)

delete [] cov1[i];

delete [] cov1;

for(int i=0;i<6;i++)

delete [] res[i];

delete [] res;

system (“pause”);

return 0;

}

Quick question. I was trying to implement this, but I don’t know how to declare the initial matrix.

I tried:

float MyMatrix[14][14];

but after I fill it and stick it in the function, I get:

error: cannot convert ‘float (*)[14]’ to ‘float**’ …

Jason,

MyMatrix has to be dynamically allocated. For example,

float **MyMatrix;

MyMatrix = new float*[14];

for(int i=0;i<14;i++)

MyMatrix[i] = new float[14];

Hope this helps!

Couldn’t you also typecast the input parameter without the need for any “new” allocated matrix?

float A[2][2] = {{1,2},{3,4}};

float Y[2][2];

MatrixInversion((float **)A, 2, (float **)Y);

Or would this not work properly?

Also one more suggestion would be to make the variables “static” instead of “global” to speed up frequent usage of the functions, but neither method will work for the recursive functions.

Thanks, that did help. FYI though, your GetMinor( float **dest, float **src, int row, int col, int order) has the source and destination switched as compared to all the function calls.

e.g. GetMinor( mat, minor, 0, i , order); //where either the “mat” and “minor” should be switched, or the corresponding in the function itself.

Took me a while to figure out why I kept getting Segmentation faults.

Jason,

Thanks a lot for pointing out the bug.

Hi Dear Janson,

Thanks for your code.

I am beginner in C++ and I don’t know how can i put your code inside of the main function, how to call MatrixInversion and how to input matrix element values.

Really, pointer-to-pointer is somehow ambiguous and I can’t understand it.

Thanks.

Sorry for confusing the main code writer.

I also want to thank chi3x10 due to main contribution.

I think your calculation of the determinant is possibly flawed. You’re calculation assumes that row and column indices in the range [0, n-1] where n is the size of the square matrix.

Generic expansion of a determinant at a row or column calls for row and column indices in the range [1, n], which means your summation should in fact be:

det += pow( -1.0, 2.0 + i ) * mat[0][i] * CalcDeterminant( minor,order-1 );

since you are expanding according to the first row of the matrix.

http://local.wasp.uwa.edu.au/~pbourke/miscellaneous/determinant/

has a very good explanation of the theory and a C implementation

Sorry, my bad…just realized that squaring the square of negative one still produces the same answer. There is nothing wrong with your determinant calculation.

🙂

Thanks very much for posting – this is very helpful!

chi3x10-

Thanks for posting this! I’m desperately trying to get my experiments under way and I’ve got enough to do without figuring this out!

One question: my compiler is setup strictly for C, no C++; I’ve translated the code most of the way (I think) but I’m sketchy on pointers. Can you direct me to anything like a C++ to C translation guide, or better yet post a C version? That’s asking a lot…

In any case, thank you!

Greg,

All you need to do is to modify memory allocation thing. Use malloc(…..) instead of new. Use free(…) instead of delete. It should work fine.

Thanks for the quick reply!

That’s what I thought (and what I did) but …

Ok, I just went back and found a bug in my translation, so now I’m not getting warnings, but I’m still segfaulting when I try to call the function.

To be clear, instead of (for example)

float *temp = new float[(order-1)*(order-1)];

float **minor = new float*[order-1];

I should use:

float *temp = malloc(sizeof(float)*((order-1)*(order-1)));

float **minor = malloc(sizeof(float)*(order-1));

Perhaps I’m having a problem initializing MyMatrix (as the above commenter).

Can I do direct assignment, i.e. (for 3×3 test case):

float **MyMatrix=malloc(sizeof(float)*3*3);

float **MyInverse=malloc(sizeof(float)*3*3);

MyMatrix[0][0]=1.2;

MyMatrix[0][1]=0.95;

//etc,etc

and then make the call as

MatrixInversion(MyMatrix, 3, MyInverse);

Thanks for your help!

Hello,

Thanks for sharing your code.

I am having a strange problem, the result being returned in Y is basically of type -1.#IND and 1.#QNAN.

I am using your code exactly as it is, and calling it as follows for a 4×4 matrix. Do you have any idea why? Thanks so much!

float **MyMatrix;

MyMatrix = new float*[4];

for(int i=0;i<4;i++)

MyMatrix[i] = new float[4];

float **Y;

Y = new float*[4];

for(int i=0;i<4;i++)

Y[i] = new float[4];

for(int i=0;i<4;i++)

for(int j=0;j<4;j++)

MyMatrix[i][j] = i * 1.0;

cout<<endl;

for (int i=0;i<4;i++)

for (int j=0;j<4;j++)

cout<<MyMatrix[i][j]<<” \n”;

MatrixInversion(MyMatrix, 4, Y);

for (int i=0;i<4;i++)

for (int j=0;j<4;j++)

cout<<Y[i][j]<<” \n”;

Tom,

Maybe the problem is that your matrix is singular. You can check if a matrix is singular by calculating the determinant of the matrix first. If determinant of a matrix is zero, the matrix is invertible.

Hope this helps

why do I get an error GetMinor undefined and CalcDeterminant undefined. trying to establish matrix determinant (n less than =to 100)

Hello! I m trying to use your code, but I get one error:

Error: Function MatrixInversion(R,3,Rinv) is not defined in current scope

This i how I declared the matrices:

double R[3][3N];

double Rinv[3][3];

And this is how I called the function:

MatrixInversion(R,3,Rinv);

I just did copy paste of your code before the main.

Thanks for the help!

hey, can anyone post here your full source code? tnx..

Hello… i used your code in my work… and it seemed to work fine.. thanks a lot…

Hi Jason – great work 🙂 … so is it ok to incorporate your code into my GPL’d FOSS project?

sure! but I hope you knows that this method would become really slow when the dimensionality of the matrix is high (like 10×10).

Hello chi3x10,

Thank you so much for your nice code. It works like a gem! I like the way you implemented the dynamically allocated variables. For many years I used the “minv” of the 40 year old Scientific Subroutine Package of IBM, which I recoded in C++. Oh well, it was time for a change!

Thanks a million!!!!

thanks for the code.

in line 80, instead of pow( -1.0, i ) you could have used (i%2==1?-1.0:1.0) which is pure c code.

Thanks for the tip 🙂

thanks !

very good implementation

Your code helped a lot in my university work. Thanks a lot!

Quick question:

I’ve gotten this code to work properly for me, but I’m having trouble using the result in the matrix Y.

For example, if I use:

for (int i=0;i<4;i++){

for (int j=0;j<4;j++){

cout << Y[i][j] << " ";

}

printf("\n");

}

the matrix comes out fine. But, if I use printf in place of cout I get jibberish numbers.

Likewise, if I declare the following assignment:

float x = Y[0][0]

and print x, I get jibberish again. Yet, doing

cout << Y[0][0];

works just fine! What am I doing wrong?

Thanks in advance.

Hi Jason,

I have found your code in google search. Thank you very much for publishing this code. It really helped me.

Since I work with large matrices, I was wondering whether you have / know another code such that is efficient for larger matrices (may be not analytic).

Thanks…

Hi, did you get an answer on this? I hate to look a gift horse in the mouth as they say, but I find this code extremely slow. A 10×10 matrix is done in about 5 seconds, but an 11×11 matrix takes over a minute. Seems very strange. Maybe I’ve done something wrong…

Hi Jason, again thanks for your excellent code! I ran out of memory when I used the routine recursively and then I realized that the allocated memory was not released correctly. So I changed it a bit as follows and again this works perfectly. Thanks a million! Art…

// release memory

for (i = 0; i < n; ++i) delete [] a[i];

delete [] a;

for (i = 0; i < n; ++i) delete [] b[i];

delete [] b;

for (i = 0; i < n; ++i) delete [] x[i];

delete [] x;

delete [] indx;

delete [] c;

Thanks Arthur.

hi jason, i am grateful for your code. You are a bomb. I am a beginner and i need more codes to improve in c++. Thanks alot.

Hi.. I am getting an error that both functions CalcDeterminant and GetMinor should have a prototype even when i have include a number of header files. Can you please help.

I am also looking for a matrix modulus code in c++ for my project.. It would be very generous of you if you could tell me where i can get it..

thank for sharing this inverse matix code 🙂

Birkenstock 37…[…]Calculation of matrix inverse in C/C++ « Jason Yu-Tseh Chi’s Notes[…]…

tv with wifi…[…]Calculation of matrix inverse in C/C++ « Jason Yu-Tseh Chi’s Notes[…]…

[…] the way of directly calculating an inverse matrix may be a way out. I modified the C++ code from https://chi3x10.wordpress.com/2008/05/28/calculate-matrix-inversion-in-c/ and share it here with those who are in the middle of C++ coding using boost library. Caution: do […]

When I run I get Declaration syntax error pointing to the main function, why?

Hello I am a beginner, can I see a sample main function plz

I have two dimensional array.

float C[50][50];

What I can enter in function

CalcDeterminant( float **mat, int order); ?

If

CalcDeterminant( C, 76130)

I got

error: cannot convert ‘float (*)[76130]’ to ‘float**’ for argument ‘1’ to ‘double CalcDeterminant(float**, int)’

If

CalcDeterminant( *C, 76130)

I got

error: cannot convert ‘float*’ to ‘float**’ for argument ‘1’ to ‘double CalcDeterminant(float**, int)’

If

CalcDeterminant( **C, 76130)

I got

error: cannot convert ‘float’ to ‘float**’ for argument ‘1’ to ‘double CalcDeterminant(float**, int)’

Could you simply put a donwloadable file, such that I do not need to remove line number after copy/paste?

@boy

Move the cursor above the source code area. A four button toolbox will appear. Click on the view source button.

plz help me in this code …

Write a C++ program that should:

Add, subtract, and multiply two matrices.

Multiply and divide matrices by scalar.

Transpose matrices.

Inverse matrices.

Calculate matrix determinant.

Compare two matrices.

Set a row or a column of the matrix to a certain value.

Calculate performance (time and memory taken by each algorithm).

please i want any one to help me in this:

Set a row or a column of the matrix to a certain value

and this too in c++:

Calculate performance (time and memory taken by each algorithm).

thank you in advance.

Sorry but no. 😦

How about a no? 😀

Undeniably believe that which you stated. Your favorite justification appeared to be on the web

the easiest thing to be aware of. I say to you, I definitely get annoyed

while people consider worries that they just do not know about.

You managed to hit the nail upon the top and also defined out the whole thing without having

side effect , people could take a signal. Will likely

be back to get more. Thanks

Just wish to say your article is as amazing.

The clarity to your submit is simply great and i can assume you’re an expert in this subject.

Well along with your permission let me to clutch your RSS feed to stay updated with impending

post. Thanks one million and please keep up the

rewarding work.

#include

#include

using namespace std;

int main ()

{

// part 1 : Storing of matrix

int matrix [3][3];

for (int i=0; i<3 ; i++)

{

cout <<"Enter row "<<i+1<< " of the 3 X 3 matrix "<<endl;

for(int j=0; j> matrix [i][j];

}

}

// part 2 : printing of entered matrix

cout <<"\n\nThe matrix you entered is ";

for (int t=0; t<3 ; t++)

{

cout <<endl;

for(int r=0; r<3 ; r++)

{

cout << matrix [t][r]<<"\t";

}

}

// part 3 : Checking singularity of the matrix

int g,h(0);

for (int t=0; t<3 ; t++)

{

g = (matrix [0][t])*((matrix[1][(t+4)%3]*matrix[2][(t+5)%3])-(matrix[1][(t+5)%3]*matrix[2][(t+4)%3]));

h=h+g;

}

if (h==0)

{

cout << "\nThe Entered matrix is a Singular \nIt does not have an INVERSE :'(";

}

else

{

// part 4 : calculating additive matrix

cout << " \n\nThe determinant of the Matrix is "<< h;

int inmatrix [3][3];

int p=0;

for (int i=0; i<3 ; i++)

{

for(int j=0; j<3 ; j++)

{

if (i==0)

{

p= (pow(-1,j))*((matrix[1][(j+4)%3]*matrix[2][(j+5)%3])-(matrix[1][(j+5)%3]*matrix[2][(j+4)%3]));

inmatrix [i][j]=p;

}

else if (i==1)

{

p= (pow(-1,j+1))*((matrix[0][(j+4)%3]*matrix[2][(j+5)%3])-(matrix[0][(j+5)%3]*matrix[2][(j+4)%3]));

inmatrix [i][j]=p;

}

else if (i==2)

{

p= (pow(-1,j))*((matrix[0][(j+4)%3]*matrix[1][(j+5)%3])-(matrix[0][(j+5)%3]*matrix[1][(j+4)%3]));

inmatrix [i][j]=p;

}

}

}

cout << "\n\nThe inverse oh the matrix is : \n\n"<<"(1/"<<h<<")";

for (int t=0; t<3 ; t++)

{

for(int r=0; r<3 ; r++)

{

cout << "\t"<<inmatrix [r][t]<<"\t";

}

cout <<endl;

}

}

return 0;

}

There are syntax errors in your code…

Running the program generates a run-time error. This error is error-related. I think a good way to solve this problem is to don’t use double pointer constructions, but only 1d arrays and use the order to index them as 2d arrays.

Edit: Problem is memory related

I really like what you guys are usually up too.

Such clever work and reporting! Keep up the excellent works guys I’ve added you guys to

my personal blogroll.

Hello, can anyone tell me how to call this function? I need to store the results in a matrix called “inverse”

Let me be a bit more clear…I need to call this function from my main function…and the matrix that I am wanting to get is called inverse…double **inverse to be exact.

i don’t understand what is order that you are passing to the functions?

Hello guys! I tried running the program but got the errors that ‘CalcDeterminant’ and ‘GetMinor’ were not declared in the scope. Please, how can I declare them.

Thank you for this. Although, could you please explain what the fuck this means:

float *temp = new float[(order-1)*(order-1)];

float **minor = new float*[order-1];

for(int i=0;i<order-1;i++)

minor[i] = temp+(i*(order-1));

Hi , thank you for the wonderful implementation. Using your code, I am trying to find inverse of matrix of order 12 and it takes almost 5 minutes to complete even with O3 optimization. Is there any way to run it faster ?

This is not a practical program for inverting a matrix for two reasons: it is much to slow and loss of accuracy due to cancellation among many large terms in brute force calculation of all those determinants.

Gauss solved the problem of accurate numeric solution of simultaneous linear equations by pivoting on the largest remaining term on the main diagonal and swapping rows to make the triangle below the main diagonal all zeros. Then Jordan has a method of making the upper triangle zero also.

Inverting a matrix is equivalent to solving the set of N linear equations in N unknowns N times, once for each basis element of possible right hand sides. A systematic way of doing this is to put the NxN matrix in the left-hand side of an N by 2N two-dimensional array and put the unit NxN matrix in the right side.

I programmed this as a C++ template

template

class InvertMatrixTemplate

{

Hello! I m trying to use your code, but I get one error:

Error: Function MatrixInversion(R,3,Rinv) is not defined in current scope

This i how I declared the matrices:

double R[3][3N];

double Rinv[3][3];

And this is how I called the function:

MatrixInversion(R,3,Rinv);

I just did copy paste of your code before the main.

Thanks for the help!

Sorry, I meant

double R[3][3];

double Rinv[3][3];

This website was… how do I say it? Relevant!!

Finally I’ve found something that helped me.

Many thanks!

thanks a lot, this small program is very useful for me!

its bug’s

using this code

double beta = 180.0;

double Aa = 1.0 / beta;

float cov[4][4] = { { 1.0f / 8.0f, 1.0f / 16.0f, 1.0f / 32.0f, 1.0f / 64.0f },

{ 1.0,1.0,1.0,1.0 },

{ 3.0 * Aa,4.0 * Aa ,5.0 * Aa ,6.0 * Aa },

{ (6.0*Aa)*(6.0*Aa),(12.0*Aa)*(12.0*Aa),(20.0*Aa)*(20.0*Aa),(30.0*Aa)*(30.0*Aa) } };

float **cov1 = new float *[4];

for (int p = 0;p<4;p++)

cov1[p] = new float[4];

for (int i = 0; i<4; i++)

for (int j = 0; j<4; j++)

cov1[i][j] = cov[i][j];

for (int i = 0; i<4; i++)

{

for (int j = 0; j<4; j++)

cout << cov1[i][j] << ' ';

cout << endl;

}

float **res = new float *[4];

for (int p = 0;p<4;p++)

res[p] = new float[4];

MatrixInversion(cov1, 4, res);

for (int i = 0; i<4; i++)

{

for (int j = 0; j<4; j++)

cout << res[i][j]<< " ";

cout << endl;

}

for (int i = 0;i<4;i++)

delete[] cov1[i];

delete[] cov1;

for (int i = 0;i<4;i++)

delete[] res[i];

delete[] res;

getchar();

return 0;

in visual studio email and I will send proof

sorry my bad it works perfect!

here is the error I made

\float cov[4][4] = { { 1.0f / 8.0f, 1.0f / 16.0f, 1.0f / 32.0f, 1.0f / 64.0f },

{ 1.0,1.0,1.0,1.0 },

{ 3.0 * Aa,4.0 * Aa ,5.0 * Aa ,6.0 * Aa },

{ 6.0*pow(Aa,2),12.0*pow(Aa,2),20.0*pow(Aa,2),30.0*pow(Aa,2) } };