LUDecomposition (LUD)

LU decomposition for ordinary matrices.

)abb package    LUD     LUDecomposition

Sy  ==> Symbol
L   ==> List
V   ==> Vector
VD  ==> Vector D
MD  ==> Matrix D
FD  ==> Fraction D
MFD ==> Matrix FD
I   ==> Integer
NNI ==> NonNegativeInteger
B   ==> Boolean
OUT ==> OutputForm

ROWREC ==> Record(Indices:L C, Entries:L D)

iter ==> "iterated"::Sy
rand ==> "random"::Sy

++ Description:
++ \axiomType{LUDecomposition} contains procedures to solve linear systems of
++ equations or to compute inverses using a LU decomposition.

LUDecomposition(D:Field) : Cat == Def where

  Cat ==> with 

    LUDecomp : MD -> Record(LU:MD, Perm:V I, Pivots:L D)
      ++ \axiom{LUDecomp(A)} computes a LU decomposition of \axiom{A}
      ++ using the algorithm of Crout. \axiom{LU} contains both triangular
      ++ matrices; \axiom{Perm} is the permutation used for partial
      ++ pivoting and \axiom{Pivots} yields the used pivots.

    LUSolve : (MD, V I, V D) -> V D
      ++ \axiom{LUSolve(LU,Perm,B)} uses a previously computed LU 
      ++ decomposition to solve a linear system with right hand side
      ++ \axiom{B}. \axiom{LU} and \axiom{Perm} are as given by
      ++ \axiom{LUDecomp}.

    LUInverse : MD -> Record(Inv:MD, Pivots:L D)
      ++ \axiom{LUInverse(A)} computes the inverse of \axiom{A} using a LU
      ++ decomposition.

  Def ==> add

    LUDecomp(AA:MD):Record(LU:MD, Perm:V I, Pivots:L D) ==
      -- LU decomposition using Crout's algorithm with partial pivoting.
      A := copy AA
      minR := minRowIndex A; maxR := maxRowIndex A
      minC := minColIndex A; maxC := maxColIndex A
      maxR^=maxC => error "LU decomposition only of square matrices"
      PermV:V I := new((maxR-minR+1)::NNI,0)
      Pivs:L D := empty

      for j in minC..maxC repeat
        for i in minR..(j-1) repeat
          s := qelt(A,i,j)
          for k in minR..(i-1) repeat
            s := s - qelt(A,i,k)*qelt(A,k,j)
          qsetelt!(A,i,j,s)

        i0:I := -1
        for i in j..maxR repeat
          s := qelt(A,i,j)
          for k in minC..(j-1) repeat
            s := s - qelt(A,i,k)*qelt(A,k,j)
          qsetelt!(A,i,j,s)
          if not(zero? s) and i0<0 then
            i0 := i            -- first non-zero pivot

        i0<0 => error "singular matrix in LUDecomp"
        if j^=i0 then
          swapRows!(A,j,i0)
        qsetelt!(PermV,j,i0)
        Pivs := cons(qelt(A,j,j),Pivs)

        if j^=maxC then
          d := 1/qelt(A,j,j)
          for k in (j+1)..maxR repeat
            qsetelt!(A,k,j,d*qelt(A,k,j))

      [A,PermV,Pivs]

    LUSolve(LU:MD,Perm:V I,XX:V D):V D ==
      -- Solves LU decomposed linear system for right hand side XX
      X := copy XX
      minR := minRowIndex LU; maxR := maxRowIndex LU
      maxIndex(X)^=maxR => error "Wrong dimensions in LUSolve"
      ii:I := -1

      for i in minR..maxR repeat             -- forward substitution
        ip := qelt(Perm,i)
        s := qelt(X,ip)
        qsetelt!(X,ip,qelt(X,i))
        if ii>=0 then
          for j in ii..(i-1) repeat
            s := s - qelt(LU,i,j)*qelt(X,j)
        else
          if not zero? s then ii := i
        qsetelt!(X,i,s)

      for i in maxR..minR by -1 repeat       -- back substitution
        s := qelt(X,i)
        for j in (i+1)..maxR repeat
          s := s - qelt(LU,i,j)*qelt(X,j)
        qsetelt!(X,i,s/qelt(LU,i,i))

      X

    LUInverse(A:MD):Record(Inv:MD, Pivots:L D) ==
      -- Inversion via LU decomposition
      Alu := LUDecomp A
      n := ncols A
      res:MD := new(n,n,0)

      for i in minRowIndex(A)..maxRowIndex(A) repeat
        v:V D := new(n,0)
        qsetelt!(v,i,1)
        res := setColumn!(res,i,LUSolve(Alu.LU,Alu.Perm,v))

      [res,Alu.Pivots]

   Compiling FriCAS source code from file 
      /var/zope2/var/LatexWiki/1676865482560581358-25px001.spad using 
      old system compiler.
   LUD abbreviates package LUDecomposition 
   processing macro definition Sy ==> Symbol 

   processing macro definition L ==> List 

   processing macro definition V ==> Vector 

   processing macro definition VD ==> Vector D 

   processing macro definition MD ==> Matrix D 

   processing macro definition FD ==> Fraction D 

   processing macro definition MFD ==> Matrix FD 

   processing macro definition I ==> Integer 

   processing macro definition NNI ==> NonNegativeInteger 

   processing macro definition B ==> Boolean 

   processing macro definition OUT ==> OutputForm 

   processing macro definition ROWREC ==> Record(Indices: L C,Entries: L D) 

   processing macro definition iter ==> ::(iterated,Sy) 

   processing macro definition rand ==> ::(random,Sy) 

   processing macro definition Cat ==> -- the constructor category 
   processing macro definition Def ==> -- the constructor capsule 
