\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} <>= )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 : Ring 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 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 @ \section{domain MATRIX Matrix} <>= )abbrev domain MATRIX Matrix ++ Author: Grabmeier, Gschnitzer, Williamson ++ Date Created: 1987 ++ Date Last Updated: July 1990 ++ 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 : Ring Row ==> Vector R Col ==> Vector R 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. 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 minRowIndex x == mnRow minColIndex x == mnCol 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 (n = 1) => 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} <>= )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 : 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 ZERO := per new(m,n,0)$Matrix(R) 0 == ZERO 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} <>= )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 : 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 commutative("*") then central ++ the elements of the Ring R, viewed as diagonal matrices, commute ++ with all matrices and, indeed, are the only matrices which commute ++ with all matrices. if R has commutative("*") 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 Rep == Matrix R minr ==> minRowIndex maxr ==> maxRowIndex minc ==> minColIndex maxc ==> maxColIndex mini ==> minIndex maxi ==> maxIndex ZERO := scalarMatrix 0 0 == ZERO ONE := scalarMatrix 1 1 == ONE 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} <>= --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. @ <<*>>= <> <> <> <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}