matrix.spad.pamphlet

\documentclass{article}
\usepackage{axiom}
\begin{document}
\title{\$SPAD/src/algebra matrix.spad}
\author{Johannes Grabmeier, Oswald Gschnitzer, Clifton J. Williamson}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{domain IMATRIX IndexedMatrix}
<<domain IMATRIX IndexedMatrix>>=
)abbrev domain IMATRIX IndexedMatrix
++ Author: Grabmeier, Gschnitzer, Williamson
++ Date Created: 1987
++ Date Last Updated: July 1990
++ Basic Operations:
++ Related Domains: Matrix, RectangularMatrix, SquareMatrix,
++ StorageEfficientMatrixOperations
++ Also See:
++ AMS Classifications:
++ Keywords: matrix, linear algebra
++ Examples:
++ References:
++ Description:
++ An \spad{IndexedMatrix} is a matrix where the minimal row and column
++ indices are parameters of the type. The domains Row and Col
++ are both IndexedVectors.
++ 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.
IndexedMatrix(R,mnRow,mnCol): Exports == Implementation where
  R : Join(SemiRng, AbelianMonoid)
  mnRow, mnCol : Integer
  Row ==> IndexedVector(R,mnCol)
  Col ==> IndexedVector(R,mnRow)
  MATLIN ==> MatrixLinearAlgebraFunctions(R,Row,Col,%)

  Exports ==> MatrixCategory(R,Row,Col)

  Implementation ==>
    InnerIndexedTwoDimensionalArray(R,mnRow,mnCol,Row,Col) add

      swapRows!(x,i1,i2) ==
        (i1 < minRowIndex(x)) or (i1 > maxRowIndex(x)) or _
           (i2 < minRowIndex(x)) or (i2 > maxRowIndex(x)) =>
             error "swapRows!: index out of range"
        i1 = i2 => x
        minRow := minRowIndex x
        xx := x pretend PrimitiveArray(PrimitiveArray(R))
        n1 := i1 - minRow; n2 := i2 - minRow
        row1 := qelt(xx,n1)
        qsetelt!(xx,n1,qelt(xx,n2))
        qsetelt!(xx,n2,row1)
        xx pretend %

      if R has CommutativeRing then

        determinant x == determinant(x)$MATLIN
        minordet x == minordet(x)$MATLIN

      if R has EuclideanDomain then

        rowEchelon x == rowEchelon(x)$MATLIN

      if R has IntegralDomain then

        rank x == rank(x)$MATLIN
        nullity x == nullity(x)$MATLIN
        nullSpace x == nullSpace(x)$MATLIN

      if R has Field then

        inverse x == inverse(x)$MATLIN

