/usr/lib/open-axiom/input/lupfact.input is in open-axiom-test 1.4.1+svn~2626-2ubuntu2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 | --Copyright The Numerical Algorithms Group Limited 1991.
-- 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]]
|