-
Matrix inversion
I have looked up matrix inversions for the past two days and still cannot properly code a program that can compute the inverse of any square matrix. I found code online that works, but only if the first numbers in each row are 1. For any other numbers, the calculation is almost correct. The only difference between the code I found and the correct results is the multiplication of each element in every row after the first by some constant.
The code I have thus far:
Code:
#include<stdio.h>
#include<stdlib.h>
int main(void)
{
int n=3;
double A[3][3]= {{1,-1,3},{2,1,2},{-2,-2,1}};
double B[n][2*n];
int i,j,k;
double M;
/*-------------Create Matrix w/identity---------------*/
for( i=0;i<n;i++)
{
for( j=n;j<2*n;j++)
{
if(j==(i+n))
{
B[i][j]=1.0;
}
else
{
B[i][j]=0.0;
}
}
}
for(i=0; i<n; i++)
{
for(j=0; j<n; j++)
{
B[i][j]=A[i][j];
}
}
/*-----------Calculate inverse----------*/
for(k=0;k<n;k++)
{
for(i=0;i<n;i++)
{
if(i!=k)
{
if(B[n-1][n-1]== 0)
{
printf("This result requires a non-singular matrix");
exit(0);
}
else if (B[k][k]==0)
k++;
M=(B[i][k]/B[k][k]);
for(j=0;j<2*n;j++)
{
B[i][j]=(B[i][j] - B[k][j]*M);
}
}
}
}
for(i=0;i<n;i++)
{
for(j=0;j<2*n;j++)
{
B[i][j]=(B[i][j]/B[i][i]);
}
}
/*-------------Display matrix--------------------*/
for(i=0;i<n;i++)
{
printf("\n");
for(j=n;j<2*n;j++)
{
printf("%.4lf\t",B[i][j]);
}
}
return 0;
}
I also tried looking for algorithms explaining the calculation, but they did not work either. I have been very frustrated over the past few days. Any help would be greatly appreciated. Thanks in advance.
-
Had you managed to look it up on Wikipedia, it would have pointed you to here, which conveniently is the algorithm you're trying to do. (If you had looked it up somewhere else, you would still have gotten the row-reduction algorithm pointed to above, just maybe on a different server somewhere, or perhaps in a book.)
You're doing the down-sweep and the up-sweep at the same time, but I think that works.
The problem would appear to be the last division -- when j=i, you're going to change B[i][i] to 1.0 and therefore fail to divide the rest of the row. Save the value before the for-loop.
-
Thank you so much. That small change made all the difference.
Thanks again.