1 files changed, 146 insertions, 0 deletions
diff --git a/src/input/lupfact.input.pamphlet b/src/input/lupfact.input.pamphlet
new file mode 100644
index 00000000..85ffa25c
--- /dev/null
+++ b/src/input/lupfact.input.pamphlet
@@ -0,0 +1,146 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/input lupfact.input}
+\author{The Axiom Team}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{License}
+<<license>>=
+--Copyright The Numerical Algorithms Group Limited 1991.
+@
+<<*>>=
+<<license>>
+
+-- This file contains some functions that compute LUP factorizations of
+-- matrices over a field.  The main function to call is lupFactor.  It
+-- accepts one argument, which should be a non-singular square matrix.
+-- If the matrix is not square, "failed" will be returned.  If the matrix
+-- is non-singular, a 'cannot coerce "failed"' error will be displayed.
+-- lupFactor returns a Union(List Matrix field,"failed") object.  Coerce
+-- this to a List Matrix field before you try to use it.  See the comment
+-- before the definition of lupFactor for the reference for the
+-- algorithm.
+ 
+)clear all
+ 
+-- state the field here
+ 
+field := Fraction Integer
+ 
+-- next computes a permutation matrix for mult on the right
+ 
+permMat: (INT, INT, INT) -> Matrix field
+ 
+permMat(dim, i, j) ==
+  m : Matrix field :=
+    diagonalMatrix [(if i = k or j = k then 0 else 1) for k in 1..dim]
+  m(i,j) := 1
+  m(j,i) := 1
+  m
+ 
+-- find first col in first row that is nonzero or returns 0
+ 
+nonZeroCol: Matrix field -> INT
+ 
+nonZeroCol(m) ==
+  foundit := false
+  col := 1
+  for i in 1..ncols(m) while not foundit repeat
+    for j in 1..nrows(m) while not foundit repeat
+      if not(m(j,i) = 0) then
+        col := i
+        foundit := true
+  col
+ 
+-- this embeds the given square matrix in a larger square matrix
+-- where the extra space is filled with 1s on the diagonal, 0 elsewhere.
+ 
+embedMatrix: (Matrix field,NNI,NNI) -> Matrix field
+ 
+embedMatrix(m, oldDim, newDim) ==
+  n := diagonalMatrix([1 for i in 1..newDim])$(Matrix(field))
+  setsubMatrix!(n,1,1,m)
+  n
+ 
+-- following implements algorithm in "The Design and Analysis of
+-- Computer Algorithms" by Aho, Hopcroft and Ullman
+ 
+lupFactorEngine: (Matrix field, INT, INT)  -> List Matrix field
+ 
+lupFactorEngine(a, m, p) ==
+  m = 1 =>
+    l : Matrix field := diagonalMatrix [1]
+    pm : Matrix field := permMat(p,1,nonZeroCol a)
+    [l,a*pm,pm]
+  m2 : NNI := m quo 2
+  b : Matrix field := subMatrix(a,1,m2,1,p)
+  c : Matrix field := subMatrix(a,m2+1,m,1,p)
+  lup := lupFactorEngine(b,m2,p)
+  l1 := lup.1
+  u1 := lup.2
+  pm1 := lup.3
+  d : Matrix field := c * (inverse(pm1) :: Matrix(field))
+  e : Matrix field := subMatrix(u1,1,m2,1,m2)
+  f : Matrix field := subMatrix(d,1,m2,1,m2)
+  g : Matrix field := d - f * (inverse(e) :: Matrix(field)) * u1
+  pmin2 : NNI := p - m2
+  g' : Matrix field := subMatrix(g,1,nrows(g),p - pmin2 + 1,p)
+  lup := lupFactorEngine(g',m2,pmin2)
+  l2 := lup.1
+  u2 := lup.2
+  pm2 := lup.3
+  pm3 := horizConcat(zero(pmin2,m2)$(Matrix field), pm2)
+  pm3 := vertConcat(horizConcat(diagonalMatrix [1 for i in 1..m2],
+    zero(m2,pmin2)$(Matrix field)),pm3)
+  h : Matrix field := u1 * (inverse(pm3) :: Matrix(field))
+  l : Matrix field := horizConcat(l1, zero(m2,m2)$(Matrix field))
+  l := vertConcat(l,horizConcat(f * (inverse(e) :: Matrix(field)), l2))
+  u : Matrix field := horizConcat(zero(m2,m2)$(Matrix field), u2)
+  u := vertConcat(h,u)
+  pm := pm3 * pm1
+  [l,u,pm]
+ 
+-- next computes floor of log base 2 of an integer
+ 
+intLog2: NNI -> NNI
+ 
+intLog2 n == if n = 1 then 0 else 1 + intLog2(n quo 2)
+ 
+-- here is the function to call
+ 
+lupFactor: Matrix field -> Union(List Matrix field,"failed")
+ 
+lupFactor m ==
+  not((r := nrows m) = ncols m) =>
+    messagePrint("Matrix must be square")$OUTFORM
+    "failed"
+  ilog := intLog2(2)
+  not(r = 2 ** ilog) =>
+    m := embedMatrix(m,r,(n := 2 ** (ilog + 1)))
+    l := lupFactorEngine(m,n,n)
+    [subMatrix(l.1,1,r,1,r),subMatrix(l.2,1,r,1,r),
+      subMatrix(l.3,1,r,1,r)]
+  lupFactorEngine(m,r,r)
+ 
+-- Example from Aho, et al.
+ 
+m : Matrix field := zero(4,4)
+for i in 4..1 by -1 repeat m(5-i,i) := i
+m
+ 
+lupFactor m
+ 
+-- Example where the dimension does not start out a power of 2
+ 
+m := [[1,2,3],[2,3,1],[3,1,2]]
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}