matstor.spad

)abbrev package MATSTOR StorageEfficientMatrixOperations
++ Author: Clifton J. Williamson
++ Date Created: 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) -> ARR
    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)
        (p = 1) => 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)

--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.