\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}
<<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}
<<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>>
 
<<package MATSTOR StorageEfficientMatrixOperations>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}