@
\section{domain MATRIX Matrix}
<<domain MATRIX Matrix>>=
)abbrev domain MATRIX Matrix
++ Author: Grabmeier, Gschnitzer, Williamson
++ Date Created: 1987
++ Date Last Updated: February 2012
++ Basic Operations:
++ Related Domains: IndexedMatrix, RectangularMatrix, SquareMatrix
++ Also See:
++ AMS Classifications:
++ Keywords: matrix, linear algebra
++ Examples:
++ References:
++ Description:
++ \spadtype{Matrix} is a matrix domain where 1-based indexing is used
++ for both rows and columns.
Matrix(R): Exports == Implementation where
  R : Join(SemiRng, AbelianMonoid)
  Row ==> Vector R
  Col ==> Vector R
  NNI   ==> NonNegativeInteger
  mnRow ==> 1
  mnCol ==> 1
  MATLIN ==> MatrixLinearAlgebraFunctions(R,Row,Col,%)
  MATSTOR ==> StorageEfficientMatrixOperations(R)

  Exports ==> MatrixCategory(R,Row,Col) with
    diagonalMatrix: Vector R -> %
      ++ \spad{diagonalMatrix(v)} returns a diagonal matrix where the elements
      ++ of v appear on the diagonal.
    kroneckerProduct : (%,%) -> %
      ++ kroneckerProduct(a,b) calculates the Kronecker product
      ++ of the matrices {\em a} and b.
      ++ This corresponds to tensor product in other algebras.
      ++ This is used by RepresentationPackage1 and if each matrix
      ++ corresponds to a group representation (repr. of
      ++ generators) of one group, then these matrices
      ++ correspond to the tensor product of the two representations.
    kroneckerSum : (%,%) -> %
      ++ kroneckerSum(a,b) calculates the Kronecker sum
      ++ of the matrices {\em a} and b.
      ++ This corresponds to Cartesian product in other algebras.
      ++ If the Kronecker product (tensorProduct) is denoted
      ++ by A (x) B
      ++ where A is an n-by-n matrix and B is an m-by-m matrix
      ++ then the Kronecker sum (cartesianProduct) is
      ++ A (+) B = A (x) Im + In (x) B
      ++ where In is an n-by-n identity matrix and
      ++ Im is an m-by-m identity matrix
      ++ As defined here kroneckerSum applies only to square
      ++ matrices but we will stretch the definition by padding out
      ++ with zeroes. I would not want to implement this in
      ++ SquareMatrix because the dimension is fixed and therefore
      ++ we would not be able to have signature (%,%) -> %
    kroneckerSum2 : (%,%) -> %
      ++ kroneckerSum(a,b) calculates the Kronecker sum
      ++ of the matrices {\em a} and b.
      ++ This corresponds to Cartesian product in other algebras.
      ++ If the Kronecker product (tensorProduct) is denoted
      ++ by A (x) B
      ++ where A is an n-by-n matrix and B is an m-by-m matrix
      ++ then the Kronecker sum (cartesianProduct) is
      ++ A (+) B = A (x) Im + In (x) B
      ++ where In is an n-by-n identity matrix and
      ++ Im is an m-by-m identity matrix
      ++ As defined here kroneckerSum applies only to square
      ++ matrices and this version enforces that by checking row
      ++ and column dimensions are the same. I would not want to
      ++ implement this in
      ++ SquareMatrix because the dimension is fixed and therefore
      ++ we would not be able to have signature (%,%) -> %
    if R has ConvertibleTo InputForm then ConvertibleTo InputForm

    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.
