|
Row echelon form and other basic operations for sparse matrices.
LU decomposition for ordinary matrices.
(C) 1995 Werner M. Seiler (axiom-xl version)
aldor #include "axiom.as"
#pile
macro Sy == Symbol
macro L == List
macro V == Vector
macro M == Matrix
macro F == Fraction
macro I == Integer
macro NNI == NonNegativeInteger
macro B == Boolean
macro OUT == OutputForm
macro ROWREC == Record(Indices:L C, Entries:L D)
macro iter == "iterated"::Sy
macro rand == "random"::Sy
-- ------------------------- --
-- SparseEchelonMatrix (SEM) --
-- ------------------------- --
+++ Description:
+++ \axiom{SparseEchelonMatrix(C,D)} implements sparse matrices whose columns
+++ are enumerated by the \axiomType{OrderedSet} \axiom{C} and whose entries
+++ belong to the \axiomType{GcdDomain} \axiom{D}. The basic operation of
+++ this domain is the computation of an row echelon form. The used algorithm
+++ tries to maintain the sparsity and is especially adapted to matrices who
+++ are already close to a row echelon form.
SparseEchelonMatrix(C:OrderedSet,D:Ring):Cat == Def where
Cat ==> CoercibleTo OUT with
shallowlyMutable
++ Matrices may be altered destructively.
finiteAggregate
++ Matrices are finite.
coerce : % -> M D
++ \axiom{coerce(A)} yields the matrix \axiom{A} in the usual matrix type.
copy : % -> %
++ \axiom{copy(A)} returns a copy of the matrix \axiom{A}.
ncols : % -> NNI
++ \axiom{ncols(A)} returns the number of columns of the matrix \axiom{A}.
nrows : % -> NNI
++ \axiom{nrows(A)} returns the number of rows of the matrix \axiom{A}.
allIndices : % -> L C
++ \axiom{allIndices(A)} returns all indices used for enumerating the
++ columns of the matrix \axiom{A}.
elimZeroCols_! : % -> ()
++ \axiom{elimZeroCols!(A)} removes columns which contain only zeros.
++ This effects basically only the value of \axiom{allIndices(A)}.
purge_! : (%,C->B) -> ()
++ \axiom{purge!(A,crit)} eliminates all columns belonging to an index
++ \axiom{c} such that \axiom{crit(c)} yields \axiom{true}.
sortedPurge_! : (%,C->B) -> ()
++ \axiom{sortedPurge!(A,crit)} is like \axiom{purge}, however, with the
++ additional assumption that \axiom{crit} respects the ordering of the
++ indices.
new : (L C, I) -> %
++ \axiom{new(inds,nrows)} generates a new matrix with \axiom{nrows}
++ rows and columns enumerated by the indices \axiom{inds}. The matrix
++ is empty, i.e. the zero matrix.
elt : (%,I,C) -> D
++ \axiom{elt(A,i,c)} returns the entry of the matrix \axiom{A} in row
++ \axiom{i} and in the column with index \axiom{c}.
setelt_! : (%,I,C,D) -> ()
++ \axiom{setelt!(A,i,c,d)} sets the entry of the matrix \axiom{A} in
++ row \axiom{i} and in the column with index \axiom{c} to the value
++ \axiom{d}.
row : (%,I) -> ROWREC
++ \axiom{row(A,i)} returns the \axiom{i}-th row of the matrix \axiom{A}.
setRow_! : (%,I,ROWREC) -> ()
++ \axiom{setRow!(A,i,ind,ent)} sets the \axiom{i}-th row of the matrix
++ \axiom{A} to the value \axiom{r}.
setRow_! : (%,I,L C,L D) -> ()
++ \axiom{setRow!(A,i,ind,ent)} sets the \axiom{i}-th row of the matrix
++ \axiom{A}. Its indices are \axiom{ind}; the entries \axiom{ent}.
deleteRow_! : (%,I) -> ()
++ \axiom{deleteRow(A,i)} deletes the \axiom{i}-th row of the matrix
++ \axiom{A}.
consRow_! : (%,ROWREC) -> ()
++ \axiom{consRow!(A,r)} inserts the row \axiom{r} at the top of the
++ matrix \axiom{A}.
appendRow_! : (%,ROWREC) -> ()
++ \axiom{appendRow!(A,r)} appends the row \axiom{r} at the end of the
++ matrix \axiom{A}.
extract : (%,I,I) -> %
++ \axiom{extract(A,i1,i2)} extracts the rows \axiom{i1} to \axiom{i2}
++ and returns them as a new matrix.
rowEchelon : % -> Record(Ech:%, Lt:M D, Pivots:L D, Rank:NNI)
++ \axiom{primitiveRowEchelon(A)} computes a row echelon form for the
++ matrix \axiom{A}. The algorithm used is fraction-free elimination.
++ It is especially adapted to matrices already close to row echelon
++ form. The transformation matrix, the used pivots and the rank of the
++ matrix are also returned.
pivot : (%,I) -> Record(Index:C, Entry:D)
++ \axiom{pivot(A,i)} returns the leading entry of the \axiom{i}-th row
++ of the matrix \axiom{A} together with its index.
pivots : % -> ROWREC
++ \axiom{pivots(A)} returns all leading entries of the matrix \axiom{A}
++ together with their indices.
* : (M D,%) -> %
++ \axiom{L*A} implements left multiplication with a usual matrix.
join : (%,%) -> %
++ \axiom{join(A,B)} vertically concats the matrices \axiom{A} and
++ \axiom{B}.
horizJoin : (%,%) -> %
++ \axiom{horizJoin(A,B)} horizontally concats the matrices \axiom{A} and
++ \axiom{B}. It is assumed that all indices of \axiom{B} are smaller than
++ those of \axiom{A}.
horizSplit : (%,C) -> Record(Left:%, Right:%)
++ \axiom{horizSplit(A,c)} splits the matrix \axiom{A} into two at the
++ column given by \axiom{c}. The first column of the right matrix is
++ enumerated by the first index less or equal to \axiom{c}.
--if D has GcdDomain then
-- setGcdMode : Sy -> Sy
-- ++ \axiom{setGcdMode(s)} sets a new value for the flag deciding on
-- ++ the method used to compute gcd's for lists. Possible values for
-- ++ \axiom{s} are \axiom{iterated} and \axiom{random}.
-- * : (M F D,%) -> %
-- ++ \axiom{L*A} implements left multiplication with a usual matrix over
-- ++ the quotient field of \axiom{D}.
-- primitiveRowEchelon : % -> Record(Ech:%, Lt:M F D, Pivots:L D, Rank:NNI)
-- ++ \axiom{primitiveRowEchelon(A)} computes a row echelon form for the
-- ++ matrix \axiom{A}. The algorithm used is fraction-free elimination.
-- ++ Every row is made primitive by division by the gcd. The algorithm
-- ++ is especially adapted to matrices already close to row echelon
-- ++ form. The transformation matrix, the used pivots and the rank of the
-- ++ matrix are also returned.
Def ==> add
import from C, D, M D
import from L ROWREC, V ROWREC, L NNI
local
minInd:I := minIndex([i for i:NNI in 1..1])
offset:I := minInd-1
emptyRec:ROWREC := [empty()$(L C), empty()$(L D)]
noChecks?:B := false --D has lazyRep -- flag for lazy representation
seed:I := 113 -- seed for random number generation
GCDmode:Sy := iter -- flag for gcd algorithm
greater(c1:C,c2:C):B == c2<c1
greater(r1:ROWREC,r2:ROWREC):B ==
empty? r1.Indices => false
empty? r2.Indices => true
first(r2.Indices) < first(r1.Indices)
-- -------------- --
-- Representation --
-- -------------- --
-- For efficency reasons most checks for correct index ranges are omitted.
macro Rep == Record(NCols:NNI, NRows:NNI, AllInds:L C, Rows:V ROWREC)
import from Rep
ncols(A:%):NNI == rep(A).NCols
nrows(A:%):NNI == rep(A).NRows
allIndices(A:%):L C == copy rep(A).AllInds
row(A:%,i:I):ROWREC ==
-- i<0 or i>rep(A).NRows => error "index out of range"
qelt(rep(A).Rows,i)
setRow_!(A:%,i:I,r:ROWREC):() ==
-- i<0 or i>rep(A).NRows => error "index out of range"
qsetelt!(rep(A).Rows,i,r)
setRow_!(A:%,i:I,inds:L C,ents:L D):() ==
-- i<0 or i>rep(A).NRows => error "index out of range"
-- #inds ~= #ents => error "improper row"
qsetelt!(rep(A).Rows, i, [inds,ents])
new(inds:L C,n:I):% ==
free emptyRec
per [#inds, n::NNI, inds, [record explode emptyRec for i in 1..n]]
elt(A:%,i:I,c:C):D ==
r := row(A,i)
pos := position(c,r.Indices)
pos<minInd => 0$D
qelt(r.Entries,pos)
setelt_!(A:%,i:I,c:C,d:D):() ==
r := row(A,i)
pos := position(c,r.Indices)
if pos>=minInd then
qsetelt!(r.Entries,pos,d)
else
j := minInd
for ind in r.Indices while c<ind repeat
j := j+1
r.Indices := insert!(c,r.Indices,j)
r.Entries := insert!(d,r.Entries,j)
qsetelt!(rep(A).Rows,i,r)
coerce(A:%):M D ==
(nc,nr,ai,ro) := explode rep A
zero? nc => error "cannot coerce matrix with zero columns"
AA:M D := new(nr,nc,0$D)
local inds:L C
local ents:L D
for r:ROWREC in entries(ro) for i in minRowIndex(AA).. repeat
inds := r.Indices
ents := r.Entries
for ind in ai for j in minColIndex(AA).. _
while not empty? inds repeat
if ind=first inds then
qsetelt!(AA, i, j, first ents)
inds := rest inds
ents := rest ents
AA
coerce(A:%):OUT ==
zero? rep(A).NCols => 0$D ::OUT
A::(M D)::OUT
copy(A:%):% ==
(Ancol,Anrow,Ainds,Arows) := explode rep A
resRows:V ROWREC := new(Anrow,emptyRec)
for l in 1..Anrow repeat
r:ROWREC := qelt(Arows,l)
qsetelt!(resRows, l, [copy r.Indices, copy r.Entries])
per [Ancol, Anrow, copy Ainds, resRows]
-- ----------------------- --
-- Basic Matrix Operations --
-- ----------------------- --
elimZeroCols_!(A:%):() ==
newInds:L C := empty()
for r in entries(rep(A).Rows) repeat
newInds := removeDuplicates! merge!(greater, newInds, r.Indices)
rep(A).AllInds := newInds
purge_!(A:%,crit:C->B):() ==
rA := rep A
newInds:L C := empty()
for c in rA.AllInds repeat
if not crit c then
newInds := cons(c,newInds)
newInds := reverse! newInds
local r:ROWREC
if #newInds ~= #rA.AllInds then
rA.AllInds := newInds
for l in 1..rA.NRows repeat
r := qelt(rA.Rows,l)
newInds:L C := empty()
newEnts:L D := empty()
for c:C in r.Indices for e:D in r.Entries repeat
if not crit c then
newInds := cons(c,newInds)
newEnts := cons(e,newEnts)
qsetelt!(rA.Rows, l, [reverse! newInds, reverse! newEnts])
sortedPurge_!(A:%,crit:C->B):() ==
rA := rep A
local r:ROWREC
if crit first rA.AllInds then
while not(empty? rA.AllInds) and crit first rA.AllInds repeat
rA.AllInds := rest rA.AllInds
for l in 1..rA.NRows repeat
r := qelt(rA.Rows,l)
while not(empty? r.Indices) and crit first r.Indices repeat
r.Indices := rest r.Indices
r.Entries := rest r.Entries
qsetelt!(rA.Rows,l,r)
deleteRow_!(A:%,i:I):() ==
rA := rep A
((rA.NRows)::I<i) => A
nr := ((rA.NRows)::I-1)::NNI
resRows:V ROWREC := new(nr,emptyRec)
for l in 1..(i-1) repeat
qsetelt!(resRows,l,qelt(rA.Rows,l))
for l in (i+1)::NNI..rA.NRows repeat
qsetelt!(resRows,l-1,qelt(rA.Rows,l))
rA.NRows := nr
rA.Rows := resRows
consRow_!(A:%,r:ROWREC):() ==
rA := rep A
rA.NRows := rA.NRows + 1
newRows:L ROWREC := cons(r,entries rA.Rows)
rA.Rows := construct newRows
newInds := setDifference(r.Indices,rA.AllInds)
if not empty? newInds then
rA.AllInds := merge(greater,rA.AllInds,sort!(greater,newInds))
appendRow_!(A:%,r:ROWREC):() ==
rA := rep A
rA.NRows := rA.NRows + 1
newRows:L ROWREC := concat(entries rA.Rows, r)
rA.Rows := construct newRows
newInds := setDifference(r.Indices,rA.AllInds)
if not empty? newInds then
rA.AllInds := merge(greater,rA.AllInds,sort!(greater,newInds))
extract(A:%,i1:I,i2:I):% ==
nr := (i2-i1+1)::NNI
resRows:V ROWREC := new(nr,emptyRec)
newInds:L C := empty()
for i in i1..i2 repeat
qsetelt!(resRows,i-i1+1,row(A,i))
newInds := removeDuplicates! merge(greater,newInds,row(A,i).Indices)
per [rep(A).NCols, nr, newInds, resRows]
join(A1:%,A2:%):% ==
rA1 := rep A1
rA2 := rep A2
newInds := removeDuplicates! merge(greater,rA1.AllInds,rA2.AllInds)
newNRows := rA1.NRows + rA2.NRows
newRows:V ROWREC := new(newNRows,emptyRec)
for l in 1..rA1.NRows repeat
qsetelt!(newRows, l, qelt(rA1.Rows,l))
for l in 1..rA2.NRows repeat
qsetelt!(newRows, rA1.NRows+l, qelt(rA2.Rows,l))
per [#newInds, newNRows, newInds, newRows]
horizJoin(A1:%,A2:%):% ==
rA1 := rep A1
rA2 := rep A2
rA1.NRows~=rA2.NRows => error "incompatible dimensions in horizJoin"
newInds := append(rA1.AllInds,rA2.AllInds)
res:% := new(newInds,rA1.NRows)
for i in 1..rA1.NRows repeat
r1 := row(A1,i)
r2 := row(A2,i)
setRow!(res, i, append(r1.Indices,r2.Indices), _
append(r1.Entries,r2.Entries))
res
horizSplit(A:%,c:C):Record(Left:%,Right:%) ==
rA := rep A
rinds:L C := allIndices A
linds:L C := empty()
while not(empty? rinds) and (first(rinds)>c) repeat
linds := cons(first(rinds),linds)
rinds := rest rinds
empty? linds => [new(linds,rA.NRows), A]
linds := reverse! linds
empty? rinds => [A, new(rinds,rA.NRows)]
LA:% := new(linds,rA.NRows)
RA:% := new(rinds,rA.NRows)
for i in 1..rA.NRows repeat
r := row(A,i)
ri:L C := r.Indices
re:L D := r.Entries
li:L C := empty()
le:L D := empty()
while not(empty? ri) and (first(ri)>c) repeat
li := cons(first(ri),li)
le := cons(first re, le)
ri := rest ri
re := rest re
if not empty? li then
li := reverse! li
le := reverse! le
setRow!(LA,i,li,le)
if not empty? ri then
setRow!(RA,i,ri,re)
[LA,RA]
-- ----------- --
-- Row Echelon --
-- ----------- --
addRows(d1:D,r1:ROWREC,d2:D,r2:ROWREC):ROWREC ==
-- Computes linear combination of two rows.
-- Local function.
empty? r1.Indices =>
one? d2 => r2
[r2.Indices, [d2*e2 for e2 in r2.Entries]]
empty? r2.Indices =>
one? d1 => r1
[r1.Indices, [d1*e1 for e1 in r1.Entries]]
resI:L C := empty()
resE:L D := empty()
local lent1:L D
local lent2:L D
if not(noChecks?) and one? d1 then
lent1 := r1.Entries
else
lent1 := [d1*e1 for e1 in r1.Entries]
if not(noChecks?) and one? d2 then
lent2 := copy r2.Entries
else
lent2 := [d2*e2 for e2 in r2.Entries]
lind2 := copy r2.Indices
for c1 in r1.Indices for e1 in lent1 repeat
while not(empty? lind2) and c1<first(lind2) repeat
resI := cons(first lind2, resI)
resE := cons(first(lent2), resE)
lind2 := rest lind2
lent2 := rest lent2
if not(empty? lind2) and first(lind2)=c1 then
sum := e1+first(lent2)
if noChecks? or not zero? sum then
resI := cons(c1,resI)
resE := cons(sum,resE)
lind2 := rest lind2
lent2 := rest lent2
else
resI := cons(c1,resI)
resE := cons(e1,resE)
resI := concat!(reverse! resI, lind2)
resE := concat!(reverse! resE, lent2)
while not(empty? resE) and zero? first resE repeat
resI := rest resI
resE := rest resE
[resI,resE]
pivot(A:%,i:I):Record(Index:C,Entry:D) ==
r := row(A,i)
empty? r.Indices => error "empty row"
[first r.Indices, first r.Entries]
pivots(A:%):ROWREC ==
resI:L C := empty()
resE:L D := empty()
for r in entries rep(A).Rows | not empty? r.Indices repeat
resI := cons(first r.Indices, resI)
resE := cons(first r.Entries, resE)
[reverse! resI, reverse! resE]
rowEchelon(AA:%):Record(Ech:%, Lt:M D, Pivots:L D, Rank:NNI) ==
A := rep copy AA
LTr:M D := diagonalMatrix [1$D for i in 1..A.NRows]
Pivs:L D := empty()
-- check pivots
for i in 1..A.NRows repeat
r := qelt(A.Rows,i)
changed?:B := false
while not(empty? r.Entries) and zero? first r.Entries repeat
r.Entries := rest r.Entries
r.Indices := rest r.Indices
changed? := true
if changed? then
qsetelt!(A.Rows,i,r)
-- sort rows by pivots (bubble sort)
notSorted?:B := true
while notSorted? repeat
notSorted? := false
oldr:ROWREC := qelt(A.Rows,1)
for i in 2..A.NRows repeat
newr:ROWREC := qelt(A.Rows,i)
if greater(newr,oldr) then
qsetelt!(A.Rows,i,oldr)
qsetelt!(A.Rows,i-1,newr)
swapRows!(LTr,i-1,i)
notSorted? := true
else
oldr := newr
-- fraction-free elimination
early?:B := false
local pivlen,pivrow,rk:NNI
for i in 1..A.NRows repeat
r:ROWREC := qelt(A.Rows,i)
if empty? r.Indices then
rk:NNI := (i-1)::NNI
early? := true
break
-- search good pivot
pivind := first r.Indices
pivlen := #r.Indices
pivrow := i
k:I := 0
for j in (i+1)..A.NRows _
while not(empty? qelt(A.Rows,j).Indices) and _
pivind=first(qelt(A.Rows,j).Indices) repeat
len:NNI := #qelt(A.Rows,j).Indices
k := k+1
if len<pivlen then
pivlen := len
pivrow := j
piv:D := first qelt(A.Rows,pivrow).Entries
Pivs := cons(piv,Pivs)
-- elimination necessary?
if k>0 then
if pivrow~=i then
pr := qelt(A.Rows,pivrow)
qsetelt!(A.Rows,pivrow,qelt(A.Rows,i))
qsetelt!(A.Rows,i,pr)
swapRows!(LTr,i,pivrow)
-- elimination (and resorting of rows)
pr := copy qelt(A.Rows,i)
pr.Indices := rest pr.Indices
pr.Entries := rest pr.Entries
for j in (i+1)..(i+k) repeat
r := copy qelt(A.Rows,i+1)
c := first r.Entries
r.Indices := rest r.Indices
r.Entries := rest r.Entries
r := addRows(piv,r,-c,pr)
for l in 1..A.NRows repeat
f := piv*qelt(LTr,i+1,l) - c*qelt(LTr,i,l)
qsetelt!(LTr,i+1,l,f)
for l in (i+2)..(2*i+k+1-j) repeat
qsetelt!(A.Rows,l-1,qelt(A.Rows,l))
swapRows!(LTr,l-1,l)
for l in (2*i+k+2-j)..A.NRows _
while greater(qelt(A.Rows,l),r) repeat
qsetelt!(A.Rows,l-1,qelt(A.Rows,l))
swapRows!(LTr,l-1,l)
qsetelt!(A.Rows,l-1,r)
if not early? then
rk:NNI := A.NRows
[per(A), LTr, Pivs, rk]
-- -------------- --
-- Multiplication --
-- -------------- --
(LL:M D) * (AA:%) : % ==
ncols(LL)~=AA.NRows => error "improper matrix dimensions"
A := rep copy AA
rlen := nrows LL
res:% := new(A.AllInds,rlen)
for c in A.AllInds repeat
tmp:V D := new(rlen,0$D)
for i in 1..A.NRows repeat
r:ROWREC := qelt(A.Rows,i)
inds:L C := r.Indices
if not(empty? inds) and first(inds)=c then
for k in 1..rlen | not zero? qelt(LL,k,i) repeat
qsetelt!(tmp, k, qelt(tmp,k) + qelt(LL,k,i)*first(r.Entries))
r.Entries := rest r.Entries
r.Indices := rest inds
qsetelt!(A.Rows,i,r)
for k in 1..rlen | not zero? qelt(tmp,k) repeat
r := qelt(res.Rows,k)
r.Indices := cons(c,r.Indices)
r.Entries := cons(qelt(tmp,k),r.Entries)
qsetelt!(res.Rows,k,r)
for k in 1..rlen repeat
r := qelt(res.Rows,k)
r.Indices := reverse! r.Indices
r.Entries := reverse! r.Entries
qsetelt!(res.Rows,k,r)
per res
-- ---------------------------- --
-- primitive form in gcd domain --
-- ---------------------------- --
--if D has Join(GcdDomain,IntegralDomain) then
-- import from Fraction D
-- setGcdMode(s:Sy):Sy ==
-- free GCDmode
-- tmp := GCDmode
-- (s=iter) or (s=rand) =>
-- GCDmode := s
-- tmp
-- error "unknown gcd mode"
-- randomGCD(le:L D):D ==
-- -- Probabilistic technique.
-- free seed
-- #le=2 => gcd(first le, second le)
-- f := first le
-- g := second le
-- l := rest rest le
-- while not empty? l repeat
-- one? first l => return 1$D
-- f := f + (1+random(seed)$I)*first(l)
-- l := rest l
-- if not empty? l then
-- one? first l => return 1$D
-- g := g + (1+random(seed)$I)*first(l)
-- l := rest l
-- h := gcd(f,g)
-- l := [h]
-- for e in le repeat
-- tmp := exquo(e,h)
-- if tmp case "failed" then
-- l := cons(e,l)
-- one?(#l) => h
-- randomGCD l
-- iteratedGCD(le:L D):D ==
-- -- Computes gcd iteratively
-- res := gcd(first le, second le)
-- l := rest rest le
-- while not(empty?(l) or one?(res)) repeat
-- res := gcd(res, first l)
-- l := rest l
-- res
-- makePrimitive(r:ROWREC):Record(GCD:D, Row:ROWREC) ==
-- -- remove common gcd of row
-- free GCDmode
-- le := r.Entries
-- one?(#le) => [first le, [r.Indices,[1$D]]]
-- local g:D
-- if GCDmode=iterated then
-- g := iteratedGCD le
-- else
-- g := randomGCD le
-- one? g => [1,r]
-- le := [exquo(e,g)::D for e in le]
-- [g, [r.Indices,le]]
-- primitiveRowEchelon(AA:%) : _
-- Record(Ech:%, Lt:M F D, Pivots:L D, Rank:NNI) ==
-- A := rep copy AA
-- LTr:M F D := diagonalMatrix [1$(F D) for i in 1..A.NRows]
-- Pivs:L D := empty()
-- -- check pivots
-- for i in 1..A.NRows repeat
-- r := qelt(A.Rows,i)
-- changed?:B := false
-- while not(empty? r.Entries) and zero? first r.Entries repeat
-- r.Entries := rest r.Entries
-- r.Indices := rest r.Indices
-- changed? := true
-- if changed? then
-- qsetelt!(A.Rows,i,r)
-- -- sort rows by pivots (bubble sort)
-- notSorted?:B := true
-- while notSorted? repeat
-- notSorted? := false
-- oldr := qelt(A.Rows,1)
-- for i in 2..A.NRows repeat
-- newr := qelt(A.Rows,i)
-- if greater(newr,oldr) then
-- qsetelt!(A.Rows,i,oldr)
-- qsetelt!(A.Rows,i-1,newr)
-- swapRows!(LTr,i-1,i)
-- notSorted? := true
-- else
-- oldr := newr
-- -- primitive fraction-free elimination
-- early?:B := false
-- local pivlen,pivrow,rk:NNI
-- for i in 1..A.NRows until finished? repeat
-- r := qelt(A.Rows,i)
-- if empty? r.Indices then
-- rk:NNI := (i-1)::NNI
-- early? := true
-- break
-- -- search good pivot
-- pivind := first r.Indices
-- pivlen := #r.Indices
-- pivrow := i
-- k:I := 0
-- for j in (i+1)..A.NRows _
-- while not(empty? qelt(A.Rows,j).Indices) and _
-- pivind=first(qelt(A.Rows,j).Indices) repeat
-- len := #qelt(A.Rows,j).Indices
-- k := k+1
-- if len<pivlen then
-- pivlen := len
-- pivrow := j
-- -- make row primitive
-- tmp := makePrimitive qelt(A.Rows,pivrow)
-- if not one? tmp.GCD then
-- qsetelt!(A.Rows,pivrow,tmp.Row)
-- q:FD := 1/tmp.GCD
-- for l in 1..A.NRows | not zero? qelt(LTr,pivrow,l) repeat
-- qsetelt!(LTr,pivrow,l,q*qelt(LTr,pivrow,l))
-- piv:D := first qelt(A.Rows,pivrow).Entries
-- Pivs := cons(piv,Pivs)
-- -- elimination necessary?
-- if k>0 then
-- if pivrow~=i then
-- pr := qelt(A.Rows,pivrow)
-- qsetelt!(A.Rows,pivrow,qelt(A.Rows,i))
-- qsetelt!(A.Rows,i,pr)
-- swapRows!(LTr,i,pivrow)
-- -- elimination (and resorting of rows)
-- pr := copy tmp.Row
-- pr.Indices := rest pr.Indices
-- pr.Entries := rest pr.Entries
-- for j in (i+1)..(i+k) repeat
-- r := copy qelt(A.Rows,i+1)
-- c := first r.Entries
-- r.Indices := rest r.Indices
-- r.Entries := rest r.Entries
-- r := addRows(piv,r,-c,pr)
-- for l in 1..A.NRows repeat
-- fd:FD := piv *$FD qelt(LTr,i+1,l) - (c*qelt(LTr,i,l))::FD
-- qsetelt!(LTr,i+1,l,fd)
-- for l in (i+2)..(2*i+k+1-j) repeat
-- qsetelt!(A.Rows,l-1,qelt(A.Rows,l))
-- swapRows!(LTr,l-1,l)
-- for l in (2*i+k+2-j)..A.NRows _
-- while greater(qelt(A.Rows,l),r) repeat
-- qsetelt!(A.Rows,l-1,qelt(A.Rows,l))
-- swapRows!(LTr,l-1,l)
-- qsetelt!(A.Rows,l-1,r)
-- if not early? then
-- rk:NNI := A.NRows
-- [per(A), LTr, Pivs, rk]
-- mult(f:F D,d:D):D ==
-- res := numer(f)*d
-- tmp := exquo(res,denom(f))
-- tmp case "failed" => error "cannot divide in mult"
-- tmp::D
-- (LL:M F D) * (AA:%) : % ==
-- ncols(LL)~=AA.NRows => error "improper matrix dimensions"
-- A := rep copy AA
-- rlen := nrows LL
-- res:% := new(A.AllInds,rlen)
-- for c in A.AllInds repeat
-- tmp:V F D := new(rlen,0$FD)
-- for i in 1..A.NRows repeat
-- r := qelt(A.Rows,i)
-- inds := r.Indices
-- if not(empty? inds) and first(inds)=c then
-- for k in 1..rlen | not zero? qelt(LL,k,i) repeat
-- qsetelt!(tmp, k, qelt(tmp,k) + qelt(LL,k,i)*first(r.Entries))
-- r.Entries := rest r.Entries
-- r.Indices := rest inds
-- qsetelt!(A.Rows,i,r)
-- for k in 1..rlen | not zero? qelt(tmp,k) repeat
-- d:Union(value:D,failed:'failed') := retractIfCan qelt(tmp,k)
-- d case failed => error "cannot divide in *"
-- r := qelt(res.Rows,k)
-- r.Indices := cons(c,r.Indices)
-- r.Entries := cons(d::D,r.Entries)
-- qsetelt!(res.Rows,k,r)
-- for k in 1..rlen repeat
-- r := qelt(res.Rows,k)
-- r.Indices := reverse! r.Indices
-- r.Entries := reverse! r.Entries
-- qsetelt!(res.Rows,k,r)
-- per res
-- --------------------- --
-- LUDecomposition (LUD) --
-- --------------------- --
+++ Description:
+++ \axiomType{LUDecomposition} contains procedures to solve linear systems of
+++ equations or to compute inverses using a LU decomposition.
LUDecomposition(D:Field) : Cat == Def where
Cat ==> with
LUDecomp : M D -> Record(LU:M D, Perm:V I, Pivots:L D)
++ \axiom{LUDecomp(A)} computes a LU decomposition of \axiom{A}
++ using the algorithm of Crout. \axiom{LU} contains both triangular
++ matrices; \axiom{Perm} is the permutation used for partial
++ pivoting and \axiom{Pivots} yields the used pivots.
LUSolve : (M D, V I, V D) -> V D
++ \axiom{LUSolve(LU,Perm,B)} uses a previously computed LU
++ decomposition to solve a linear system with right hand side
++ \axiom{B}. \axiom{LU} and \axiom{Perm} are as given by
++ \axiom{LUDecomp}.
LUInverse : M D -> Record(Inv:M D, Pivots:L D)
++ \axiom{LUInverse(A)} computes the inverse of \axiom{A} using a LU
++ decomposition.
Def ==> add
import from D
LUDecomp(AA:M D):Record(LU:M D, Perm:V I, Pivots:L D) ==
-- LU decomposition using Crout's algorithm with partial pivoting.
A := copy AA
minR := minRowIndex A; maxR := maxRowIndex A
minC := minColIndex A; maxC := maxColIndex A
maxR~=maxC => error "LU decomposition only of square matrices"
PermV:V I := new((maxR-minR+1)::NNI,0)
Pivs:L D := empty()
for j in minC..maxC repeat
for i in minR..(j-1) repeat
s := qelt(A,i,j)
for k in minR..(i-1) repeat
s := s - qelt(A,i,k)*qelt(A,k,j)
qsetelt!(A,i,j,s)
i0:I := -1
for i in j..maxR repeat
s := qelt(A,i,j)
for k in minC..(j-1) repeat
s := s - qelt(A,i,k)*qelt(A,k,j)
qsetelt!(A,i,j,s)
if not(zero? s) and i0<0 then
i0 := i -- first non-zero pivot
i0<0 => error "singular matrix in LUDecomp"
if j~=i0 then
swapRows!(A,j,i0)
qsetelt!(PermV,j,i0)
Pivs := cons(qelt(A,j,j),Pivs)
if j~=maxC then
d := 1/qelt(A,j,j)
for k in (j+1)..maxR repeat
qsetelt!(A,k,j,d*qelt(A,k,j))
[A,PermV,Pivs]
LUSolve(LU:M D,Perm:V I,XX:V D):V D ==
-- Solves LU decomposed linear system for right hand side XX
X := copy XX
minR := minRowIndex LU; maxR := maxRowIndex LU
maxIndex(X)~=maxR => error "Wrong dimensions in LUSolve"
ii:I := -1
for i in minR..maxR repeat -- forward substitution
ip := qelt(Perm,i)
s := qelt(X,ip)
qsetelt!(X,ip,qelt(X,i))
if ii>=0 then
for j in ii..(i-1) repeat
s := s - qelt(LU,i,j)*qelt(X,j)
else
if not zero? s then ii := i
qsetelt!(X,i,s)
for i in maxR..minR by -1 repeat -- back substitution
s := qelt(X,i)
for j in (i+1)..maxR repeat
s := s - qelt(LU,i,j)*qelt(X,j)
qsetelt!(X,i,s/qelt(LU,i,i))
X
LUInverse(A:M D):Record(Inv:M D, Pivots:L D) ==
-- Inversion via LU decomposition
Alu := LUDecomp A
n := ncols A
res:M D := new(n,n,0)
for i in minRowIndex(A)..maxRowIndex(A) repeat
v:V D := new(n,0)
qsetelt!(v,i,1)
res := setColumn!(res,i,LUSolve(Alu.LU,Alu.Perm,v))
[res,Alu.Pivots]
aldor Compiling FriCAS source code from file
/var/zope2/var/LatexWiki/2184467348127326179-25px001.as using
AXIOM-XL compiler and options
-O -Fasy -Fao -Flsp -laxiom -Mno-AXL_W_WillObsolete -DAxiom -Y $AXIOM/algebra
Use the system command )set compiler args to change these
options.
#1 (Warning) Deprecated message prefix: use `ALDOR_' instead of `_AXL'
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 57:
elimZeroCols_! : % -> ()
................^
[L57 C17] #2 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 61:
purge_! : (%,C->B) -> ()
.........^
[L61 C10] #3 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 65:
sortedPurge_! : (%,C->B) -> ()
...............^
[L65 C16] #4 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 79:
setelt_! : (%,I,C,D) -> ()
..........^
[L79 C11] #5 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 87:
setRow_! : (%,I,ROWREC) -> ()
..........^
[L87 C11] #6 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 91:
setRow_! : (%,I,L C,L D) -> ()
..........^
[L91 C11] #7 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 95:
deleteRow_! : (%,I) -> ()
.............^
[L95 C14] #8 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 99:
consRow_! : (%,ROWREC) -> ()
...........^
[L99 C12] #9 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 103:
appendRow_! : (%,ROWREC) -> ()
.............^
[L103 C14] #10 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 210:
setRow_!(A:%,i:I,r:ROWREC):() ==
..........^
[L210 C11] #11 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 215:
setRow_!(A:%,i:I,inds:L C,ents:L D):() ==
..........^
[L215 C11] #12 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 233:
setelt_!(A:%,i:I,c:C,d:D):() ==
..........^
[L233 C11] #13 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 274:
r:ROWREC := qelt(Arows,l)
.........................^
[L274 C26] #20 (Error) (After Macro Expansion) Argument 1 of `qelt' did not match any possible parameter type.
The rejected type is Vector(Record(Indices: List(C), Entries: List(D))).
Expected one of:
-- List(C)
-- List(D)
-- List(Record(Indices: List(C), Entries: List(D)))
-- List(NonNegativeInteger)
Expanded expression was: Arows
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 275:
qsetelt!(resRows, l, [copy r.Indices, copy r.Entries])
.................^
[L275 C18] #21 (Error) (After Macro Expansion) Argument 1 of `qsetelt!' did not match any possible parameter type.
The rejected type is Vector(Record(Indices: List(C), Entries: List(D))).
Expected one of:
-- List(C)
-- List(D)
-- List(Record(Indices: List(C), Entries: List(D)))
-- List(NonNegativeInteger)
Expanded expression was: resRows
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 283:
elimZeroCols_!(A:%):() ==
................^
[L283 C17] #14 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 290:
purge_!(A:%,crit:C->B):() ==
.........^
[L290 C10] #15 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 301:
r := qelt(rA.Rows,l)
......................^
[L301 C23] #22 (Error) (After Macro Expansion) Argument 1 of `qelt' did not match any possible parameter type.
The rejected type is Rows: Vector(Record(Indices: List(C), Entries: List(D))).
Expected one of:
-- List(C)
-- List(D)
-- List(Record(Indices: List(C), Entries: List(D)))
-- List(NonNegativeInteger)
Expanded expression was: rA(Rows)
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 308:
qsetelt!(rA.Rows, l, [reverse! newInds, reverse! newEnts])
.....................^
[L308 C22] #23 (Error) (After Macro Expansion) Argument 1 of `qsetelt!' did not match any possible parameter type.
The rejected type is Rows: Vector(Record(Indices: List(C), Entries: List(D))).
Expected one of:
-- List(C)
-- List(D)
-- List(Record(Indices: List(C), Entries: List(D)))
-- List(NonNegativeInteger)
Expanded expression was: rA(Rows)
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 311:
sortedPurge_!(A:%,crit:C->B):() ==
...............^
[L311 C16] #16 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 318:
r := qelt(rA.Rows,l)
......................^
[L318 C23] #24 (Error) (After Macro Expansion) Argument 1 of `qelt' did not match any possible parameter type.
The rejected type is Rows: Vector(Record(Indices: List(C), Entries: List(D))).
Expected one of:
-- List(C)
-- List(D)
-- List(Record(Indices: List(C), Entries: List(D)))
-- List(NonNegativeInteger)
Expanded expression was: rA(Rows)
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 322:
qsetelt!(rA.Rows,l,r)
.....................^
[L322 C22] #25 (Error) (After Macro Expansion) Argument 1 of `qsetelt!' did not match any possible parameter type.
The rejected type is Rows: Vector(Record(Indices: List(C), Entries: List(D))).
Expected one of:
-- List(C)
-- List(D)
-- List(Record(Indices: List(C), Entries: List(D)))
-- List(NonNegativeInteger)
Expanded expression was: rA(Rows)
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 325:
deleteRow_!(A:%,i:I):() ==
.............^
[L325 C14] #17 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 333:
qsetelt!(resRows,l-1,qelt(rA.Rows,l))
....................................^
[L333 C37] #26 (Error) (After Macro Expansion) Argument 1 of `qelt' did not match any possible parameter type.
The rejected type is Rows: Vector(Record(Indices: List(C), Entries: List(D))).
Expected one of:
-- List(C)
-- List(D)
-- List(Record(Indices: List(C), Entries: List(D)))
-- List(NonNegativeInteger)
Expanded expression was: rA(Rows)
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 338:
consRow_!(A:%,r:ROWREC):() ==
...........^
[L338 C12] #18 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 348:
appendRow_!(A:%,r:ROWREC):() ==
.............^
[L348 C14] #19 (Warning) Escape character ignored. Do you mean '__'?
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 375:
qsetelt!(newRows, l, qelt(rA1.Rows,l))
.....................................^
[L375 C38] #27 (Error) (After Macro Expansion) Argument 1 of `qelt' did not match any possible parameter type.
The rejected type is Rows: Vector(Record(Indices: List(C), Entries: List(D))).
Expected one of:
-- List(C)
-- List(D)
-- List(Record(Indices: List(C), Entries: List(D)))
-- List(NonNegativeInteger)
Expanded expression was: rA1(Rows)
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 377:
qsetelt!(newRows, rA1.NRows+l, qelt(rA2.Rows,l))
...............................................^
[L377 C48] #28 (Error) (After Macro Expansion) Argument 1 of `qelt' did not match any possible parameter type.
The rejected type is Rows: Vector(Record(Indices: List(C), Entries: List(D))).
Expected one of:
-- List(C)
-- List(D)
-- List(Record(Indices: List(C), Entries: List(D)))
-- List(NonNegativeInteger)
Expanded expression was: rA2(Rows)
"/var/zope2/var/LatexWiki/2184467348127326179-25px001.as", line 386:
res:% := new(newInds,rA1.NRows)
...................^
[L386 C20] #29 (Error) (After Macro Expansion) Argument 1 of `new' did not match any possible parameter type.
The rejected type is List(C).
Expected one of:
-- SingleInteger
-- NonNegativeInteger
Expanded expression was: newInds
[L386 C20] #30 (Fatal Error) (After Macro Expansion) Too many errors (use `-M emax=n' or `-M no-emax' to change the limit).
Expanded expression was: ()
The )library system command was not called after compilation.
|