path: root/src/algebra/matcat.spad.pamphlet

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra matcat.spad}
\author{Johannes Grabmeier, Oswald Gschnitzer, Clifton J. Williamson}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{category MATCAT MatrixCategory}
<<category MATCAT MatrixCategory>>=
)abbrev category MATCAT MatrixCategory
++ Authors: Grabmeier, Gschnitzer, Williamson
++ Date Created: 1987
++ Date Last Updated: July 1990
++ Basic Operations:
++ Related Domains: Matrix(R)
++ Also See:
++ AMS Classifications:
++ Keywords: matrix, linear algebra
++ Examples:
++ References:
++ Description:
++   \spadtype{MatrixCategory} is a general matrix category which allows
++   different representations and indexing schemes.  Rows and
++   columns may be extracted with rows returned as objects of
++   type Row and colums returned as objects of type Col.
++   A domain belonging to this category will be shallowly mutable.
++   The index of the 'first' row may be obtained by calling the
++   function \spadfun{minRowIndex}.  The index of the 'first' column may
++   be obtained by calling the function \spadfun{minColIndex}.  The index of
++   the first element of a Row is the same as the index of the
++   first column in a matrix and vice versa.
MatrixCategory(R,Row,Col): Category == Definition where
  R   : Ring
  Row : FiniteLinearAggregate R
  Col : FiniteLinearAggregate R
  macro NNI == NonNegativeInteger
  macro I == Integer

  Definition ==> TwoDimensionalArrayCategory(R,Row,Col) with
--% Predicates
     square?  : % -> Boolean
       ++ \spad{square?(m)} returns true if m is a square matrix
       ++ (i.e. if m has the same number of rows as columns) and false otherwise.
     diagonal?: % -> Boolean
       ++ \spad{diagonal?(m)} returns true if the matrix m is square and
       ++ diagonal (i.e. all entries of m not on the diagonal are zero) and
       ++ false otherwise.
     symmetric?: % -> Boolean
       ++ \spad{symmetric?(m)} returns true if the matrix m is square and
       ++ symmetric (i.e. \spad{m[i,j] = m[j,i]} for all i and j) and false
       ++ otherwise.
     antisymmetric?: % -> Boolean
       ++ \spad{antisymmetric?(m)} returns true if the matrix m is square and
       ++ antisymmetric (i.e. \spad{m[i,j] = -m[j,i]} for all i and j) and false
       ++ otherwise.

--% Creation

     zero: (NonNegativeInteger,NonNegativeInteger) -> %
       ++ \spad{zero(m,n)} returns an m-by-n zero matrix.
     matrix: List List R -> %
       ++ \spad{matrix(l)} converts the list of lists l to a matrix, where the
       ++ list of lists is viewed as a list of the rows of the matrix.
     matrix: (NNI, NNI, (I,I) -> R) -> %
       ++ \spad{matrix(n,m,f)} construcys and \spad{n * m} matrix with
       ++ the \spad{(i,j)} entry equal to \spad{f(i,j)}.
     scalarMatrix: (NonNegativeInteger,R) -> %
       ++ \spad{scalarMatrix(n,r)} returns an n-by-n matrix with r's on the
       ++ diagonal and zeroes elsewhere.
     diagonalMatrix: List R -> %
       ++ \spad{diagonalMatrix(l)} returns a diagonal matrix with the elements
       ++ of l on the diagonal.
     diagonalMatrix: List % -> %
       ++ \spad{diagonalMatrix([m1,...,mk])} creates a block diagonal matrix
       ++ M with block matrices {\em m1},...,{\em mk} down the diagonal,
       ++ with 0 block matrices elsewhere.
       ++ More precisly: if \spad{ri := nrows mi}, \spad{ci := ncols mi},
       ++ then m is an (r1+..+rk) by (c1+..+ck) - matrix  with entries
       ++ \spad{m.i.j = ml.(i-r1-..-r(l-1)).(j-n1-..-n(l-1))}, if
       ++ \spad{(r1+..+r(l-1)) < i <= r1+..+rl} and
       ++ \spad{(c1+..+c(l-1)) < i <= c1+..+cl},
       ++ \spad{m.i.j} = 0  otherwise.
     coerce: Col -> %
       ++ \spad{coerce(col)} converts the column col to a column matrix.
     transpose: Row -> %
       ++ \spad{transpose(r)} converts the row r to a row matrix.

