Edit detail for Aldor Jet Sparse revision 1 of 1

changed:
-
Row echelon form and other basic operations for sparse matrices.

LU decomposition for ordinary matrices.

(C) 1995  Werner M. Seiler (axiom-xl version)

\begin{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]
\end{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]

   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.