Matrix inversion

[Jim Hurley]

Michael,

Thank you for these routines. I think if I can get them to work I 
can find the matrix inverse. My source for numerical work is : 
Numerical analysis, by Richard Burden. It gives only  algorithms and 
not program routines.

Unfortunately I am having trouble with your handlers. I tried a 
simple set of linear equations:

x + 2y = 1
x + y = 0

With the solution: x = -1 and y = 1.

But the routine below, LUBKSB, yields  x = -2 and y = 1.

I  used

aMatrix = 1,2
                 1,1

  with  n = 2

and

rhsMatrix =1
                   0

Any ideas?

Jim


At 2:40 PM +1000 3/31/03, Michael J. Lew wrote:
>Here are a couple of matrix handlers that I translated from Fortran 
>routines published in the excellent Numerical Recipes book by Press 
>et al. Not enough entirely to invert matrices (that wasn't what I 
>needed them for), but the small additional routine is fully 
>documented in a downloadable pdf from 
><http://lib-www.lanl.gov/numerical/bookfpdf.html>.
>
>For anyone who does numerical work I highly recommend the Numerical 
>Recipes books.
>
>#--****** LUBKSB ******
># Function Revolutionised from Numerical Recipes by MJL
># --Input is aMatrix, an n by n square matrix and rhsMatrix, an n by 1 vector
># also uses the global indx which is built by LUDCMP
># Returns an n by 1 vector solution
># For use with LUDCMP to solve linear equations or invert a matrix
># Note that the input and output parameters are not fully equivalent 
>to the originals
>
>function LUBKSB aMatrix,n,rhsMatrix
>   global indx
>   put aMatrix into a
>   put rhsMatrix into b
>   set the numberformat to "0.#########################"
>   put 0 into ii
>   repeat with i =1 to n
>      put indx[i] into ip
>     put b[ip] into mySum
>     put b[i] into b[ip]
>      if ii<>0 then
>       repeat with j = ii to i-1
>         put  mySum-a[i,j]*b[j] into mySum
>        end repeat --j
>     else
>       if mySum <> 0 then put i into ii
>      end if
>     put mySum into b[i]
>    end repeat --i
>   repeat with i = n down to 1
>     put b[i] into mySum
>     if i<n then
>       repeat with j = i+1 to n
>         put mySum-a[i,j]*b[j] into mySum
>       end repeat --j
>      end if
>     put mySum/a[i,i] into b[i]
>    end repeat --i
>   return b
>end LUBKSB
>
>
>
>#--***** LUDCMP *****
># Function Revolutionised by MJL from Numerical Recipes subroutine ludcmp
># input is an n by n square matrix, aMatrix
># returns aMatrix modified to be an LU decomposed version of itself
># also returns (a global) indx which is a vector of the row perms 
>affected by partial pivoting (!)
># For use in conjunction with LUBKSB to solve linear equations or 
>invert a matrix
># Note that the input and output parameters are not fully equivalent 
>to the originals
>
>function LUDCMP amatrix,n
>   global indx
>   set the numberformat to "0.#########################"
>   put 10^-22 into tiny
>   put aMatrix into a
>   put 1 into d
>   repeat with i = 1 to n
>     put 0 into big
>     repeat with j=1 to n
>       if abs(a[i,j]) >big then put a[i,j] into big
>     end repeat --j
>     if big =0 then
>       answer "Singular matrix in LUDCMP. Procedure failed"
>       exit to top
>     end if
>     put 1/big into vv[i]
>   end repeat --i
>   repeat with j=1 to n
>     repeat with i=1 to j-1
>       put a[i,j] into mySum
>       repeat with k = 1 to i-1
>         put mySum-a[i,k]*a[k,j] into mySum
>       end repeat --k
>       put mySum into a[i,j]
>     end repeat --i
>     put 0 into big
>     repeat with i=j to n
>       put a[i,j] into mySum
>       repeat with k=1 to j-1
>         put mySum-a[i,k]*a[k,j] into mySum
>       end repeat --k
>       put mySum into a[i,j]
>       put abs(mySum)*vv[i] into dum
>       if dum>big then
>         put dum into big
>         put i into imax
>       end if
>     end repeat --i
>     if j<>imax then
>       repeat with k=1 to n
>         put a[imax,k] into dum
>         put a[j,k] into a[imax,k]
>         put dum into a[j,k]
>       end repeat --k
>       put -d into d
>       put vv[j] into vv[imax]
>     end if
>     put imax into indx[j]
>     if a[j,j] =0 then put tiny into a[j,j]
>     if j<>n then
>       put 1/(a[j,j]) into dum
>       repeat with i=j+1 to n
>         put dum*a[i,j] into a[i,j]
>       end repeat --i
>     end if
>   end repeat --j
>   return a
>end LUDCMP
>
>
>
>--
>Michael J. Lew
>
>Senior Lecturer
>Department of Pharmacology
>The University of Melbourne
>Parkville 3010
>Victoria
>Australia
>
>Phone +613 8344 8304
>
>**
>New email address: michaell at unimelb.edu.au
>**