\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}