-- matrix: Vector Vector R -> %
-- ++ \spad{matrix(v)} converts the vector of vectors v to a matrix, where
-- ++ the vector of vectors is viewed as a vector of the rows of the
-- ++ matrix
-- diagonalMatrix: Vector % -> %
-- ++ \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.
-- vectorOfVectors: % -> Vector Vector R
-- ++ \spad{vectorOfVectors(m)} returns the rows of the matrix m as a
-- ++ vector of vectors

  Implementation ==>
   InnerIndexedTwoDimensionalArray(R,mnRow,mnCol,Row,Col) add
    minr ==> minRowIndex
    maxr ==> maxRowIndex
    minc ==> minColIndex
    maxc ==> maxColIndex
    mini ==> minIndex
    maxi ==> maxIndex
    Qelt2 ==> QAREF2O$Lisp
    Qsetelt2 ==> QSETAREF2O$Lisp

    minRowIndex x == mnRow
    minColIndex x == mnCol

    -- qelt and qsetelt! are logically unnecessary, but good for
    -- performance
    qelt(m, i, j) == Qelt2(m, i, j, mnRow@Integer, mnCol@Integer)
    qsetelt!(m,i,j,r) == Qsetelt2(m, i, j, r, mnRow@Integer, mnCol@Integer)

    swapRows!(x,i1,i2) ==
        (i1 < minRowIndex(x)) or (i1 > maxRowIndex(x)) or _
           (i2 < minRowIndex(x)) or (i2 > maxRowIndex(x)) =>
             error "swapRows!: index out of range"
        i1 = i2 => x
        minRow := minRowIndex x
        xx := x pretend PrimitiveArray(PrimitiveArray(R))
        n1 := i1 - minRow; n2 := i2 - minRow
        row1 := qelt(xx,n1)
        qsetelt!(xx,n1,qelt(xx,n2))
        qsetelt!(xx,n2,row1)
        xx pretend %

    kroneckerProduct(a: %, b: %) ==
       n      : NNI := nrows a
       m      : NNI := nrows b
       nc     : NNI := ncols a
       mc     : NNI := ncols b
       c      : %   := zero(n * m, nc * mc)
       indexr : NNI := 1                             --   row index
       for i in 1..n repeat
          for k in 1..m repeat
             indexc : NNI := 1                       --   column index
             for j in 1..nc repeat
                for l in 1..mc repeat
                   c(indexr,indexc) := a(i,j) * b(k,l)
                   indexc := indexc + 1
             indexr := indexr + 1
       c

    kroneckerSum(a:%, b:%) ==
       n : NNI := nrows a
       m : NNI := nrows b
       nc : NNI := ncols a
       mc : NNI := ncols b
       c : % := zero(n * m, nc * mc)
       indexr : NNI := 1                             --   row index
       for i in 1..n repeat
          for k in 1..m repeat
             indexc : NNI := 1                       --   column index
             for j in 1..nc repeat
                for l in 1..mc repeat
                   if i = j then
                     c(indexr,indexc) := c(indexr,indexc) + b(k,l)
                   if k = l then
                     c(indexr,indexc) := c(indexr,indexc) + a(i,j)
                   indexc := indexc + 1
             indexr := indexr + 1
       c

    -- local function used by kroneckerSum2 to calculate row or
    -- column indexes from indexes of individual operands.
    expd(dima:NNI,ai:NNI,bi:NNI):NNI ==
      (subtractIfCan(bi,1::NNI)::NNI)*dima +ai

    kroneckerSum2(a:%, b:%) ==
       n : NNI := nrows a
       m : NNI := nrows b
       if n ~= ncols a then error "kroneckerSum requires square matrix"
       if m ~= ncols b then error "kroneckerSum requires square matrix"
       c : % := zero(n * m, n * m)
       for i in 1..n repeat
         for k in 1..m repeat
           for j in 1..n repeat
             c(expd(m,k,i),expd(m,k,j)) := c(expd(m,k,i),expd(m,k,j)) + a(i,j)             
           for l in 1..m repeat
             c(expd(m,k,i),expd(m,l,i)) := c(expd(m,k,i),expd(m,l,i)) + b(k,l)
       c

    if R has CommutativeRing then

        determinant x == determinant(x)$MATLIN
        minordet x == minordet(x)$MATLIN

    if R has EuclideanDomain then

        rowEchelon x == rowEchelon(x)$MATLIN

    if R has IntegralDomain then

        rank x == rank(x)$MATLIN
        nullity x == nullity(x)$MATLIN
        nullSpace x == nullSpace(x)$MATLIN

    if R has Field then

        inverse x == inverse(x)$MATLIN

