Calculation of matrix inverse in C/C++

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.
GetMinor:
Get minor matrix.

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

Advertisements

71 Responses to Calculation of matrix inverse in C/C++

  1. Fab says:

    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 ?

  2. chi3x10 says:

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

  3. Vessi says:

    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;

    }

  4. Jason says:

    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**’ …

  5. chi3x10 says:

    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!

    • Dennis Dunn says:

      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.

  6. Jason says:

    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.

  7. chi3x10 says:

    Jason,
    Thanks a lot for pointing out the bug.

  8. Hamid says:

    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.

  9. Hamid says:

    Sorry for confusing the main code writer.
    I also want to thank chi3x10 due to main contribution.

  10. Cixelsid says:

    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

  11. Cixelsid says:

    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.

  12. Shannon says:

    Thanks very much for posting – this is very helpful!

  13. Greg says:

    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!

  14. chi3x10 says:

    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.

  15. Greg says:

    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!

  16. Tom says:

    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”;

  17. chi3x10 says:

    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

  18. reed says:

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

    • Mario says:

      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!

  19. dcrow says:

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

  20. Naz says:

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

  21. chris says:

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

    • chi3x10 says:

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

  22. Arthur says:

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

  23. Ali says:

    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.

  24. omid says:

    thanks !

    very good implementation

  25. Andres says:

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

  26. Jack says:

    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.

  27. Meital says:

    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…

    • Ken says:

      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…

  28. Arthur Grunwald says:

    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;

  29. aguheva thompson junior says:

    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.

  30. anshul says:

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

  31. toto says:

    thank for sharing this inverse matix code 🙂

  32. Birkenstock 37…

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

  33. tv with wifi says:

    tv with wifi…

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

  34. […] 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 […]

  35. jaytee says:

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

  36. jaytee says:

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

  37. 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)’

  38. boy says:

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

    • chi3x10 says:

      @boy

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

  39. shimaa subhy says:

    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).

  40. samar says:

    and this too in c++:
     Calculate performance (time and memory taken by each algorithm).
    thank you in advance.

  41. Visnyei says:

    How about a no? 😀

  42. 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

  43. youtube says:

    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.

  44. DoNNy says:

    #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;
    }

  45. Allard says:

    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.

  46. 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.

  47. Jawad Rehman says:

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

    • Jawad Rehman says:

      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.

  48. abdelrhman says:

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

  49. JVN says:

    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.

  50. Larrimus says:

    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));

  51. 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 ?

  52. Graham Frye says:

    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
    {

  53. Mario says:

    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!

  54. This website was… how do I say it? Relevant!!
    Finally I’ve found something that helped me.
    Many thanks!

  55. jim says:

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

  56. ollp says:

    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

    • ollp says:

      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) } };

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: