\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra matstor.spad} \author{Clifton J. Williamson} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package MATSTOR StorageEfficientMatrixOperations} <>= )abbrev package MATSTOR StorageEfficientMatrixOperations ++ Author: Clifton J. Williamson ++ Date Created: 18 July 1990 ++ Date Last Updated: 18 July 1990 ++ Basic Operations: ++ Related Domains: Matrix(R) ++ Also See: ++ AMS Classifications: ++ Keywords: matrix, linear algebra ++ Examples: ++ References: ++ Description: ++ This package provides standard arithmetic operations on matrices. ++ The functions in this package store the results of computations ++ in existing matrices, rather than creating new matrices. This ++ package works only for matrices of type Matrix and uses the ++ internal representation of this type. StorageEfficientMatrixOperations(R): Exports == Implementation where R : Ring M ==> Matrix R NNI ==> NonNegativeInteger ARR ==> PrimitiveArray R REP ==> PrimitiveArray PrimitiveArray R Exports ==> with copy! : (M,M) -> M ++ \spad{copy!(c,a)} copies the matrix \spad{a} into the matrix c. ++ Error: if \spad{a} and c do not have the same ++ dimensions. plus! : (M,M,M) -> M ++ \spad{plus!(c,a,b)} computes the matrix sum \spad{a + b} and stores the ++ result in the matrix c. ++ Error: if \spad{a}, b, and c do not have the same dimensions. minus! : (M,M) -> M ++ \spad{minus!(c,a)} computes \spad{-a} and stores the result in the ++ matrix c. ++ Error: if a and c do not have the same dimensions. minus! : (M,M,M) -> M ++ \spad{!minus!(c,a,b)} computes the matrix difference \spad{a - b} ++ and stores the result in the matrix c. ++ Error: if \spad{a}, b, and c do not have the same dimensions. leftScalarTimes! : (M,R,M) -> M ++ \spad{leftScalarTimes!(c,r,a)} computes the scalar product ++ \spad{r * a} and stores the result in the matrix c. ++ Error: if \spad{a} and c do not have the same dimensions. rightScalarTimes! : (M,M,R) -> M ++ \spad{rightScalarTimes!(c,a,r)} computes the scalar product ++ \spad{a * r} and stores the result in the matrix c. ++ Error: if \spad{a} and c do not have the same dimensions. times! : (M,M,M) -> M ++ \spad{times!(c,a,b)} computes the matrix product \spad{a * b} ++ and stores the result in the matrix c. ++ Error: if \spad{a}, b, and c do not have ++ compatible dimensions. power! : (M,M,M,M,NNI) -> M ++ \spad{power!(a,b,c,m,n)} computes m ** n and stores the result in ++ \spad{a}. The matrices b and c are used to store intermediate results. ++ Error: if \spad{a}, b, c, and m are not square ++ and of the same dimensions. ** : (M,NNI) -> M ++ \spad{x ** n} computes the n-th power ++ of a square matrix. The power n is assumed greater than 1. Implementation ==> add rep : M -> REP rep m == m pretend REP copy!(c,a) == m := nrows a; n := ncols a not((nrows c) = m and (ncols c) = n) => error "copy!: matrices of incompatible dimensions" aa := rep a; cc := rep c for i in 0..(m-1) repeat aRow := qelt(aa,i); cRow := qelt(cc,i) for j in 0..(n-1) repeat qsetelt!(cRow,j,qelt(aRow,j)) c plus!(c,a,b) == m := nrows a; n := ncols a not((nrows b) = m and (ncols b) = n) => error "plus!: matrices of incompatible dimensions" not((nrows c) = m and (ncols c) = n) => error "plus!: matrices of incompatible dimensions" aa := rep a; bb := rep b; cc := rep c for i in 0..(m-1) repeat aRow := qelt(aa,i); bRow := qelt(bb,i); cRow := qelt(cc,i) for j in 0..(n-1) repeat qsetelt!(cRow,j,qelt(aRow,j) + qelt(bRow,j)) c minus!(c,a) == m := nrows a; n := ncols a not((nrows c) = m and (ncols c) = n) => error "minus!: matrices of incompatible dimensions" aa := rep a; cc := rep c for i in 0..(m-1) repeat aRow := qelt(aa,i); cRow := qelt(cc,i) for j in 0..(n-1) repeat qsetelt!(cRow,j,-qelt(aRow,j)) c minus!(c,a,b) == m := nrows a; n := ncols a not((nrows b) = m and (ncols b) = n) => error "minus!: matrices of incompatible dimensions" not((nrows c) = m and (ncols c) = n) => error "minus!: matrices of incompatible dimensions" aa := rep a; bb := rep b; cc := rep c for i in 0..(m-1) repeat aRow := qelt(aa,i); bRow := qelt(bb,i); cRow := qelt(cc,i) for j in 0..(n-1) repeat qsetelt!(cRow,j,qelt(aRow,j) - qelt(bRow,j)) c leftScalarTimes!(c,r,a) == m := nrows a; n := ncols a not((nrows c) = m and (ncols c) = n) => error "leftScalarTimes!: matrices of incompatible dimensions" aa := rep a; cc := rep c for i in 0..(m-1) repeat aRow := qelt(aa,i); cRow := qelt(cc,i) for j in 0..(n-1) repeat qsetelt!(cRow,j,r * qelt(aRow,j)) c rightScalarTimes!(c,a,r) == m := nrows a; n := ncols a not((nrows c) = m and (ncols c) = n) => error "rightScalarTimes!: matrices of incompatible dimensions" aa := rep a; cc := rep c for i in 0..(m-1) repeat aRow := qelt(aa,i); cRow := qelt(cc,i) for j in 0..(n-1) repeat qsetelt!(cRow,j,qelt(aRow,j) * r) c copyCol!: (ARR,REP,Integer,Integer) -> Void copyCol!(bCol,bb,j,n1) == for i in 0..n1 repeat qsetelt!(bCol,i,qelt(qelt(bb,i),j)) times!(c,a,b) == m := nrows a; n := ncols a; p := ncols b not((nrows b) = n and (nrows c) = m and (ncols c) = p) => error "times!: matrices of incompatible dimensions" aa := rep a; bb := rep b; cc := rep c bCol : ARR := new(n,0) m1 := (m :: Integer) - 1; n1 := (n :: Integer) - 1 for j in 0..(p-1) repeat copyCol!(bCol,bb,j,n1) for i in 0..m1 repeat aRow := qelt(aa,i); cRow := qelt(cc,i) sum : R := 0 for k in 0..n1 repeat sum := sum + qelt(aRow,k) * qelt(bCol,k) qsetelt!(cRow,j,sum) c power!(a,b,c,m,p) == mm := nrows a; nn := ncols a not(mm = nn) => error "power!: matrix must be square" not((nrows b) = mm and (ncols b) = nn) => error "power!: matrices of incompatible dimensions" not((nrows c) = mm and (ncols c) = nn) => error "power!: matrices of incompatible dimensions" not((nrows m) = mm and (ncols m) = nn) => error "power!: matrices of incompatible dimensions" flag := false copy!(b,m) repeat if odd? p then flag => times!(c,b,a) copy!(a,c) flag := true copy!(a,b) one? p => return a p := p quo 2 times!(c,b,b) copy!(b,c) m ** n == not square? m => error "**: matrix must be square" a := copy m; b := copy m; c := copy m power!(a,b,c,m,n) @ \section{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. @ <<*>>= <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}