Matrix inversion

[Michael J. Lew]

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
**