-- matrix(v: Vector Vector R) ==
-- (rows := # v) = 0 => new(0,0,0)
-- -- error check: this is a top level function
-- cols := # v.mini(v)
-- for k in (mini(v) + 1)..maxi(v) repeat
-- cols ~= # v.k => error "matrix: rows of different lengths"
-- ans := new(rows,cols,0)
-- for i in minr(ans)..maxr(ans) for k in mini(v)..maxi(v) repeat
-- vv := v.k
-- for j in minc(ans)..maxc(ans) for l in mini(vv)..maxi(vv) repeat
-- ans(i,j) := vv.l
-- ans

    diagonalMatrix(v: Vector R) ==
      n := #v; ans := zero(n,n)
      for i in minr(ans)..maxr(ans) for j in minc(ans)..maxc(ans) _
          for k in mini(v)..maxi(v) repeat qsetelt!(ans,i,j,qelt(v,k))
      ans

-- diagonalMatrix(vec: Vector %) ==
-- rows : NonNegativeInteger := 0
-- cols : NonNegativeInteger := 0
-- for r in mini(vec)..maxi(vec) repeat
-- mat := vec.r
-- rows := rows + nrows mat; cols := cols + ncols mat
-- ans := zero(rows,cols)
-- loR := minr ans; loC := minc ans
-- for r in mini(vec)..maxi(vec) repeat
-- mat := vec.r
-- 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
-- ans(i,j) := mat(k,l)
-- loR := hiR + 1; loC := hiC + 1
-- ans

-- vectorOfVectors x ==
-- vv : Vector Vector R := new(nrows x,0)
-- cols := ncols x
-- for k in mini(vv)..maxi(vv) repeat
-- vv.k := new(cols,0)
-- for i in minr(x)..maxr(x) for k in mini(vv)..maxi(vv) repeat
-- v := vv.k
-- for j in minc(x)..maxc(x) for l in mini(v)..maxi(v) repeat
-- v.l := x(i,j)
-- vv

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

@

\section{domain RMATRIX RectangularMatrix}
<<domain RMATRIX RectangularMatrix>>=
)abbrev domain RMATRIX RectangularMatrix
++ Author: Grabmeier, Gschnitzer, Williamson
++ Date Created: 1987
++ Date Last Updated: July 1990
++ Basic Operations:
++ Related Domains: IndexedMatrix, Matrix, SquareMatrix
++ Also See:
++ AMS Classifications:
++ Keywords: matrix, linear algebra
++ Examples:
++ References:
++ Description:
++ \spadtype{RectangularMatrix} is a matrix domain where the number of rows
++ and the number of columns are parameters of the domain.
RectangularMatrix(m,n,R): Exports == Implementation where
  m,n : NonNegativeInteger
  R : Join(SemiRng, AbelianMonoid)
  Row ==> DirectProduct(n,R)
  Col ==> DirectProduct(m,R)
  Exports ==> Join(RectangularMatrixCategory(m,n,R,Row,Col),_
                   CoercibleTo Matrix R) with

    if R has Field then VectorSpace R

    if R has ConvertibleTo InputForm then ConvertibleTo InputForm

    rectangularMatrix: Matrix R -> %
      ++ \spad{rectangularMatrix(m)} converts a matrix of type \spadtype{Matrix}
      ++ to a matrix of type \spad{RectangularMatrix}.
    coerce: % -> Matrix R
      ++ \spad{coerce(m)} converts a matrix of type \spadtype{RectangularMatrix}
      ++ to a matrix of type \spad{Matrix}.

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

    ZERO := new(m,n,0)$Matrix(R) pretend %
    0 == ZERO

    coerce(x:%):OutputForm == coerce(x pretend Matrix R)$Matrix(R)

    matrix(l: List List R) ==
      -- error check: this is a top level function
      #l ~= m => error "matrix: wrong number of rows"
      for ll in l repeat
        #ll ~= n => error "matrix: wrong number of columns"
      ans : Matrix R := new(m,n,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 pretend %

    row(x,i) == directProduct row(x pretend Matrix(R),i)
    column(x,j) == directProduct column(x pretend Matrix(R),j)

    coerce(x:%):Matrix(R) == copy(x pretend Matrix(R))

    rectangularMatrix x ==
      (nrows(x) ~= m) or (ncols(x) ~= n) =>
        error "rectangularMatrix: matrix of bad dimensions"
      copy(x) pretend %

    if R has EuclideanDomain then

      rowEchelon x == rowEchelon(x pretend Matrix(R)) pretend %

      columnSpace x ==
        [directProduct c for c in columnSpace(x pretend Matrix R)]

    if R has IntegralDomain then

      rank x == rank(x pretend Matrix(R))
      nullity x == nullity(x pretend Matrix(R))
      nullSpace x ==
        [directProduct c for c in nullSpace(x pretend Matrix(R))]

    if R has Field then

      dimension() == (m * n) :: CardinalNumber

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

@
\section{domain SQMATRIX SquareMatrix}
<<domain SQMATRIX SquareMatrix>>=
)abbrev domain SQMATRIX SquareMatrix
++ Author: Grabmeier, Gschnitzer, Williamson
++ Date Created: 1987
++ Date Last Updated: July 1990
++ Basic Operations:
++ Related Domains: IndexedMatrix, Matrix, RectangularMatrix
++ Also See:
++ AMS Classifications:
++ Keywords: matrix, linear algebra
++ Examples:
++ References:
++ Description:
++ \spadtype{SquareMatrix} is a matrix domain of square matrices, where the
++ number of rows (= number of columns) is a parameter of the type.
SquareMatrix(ndim,R): Exports == Implementation where
  ndim : NonNegativeInteger
  R : Join(SemiRng, AbelianMonoid)
  Row ==> DirectProduct(ndim,R)
  Col ==> DirectProduct(ndim,R)
  MATLIN ==> MatrixLinearAlgebraFunctions(R,Row,Col,%)

  Exports ==> Join(SquareMatrixCategory(ndim,R,Row,Col),_
                   CoercibleTo Matrix R) with

    transpose: % -> %
      ++ \spad{transpose(m)} returns the transpose of the matrix m.
    squareMatrix: Matrix R -> %
      ++ \spad{squareMatrix(m)} converts a matrix of type \spadtype{Matrix}
      ++ to a matrix of type \spadtype{SquareMatrix}.
    coerce: % -> Matrix R
      ++ \spad{coerce(m)} converts a matrix of type \spadtype{SquareMatrix}
      ++ to a matrix of type \spadtype{Matrix}.
