0

I wrote a inversion function for an n*n square matrix.

void inverseMatrix(int n, float **matrix)
{
    float ratio,a;
    int i, j, k;

   for(i = 0; i < n; i++)
   {
      for(j = n; j < 2*n; j++)
      {
         if(i==(j-n))
            matrix[i][j] = 1.0;
         else
            matrix[i][j] = 0.0;
      }
  }

  for(i = 0; i < n; i++)
  {
      for(j = 0; j < n; j++)
      {
          if(i!=j)
          {
              ratio = matrix[j][i]/matrix[i][i];
              for(k = 0; k < 2*n; k++)
              {
                  matrix[j][k] -= ratio * matrix[i][k];
              }
          }
      }
  }

  for(i = 0; i < n; i++)
  {
      a = matrix[i][i];
      for(j = 0; j < 2*n; j++)
      {
          matrix[i][j] /= a;
      }
  }

//return matrix;
}

This works fine in almost all cases but is failing in some cases like the ones shown here:

1 1 1 0       1 1 1 0
1 1 2 0       1 1 2 0  
1 2 0 1       1 2 1 0
1 2 0 2       1 2 0 2

What could be the case I am overlooking?

Thanks!

re3el
  • 735
  • 2
  • 12
  • 28
  • 3
    Can you be more specific, when you say, "fail," do you mean segfaults, or just returns the wrong result? – iluvcapra Aug 17 '15 at 05:32
  • 1
    You're not handling the divide by zero case whenever you have a zero on the diagonal. – Paul R Aug 17 '15 at 05:32
  • @PaulR: How can I solve that? Thanks – re3el Aug 17 '15 at 05:33
  • 2
    You probably need a more robust algorithm. – Paul R Aug 17 '15 at 05:39
  • @iluvcapra: No problems with compilation. I am getting garbage values for some matrices. As mentioned by PaulR, I overlooked the divide by zero case. Any help for solving it would be thankful. – re3el Aug 17 '15 at 05:52
  • Why do you reinvent the wheel? http://arma.sourceforge.net/ http://eigen.tuxfamily.org/index.php?title=Main_Page – user1436187 Aug 17 '15 at 05:57
  • @user1436187: I don't want to use library functions as I have some other defined data types for matrices. I need the code to modify it as per my requirements. – re3el Aug 17 '15 at 05:59
  • @re3el you need to check if the determinant of the matrix is zero or not in first place for confirming if the inverse of matrix exists.and you will have to use recursion if you want to calculate determinant. – nikhil mehta Aug 17 '15 at 06:00
  • @nikhilmehta determinants are usually computed by using row operations (which have a predictable effect on the determinant) to transform into a form from which the determinant is trivial to calculate. You seem to be thinking in terms of cofactor expansion -- but that gets prohibitively expensive (in general) as the dimension increases. – John Coleman Aug 17 '15 at 16:51

2 Answers2

0

see http://www.sourcecodesworld.com/source/show.asp?ScriptID=1086.

uses the Gauss Jordan algorithm

#include<stdio.h>
#include<stdlib.h>

int main()
{
    float **A,**I,temp;
    int i,j,k,matsize;

    printf("Enter the size of the matrix(i.e. value of 'n' as size is
nXn):");
    scanf("%d",&matsize);

    A=(float **)malloc(matsize*sizeof(float *));            //allocate memory
dynamically for matrix A(matsize X matsize)
    for(i=0;i<matsize;i++)
        A[i]=(float *)malloc(matsize*sizeof(float));

    I=(float **)malloc(matsize*sizeof(float *));            //memory allocation for
indentity matrix I(matsize X matsize)
    for(i=0;i<matsize;i++)
        I[i]=(float *)malloc(matsize*sizeof(float));

    printf("Enter the matrix: ");                           // ask the user for matrix A
    for(i=0;i<matsize;i++)
        for(j=0;j<matsize;j++)
            scanf("%f",&A[i][j]);

    for(i=0;i<matsize;i++)                                  //automatically initialize the unit matrix, e.g.
        for(j=0;j<matsize;j++)                              //  -       -
            if(i==j)                                        // | 1  0  0 |
                I[i][j]=1;                                  // | 0  1  0 |
            else                                            // | 0  0  1 |
                I[i][j]=0;                                  //  -       -
/*---------------LoGiC starts here------------------*/      //procedure // to make the matrix A to unit matrix

    for(k=0;k<matsize;k++)                                  //by some row operations,and the same row operations of
    {                                                       //Unit mat. I gives the inverse of matrix A
        temp=A[k][k];                   //'temp'  
        // stores the A[k][k] value so that A[k][k]  will not change
        for(j=0;j<matsize;j++)      //during the operation //A[i] //[j]/=A[k][k]  when i=j=k
        {
            A[k][j]/=temp;                                  //it performs // the following row operations to make A to unit matrix
            I[k][j]/=temp;                                  //R0=R0/A[0][0],similarly for I also
R0=R0/A[0][0]
        }                                                   //R1=R1-R0*A[1][0] similarly for I
        for(i=0;i<matsize;i++)                              //R2=R2-R0*A[2][0]      ,,
        {
            temp=A[i][k];                       //R1=R1/A[1][1]
            for(j=0;j<matsize;j++)             //R0=R0-R1*A[0][1]
            {                                   //R2=R2-R1*A[2][1]
                if(i==k)
                    break;                      //R2=R2/A[2][2]
                A[i][j]-=A[k][j]*temp;          //R0=R0-R2*A[0][2]
                I[i][j]-=I[k][j]*temp;          //R1=R1-R2*A[1][2]
            }
        }
    }
/*---------------LoGiC ends here--------------------*/
    printf("The inverse of the matrix is: ");               //Print the //matrix I that now contains the inverse of mat. A
    for(i=0;i<matsize;i++)
    {
        for(j=0;j<matsize;j++)
            printf("%f  ",I[i][j]);
        printf(" ");
    }
    return 0;
}
ralf htp
  • 9,149
  • 4
  • 22
  • 34
  • 1
    Please don't encourage the bad habit of [casting the result of malloc in C](http://stackoverflow.com/questions/605845/do-i-cast-the-result-of-malloc). – Paul R Aug 18 '15 at 16:49
  • Encouraging fragmented pointer-to-pointer lookup tables is far worse practice that casting the result of malloc. Prevent the whole program from being slow, burdensome and unsafe vs sate some purists on an internet site... – Lundin Sep 22 '15 at 11:10
  • A good example of algorithm with such an ugly code... – Peter May 06 '23 at 23:50
