Matrix Inversion Subroutine

[Vipro]

Hello:
I have found this matrix inversion subroutine on MS Fortran and I need to
use it but I don't have clear what the variables represent. Please somebody
can help me.

Thanks

Rene Suarez. Please reply to (rene...@hotmail.com)

SUBROUTINE gaussj(a,n,np,b,m,mp)
INTEGER m,mp,n,np,NMAX
REAL a(np,np),b(np,mp)
PARAMETER (NMAX=50)
INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),indxr(NMAX),ipiv(NMAX)
REAL big,dum,pivinv
do 11 j=1,n
ipiv(j)=0
11 continue
do 22 i=1,n
big=0.
do 13 j=1,n
if(ipiv(j).ne.1)then
do 12 k=1,n
if (ipiv(k).eq.0) then
if (abs(a(j,k)).ge.big)then
big=abs(a(j,k))
irow=j
icol=k
endif

else if (ipiv(k).gt.1) then
pause 'singular matrix in gaussj'
endif
12 continue
endif
13 continue
ipiv(icol)=ipiv(icol)+1
if (irow.ne.icol) then
do 14 l=1,n
dum=a(irow,l)
a(irow,l)=a(icol,l)
a(icol,l)=dum
14 continue
do 15 l=1,m
dum=b(irow,l)
b(irow,l)=b(icol,l)
b(icol,l)=dum
15 continue
endif
indxr(i)=irow
indxc(i)=icol
if (a(icol,icol).eq.0.) pause 'singular matrix in gaussj'

pivinv=1./a(icol,icol)
a(icol,icol)=1.
do 16 l=1,n
a(icol,l)=a(icol,l)*pivinv
16 continue
do 17 l=1,m
b(icol,l)=b(icol,l)*pivinv
17 continue
do 21 ll=1,n
if(ll.ne.icol)then
dum=a(ll,icol)
a(ll,icol)=0.
do 18 l=1,n
a(ll,l)=a(ll,l)-a(icol,l)*dum
18 continue
do 19 l=1,m
b(ll,l)=b(ll,l)-b(icol,l)*dum
19 continue
endif
21 continue
22 continue
do 24 l=n,1,-1
if(indxr(l).ne.indxc(l))then

do 23 k=1,n
dum=a(k,indxr(l))
a(k,indxr(l))=a(k,indxc(l))
a(k,indxc(l))=dum
23 continue
endif
24 continue
return
END

[Jerry A Green]

Hello Rene:

Let's give it a go. First of all I am going to refer to all of the
variables by their uppercase names.

Then the primary variables that you need to be concerned with are
A, N, NP, B, M, and MP

A is the input matrix. The max dimension of A is NP by NP (This is also
how A should be dimensioned in the calling program).

N is the number of rows (and columns) of the matrix that you want to
Invert (note that N must be <= to NP)

B is the output matrix. The max dimension of B is MP by MP (This is also
how B should be dimensioned in the calling program).

M is the number of rows (and columns) of the Inverted matrix.

Jerry . . .

--
Jerry A Green mailto:j...@cs-software.com
Custom Solutions http://www.cs-software.com/
209 Bayberry Run
Summerville, SC 29485-8778 Your source for discounted
Voice: (843) 871 9081 Fortran compilers and
Fax: (843) 873 8626 related software

[Bill Leen]

This is from the book "Numerical Recipes" by Press, Flannery et al.
On page 28, the description of input and output is:
"A is an input matrix of N by N elements, stored in an array of
physical dimensions NP by NP. B is is an input matrix of N by M
containing the M right-hand side vectors, stored in an array of physical
dimensions NP by MP. On output, A is replaced by its matrix inverse,
and B is replaced by the corresponding set of solution vectors."

Also see http://cfata2.harvard.edu/nr/ , the Numerical Recipes home
page.

Bill Leen