-- symdecomp : % -> Record(sym:%,antisym:%)
-- ++ \spad{symdecomp(m)} decomposes the matrix m as a sum of a symmetric
-- ++ matrix \spad{m1} and an antisymmetric matrix \spad{m2}. The object
-- ++ returned is the Record \spad{[m1,m2]}
-- if R has CommutativeStar then
-- minorsVect: -> Vector(Union(R,"uncomputed")) --range: 1..2^n-1
-- ++ \spad{minorsVect(m)} returns a vector of the minors of the matrix m
    if R has CommutativeStar and R has unitsKnown then unitsKnown
      ++ the invertible matrices are simply the matrices whose determinants
      ++ are units in the Ring R.
    if R has ConvertibleTo InputForm then ConvertibleTo InputForm

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

    ZERO := scalarMatrix 0
    0 == ZERO

    if R has Monoid then

        ONE := scalarMatrix 1
        1 == ONE

    if R has Ring then

        characteristic() == characteristic()$R

    matrix(l: List List R) ==
      -- error check: this is a top level function
      #l ~= ndim => error "matrix: wrong number of rows"
      for ll in l repeat
        #ll ~= ndim => error "matrix: wrong number of columns"
      ans : Matrix R := new(ndim,ndim,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 pretend %

    row(x,i) == directProduct row(x pretend Matrix(R),i)
    column(x,j) == directProduct column(x pretend Matrix(R),j)
    coerce(x:%):OutputForm == coerce(x pretend Matrix R)$Matrix(R)

    scalarMatrix r == scalarMatrix(ndim,r)$Matrix(R) pretend %

    diagonalMatrix l ==
      #l ~= ndim =>
        error "diagonalMatrix: wrong number of entries in list"
      diagonalMatrix(l)$Matrix(R) pretend %

    coerce(x:%):Matrix(R) == copy(x pretend Matrix(R))

    squareMatrix x ==
      (nrows(x) ~= ndim) or (ncols(x) ~= ndim) =>
        error "squareMatrix: matrix of bad dimensions"
      copy(x) pretend %

    x:% * v:Col ==
      directProduct((x pretend Matrix(R)) * (v :: Vector(R)))

    v:Row * x:% ==
      directProduct((v :: Vector(R)) * (x pretend Matrix(R)))

    if R has CommutativeRing then

      determinant x == determinant(x pretend Matrix(R))
      minordet x == minordet(x pretend Matrix(R))
      Pfaffian x == Pfaffian(x pretend Matrix(R))

    if R has EuclideanDomain then

      rowEchelon x == rowEchelon(x pretend Matrix(R)) pretend %

    if R has IntegralDomain then

      rank x == rank(x pretend Matrix(R))
      nullity x == nullity(x pretend Matrix(R))
      nullSpace x ==
        [directProduct c for c in nullSpace(x pretend Matrix(R))]

    if R has Field then

      dimension() == (m * n) :: CardinalNumber

      inverse x ==
        (u := inverse(x pretend Matrix(R))) case "failed" => "failed"
        (u :: Matrix(R)) pretend %

      x:% ^ n:Integer ==
        ((x pretend Matrix(R)) ^ n) pretend %

      recip x == inverse x

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


@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--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>>

<<domain IMATRIX IndexedMatrix>>
<<domain MATRIX Matrix>>
<<domain RMATRIX RectangularMatrix>>
<<domain SQMATRIX SquareMatrix>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}