1)abbrev domain IMATRIX IndexedMatrix
2++ Author: Grabmeier, Gschnitzer, Williamson
3++ Date Created: 1987
4++ Basic Operations:
5++ Related Domains: Matrix, RectangularMatrix, SquareMatrix,
6++   StorageEfficientMatrixOperations
7++ Also See:
8++ AMS Classifications:
9++ Keywords: matrix, linear algebra
10++ Examples:
11++ References:
12++ Description:
13++   An \spad{IndexedMatrix} is a matrix where the minimal row and column
14++   indices are parameters of the type.  The domains Row and Col
15++   are both IndexedVectors.
16++   The index of the 'first' row may be obtained by calling the
17++   function \spadfun{minRowIndex}.  The index of the 'first' column may
18++   be obtained by calling the function \spadfun{minColIndex}.  The index of
19++   the first element of a 'Row' is the same as the index of the
20++   first column in a matrix and vice versa.
21IndexedMatrix(R, mnRow, mnCol) : Exports == Implementation where
22  R : AbelianMonoid
23  mnRow, mnCol : Integer
24  Row ==> IndexedVector(R, mnCol)
25  Col ==> IndexedVector(R, mnRow)
26  MATLIN ==> MatrixLinearAlgebraFunctions(R, Row, Col, %)
27
28  Exports ==> MatrixCategory(R, Row, Col)
29
30  Implementation ==>
31    InnerIndexedTwoDimensionalArray(R, mnRow, mnCol, Row, Col) add
32
33      Qelt2 ==> QAREF2O$Lisp
34      Qsetelt2 ==> QSETAREF2O$Lisp
35
36      swapRows!(x, i1, i2) ==
37        (i1 < minRowIndex(x)) or (i1 > maxRowIndex(x)) or _
38           (i2 < minRowIndex(x)) or (i2 > maxRowIndex(x)) =>
39             error "swapRows!: index out of range"
40        i1 = i2 => x
41        ro := mnRow
42        co := mnCol
43        for j in co..maxColIndex(x) repeat
44            t1 : R := Qelt2(x, i1, j, ro, co)
45            t2 : R := Qelt2(x, i2, j, ro, co)
46            Qsetelt2(x, i1, j, t2, ro, co)
47            Qsetelt2(x, i2, j, t1, ro, co)
48        x
49
50      if R has CommutativeRing then
51
52        determinant x == determinant(x)$MATLIN
53        minordet    x == minordet(x)$MATLIN
54
55      if R has EuclideanDomain then
56
57        rowEchelon  x == rowEchelon(x)$MATLIN
58
59      if R has IntegralDomain then
60
61        rank        x == rank(x)$MATLIN
62        nullity     x == nullity(x)$MATLIN
63        nullSpace   x == nullSpace(x)$MATLIN
64
65      if R has Field then
66
67        inverse     x == inverse(x)$MATLIN
68
69)abbrev domain MATRIX Matrix
70++ Author: Grabmeier, Gschnitzer, Williamson
71++ Date Created: 1987
72++ Basic Operations:
73++ Related Domains: IndexedMatrix, RectangularMatrix, SquareMatrix
74++ Also See:
75++ AMS Classifications:
76++ Keywords: matrix, linear algebra
77++ Examples:
78++ References:
79++ Description:
80++   \spadtype{Matrix} is a matrix domain where 1-based indexing is used
81++   for both rows and columns.
82Matrix(R) : Exports == Implementation where
83  R : AbelianMonoid
84  Row ==> Vector R
85  Col ==> Vector R
86  mnRow ==> 1
87  mnCol ==> 1
88  MATLIN ==> MatrixLinearAlgebraFunctions(R, Row, Col, %)
89
90  Exports ==> MatrixCategory(R, Row, Col) with
91    diagonalMatrix : Vector R -> %
92      ++ \spad{diagonalMatrix(v)} returns a diagonal matrix where the elements
93      ++ of v appear on the diagonal.
94
95    if R has ConvertibleTo InputForm then ConvertibleTo InputForm
96
97    if R has IntegralDomain then
98      invertIfCan : % -> Union(%,"failed")
99        ++ \spad{invertIfCan(m)} returns the inverse of the matrix m.
100        ++ If the matrix is not invertible, "failed" is returned.
101        ++ Error: if the matrix is not square.
102--     matrix: Vector Vector R -> %
103--       ++ \spad{matrix(v)} converts the vector of vectors v to a matrix, where
104--       ++ the vector of vectors is viewed as a vector of the rows of the
105--       ++ matrix
106--     diagonalMatrix: Vector % -> %
107--       ++ \spad{diagonalMatrix([m1, ..., mk])} creates a block diagonal matrix
108--       ++ M with block matrices {\em m1}, ..., {\em mk} down the diagonal,
109--       ++ with 0 block matrices elsewhere.
110--     vectorOfVectors: % -> Vector Vector R
111--       ++ \spad{vectorOfVectors(m)} returns the rows of the matrix m as a
112--       ++ vector of vectors
113
114  Implementation ==>
115   InnerIndexedTwoDimensionalArray(R, mnRow, mnCol, Row, Col) add
116    minr ==> minRowIndex
117    maxr ==> maxRowIndex
118    minc ==> minColIndex
119    maxc ==> maxColIndex
120    mini ==> minIndex
121    maxi ==> maxIndex
122    Qelt2 ==> QAREF2O$Lisp
123    Qsetelt2 ==> QSETAREF2O$Lisp
124
125    import from List(List(R))
126
127    minRowIndex x == mnRow
128    minColIndex x == mnCol
129
130    -- qelt, qsetelt!, swapRows! and copy are logically unnecessary,
131    -- but good for performance
132
133    qelt(m, i, j) == Qelt2(m, i, j, mnRow@Integer, mnCol@Integer)
134    qsetelt!(m, i, j, r) == Qsetelt2(m, i, j, r, mnRow@Integer, mnCol@Integer)
135
136    swapRows!(x, i1, i2) ==
137        (i1 < minRowIndex(x)) or (i1 > maxRowIndex(x)) or _
138           (i2 < minRowIndex(x)) or (i2 > maxRowIndex(x)) =>
139             error "swapRows!: index out of range"
140        i1 = i2 => x
141        for j in minColIndex(x)..maxColIndex(x) repeat
142            t1 : R := Qelt2(x, i1, j, mnRow@Integer, mnCol@Integer)
143            t2 : R := Qelt2(x, i2, j, mnRow@Integer, mnCol@Integer)
144            Qsetelt2(x, i1, j, t2, mnRow@Integer, mnCol@Integer)
145            Qsetelt2(x, i2, j, t1, mnRow@Integer, mnCol@Integer)
146        x
147
148    copy m ==
149        ans : % := MAKE_MATRIX(nrows m, ncols m)$Lisp
150        for i in minRowIndex(m)..maxRowIndex(m) repeat
151            for j in minColIndex(m)..maxColIndex(m) repeat
152                qsetelt!(ans, i, j, qelt(m, i, j))
153        ans
154
155    if R has CommutativeRing then
156
157        determinant x == determinant(x)$MATLIN
158        minordet    x == minordet(x)$MATLIN
159
160    if R has EuclideanDomain then
161
162        rowEchelon  x == rowEchelon(x)$MATLIN
163
164    if R has IntegralDomain then
165
166        rank        x == rank(x)$MATLIN
167        nullity     x == nullity(x)$MATLIN
168        nullSpace   x == nullSpace(x)$MATLIN
169
170    if R has Field then
171
172        inverse     x == inverse(x)$MATLIN
173
174    if R has IntegralDomain then
175
176        invertIfCan(x) == invertIfCan(x)$MATLIN
177
178--     matrix(v: Vector Vector R) ==
179--       (rows := # v) = 0 => new(0, 0, 0)
180--       -- error check: this is a top level function
181--       cols := # v.mini(v)
182--       for k in (mini(v) + 1)..maxi(v) repeat
183--         cols ~= # v.k => error "matrix: rows of different lengths"
184--       ans := new(rows, cols, 0)
185--       for i in minr(ans)..maxr(ans) for k in mini(v)..maxi(v) repeat
186--         vv := v.k
187--         for j in minc(ans)..maxc(ans) for l in mini(vv)..maxi(vv) repeat
188--           ans(i, j) := vv.l
189--       ans
190
191    diagonalMatrix(v : Vector R) ==
192      n := #v; ans := zero(n, n)
193      for i in minr(ans)..maxr(ans) for j in minc(ans)..maxc(ans) _
194          for k in mini(v)..maxi(v) repeat qsetelt!(ans, i, j, qelt(v, k))
195      ans
196
197--     diagonalMatrix(vec: Vector %) ==
198--       rows : NonNegativeInteger := 0
199--       cols : NonNegativeInteger := 0
200--       for r in mini(vec)..maxi(vec) repeat
201--         mat := vec.r
202--         rows := rows + nrows mat; cols := cols + ncols mat
203--       ans := zero(rows, cols)
204--       loR := minr ans; loC := minc ans
205--       for r in mini(vec)..maxi(vec) repeat
206--         mat := vec.r
207--         hiR := loR + nrows(mat) - 1; hiC := loC + nrows(mat) - 1
208--         for i in loR..hiR for k in minr(mat)..maxr(mat) repeat
209--           for j in loC..hiC for l in minc(mat)..maxc(mat) repeat
210--             ans(i, j) := mat(k, l)
211--         loR := hiR + 1; loC := hiC + 1
212--       ans
213
214--     vectorOfVectors x ==
215--       vv : Vector Vector R := new(nrows x, 0)
216--       cols := ncols x
217--       for k in mini(vv)..maxi(vv) repeat
218--         vv.k := new(cols, 0)
219--       for i in minr(x)..maxr(x) for k in mini(vv)..maxi(vv) repeat
220--         v := vv.k
221--         for j in minc(x)..maxc(x) for l in mini(v)..maxi(v) repeat
222--           v.l := x(i, j)
223--       vv
224
225    if R has ConvertibleTo InputForm then
226      convert(x : %) : InputForm ==
227         convert [convert('matrix)@InputForm,
228                  convert listOfLists x]$List(InputForm)
229
230)abbrev domain RMATRIX RectangularMatrix
231++ Author: Grabmeier, Gschnitzer, Williamson
232++ Date Created: 1987
233++ Basic Operations:
234++ Related Domains: IndexedMatrix, Matrix, SquareMatrix
235++ Also See:
236++ AMS Classifications:
237++ Keywords: matrix, linear algebra
238++ Examples:
239++ References:
240++ Description:
241++   \spadtype{RectangularMatrix} is a matrix domain where the number of rows
242++   and the number of columns are parameters of the domain.
243RectangularMatrix(m, n, R) : Exports == Implementation where
244  m, n : NonNegativeInteger
245  R   : Join(SemiRng, AbelianMonoid)
246  Row ==> DirectProduct(n, R)
247  Col ==> DirectProduct(m, R)
248  Exports ==> Join(RectangularMatrixCategory(m, n, R, Row, Col), _
249                   CoercibleTo Matrix R) with
250
251    if R has ConvertibleTo InputForm then ConvertibleTo InputForm
252
253    rectangularMatrix : Matrix R -> %
254      ++ \spad{rectangularMatrix(m)} converts a matrix of type \spadtype{Matrix}
255      ++ to a matrix of type \spad{RectangularMatrix}.
256    coerce : % -> Matrix R
257      ++ \spad{coerce(m)} converts a matrix of type \spadtype{RectangularMatrix}
258      ++ to a matrix of type \spad{Matrix}.
259
260  Implementation ==> Matrix R add
261    minr ==> minRowIndex
262    maxr ==> maxRowIndex
263    minc ==> minColIndex
264    maxc ==> maxColIndex
265    mini ==> minIndex
266    maxi ==> maxIndex
267
268    ZERO := new(m, n, 0)$Matrix(R) pretend %
269    0    == ZERO
270
271    coerce(x : %) : OutputForm == coerce(x pretend Matrix R)$Matrix(R)
272
273    matrix(l : List List R) ==
274      -- error check: this is a top level function
275      #l ~= m => error "matrix: wrong number of rows"
276      for ll in l repeat
277        #ll ~= n => error "matrix: wrong number of columns"
278      ans : Matrix R := new(m, n, 0)
279      for i in minr(ans)..maxr(ans) for ll in l repeat
280        for j in minc(ans)..maxc(ans) for r in ll repeat
281          qsetelt!(ans, i, j, r)
282      ans pretend %
283
284    row(x, i)    == directProduct row(x pretend Matrix(R), i)
285    column(x, j) == directProduct column(x pretend Matrix(R), j)
286
287    coerce(x : %) : Matrix(R) == copy(x pretend Matrix(R))
288
289    rectangularMatrix x ==
290      (nrows(x) ~= m) or (ncols(x) ~= n) =>
291        error "rectangularMatrix: matrix of bad dimensions"
292      copy(x) pretend %
293
294    if R has EuclideanDomain then
295
296      rowEchelon x == rowEchelon(x pretend Matrix(R)) pretend %
297
298      columnSpace x ==
299        [directProduct c for c in columnSpace(x pretend Matrix R)]
300
301    if R has IntegralDomain then
302
303      rank x    == rank(x pretend Matrix(R))
304      nullity x == nullity(x pretend Matrix(R))
305      nullSpace x ==
306        [directProduct c for c in nullSpace(x pretend Matrix(R))]
307
308    if R has ConvertibleTo InputForm then
309      convert(x : %) : InputForm ==
310         convert [convert('rectangularMatrix)@InputForm,
311                  convert(x::Matrix(R))]$List(InputForm)
312
313)abbrev domain SQMATRIX SquareMatrix
314++ Author: Grabmeier, Gschnitzer, Williamson
315++ Date Created: 1987
316++ Basic Operations:
317++ Related Domains: IndexedMatrix, Matrix, RectangularMatrix
318++ Also See:
319++ AMS Classifications:
320++ Keywords: matrix, linear algebra
321++ Examples:
322++ References:
323++ Description:
324++   \spadtype{SquareMatrix} is a matrix domain of square matrices, where the
325++   number of rows (= number of columns) is a parameter of the type.
326SquareMatrix(ndim, R) : Exports == Implementation where
327  ndim : NonNegativeInteger
328  R    : Join(SemiRng, AbelianMonoid)
329  Row ==> DirectProduct(ndim, R)
330  Col ==> DirectProduct(ndim, R)
331  MATLIN ==> MatrixLinearAlgebraFunctions(R, Row, Col, %)
332
333  Exports ==> Join(SquareMatrixCategory(ndim, R, Row, Col), _
334                   CoercibleTo Matrix R) with
335
336    transpose : % -> %
337      ++ \spad{transpose(m)} returns the transpose of the matrix m.
338    squareMatrix : Matrix R -> %
339      ++ \spad{squareMatrix(m)} converts a matrix of type \spadtype{Matrix}
340      ++ to a matrix of type \spadtype{SquareMatrix}.
341    coerce : % -> Matrix R
342      ++ \spad{coerce(m)} converts a matrix of type \spadtype{SquareMatrix}
343      ++ to a matrix of type \spadtype{Matrix}.
344--  symdecomp : % -> Record(sym: %, antisym: %)
345--    ++ \spad{symdecomp(m)} decomposes the matrix m as a sum of a symmetric
346--    ++ matrix \spad{m1} and an antisymmetric matrix \spad{m2}. The object
347--    ++ returned is the Record \spad{[m1, m2]}
348--  if R has CommutativeStar then
349--    minorsVect: -> Vector(Union(R,"uncomputed")) --range: 1..2^n-1
350--      ++ \spad{minorsVect(m)} returns a vector of the minors of the matrix m
351    if R has CommutativeStar and R has unitsKnown then unitsKnown
352      ++ the invertible matrices are simply the matrices whose determinants
353      ++ are units in the Ring R.
354    if R has ConvertibleTo InputForm then ConvertibleTo InputForm
355
356  Implementation ==> Matrix R add
357    minr ==> minRowIndex
358    maxr ==> maxRowIndex
359    minc ==> minColIndex
360    maxc ==> maxColIndex
361    mini ==> minIndex
362    maxi ==> maxIndex
363
364    ZERO := scalarMatrix 0
365    0    == ZERO
366
367    if R has Monoid then
368
369        ONE  := scalarMatrix 1
370        1    == ONE
371
372    if R has Ring then
373
374        characteristic() == characteristic()$R
375
376    matrix(l : List List R) ==
377      -- error check: this is a top level function
378      #l ~= ndim => error "matrix: wrong number of rows"
379      for ll in l repeat
380        #ll ~= ndim => error "matrix: wrong number of columns"
381      ans : Matrix R := new(ndim, ndim, 0)
382      for i in minr(ans)..maxr(ans) for ll in l repeat
383        for j in minc(ans)..maxc(ans) for r in ll repeat
384          qsetelt!(ans, i, j, r)
385      ans pretend %
386
387    row(x, i)    == directProduct row(x pretend Matrix(R), i)
388    column(x, j) == directProduct column(x pretend Matrix(R), j)
389    coerce(x : %) : OutputForm == coerce(x pretend Matrix R)$Matrix(R)
390
391    scalarMatrix r == scalarMatrix(ndim, r)$Matrix(R) pretend %
392
393    diagonalMatrix l ==
394      #l ~= ndim =>
395        error "diagonalMatrix: wrong number of entries in list"
396      diagonalMatrix(l)$Matrix(R) pretend %
397
398    coerce(x : %) : Matrix(R) == copy(x pretend Matrix(R))
399
400    squareMatrix x ==
401      (nrows(x) ~= ndim) or (ncols(x) ~= ndim) =>
402        error "squareMatrix: matrix of bad dimensions"
403      copy(x) pretend %
404
405    x : % * v : Col ==
406      directProduct((x pretend Matrix(R)) * (v :: Vector(R)))
407
408    v : Row * x : % ==
409      directProduct((v :: Vector(R)) * (x pretend Matrix(R)))
410
411    if R has CommutativeRing then
412
413      determinant x == determinant(x pretend Matrix(R))
414      minordet x    == minordet(x pretend Matrix(R))
415      Pfaffian x    == Pfaffian(x pretend Matrix(R))
416
417    if R has EuclideanDomain then
418
419      rowEchelon x == rowEchelon(x pretend Matrix(R)) pretend %
420
421      columnSpace x ==
422          [directProduct c for c in columnSpace(x pretend Matrix R)]
423
424    if R has IntegralDomain then
425
426      rank x    == rank(x pretend Matrix(R))
427      nullity x == nullity(x pretend Matrix(R))
428      nullSpace x ==
429        [directProduct c for c in nullSpace(x pretend Matrix(R))]
430      recip(x) ==
431          (u := invertIfCan(x pretend Matrix(R))) case "failed" => "failed"
432          (u :: Matrix(R)) pretend %
433
434    if R has Field then
435
436      -- dimension() == (m * n) :: CardinalNumber
437
438      inverse x ==
439        (u := inverse(x pretend Matrix(R))) case "failed" => "failed"
440        (u :: Matrix(R)) pretend %
441
442      x : % ^ n : Integer ==
443        ((x pretend Matrix(R)) ^ n) pretend %
444
445      recip x == inverse x
446
447    if R has ConvertibleTo InputForm then
448      convert(x : %) : InputForm ==
449         convert [convert('squareMatrix)@InputForm,
450                  convert(x::Matrix(R))]$List(InputForm)
451
452
453--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
454--All rights reserved.
455--
456--Redistribution and use in source and binary forms, with or without
457--modification, are permitted provided that the following conditions are
458--met:
459--
460--    - Redistributions of source code must retain the above copyright
461--      notice, this list of conditions and the following disclaimer.
462--
463--    - Redistributions in binary form must reproduce the above copyright
464--      notice, this list of conditions and the following disclaimer in
465--      the documentation and/or other materials provided with the
466--      distribution.
467--
468--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
469--      names of its contributors may be used to endorse or promote products
470--      derived from this software without specific prior written permission.
471--
472--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
473--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
474--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
475--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
476--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
477--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
478--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
479--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
480--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
481--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
482--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
483