------------------------------------------------------------------------
   initializing NRLIB LUD for LUDecomposition 
   compiling into NRLIB LUD 
   compiling exported LUDecomp : Matrix D -> Record(LU: Matrix D,Perm: Vector Integer,Pivots: List D)
Time: 0.14 SEC.

   compiling exported LUSolve : (Matrix D,Vector Integer,Vector D) -> Vector D
Time: 0.12 SEC.

   compiling exported LUInverse : Matrix D -> Record(Inv: Matrix D,Pivots: List D)
Time: 0.03 SEC.

(time taken in buildFunctor:  1)

;;;     ***       |LUDecomposition| REDEFINED

;;;     ***       |LUDecomposition| REDEFINED
Time: 0.01 SEC.

   Warnings: 
      [1] LUDecomp:  i0 has no value

   Cumulative Statistics for Constructor LUDecomposition
      Time: 0.30 seconds

   finalizing NRLIB LUD 
   Processing LUDecomposition for Browser database:
--------(LUDecomp ((Record (: LU MD) (: Perm (V I)) (: Pivots (L D))) MD))---------
--------(LUSolve ((V D) MD (V I) (V D)))---------
--------(LUInverse ((Record (: Inv MD) (: Pivots (L D))) MD))---------
--------constructor---------
------------------------------------------------------------------------
   LUDecomposition is now explicitly exposed in frame initial 
   LUDecomposition will be automatically loaded when needed from 
      /var/zope2/var/LatexWiki/LUD.NRLIB/code

A:=matrix [[subscript('a,[10*i+j]) for i in 1..3] for j in 1..3]

(1)

diagProduct(x) == reduce(*,[x(i,i) for i in 1..nrows(x)])

B:=LUDecomp A;

B.LU

(2)

B.Perm

(3)

B.Pivots

(4)

diagProduct(B.LU)=determinant A

Compiling function diagProduct with type Matrix Fraction Polynomial 
      Integer -> Fraction Polynomial Integer

(5)

%::Boolean

(6)