Edit detail for Jet LUDecomposition revision 1 of 3

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)
****** comp fails at level 3 with expression: ******
error in function LUDecomp 

(SEQ
 (LET A
   (|copy| AA))
 (LET |minR|
   (|minRowIndex| A))
 (LET |maxR|
   (|maxRowIndex| A))
 (LET |minC|
   (|minColIndex| A))
 (LET |maxC|
   (|maxColIndex| A))
 (LET #:G682
   | << |
   (^= |maxR| |maxC|)
   | >> |)
 (|exit| 1
  (IF #:G682
      (|error| "LU decomposition only of square matrices")
      (SEQ
       (LET (|:| |PermV| (|Vector| (|Integer|)))
         (|new| (|::| (+ (- |maxR| |minR|) 1) (|NonNegativeInteger|)) 0))
       (LET (|:| |Pivs| (|List| D))
         |empty|)
       (REPEAT (STEP |j| |minC| 1 |maxC|)
               (SEQ
                (REPEAT (STEP |i| |minR| 1 (- |j| 1))
                        (SEQ
                         (LET |s|
                           (|qelt| A |i| |j|))
                         (REPEAT (STEP |k| |minR| 1 (- |i| 1))
                                 (LET |s|
                                   (- |s|
                                      (* (|qelt| A |i| |k|)
                                         (|qelt| A |k| |j|)))))
                         (|exit| 1 (|qsetelt!| A |i| |j| |s|))))
                (LET (|:| |i0| (|Integer|))
                  (- 1))
                (REPEAT (STEP |i| |j| 1 |maxR|)
                        (SEQ
                         (LET |s|
                           (|qelt| A |i| |j|))
                         (REPEAT (STEP |k| |minC| 1 (- |j| 1))
                                 (LET |s|
                                   (- |s|
                                      (* (|qelt| A |i| |k|)
                                         (|qelt| A |k| |j|)))))
                         (|qsetelt!| A |i| |j| |s|)
                         (LET #:G683
                           (|zero?| |s|))
                         (|exit| 1
                          (IF #:G683
                              |noBranch|
                              (IF (< |i0| 0)
                                  (LET |i0|
                                    |i|)
                                  |noBranch|)))))
                (|exit| 1
                 (IF (< |i0| 0)
                     (|error| "singular matrix in LUDecomp")
                     (SEQ
                      (SEQ
                       (LET #:G684
                         (^= |j| |i0|))
                       (|exit| 1
                        (IF #:G684
                            (|swapRows!| A |j| |i0|)
                            |noBranch|)))
                      (|qsetelt!| |PermV| |j| |i0|)
                      (LET |Pivs|
                        (|cons| (|qelt| A |j| |j|) |Pivs|))
                      (LET #:G685
                        (^= |j| |maxC|))
                      (|exit| 1
                       (IF #:G685
                           (SEQ
                            (LET |d|
                              (/ 1 (|qelt| A |j| |j|)))
                            (|exit| 1
                             (REPEAT (STEP |k| (+ |j| 1) 1 |maxR|)
                                     (|qsetelt!| A |k| |j|
                                      (* |d| (|qelt| A |k| |j|))))))
                           |noBranch|)))))))
       (|exit| 1 (|construct| A |PermV| |Pivs|))))))
****** level 3  ******
$x:= (^= maxR maxC)
$m:= $EmptyMode
$f:=
((((|maxC| #) (|minC| #) (|maxR| #) (|minR| #) ...)))

   >> Apparent user error:
   NoValueMode
    is an unknown mode

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;

   There are no library operations named LUDecomp
      Use HyperDoc Browse or issue
                              )what op LUDecomp
      to learn if there is any operation containing " LUDecomp " in its
      name.

   Cannot find a definition or applicable library operation named
      LUDecomp with argument type(s)
                         Matrix(Polynomial(Integer))

      Perhaps you should use "@" to indicate the required return type,
      or "$" to specify which version of the function you need.
B.LU

   There are no library operations named B
      Use HyperDoc Browse or issue
                                 )what op B
      to learn if there is any operation containing " B " in its name.

   Cannot find a definition or applicable library operation named B
      with argument type(s)
                                Variable(LU)

      Perhaps you should use "@" to indicate the required return type,
      or "$" to specify which version of the function you need.
B.Perm

   There are no library operations named B
      Use HyperDoc Browse or issue
                                 )what op B
      to learn if there is any operation containing " B " in its name.

   Cannot find a definition or applicable library operation named B
      with argument type(s)
                               Variable(Perm)

      Perhaps you should use "@" to indicate the required return type,
      or "$" to specify which version of the function you need.
B.Pivots

   There are no library operations named B
      Use HyperDoc Browse or issue
                                 )what op B
      to learn if there is any operation containing " B " in its name.

   Cannot find a definition or applicable library operation named B
      with argument type(s)
                              Variable(Pivots)

      Perhaps you should use "@" to indicate the required return type,
      or "$" to specify which version of the function you need.
diagProduct(B.LU)=determinant A

   There are no library operations named B
      Use HyperDoc Browse or issue
                                 )what op B
      to learn if there is any operation containing " B " in its name.

   Cannot find a definition or applicable library operation named B
      with argument type(s)
                                Variable(LU)

      Perhaps you should use "@" to indicate the required return type,
      or "$" to specify which version of the function you need.
%::Boolean

   Cannot convert from type Void to Boolean for value
   "()"