Becareful about the missing parts in the program due to HTML syntax messing up with c code , use the source code download at the bottom of the page

Matrix Inversion is an essential part in mesh analysis.
The algorithm used is Shipley-Coleman
Theoretically, with this algorithm, you can find inverse matrix is with an indefinite size.
So much for words..
here is how the program works.
This program is written in PocketC PalmOS, but you can slightly play on it and
use it for Java, C, C++.
For a matrix A size 5x5 (nxn), let's assume the values will be like:
size_n=5
 

So in Palm you open a Memo with a header :"\matrix"
ex:
 

Listing1----------------Beginning-------------------------------- /matrix 5 0.34 0.54 -3.2 4.0 7.0 -90.123 0.0045 -5.689 43.12 10.23 -1.023 0.987 43.231 -112.34 23.12 -0.3 42.23 6.7 -3.12 1.23 12.56 32.98 4.100023 12.54 Listing1----------------End---------------------------------------

/matrix size_n a11 a12 a13 a14 a15 a21 a22 a23 a24 a25 a31 a32 a33 a34 a35 a41 a42 a43 a44 a45 a51 a52 a53 a54 a55

Algorithm:
For algorithm, I consulted my 2nd year Programming Notes by H. Yin (UMIST-EE)
algorithm is quite straight-forward. You walk on the diagonal, and everytime you move one diagonal
element forward , you carry out a set of operations.
Ok, here is the way to do it.
for a matrix n*n
you begin with l=1 and carry out till n
(!= : not equal to)
1) a'_LL=1/a_LL
2)a'_ML=a_MLa'_LL                            for all M!=L (Column elements)
3)a'_MP=a_MP-a'_MLa_LP                   for all M!=L ad P!=L (other elements)
4)a'_LP=-a'_LLa_LP                             for all P!=L (Row elements)

Seems a bit confusing? Let's look at the program

Program Cycle:
1) Look for a memo with a heading : "/matrix"
mmfind("/matrix");
2)Skip the first line
mmgetl();
3) Read the size value from the second line and put it in a integer var k
k=(int)mmgetl();
4) Until you reach the end of the memo, read the values and put them in a float matrix a
Note: a is defined in "inverse" file as :
float a[100];
(so you can deal matrices up to 10x10)
5) Call the subroutine invmatrix (k)
invmatrix(k)
6)close memo
mmclose();

here is the listing:
------------Listing--prog.c---------------Beginning---------------
//prog.c
include "inverse"

main()
{
int g,k;
g=0;
mmfind("/matrix2");
mmgetl();
puts("\n");
k=(int)mmgetl();
while(!mmeof()){
  a[g++]=(float)mmgetl();
 
}
invmatrix(k);
mmclose();
}
---------Listing--prog.c------------------End-----------------------

And here is the inverse subroutine

---------Listing--inverse-----------------Beginning------------------
/$inverse

float a[100];
 

invmatrix(int n){
int l,m,p,b;
int pos,p_old;
int p_old2;
for(l=0;l<n;l++){
  p_old2=l*n+l;
  a[p_old2]=1/a[p_old2];
  for(m=0;m<n;m++)
   if(m!=l){
     p_old=m*n+l;
     a[p_old]=a[p_old]*a[p_old2];
     }

for(m=0;m<n;m++){
 for(p=0;p<n;p++){
  if(m!=l){
  if(p!=l){
   pos=m*n+p;
   p_old=m*n+l;
   p_old2=l*n+p;
   a[pos]=a[pos]-(a[p_old]*a[p_old2]);
   }
  }
 }
}

 p_old=l*n+l;
for(p=0;p<n;p++)
 if(p!=l){
   pos=l*n+p;
   a[pos]=-a[p_old]*a[pos];
 }
}
pos=n*n;
for(b=0;b<pos;b++)
 puts(format(a[b],7)+"\n");
}
--------------------Listing--inverse--------End---------------------------------------------

invmatrix file:
first defines a float array a[100]
because PocketC doesn't have a[n][n] notation for floats, operations done
in a single dimension array
a_MP=m*n + p
 matrix A size n*n
Afterwards the algorithm carrys on.

Source Code - prog.c

Source Code - inverse.h

 If you have any problems, drop me a line

sanliuk@yahoo.co.uk