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
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
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}
|