--% Creation of new matrices from old

     transpose: % -> %
       ++ \spad{transpose(m)} returns the transpose of the matrix m.
     squareTop: % -> %
       ++ \spad{squareTop(m)} returns an n-by-n matrix consisting of the first
       ++ n rows of the m-by-n matrix m. Error: if
       ++ \spad{m < n}.
     horizConcat: (%,%) -> %
       ++ \spad{horizConcat(x,y)} horizontally concatenates two matrices with
       ++ an equal number of rows. The entries of y appear to the right
       ++ of the entries of x.  Error: if the matrices
       ++ do not have the same number of rows.
     vertConcat: (%,%) -> %
       ++ \spad{vertConcat(x,y)} vertically concatenates two matrices with an
       ++ equal number of columns. The entries of y appear below
       ++ of the entries of x.  Error: if the matrices
       ++ do not have the same number of columns.

--% Part extractions/assignments

     listOfLists: % -> List List R
       ++ \spad{listOfLists(m)} returns the rows of the matrix m as a list
       ++ of lists.
     elt: (%,List Integer,List Integer) -> %
       ++ \spad{elt(x,rowList,colList)} returns an m-by-n matrix consisting
       ++ of elements of x, where \spad{m = # rowList} and \spad{n = # colList}.
       ++ If \spad{rowList = [i<1>,i<2>,...,i<m>]} and \spad{colList =
       ++ [j<1>,j<2>,...,j<n>]}, then the \spad{(k,l)}th entry of
       ++ \spad{elt(x,rowList,colList)} is \spad{x(i<k>,j<l>)}.
     setelt: (%,List Integer,List Integer, %) -> %
       ++ \spad{setelt(x,rowList,colList,y)} destructively alters the matrix x.
       ++ If y is \spad{m}-by-\spad{n}, \spad{rowList = [i<1>,i<2>,...,i<m>]}
       ++ and \spad{colList = [j<1>,j<2>,...,j<n>]}, then \spad{x(i<k>,j<l>)}
       ++ is set to \spad{y(k,l)} for \spad{k = 1,...,m} and \spad{l = 1,...,n}.
     swapRows!: (%,Integer,Integer) -> %
       ++ \spad{swapRows!(m,i,j)} interchanges the \spad{i}th and \spad{j}th
       ++ rows of m. This destructively alters the matrix.
     swapColumns!: (%,Integer,Integer) -> %
       ++ \spad{swapColumns!(m,i,j)} interchanges the \spad{i}th and \spad{j}th
       ++ columns of m. This destructively alters the matrix.
     subMatrix: (%,Integer,Integer,Integer,Integer) -> %
       ++ \spad{subMatrix(x,i1,i2,j1,j2)} extracts the submatrix
       ++ \spad{[x(i,j)]} where the index i ranges from \spad{i1} to \spad{i2}
       ++ and the index j ranges from \spad{j1} to \spad{j2}.
     setsubMatrix!: (%,Integer,Integer,%) -> %
       ++ \spad{setsubMatrix(x,i1,j1,y)} destructively alters the
       ++ matrix x. Here \spad{x(i,j)} is set to \spad{y(i-i1+1,j-j1+1)} for
       ++ \spad{i = i1,...,i1-1+nrows y} and \spad{j = j1,...,j1-1+ncols y}.

--% Arithmetic

     +: (%,%) -> %
       ++ \spad{x + y} is the sum of the matrices x and y.
       ++ Error: if the dimensions are incompatible.
     -: (%,%) -> %
       ++ \spad{x - y} is the difference of the matrices x and y.
       ++ Error: if the dimensions are incompatible.
     -:  %    -> %
       ++ \spad{-x} returns the negative of the matrix x.
     *: (%,%) -> %
       ++ \spad{x * y} is the product of the matrices x and y.
       ++ Error: if the dimensions are incompatible.
     *: (R,%) -> %
       ++ \spad{r*x} is the left scalar multiple of the scalar r and the
       ++ matrix x.
     *: (%,R) -> %
       ++ \spad{x * r} is the right scalar multiple of the scalar r and the
       ++ matrix x.
     *: (Integer,%) -> %
       ++ \spad{n * x} is an integer multiple.
     *: (%,Col) -> Col
       ++ \spad{x * c} is the product of the matrix x and the column vector c.
       ++ Error: if the dimensions are incompatible.
     *: (Row,%) -> Row
       ++ \spad{r * x} is the product of the row vector r and the matrix x.
       ++ Error: if the dimensions are incompatible.
     **: (%,NonNegativeInteger) -> %
       ++ \spad{x ** n} computes a non-negative integral power of the matrix x.
       ++ Error: if the matrix is not square.
     if R has IntegralDomain then
       exquo: (%,R) -> Union(%,"failed")
         ++ \spad{exquo(m,r)} computes the exact quotient of the elements
         ++ of m by r, returning \axiom{"failed"} if this is not possible.
     if R has Field then
       /: (%,R) -> %
         ++ \spad{m/r} divides the elements of m by r. Error: if \spad{r = 0}.

--% Linear algebra

     if R has EuclideanDomain then
       rowEchelon: % -> %
         ++ \spad{rowEchelon(m)} returns the row echelon form of the matrix m.
     if R has IntegralDomain then
       rank: % -> NonNegativeInteger
         ++ \spad{rank(m)} returns the rank of the matrix m.
       nullity: % -> NonNegativeInteger
         ++ \spad{nullity(m)} returns the nullity of the matrix m. This is
         ++ the dimension of the null space of the matrix m.
       nullSpace: % -> List Col
         ++ \spad{nullSpace(m)} returns a basis for the null space of
         ++ the matrix m.
     if R has commutative("*") then
       determinant: % -> R
         ++ \spad{determinant(m)} returns the determinant of the matrix m.
         ++ Error: if the matrix is not square.
       minordet: % -> R
         ++ \spad{minordet(m)} computes the determinant of the matrix m using
         ++ minors. Error: if the matrix is not square.
     if R has Field then
       inverse: % -> Union(%,"failed")
         ++ \spad{inverse(m)} returns the inverse of the matrix m.
         ++ If the matrix is not invertible, "failed" is returned.
         ++ Error: if the matrix is not square.
       **: (%,Integer) -> %
         ++ \spad{m**n} computes an integral power of the matrix m.
         ++ Error: if matrix is not square or if the matrix
         ++ is square but not invertible.

   add
     minr ==> minRowIndex
     maxr ==> maxRowIndex
     minc ==> minColIndex
     maxc ==> maxColIndex
     mini ==> minIndex
     maxi ==> maxIndex

--% Predicates

     square? x == nrows x = ncols x

     diagonal? x ==
       not square? x => false
       for i in minr x .. maxr x repeat
         for j in minc x .. maxc x | (j - minc x) ~= (i - minr x) repeat
           not zero? qelt(x, i, j) => return false
       true

     symmetric? x ==
       (nRows := nrows x) ~= ncols x => false
       mr := minRowIndex x; mc := minColIndex x
       for i in 0..(nRows - 1) repeat
         for j in (i + 1)..(nRows - 1) repeat
           qelt(x,mr + i,mc + j) ~= qelt(x,mr + j,mc + i) => return false
       true

     antisymmetric? x ==
       (nRows := nrows x) ~= ncols x => false
       mr := minRowIndex x; mc := minColIndex x
       for i in 0..(nRows - 1) repeat
         for j in i..(nRows - 1) repeat
           qelt(x,mr + i,mc + j) ~= -qelt(x,mr + j,mc + i) =>
             return false
       true

--% Creation of matrices

     zero(rows,cols) == new(rows,cols,0)

     matrix(l: List List R) ==
       null l => new(0,0,0)
       -- error check: this is a top level function
       rows : NonNegativeInteger := 1; cols := # first l
       cols = 0 => error "matrices with zero columns are not supported"
       for ll in rest l repeat
         cols ~= # ll => error "matrix: rows of different lengths"
         rows := rows + 1
       ans := new(rows,cols,0)
       for i in minr(ans)..maxr(ans) for ll in l repeat
         for j in minc(ans)..maxc(ans) for r in ll repeat
           qsetelt!(ans,i,j,r)
       ans

     matrix(n,m,f) ==
       mat := new(n,m,0)
       for i in minr mat..maxr mat repeat
         for j in minc mat..maxc mat repeat
           qsetelt!(mat,i,j,f(i,j))
       mat

     scalarMatrix(n,r) ==
       ans := zero(n,n)
       for i in minr(ans)..maxr(ans) for j in minc(ans)..maxc(ans) repeat
         qsetelt!(ans,i,j,r)
       ans

     diagonalMatrix(l: List R) ==
       n := #l; ans := zero(n,n)
       for i in minr(ans)..maxr(ans) for j in minc(ans)..maxc(ans) _
           for r in l repeat qsetelt!(ans,i,j,r)
       ans

     diagonalMatrix(list: List %) ==
       rows : NonNegativeInteger := 0
       cols : NonNegativeInteger := 0
       for mat in list repeat
         rows := rows + nrows mat
         cols := cols + ncols mat
       ans := zero(rows,cols)
       loR := minr ans; loC := minc ans
       for mat in list repeat
         hiR := loR + nrows(mat) - 1; hiC := loC + nrows(mat) - 1
         for i in loR..hiR for k in minr(mat)..maxr(mat) repeat
           for j in loC..hiC for l in minc(mat)..maxc(mat) repeat
             qsetelt!(ans,i,j,qelt(mat,k,l))
         loR := hiR + 1; loC := hiC + 1
       ans

     coerce(v:Col) ==
       x := new(#v,1,0)
       one := minc(x)
       for i in minr(x)..maxr(x) for k in mini(v)..maxi(v) repeat
         qsetelt!(x,i,one,qelt(v,k))
       x

     transpose(v:Row) ==
       x := new(1,#v,0)
       one := minr(x)
       for j in minc(x)..maxc(x) for k in mini(v)..maxi(v) repeat
         qsetelt!(x,one,j,qelt(v,k))
       x

     transpose(x:%) ==
       ans := new(ncols x,nrows x,0)
       for i in minr(ans)..maxr(ans) repeat
         for j in minc(ans)..maxc(ans) repeat
           qsetelt!(ans,i,j,qelt(x,j,i))
       ans

     squareTop x ==
       nrows x < (cols := ncols x) =>
         error "squareTop: number of columns exceeds number of rows"
       ans := new(cols,cols,0)
       for i in minr(x)..(minr(x) + cols - 1) repeat
         for j in minc(x)..maxc(x) repeat
           qsetelt!(ans,i,j,qelt(x,i,j))
       ans

     horizConcat(x,y) ==
       (rows := nrows x) ~= nrows y =>
         error "HConcat: matrices must have same number of rows"
       ans := new(rows,(cols := ncols x) + ncols y,0)
       for i in minr(x)..maxr(x) repeat
         for j in minc(x)..maxc(x) repeat
           qsetelt!(ans,i,j,qelt(x,i,j))
       for i in minr(y)..maxr(y) repeat
         for j in minc(y)..maxc(y) repeat
           qsetelt!(ans,i,j + cols,qelt(y,i,j))
       ans

     vertConcat(x,y) ==
       (cols := ncols x) ~= ncols y =>
         error "HConcat: matrices must have same number of columns"
       ans := new((rows := nrows x) + nrows y,cols,0)
       for i in minr(x)..maxr(x) repeat
         for j in minc(x)..maxc(x) repeat
           qsetelt!(ans,i,j,qelt(x,i,j))
       for i in minr(y)..maxr(y) repeat
         for j in minc(y)..maxc(y) repeat
           qsetelt!(ans,i + rows,j,qelt(y,i,j))
       ans

--% Part extraction/assignment

     listOfLists x ==
       ll : List List R := nil()
       for i in maxr(x)..minr(x) by -1 repeat
         l : List R := nil()
         for j in maxc(x)..minc(x) by -1 repeat
           l := cons(qelt(x,i,j),l)
         ll := cons(l,ll)
       ll

     swapRows!(x,i1,i2) ==
       (i1 < minr(x)) or (i1 > maxr(x)) or (i2 < minr(x)) or _
           (i2 > maxr(x)) => error "swapRows!: index out of range"
       i1 = i2 => x
       for j in minc(x)..maxc(x) repeat
         r := qelt(x,i1,j)
         qsetelt!(x,i1,j,qelt(x,i2,j))
         qsetelt!(x,i2,j,r)
       x

     swapColumns!(x,j1,j2) ==
       (j1 < minc(x)) or (j1 > maxc(x)) or (j2 < minc(x)) or _
           (j2 > maxc(x)) => error "swapColumns!: index out of range"
       j1 = j2 => x
       for i in minr(x)..maxr(x) repeat
         r := qelt(x,i,j1)
         qsetelt!(x,i,j1,qelt(x,i,j2))
         qsetelt!(x,i,j2,r)
       x

     elt(x:%,rowList:List Integer,colList:List Integer) ==
       for ei in rowList repeat
         (ei < minr(x)) or (ei > maxr(x)) =>
           error "elt: index out of range"
       for ej in colList repeat
         (ej < minc(x)) or (ej > maxc(x)) =>
           error "elt: index out of range"
       y := new(# rowList,# colList,0)
       for ei in rowList for i in minr(y)..maxr(y) repeat
         for ej in colList for j in minc(y)..maxc(y) repeat
           qsetelt!(y,i,j,qelt(x,ei,ej))
       y

     setelt(x:%,rowList:List Integer,colList:List Integer,y:%) ==
       for ei in rowList repeat
         (ei < minr(x)) or (ei > maxr(x)) =>
           error "setelt: index out of range"
       for ej in colList repeat
         (ej < minc(x)) or (ej > maxc(x)) =>
           error "setelt: index out of range"
       ((# rowList) ~= (nrows y)) or ((# colList) ~= (ncols y)) =>
         error "setelt: matrix has bad dimensions"
       for ei in rowList for i in minr(y)..maxr(y) repeat
         for ej in colList for j in minc(y)..maxc(y) repeat
           qsetelt!(x,ei,ej,qelt(y,i,j))
       y

     subMatrix(x,i1,i2,j1,j2) ==
       (i2 < i1) => error "subMatrix: bad row indices"
       (j2 < j1) => error "subMatrix: bad column indices"
       (i1 < minr(x)) or (i2 > maxr(x)) =>
         error "subMatrix: index out of range"
       (j1 < minc(x)) or (j2 > maxc(x)) =>
         error "subMatrix: index out of range"
       rows := (i2 - i1 + 1) pretend NonNegativeInteger
       cols := (j2 - j1 + 1) pretend NonNegativeInteger
       y := new(rows,cols,0)
       for i in minr(y)..maxr(y) for k in i1..i2 repeat
         for j in minc(y)..maxc(y) for l in j1..j2 repeat
           qsetelt!(y,i,j,qelt(x,k,l))
       y

     setsubMatrix!(x,i1,j1,y) ==
       i2 := i1 + nrows(y) -1
       j2 := j1 + ncols(y) -1
       (i1 < minr(x)) or (i2 > maxr(x)) =>
         error "setsubMatrix!: inserted matrix too big, use subMatrix to restrict it"
       (j1 < minc(x)) or (j2 > maxc(x)) =>
         error "setsubMatrix!: inserted matrix too big, use subMatrix to restrict it"
       for i in minr(y)..maxr(y) for k in i1..i2 repeat
         for j in minc(y)..maxc(y) for l in j1..j2 repeat
           qsetelt!(x,k,l,qelt(y,i,j))
       x

--% Arithmetic

     x + y ==
       ((r := nrows x) ~= nrows y) or ((c := ncols x) ~= ncols y) =>
         error "can't add matrices of different dimensions"
       ans := new(r,c,0)
       for i in minr(x)..maxr(x) repeat
         for j in minc(x)..maxc(x) repeat
           qsetelt!(ans,i,j,qelt(x,i,j) + qelt(y,i,j))
       ans

     x - y ==
       ((r := nrows x) ~= nrows y) or ((c := ncols x) ~= ncols y) =>
         error "can't subtract matrices of different dimensions"
       ans := new(r,c,0)
       for i in minr(x)..maxr(x) repeat
         for j in minc(x)..maxc(x) repeat
           qsetelt!(ans,i,j,qelt(x,i,j) - qelt(y,i,j))
       ans

     - x == map(- #1,x)

     a:R * x:% == map(a * #1,x)
     x:% * a:R == map(#1 * a,x)
     m:Integer * x:% == map(m * #1,x)

     x:% * y:% ==
       (ncols x ~= nrows y) =>
         error "can't multiply matrices of incompatible dimensions"
       ans := new(nrows x,ncols y,0)
       for i in minr(x)..maxr(x) repeat
         for j in minc(y)..maxc(y) repeat
           entry :=
             sum : R := 0
             for k in minr(y)..maxr(y) for l in minc(x)..maxc(x) repeat
               sum := sum + qelt(x,i,l) * qelt(y,k,j)
             sum
           qsetelt!(ans,i,j,entry)
       ans

     positivePower:(%,Integer) -> %
     positivePower(x,n) ==
       one? n => x
       odd? n => x * positivePower(x,n - 1)
       y := positivePower(x,n quo 2)
       y * y

     x:% ** n:NonNegativeInteger ==
       not((nn:= nrows x) = ncols x) => error "**: matrix must be square"
       zero? n => scalarMatrix(nn,1)
       positivePower(x,n)

     --if R has ConvertibleTo InputForm then
       --convert(x:%):InputForm ==
         --convert [convert('matrix)@InputForm,
                  --convert listOfLists x]$List(InputForm)

     if Col has shallowlyMutable then

       x:% * v:Col ==
         ncols(x) ~= #v =>
           error "can't multiply matrix A and vector v if #cols A ~= #v"
         w : Col := new(nrows x,0)
         for i in minr(x)..maxr(x) for k in mini(w)..maxi(w) repeat
           w.k :=
             sum : R := 0
             for j in minc(x)..maxc(x) for l in mini(v)..maxi(v) repeat
               sum := sum + qelt(x,i,j) * v(l)
             sum
         w

     if Row has shallowlyMutable then

       v:Row * x:% ==
         nrows(x) ~= #v =>
           error "can't multiply vector v and matrix A if #rows A ~= #v"
         w : Row := new(ncols x,0)
         for j in minc(x)..maxc(x) for k in mini(w)..maxi(w) repeat
           w.k :=
             sum : R := 0
             for i in minr(x)..maxr(x) for l in mini(v)..maxi(v) repeat
               sum := sum + qelt(x,i,j) * v(l)
             sum
         w

     if R has IntegralDomain then
       x exquo a ==
         ans := new(nrows x,ncols x,0)
         for i in minr(x)..maxr(x) repeat
           for j in minc(x)..maxc(x) repeat
             entry :=
               (r := (qelt(x,i,j) exquo a)) case "failed" =>
                 return "failed"
               r :: R
             qsetelt!(ans,i,j,entry)
         ans

     if R has Field then
       x / r == map(#1 / r,x)

       x:% ** n:Integer ==
         not((nn:= nrows x) = ncols x) => error "**: matrix must be square"
         zero? n => scalarMatrix(nn,1)
         positive? n => positivePower(x,n)
         (xInv := inverse x) case "failed" =>
           error "**: matrix must be invertible"
         positivePower(xInv :: %,-n)

@
\section{category RMATCAT RectangularMatrixCategory}
<<category RMATCAT RectangularMatrixCategory>>=
)abbrev category RMATCAT RectangularMatrixCategory
++ Authors: Grabmeier, Gschnitzer, Williamson
++ Date Created: 1987
++ Date Last Updated: July 1990
++ Basic Operations:
++ Related Domains: RectangularMatrix(m,n,R)
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Examples:
++ References:
++ Description:
++   \spadtype{RectangularMatrixCategory} is a category of matrices of fixed
++   dimensions.  The dimensions of the matrix will be parameters of the
++   domain.  Domains in this category will be R-modules and will be
++   non-mutable.
RectangularMatrixCategory(m,n,R,Row,Col): Category == Definition where
  m,n : NonNegativeInteger
  R   : Ring
  Row : DirectProductCategory(n,R)
  Col : DirectProductCategory(m,R)

  Definition == Join(MatrixCategory(R,Row,Col),BiModule(R,R)) with
    if R has CommutativeRing then Module(R)
--% Arithmetic

    if R has IntegralDomain then
      exquo: (%,R) -> Union(%,"failed")
        ++ \spad{exquo(m,r)} computes the exact quotient of the elements
        ++ of m by r, returning \axiom{"failed"} if this is not possible.
    if R has Field then
      /: (%,R) -> %
        ++ \spad{m/r} divides the elements of m by r. Error: if \spad{r = 0}.

--% Linear algebra

    if R has EuclideanDomain then
      rowEchelon: % -> %
        ++ \spad{rowEchelon(m)} returns the row echelon form of the matrix m.
    if R has IntegralDomain then
      rank: % -> NonNegativeInteger
        ++ \spad{rank(m)} returns the rank of the matrix m.
      nullity: % -> NonNegativeInteger
        ++ \spad{nullity(m)} returns the nullity of the matrix m. This is
        ++ the dimension of the null space of the matrix m.
      nullSpace: % -> List Col
        ++ \spad{nullSpace(m)}+ returns a basis for the null space of
        ++ the matrix m.
   add
     nrows x == m
     ncols x == n
     square? x == m = n

     diagonal? x ==
       not square? x => false
       for i in minRowIndex x .. maxRowIndex x repeat
         for j in minColIndex x .. maxColIndex x
           | (j - minColIndex x) ~= (i - minRowIndex x) repeat
             not zero? qelt(x, i, j) => return false
       true

     symmetric? x ==
       m ~= n => false
       mr := minRowIndex x; mc := minColIndex x
       for i in 0..(n - 1) repeat
         for j in (i + 1)..(n - 1) repeat
           qelt(x,mr + i,mc + j) ~= qelt(x,mr + j,mc + i) => return false
       true

     antisymmetric? x ==
       m ~= n => false
       mr := minRowIndex x; mc := minColIndex x
       for i in 0..(n - 1) repeat
         for j in i..(n - 1) repeat
           qelt(x,mr + i,mc + j) ~= -qelt(x,mr + j,mc + i) =>
             return false
       true

@
\section{category SMATCAT SquareMatrixCategory}
<<category SMATCAT SquareMatrixCategory>>=
)abbrev category SMATCAT SquareMatrixCategory
++ Authors: Grabmeier, Gschnitzer, Williamson
++ Date Created: 1987
++ Date Last Updated: July 1990
++ Basic Operations:
++ Related Domains: SquareMatrix(ndim,R)
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Examples:
++ References:
++ Description:
++   \spadtype{SquareMatrixCategory} is a general square matrix category which
++   allows different representations and indexing schemes.  Rows and
++   columns may be extracted with rows returned as objects of
++   type Row and colums returned as objects of type Col.
SquareMatrixCategory(ndim,R,Row,Col): Category == Definition where
  ndim : NonNegativeInteger
  R    : Ring
  Row  : DirectProductCategory(ndim,R)
  Col  : DirectProductCategory(ndim,R)
  I ==> Integer

  Definition ==> Join(DifferentialExtension R, BiModule(R, R),_
                      RectangularMatrixCategory(ndim,ndim,R,Row,Col),_
                      FullyRetractableTo R,_
                      FullyLinearlyExplicitRingOver R) with
    if R has CommutativeRing then Module(R)
    scalarMatrix: R -> %
      ++ \spad{scalarMatrix(r)} returns an n-by-n matrix with r's on the
      ++ diagonal and zeroes elsewhere.
    diagonalMatrix: List R -> %
      ++ \spad{diagonalMatrix(l)} returns a diagonal matrix with the elements
      ++ of l on the diagonal.
    diagonal: % -> Row
      ++ \spad{diagonal(m)} returns a row consisting of the elements on the
      ++ diagonal of the matrix m.
    trace: % -> R
      ++ \spad{trace(m)} returns the trace of the matrix m. this is the sum
      ++ of the elements on the diagonal of the matrix m.
    diagonalProduct: % -> R
      ++ \spad{diagonalProduct(m)} returns the product of the elements on the
      ++ diagonal of the matrix m.
    *: (%,Col) -> Col
      ++ \spad{x * c} is the product of the matrix x and the column vector c.
      ++ Error: if the dimensions are incompatible.
    *: (Row,%) -> Row
      ++ \spad{r * x} is the product of the row vector r and the matrix x.
      ++ Error: if the dimensions are incompatible.

--% Linear algebra

    if R has commutative("*") then
      Algebra R
      determinant: % -> R
        ++ \spad{determinant(m)} returns the determinant of the matrix m.
      minordet: % -> R
        ++ \spad{minordet(m)} computes the determinant of the matrix m
        ++ using minors.
    if R has Field then
      inverse: % -> Union(%,"failed")
        ++ \spad{inverse(m)} returns the inverse of the matrix m, if that
        ++ matrix is invertible and returns "failed" otherwise.
      **: (%,Integer) -> %
        ++ \spad{m**n} computes an integral power of the matrix m.
        ++ Error: if the matrix is not invertible.

   add
    minr ==> minRowIndex
    maxr ==> maxRowIndex
    minc ==> minColIndex
    maxc ==> maxColIndex
    mini ==> minIndex
    maxi ==> maxIndex

    positivePower:(%,Integer) -> %
    positivePower(x,n) ==
      one? n => x
      odd? n => x * positivePower(x,n - 1)
      y := positivePower(x,n quo 2)
      y * y

    x:% ** n:NonNegativeInteger ==
      zero? n => scalarMatrix 1
      positivePower(x,n)

    coerce(r:R) == scalarMatrix r

    equation2R: Vector % -> Matrix R

    differentiate(x:%,d:R -> R) == map(d,x)

    diagonal x ==
      v:Vector(R) := new(ndim,0)
      for i in minr x .. maxr x
        for j in minc x .. maxc x
          for k in minIndex v .. maxIndex v repeat
            qsetelt!(v, k, qelt(x, i, j))
      directProduct v

    retract(x:%):R ==
      diagonal? x => retract diagonal x
      error "Not retractable"

    retractIfCan(x:%):Union(R, "failed") ==
      diagonal? x => retractIfCan diagonal x
      "failed"

    equation2R v ==
      ans:Matrix(Col) := new(ndim,#v,0)
      for i in minr ans .. maxr ans repeat
        for j in minc ans .. maxc ans repeat
          qsetelt!(ans, i, j, column(qelt(v, j), i))
      reducedSystem ans

    reducedSystem(x:Matrix %):Matrix(R) ==
      empty? x => new(0,0,0)
      reduce(vertConcat, [equation2R row(x, i)
                               for i in minr x .. maxr x])$List(Matrix R)

    reducedSystem(m:Matrix %, v:Vector %):
     Record(mat:Matrix R, vec:Vector R) ==
      vh:Vector(R) :=
        empty? v => new(0,0)
        rh := reducedSystem(v::Matrix %)@Matrix(R)
        column(rh, minColIndex rh)
      [reducedSystem(m)@Matrix(R), vh]

    trace x ==
      tr : R := 0
      for i in minr(x)..maxr(x) for j in minc(x)..maxc(x) repeat
        tr := tr + x(i,j)
      tr

    diagonalProduct x ==
      pr : R := 1
      for i in minr(x)..maxr(x) for j in minc(x)..maxc(x) repeat
        pr := pr * x(i,j)
      pr

    if R has Field then

      x:% ** n:Integer ==
        zero? n => scalarMatrix 1
        positive? n => positivePower(x,n)
        (xInv := inverse x) case "failed" =>
          error "**: matrix must be invertible"
        positivePower(xInv :: %,-n)

@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--Copyright (C) 2007-2013, Gabriel Dos Reis.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>

<<category MATCAT MatrixCategory>>
<<category RMATCAT RectangularMatrixCategory>>
<<category SMATCAT SquareMatrixCategory>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}