$SPAD/src/input lupfact.input

Timothy Daly December 3, 2016

Abstract

Contents

2

1 License

— license —

--Copyright The Numerical Algorithms Group Limited 1991.

———-

— * —

)set break resume
)spool lupfact.output
)set message test on
)set message auto off
)clear all

———-

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 nonsingular 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.

State the field here — * —

--S 1 of 18
field := Fraction Integer
--R
--R
--R   (1)  Fraction(Integer)
--R                                                                 Type: Domain
--E 1

———-

Next computes a permutation matrix for mult on the right — * —

--S 2 of 18
permMat: (INT, INT, INT) -> Matrix field
--R
--R                                                                   Type: Void
--E 2
--S 3 of 18
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
--R
--R                                                                   Type: Void
--E 3

———-

Find first col in first row that is nonzero or returns 0 — * —

--S 4 of 18
nonZeroCol: Matrix field -> INT
--R
--R                                                                   Type: Void
--E 4
--S 5 of 18
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
--R
--R                                                                   Type: Void
--E 5

———-

This embeds the given square matrix in a larger square matrix where the extra space is filled with 1s on the diagonal, 0 elsewhere.

— * —

--S 6 of 18
embedMatrix: (Matrix field,NNI,NNI) -> Matrix field
--R
--R                                                                   Type: Void
--E 6
--S 7 of 18
embedMatrix(m, oldDim, newDim) ==
  n := diagonalMatrix([1 for i in 1..newDim])$(Matrix(field))
  setsubMatrix!(n,1,1,m)
  n
--R
--R                                                                   Type: Void

--E 7

———-

The following implements algorithm in ”The Design and Analysis of Computer Algorithms” by Aho, Hopcroft and Ullman

— * —

--S 8 of 18
lupFactorEngine: (Matrix field, INT, INT)  -> List Matrix field
--R
--R                                                                   Type: Void
--E 8
--S 9 of 18
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]
--R
--R                                                                   Type: Void
--E 9

———-

Next computes floor of log base 2 of an integer — * —

--S 10 of 18
intLog2: NNI -> NNI
--R
--R                                                                   Type: Void
--E 10
--S 11 of 18
intLog2 n == if n = 1 then 0 else 1 + intLog2(n quo 2)
--R
--R                                                                   Type: Void
--E 11

———-

Here is the function to call — * —

--S 12 of 18
lupFactor: Matrix field -> Union(List Matrix field,"failed")
--R
--R                                                                   Type: Void
--E 12
--S 13 of 18
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)
--R
--R                                                                   Type: Void
--E 13

———-

Example from Aho, et al.

— * —

--S 14 of 18

m : Matrix field := zero(4,4)
--R
--R
--R         +0  0  0  0+
--R         |          |
--R         |0  0  0  0|
--R   (14)  |          |
--R         |0  0  0  0|
--R         |          |
--R         +0  0  0  0+
--R                                              Type: Matrix(Fraction(Integer))
--E 14
--S 15 of 18
for i in 4..1 by -1 repeat m(5-i,i) := i
--R
--R                                                                   Type: Void
--E 15
--S 16 of 18
m
--R
--R
--R         +0  0  0  4+
--R         |          |
--R         |0  0  3  0|
--R   (16)  |          |
--R         |0  2  0  0|
--R         |          |
--R         +1  0  0  0+
--R                                              Type: Matrix(Fraction(Integer))
--E 16
--S 17 of 18
lupFactor m
--R
--R   Compiling function intLog2 with type NonNegativeInteger ->
--R      NonNegativeInteger
--R   Compiling function embedMatrix with type (Matrix(Fraction(Integer)),
--R      NonNegativeInteger,NonNegativeInteger) -> Matrix(Fraction(Integer
--R      ))
--R   Compiling function nonZeroCol with type Matrix(Fraction(Integer))
--R       -> Integer
--R   Compiling function permMat with type (Integer,Integer,Integer) ->
--R      Matrix(Fraction(Integer))
--R   Compiling function lupFactorEngine with type (Matrix(Fraction(
--R      Integer)),Integer,Integer) -> List(Matrix(Fraction(Integer)))
--R   Compiling function lupFactor with type Matrix(Fraction(Integer)) ->
--R      Union(List(Matrix(Fraction(Integer))),"failed")
--I   Compiling function G7376 with type Integer -> Boolean

--R
--R          +1  0  0  0+ +4  0  0  0+ +0  0  0  1+
--R          |          | |          | |          |
--R          |0  1  0  0| |0  3  0  0| |0  0  1  0|
--R   (17)  [|          |,|          |,|          |]
--R          |0  0  1  0| |0  0  2  0| |0  1  0  0|
--R          |          | |          | |          |
--R          +0  0  0  1+ +0  0  0  1+ +1  0  0  0+
--R                             Type: Union(List(Matrix(Fraction(Integer))),...)
--E 17

———-

Example where the dimension does not start out a power of 2 — * —

--S 18 of 18
m := [[1,2,3],[2,3,1],[3,1,2]]
--R
--R
--R         +1  2  3+
--R         |       |
--R   (18)  |2  3  1|
--R         |       |
--R         +3  1  2+
--R                                              Type: Matrix(Fraction(Integer))
--E 18
)spool
)lisp (bye)

———-

References^{[1] nothing}

9