0

Diagonal elems have to be scaled to 1s first before zeroing lower triangle elems (your second nested loop). It turns out that diagonal contains 0 => no inverse exists OR we get a row echelon form. With a reversed iteration on diagonal we can reduce it and get the inverse. The code is far from optimal. When you have FPU, double could be a better choice than float for such numerical computations.

Note that zeroing submatrices, row swaps etc. could be replaced by far more optimal solutions. Matrix is a custom type and IsFloat0 is a custom function but all should be clear of naming and context. Enjoy code:

const uint sz = 4;
Matrix< double > mx;
mx.Resize( 2 * sz, sz );
mx.Zero();
for( uint rdx = 0; rdx < mx.NumRow(); ++rdx )
{
        mx( rdx, rdx + mx.NumRow() ) = 1.0; // eye
}
mx( 0, 0 ) = 1.0; mx( 0, 1 ) = 1.0; mx( 0, 2 ) = 1.0; mx( 0, 3 ) = 0.0;
mx( 1, 0 ) = 1.0; mx( 1, 1 ) = 1.0; mx( 1, 2 ) = 2.0; mx( 1, 3 ) = 0.0;
mx( 2, 0 ) = 1.0; mx( 2, 1 ) = 2.0; mx( 2, 2 ) = 0.0; mx( 2, 3 ) = 1.0;
mx( 3, 0 ) = 1.0; mx( 3, 1 ) = 2.0; mx( 3, 2 ) = 0.0; mx( 3, 3 ) = 2.0;

// pivot iteration
uint idx;
for( idx = 0; idx < sz; ++idx )
{
    // search for non-0 pivot
    uint sdx = sz;
    for( uint rdx = idx; rdx < sz; ++rdx )
    {
        if( !Util::IsFloat0( mx( rdx, idx ) ) ) { sdx = rdx; rdx = sz - 1; }
    }
    if( sdx < sz )
    {
        // swap rows
        if( idx != sdx )
        {
            for( uint cdx = 0; cdx < ( sz << 1 ); ++cdx )
            {
                double swp;
                swp = mx( idx, cdx );
                mx( idx, cdx ) = mx( sdx, cdx );
                mx( sdx, cdx ) = swp;
            }
        }
        // 1 pivot and 0 col
        {
            double sc = 1.0 / mx( idx, idx );
            for( uint cdx = 0; cdx < ( sz << 1 ); ++cdx )
            {
                mx( idx, cdx ) *= sc; // 1
            }

            for( uint rdx = 1 + idx; rdx < sz; ++rdx )
            {
                double sd = mx( rdx, idx );
                for( uint cdx = 0; cdx < ( sz << 1 ); ++cdx )
                {
                    mx( rdx, cdx ) -= sd * mx( idx, cdx ); // 0
                }
            }
        }
    }
    else { idx = sz; }
}
if( sz < idx ) { mx.Zero(); }
else
{
    for( idx = 0; idx < sz; ++idx )
    {
        uint ydx = sz - 1 - idx;
        for( uint rdx = 0; rdx < ydx; ++rdx )
        {
            double sc = mx( rdx, ydx );
            for( uint cdx = 0; cdx < ( sz << 1 ); ++cdx )
            {
                mx( rdx, cdx ) -= sc * mx( ydx, cdx ); // 0
            }
        }
    }
}
renonsz
  • 571
  • 1
  • 4
  • 17