hi all , I found the followin interesting programe on one of the sites that I browsed, it is used to invert a matrix using two different methods, I have a matrix
given by
double a[3][3]={{1.,2.,6.},{3 .,4., 5.},{3.,4.,7.}} ;
and I am trying to find it is inverse using the codes shown below, I did know how to fit the my matrix shown above into the programs , what I understood from the code comments is that a is the matrix that is to be inverted and x is where the result will be stored and n is the order , in my case it is 3, but I was not able to figure out the index, anybody knows what it refers to?
void migs (a,n,x,indx)
given by
double a[3][3]={{1.,2.,6.},{3 .,4., 5.},{3.,4.,7.}} ;
and I am trying to find it is inverse using the codes shown below, I did know how to fit the my matrix shown above into the programs , what I understood from the code comments is that a is the matrix that is to be inverted and x is where the result will be stored and n is the order , in my case it is 3, but I was not able to figure out the index, anybody knows what it refers to?
void migs (a,n,x,indx)
Code:
/* Updated 10/24/2001. */
/************************* Program 4.4 ****************************/
/* */
/************************************************************************/
/* Please Note: */
/* */
/* (1) This computer program is written by Tao Pang in conjunction with */
/* his book, "An Introduction to Computational Physics," published */
/* by Cambridge University Press in 1997. */
/* */
/* (2) No warranties, express or implied, are made for this program. */
/* */
/************************************************************************/
#include <stdio.h>
#include <math.h>
#define NMAX 100
void migs (a,n,x,indx)
/* Function to invert matrix a[][] with the inverse stored
in x[][] in the output. Copyright (c) Tao Pang 2001. */
int n;
int indx[];
double a[NMAX][NMAX],x[NMAX][NMAX];
{
int i,j,k;
double b[NMAX][NMAX];
void elgs();
if (n > NMAX)
{
printf("The matrix dimension is too large.\n");
exit(1);
}
for (i = 0; i < n; ++i)
{
for (j = 0; j < n; ++j)
{
b[i][j] = 0;
}
}
for (i = 0; i < n; ++i)
{
b[i][i] = 1;
}
elgs (a,n,indx);
for (i = 0; i < n-1; ++i)
{
for (j = i+1; j < n; ++j)
{
for (k = 0; k < n; ++k)
{
b[indx[j]][k] = b[indx[j]][k]-a[indx[j]][i]*b[indx[i]][k];
}
}
}
for (i = 0; i < n; ++i)
{
x[n-1][i] = b[indx[n-1]][i]/a[indx[n-1]][n-1];
for (j = n-2; j >= 0; j = j-1)
{
x[j][i] = b[indx[j]][i];
for (k = j+1; k < n; ++k)
{
x[j][i] = x[j][i]-a[indx[j]][k]*x[k][i];
}
x[j][i] = x[j][i]/a[indx[j]][j];
}
}
}
void elgs (a,n,indx)
/* Function to perform the partial-pivoting Gaussian elimination.
a[][] is the original matrix in the input and transformed
matrix plus the pivoting element ratios below the diagonal
in the output. indx[] records the pivoting order.
Copyright (c) Tao Pang 2001. */
int n;
int indx[];
double a[NMAX][NMAX];
{
int i, j, k, itmp;
double c1, pi, pi1, pj;
double c[NMAX];
if (n > NMAX)
{
printf("The matrix dimension is too large.\n");
exit(1);
}
/* Initialize the index */
for (i = 0; i < n; ++i)
{
indx[i] = i;
}
/* Find the rescaling factors, one from each row */
for (i = 0; i < n; ++i)
{
c1 = 0;
for (j = 0; j < n; ++j)
{
if (fabs(a[i][j]) > c1) c1 = fabs(a[i][j]);
}
c[i] = c1;
}
/* Search the pivoting (largest) element from each column */
for (j = 0; j < n-1; ++j)
{
pi1 = 0;
for (i = j; i < n; ++i)
{
pi = fabs(a[indx[i]][j])/c[indx[i]];
if (pi > pi1)
{
pi1 = pi;
k = i;
}
}
/* Interchange the rows via indx[] to record pivoting order */
itmp = indx[j];
indx[j] = indx[k];
indx[k] = itmp;
for (i = j+1; i < n; ++i)
{
pj = a[indx[i]][j]/a[indx[j]][j];
/* Record pivoting ratios below the diagonal */
a[indx[i]][j] = pj;
/* Modify other elements accordingly */
for (k = j+1; k < n; ++k)
{
a[indx[i]][k] = a[indx[i]][k]-pj*a[indx[j]][k];
}
}
}
}
Comment