|
SparseEchelonMatrix (SEM)
Row echelon form and other basic operations for sparse matrices.
spad )abb domain SEM SparseEchelonMatrix
Sy ==> Symbol
L ==> List
V ==> Vector
VD ==> Vector D
MD ==> Matrix D
FD ==> Fraction D
MFD ==> Matrix FD
I ==> Integer
NNI ==> NonNegativeInteger
B ==> Boolean
OUT ==> OutputForm
ROWREC ==> Record(Indices:L C, Entries:L D)
iter ==> "iterated"::Sy
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 : $ -> MD
++ \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_! : $ -> Void
++ \axiom{elimZeroCols!(A)} removes columns which contain only zeros.
++ This effects basically only the value of \axiom{allIndices(A)}.
purge_! : ($,C->B) -> Void
++ \axiom{purge!(A,crit)} eliminates all columns belonging to an index
++ \axiom{c} such that \axiom{crit(c)} yields \axiom{true}.
sortedPurge_! : ($,C->B) -> Void
++ \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) -> Void
++ \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) -> Void
++ \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) -> Void
++ \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) -> Void
++ \axiom{deleteRow(A,i)} deletes the \axiom{i}-th row of the matrix
++ \axiom{A}.
consRow_! : ($,ROWREC) -> Void
++ \axiom{consRow!(A,r)} inserts the row \axiom{r} at the top of the
++ matrix \axiom{A}.
appendRow_! : ($,ROWREC) -> Void
++ \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:MD, 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.
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}.
primitiveRowEchelon : $ -> Record(Ech:$, Lt:MFD, 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.
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.
"*" : (MD,$) -> $
++ \axiom{L*A} implements left multiplication with a usual matrix.
if D has IntegralDomain then
"*" : (MFD,$) -> $
++ \axiom{L*A} implements left multiplication with a usual matrix over
++ the quotient field of \axiom{D}.
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}.
Def ==> add
minInd:I := minIndex([i for i in 1..1])
offset:I := minInd-1
emptyRec:ROWREC := [empty,empty]
noChecks?:B := D has lazyRep -- flag for lazy representation
seed:I := 113 -- seed for random number generation
GCDmode:Sy := iter -- flag for gcd algorithm
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.
Rep := Record(NCols:NNI, NRows:NNI, AllInds:L C, Rows:V ROWREC)
ncols(A:$):NNI == A.NCols
nrows(A:$):NNI == A.NRows
allIndices(A:$):L C == copy A.AllInds
row(A:$,i:I):ROWREC ==
-- i<0 or i>A.NRows => error "index out of range"
qelt(A.Rows,i)
setRow_!(A:$,i:I,r:ROWREC):Void ==
-- i<0 or i>A.NRows => error "index out of range"
qsetelt!(A.Rows,i,r)
void
setRow_!(A:$,i:I,inds:L C,ents:L D):Void ==
-- i<0 or i>A.NRows => error "index out of range"
-- #inds ^= #ents => error "improper row"
qsetelt!(A.Rows, i, [inds,ents])
void
new(inds:L C,n:I):$ ==
[#inds, n::NNI, inds, [copy 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):Void ==
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!(A.Rows,i,r)
void
coerce(A:$):MD ==
zero? A.NCols => error "cannot coerce matrix with zero columns"
AA:MD := new(A.NRows,A.NCols,0$D)
for r in entries(A.Rows) for i in minRowIndex(AA).. repeat
inds := r.Indices
ents := r.Entries
for ind in A.AllInds 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? A.NCols => 0$D ::OUT
A::MD::OUT
copy(A:$):$ ==
resRows:V ROWREC := new(A.NRows,emptyRec)
for l in 1..A.NRows repeat
r := qelt(A.Rows,l)
qsetelt!(resRows, l, [copy r.Indices, copy r.Entries])
[A.NCols, A.NRows, copy A.AllInds, resRows]
-- ----------------------- --
-- Basic Matrix Operations --
-- ----------------------- --
elimZeroCols_!(A:$):Void ==
newInds:L C := empty
for r in entries(A.Rows) repeat
newInds := removeDuplicates! merge!(#2<#1, newInds, r.Indices)
A.AllInds := newInds
void
purge_!(A:$,crit:C->B):Void ==
newInds:L C := empty
for c in A.AllInds repeat
if not crit c then
newInds := cons(c,newInds)
newInds := reverse! newInds
if #newInds ^= #A.AllInds then
A.AllInds := newInds
for l in 1..A.NRows repeat
r := qelt(A.Rows,l)
newInds:L C := empty
newEnts:L D := empty
for c in r.Indices for e in r.Entries repeat
if not crit c then
newInds := cons(c,newInds)
newEnts := cons(e,newEnts)
qsetelt!(A.Rows, l, [reverse! newInds, reverse! newEnts])
void
sortedPurge_!(A:$,crit:C->B):Void ==
if crit first A.AllInds then
while not(empty? A.AllInds) and crit first A.AllInds repeat
A.AllInds := rest A.AllInds
for l in 1..A.NRows repeat
r := qelt(A.Rows,l)
while not(empty? r.Indices) and crit first r.Indices repeat
r.Indices := rest r.Indices
r.Entries := rest r.Entries
qsetelt!(A.Rows,l,r)
void
deleteRow_!(A:$,i:I):Void ==
i>A.NRows => A
nr := (A.NRows-1)::NNI
resRows:V ROWREC := new(nr,emptyRec)
for l in 1..(i-1) repeat
qsetelt!(resRows,l,qelt(A.Rows,l))
for l in (i+1)..A.NRows repeat
qsetelt!(resRows,l-1,qelt(A.Rows,l))
A.NRows := nr
A.Rows := resRows
void
consRow_!(A:$,r:ROWREC):Void ==
A.NRows := A.NRows + 1
newRows:L ROWREC := cons(r,entries A.Rows)
A.Rows := construct newRows
newInds := setDifference(r.Indices,A.AllInds)
if not empty? newInds then
A.AllInds := merge(#2<#1,A.AllInds,sort!(#2<#1,newInds))
void
appendRow_!(A:$,r:ROWREC):Void ==
A.NRows := A.NRows + 1
newRows:L ROWREC := concat(entries A.Rows, r)
A.Rows := construct newRows
newInds := setDifference(r.Indices,A.AllInds)
if not empty? newInds then
A.AllInds := merge(#2<#1,A.AllInds,sort!(#2<#1,newInds))
void
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(#2<#1,newInds,row(A,i).Indices)
[A.NCols, nr, newInds, resRows]
join(A1:$,A2:$):$ ==
newInds := removeDuplicates! merge(#2<#1,A1.AllInds,A2.AllInds)
newNRows := A1.NRows + A2.NRows
newRows:V ROWREC := new(newNRows,emptyRec)
for l in 1..A1.NRows repeat
qsetelt!(newRows, l, qelt(A1.Rows,l))
for l in 1..A2.NRows repeat
qsetelt!(newRows, A1.NRows+l, qelt(A2.Rows,l))
[#newInds, newNRows, newInds, newRows]
horizJoin(A1:$,A2:$):$ ==
A1.NRows^=A2.NRows => error "incompatible dimensions in horizJoin"
newInds := append(A1.AllInds,A2.AllInds)
res:$ := new(newInds,A1.NRows)
for i in 1..A1.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:$) ==
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,A.NRows), A]
linds := reverse! linds
empty? rinds => [A, new(rinds,A.NRows)]
LA:$ := new(linds,A.NRows)
RA:$ := new(rinds,A.NRows)
for i in 1..A.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
lent1:L D
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 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:MD, Pivots:L D, Rank:NNI) ==
A := copy AA
LTr:MD := 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)
sorted?:B := false
until sorted? repeat
sorted? := true
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)
sorted? := false
else
oldr := newr
-- fraction-free elimination
finished?:B := false
pivlen,pivrow,rk:NNI
for i in 1..A.NRows until finished? repeat
r := qelt(A.Rows,i)
finished? := empty? r.Indices
if finished? then
rk:NNI := (i-1)::NNI
else -- 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
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 finished? then
rk:NNI := A.NRows
[A, LTr, Pivs, rk]
if D has GcdDomain then
setGcdMode(s:Sy):Sy ==
tmp := GCDmode
(s=iter) or (s=rand) =>
GCDmode := s
tmp
error "unknown gcd mode"
randomGCD(le:L D):D ==
-- Probabilistic technique.
#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(113)$I)*first(l)
l := rest l
if not empty? l then
one? first l => return 1$D
g := g + (1+random(113)$I)*first(l)
l := rest l
h := gcd(f,g)
l := [h]
for e in le repeat
tmp := e exquo 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
le := r.Entries
one?(#le) => [first le, [r.Indices,[1$D]]]
g:D
if GCDmode=iterated then
g := iteratedGCD le
else
g := randomGCD le
one? g => [1,r]
le := [(e exquo g)::D for e in le]
[g, [r.Indices,le]]
primitiveRowEchelon(AA:$) : _
Record(Ech:$, Lt:MFD, Pivots:L D, Rank:NNI) ==
A := copy AA
LTr:MFD := diagonalMatrix [1$FD 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)
sorted?:B := false
until sorted? repeat
sorted? := true
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)
sorted? := false
else
oldr := newr
-- primitive fraction-free elimination
finished?:B := false
pivlen,pivrow,rk:NNI
for i in 1..A.NRows until finished? repeat
r := qelt(A.Rows,i)
finished? := empty? r.Indices
if finished? then
rk:NNI := (i-1)::NNI
else -- 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 finished? then
rk:NNI := A.NRows
[A, LTr, Pivs, rk]
-- -------------- --
-- Multiplication --
-- -------------- --
L:MD * AA:$ ==
ncols(L)^=AA.NRows => error "improper matrix dimensions"
A := copy AA
rlen := nrows L
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 := qelt(A.Rows,i)
inds := r.Indices
if not(empty? inds) and first(inds)=c then
for k in 1..rlen | not zero? qelt(L,k,i) repeat
qsetelt!(tmp, k, qelt(tmp,k) + qelt(L,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)
res
if D has IntegralDomain then
mult(f:FD,d:D):D ==
res := numer(f)*d
tmp := res exquo denom(f)
tmp case "failed" => error "cannot divide in mult"
tmp::D
L:MFD * AA:$ ==
ncols(L)^=AA.NRows => error "improper matrix dimensions"
A := copy AA
rlen := nrows L
res:$ := new(A.AllInds,rlen)
for c in A.AllInds repeat
tmp:V FD := 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(L,k,i) repeat
qsetelt!(tmp, k, qelt(tmp,k) + qelt(L,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(D,"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)
res
spad Compiling FriCAS source code from file
/var/zope2/var/LatexWiki/8860614091544751208-25px001.spad using
old system compiler.
SEM abbreviates domain SparseEchelonMatrix
processing macro definition Sy ==> Symbol
processing macro definition L ==> List
processing macro definition V ==> Vector
processing macro definition VD ==> Vector D
processing macro definition MD ==> Matrix D
processing macro definition FD ==> Fraction D
processing macro definition MFD ==> Matrix FD
processing macro definition I ==> Integer
processing macro definition NNI ==> NonNegativeInteger
processing macro definition B ==> Boolean
processing macro definition OUT ==> OutputForm
processing macro definition ROWREC ==> Record(Indices: L C,Entries: L D)
processing macro definition iter ==> ::(iterated,Sy)
processing macro definition rand ==> ::(random,Sy)
processing macro definition Cat ==> -- the constructor category
processing macro definition Def ==> -- the constructor capsule
------------------------------------------------------------------------
initializing NRLIB SEM for SparseEchelonMatrix
compiling into NRLIB SEM
augmenting D: (ATTRIBUTE D lazyRep)
compiling local greater : (Record(Indices: List C,Entries: List D),Record(Indices: List C,Entries: List D)) -> Boolean
Time: 0.21 SEC.
compiling exported ncols : $ -> NonNegativeInteger
SEM;ncols;$Nni;2 is replaced by QVELTA0
Time: 0.01 SEC.
compiling exported nrows : $ -> NonNegativeInteger
SEM;nrows;$Nni;3 is replaced by QVELTA1
Time: 0 SEC.
compiling exported allIndices : $ -> List C
Time: 0 SEC.
compiling exported row : ($,Integer) -> Record(Indices: List C,Entries: List D)
Time: 0.04 SEC.
compiling exported setRow! : ($,Integer,Record(Indices: List C,Entries: List D)) -> Void
Time: 0.06 SEC.
compiling exported setRow! : ($,Integer,List C,List D) -> Void
Time: 0.04 SEC.
compiling exported new : (List C,Integer) -> $
Time: 0.14 SEC.
compiling exported elt : ($,Integer,C) -> D
Time: 0.01 SEC.
compiling exported setelt! : ($,Integer,C,D) -> Void
Time: 0.06 SEC.
compiling exported coerce : $ -> Matrix D
Time: 0.19 SEC.
compiling exported coerce : $ -> OutputForm
Time: 0.04 SEC.
compiling exported copy : $ -> $
Time: 0.05 SEC.
compiling exported elimZeroCols! : $ -> Void
Time: 0.20 SEC.
compiling exported purge! : ($,C -> Boolean) -> Void
Time: 0.08 SEC.
compiling exported sortedPurge! : ($,C -> Boolean) -> Void
Time: 0.07 SEC.
compiling exported deleteRow! : ($,Integer) -> Void
Time: 0.07 SEC.
compiling exported consRow! : ($,Record(Indices: List C,Entries: List D)) -> Void
Time: 0.11 SEC.
compiling exported appendRow! : ($,Record(Indices: List C,Entries: List D)) -> Void
Time: 0.21 SEC.
compiling exported extract : ($,Integer,Integer) -> $
Time: 0.05 SEC.
compiling exported join : ($,$) -> $
Time: 0.06 SEC.
compiling exported horizJoin : ($,$) -> $
Time: 0.04 SEC.
compiling exported horizSplit : ($,C) -> Record(Left: $,Right: $)
Time: 0.05 SEC.
compiling local addRows : (D,Record(Indices: List C,Entries: List D),D,Record(Indices: List C,Entries: List D)) -> Record(Indices: List C,Entries: List D)
Time: 0.09 SEC.
compiling exported pivot : ($,Integer) -> Record(Index: C,Entry: D)
Time: 0.11 SEC.
compiling exported pivots : $ -> Record(Indices: List C,Entries: List D)
Time: 0.12 SEC.
compiling exported rowEchelon : $ -> Record(Ech: $,Lt: Matrix D,Pivots: List D,Rank: NonNegativeInteger)
Time: 1.50 SEC.
****** Domain: D already in scope
augmenting D: (GcdDomain)
augmenting $: (SIGNATURE $ setGcdMode ((Symbol) (Symbol)))
augmenting $: (SIGNATURE $ primitiveRowEchelon ((Record (: Ech $) (: Lt (Matrix (Fraction D))) (: Pivots (List D)) (: Rank (NonNegativeInteger))) $))
compiling exported setGcdMode : Symbol -> Symbol
Time: 0.05 SEC.
compiling local randomGCD : List D -> D
Time: 0.07 SEC.
compiling local iteratedGCD : List D -> D
Time: 0.01 SEC.
compiling local makePrimitive : Record(Indices: List C,Entries: List D) -> Record(GCD: D,Row: Record(Indices: List C,Entries: List D))
Time: 0.03 SEC.
compiling exported primitiveRowEchelon : $ -> Record(Ech: $,Lt: Matrix Fraction D,Pivots: List D,Rank: NonNegativeInteger)
Time: 2.09 SEC.
compiling exported * : (Matrix D,$) -> $
Time: 0.36 SEC.
****** Domain: D already in scope
augmenting D: (IntegralDomain)
augmenting $: (SIGNATURE $ * ($ (Matrix (Fraction D)) $))
compiling local mult : (Fraction D,D) -> D
Time: 0.15 SEC.
compiling exported * : (Matrix Fraction D,$) -> $
Time: 0.34 SEC.
****** Domain: D already in scope
augmenting D: (GcdDomain)
augmenting $: (SIGNATURE $ setGcdMode ((Symbol) (Symbol)))
augmenting $: (SIGNATURE $ primitiveRowEchelon ((Record (: Ech $) (: Lt (Matrix (Fraction D))) (: Pivots (List D)) (: Rank (NonNegativeInteger))) $))
****** Domain: D already in scope
augmenting D: (IntegralDomain)
augmenting $: (SIGNATURE $ * ($ (Matrix (Fraction D)) $))
(time taken in buildFunctor: 1)
;;; *** |SparseEchelonMatrix| REDEFINED
;;; *** |SparseEchelonMatrix| REDEFINED
Time: 0.01 SEC.
Warnings:
[1] coerce: not known that (Ring) is of mode (CATEGORY D (ATTRIBUTE lazyRep))
[2] purge!: newInds has no value
[3] purge!: newEnts has no value
[4] addRows: resI has no value
[5] addRows: resE has no value
[6] addRows: lent2 has no value
[7] rowEchelon: pivrow has no value
[8] rowEchelon: rk has no value
[9] makePrimitive: iterated has no value
[10] primitiveRowEchelon: pivrow has no value
[11] primitiveRowEchelon: rk has no value
Cumulative Statistics for Constructor SparseEchelonMatrix
Time: 6.72 seconds
finalizing NRLIB SEM
Processing SparseEchelonMatrix for Browser database:
--------(shallowlyMutable (attribute))---------
--------(finiteAggregate (attribute))---------
--------(coerce (MD $))---------
--------(copy ($ $))---------
--------(ncols (NNI $))---------
--------(nrows (NNI $))---------
--------(allIndices ((L C) $))---------
--------(elimZeroCols! ((Void) $))---------
--------(purge! ((Void) $ (Mapping B C)))---------
--------(sortedPurge! ((Void) $ (Mapping B C)))---------
--------(new ($ (L C) I))---------
--------(elt (D $ I C))---------
--------(setelt! ((Void) $ I C D))---------
--------(row (ROWREC $ I))---------
--------(setRow! ((Void) $ I ROWREC))---------
--------(setRow! ((Void) $ I (L C) (L D)))---------
--------(deleteRow! ((Void) $ I))---------
--------(consRow! ((Void) $ ROWREC))---------
--------(appendRow! ((Void) $ ROWREC))---------
--------(extract ($ $ I I))---------
--------(rowEchelon ((Record (: Ech $) (: Lt MD) (: Pivots (L D)) (: Rank NNI)) $))---------
--------(setGcdMode (Sy Sy))---------
--------(primitiveRowEchelon ((Record (: Ech $) (: Lt MFD) (: Pivots (L D)) (: Rank NNI)) $))---------
--------(pivot ((Record (: Index C) (: Entry D)) $ I))---------
--------(pivots (ROWREC $))---------
--------(* ($ MD $))---------
--------(* ($ MFD $))---------
--------(join ($ $ $))---------
--------(horizJoin ($ $ $))---------
--------(horizSplit ((Record (: Left $) (: Right $)) $ C))---------
--------constructor---------
------------------------------------------------------------------------
SparseEchelonMatrix is now explicitly exposed in frame initial
SparseEchelonMatrix will be automatically loaded when needed from
/var/zope2/var/LatexWiki/SEM.NRLIB/code
|