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

\documentclass{article}
\usepackage{open-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 MATRIX Matrix}
<<domain MATRIX Matrix>>=
)abbrev domain MATRIX Matrix
++ Author: Grabmeier, Gschnitzer, Williamson
++ Date Created: 1987
++ Date Last Updated: July 1990
++ Basic Operations:
++ Related Domains: 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 : Ring
  Row ==> Vector R
  Col ==> Vector R
  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.

    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 == InnerTwoDimensionalArray(R,Row,Col) add
    minr ==> minRowIndex
    maxr ==> maxRowIndex
    minc ==> minColIndex
    maxc ==> maxColIndex
    mini ==> minIndex
    maxi ==> maxIndex
 
    minRowIndex x == 1
    minColIndex x == 1
 
    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 $
 
    positivePower:($,Integer,NonNegativeInteger) -> $
    positivePower(x,n,nn) ==
      one? n => x
      -- no need to allocate space for 3 additional matrices
      n = 2 => x * x
      n = 3 => x * x * x
      n = 4 => (y := x * x; y * y)
      a := new(nn,nn,0) pretend Matrix(R)
      b := new(nn,nn,0) pretend Matrix(R)
      c := new(nn,nn,0) pretend Matrix(R)
      xx := x pretend Matrix(R)
      power!(a,b,c,xx,n :: NonNegativeInteger)$MATSTOR pretend $
 
    x:$ ** n:NonNegativeInteger ==
      not((nn := nrows x) = ncols x) =>
        error "**: matrix must be square"
      zero? n => scalarMatrix(nn,1)
      positivePower(x,n,nn)
 
    if R has commutative("*") 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
 
        x:$ ** n:Integer ==
          nn := nrows x
          not(nn = ncols x) =>
            error "**: matrix must be square"
          zero? n => scalarMatrix(nn,1)
          positive? n => positivePower(x,n,nn)
          (xInv := inverse x) case "failed" =>
            error "**: matrix must be invertible"
          positivePower(xInv :: $,-n,nn)
 
--     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: 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   : Ring
  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}.
 
  Implementation ==> Matrix R add
    minr ==> minRowIndex
    maxr ==> maxRowIndex
    minc ==> minColIndex
    maxc ==> maxColIndex
    mini ==> minIndex
    maxi ==> maxIndex
 
    0 == per new(m,n,0)$Matrix(R)
 
    coerce(x:$):OutputForm == rep(x)::OutputForm

    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)
      per ans
 
    row(x,i)    == directProduct row(rep x,i)
    column(x,j) == directProduct column(rep x,j)
 
    coerce(x:$):Matrix(R) == copy rep x
 
    rectangularMatrix x ==
      (nrows(x) ~= m) or (ncols(x) ~= n) =>
        error "rectangularMatrix: matrix of bad dimensions"
      per copy(x)
 
    if R has EuclideanDomain then
 
      rowEchelon x == per rowEchelon(rep x)
 
    if R has IntegralDomain then
 
      rank x    == rank rep x
      nullity x == nullity rep x
      nullSpace x ==
        [directProduct c for c in nullSpace rep x]
 
    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: 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    : Ring
  Row ==> DirectProduct(ndim,R)
  Col ==> DirectProduct(ndim,R)
  MATLIN ==> MatrixLinearAlgebraFunctions(R,Row,Col,$)
 
  Exports ==> Join(SquareMatrixCategory(ndim,R,Row,Col),_
                   CoercibleTo Matrix R) with

    new: R -> %
      ++ \spad{new(c)} constructs a new \spadtype{SquareMatrix}
      ++ object of dimension  \spad{ndim} with initial entries equal
      ++ to \spad{c}. 
    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}.
--  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 commutative("*") 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 ConvertibleTo InputForm then ConvertibleTo InputForm
 
  Implementation ==> Matrix R add
    Rep == Matrix R
    minr ==> minRowIndex
    maxr ==> maxRowIndex
    minc ==> minColIndex
    maxc ==> maxColIndex
    mini ==> minIndex
    maxi ==> maxIndex
 
    0 == scalarMatrix 0
    1 == scalarMatrix 1

    characteristic == characteristic$R

    new c == per new(ndim,ndim,c)$Rep
 
    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 := new(ndim,ndim,0)$Rep
      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)
      per ans
 
    row(x,i)    == directProduct row(rep x,i)
    column(x,j) == directProduct column(rep x,j)
    coerce(x:$):OutputForm == rep(x)::OutputForm
 
    scalarMatrix r == per scalarMatrix(ndim,r)$Matrix(R)
 
    diagonalMatrix l ==
      #l ~= ndim =>
        error "diagonalMatrix: wrong number of entries in list"
      per diagonalMatrix(l)$Matrix(R)
 
    coerce(x: %): Matrix(R) == copy rep x
 
    squareMatrix x ==
      (nrows(x) ~= ndim) or (ncols(x) ~= ndim) =>
        error "squareMatrix: matrix of bad dimensions"
      per copy x
 
    x:% * v:Col ==
      directProduct(rep(x) * (v :: Vector(R)))
 
    v:Row * x:$ ==
      directProduct((v :: Vector(R)) * rep(x))
 
    x:$ ** n:NonNegativeInteger ==
      per(rep(x) ** n)
 
    if R has commutative("*") then
 
      determinant x == determinant rep x
      minordet x    == minordet rep x
 
    if R has EuclideanDomain then
 
      rowEchelon x == per rowEchelon rep x
 
    if R has IntegralDomain then
 
      rank x    == rank rep x
      nullity x == nullity rep x
      nullSpace x ==
        [directProduct c for c in nullSpace rep x]
 
    if R has Field then
 
      dimension == (m * n) :: CardinalNumber
 
      inverse x ==
        (u := inverse rep x) case "failed" => "failed"
        per(u :: Matrix(R))
 
      x:$ ** n:Integer ==
        per(rep(x) ** n)
 
      recip x == inverse x
 
    if R has ConvertibleTo InputForm then
      convert(x:$):InputForm ==
         convert [convert('squareMatrix)@InputForm,
                  convert(rep x)@InputForm]$List(InputForm)


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

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