====================================================================
Matrix Category
====================================================================

**Predicates**

``square?(m)`` returns true if m is a square matrix
(if m has the same number of rows as columns) and false otherwise. ::

  square matrix [[j**i for i in 0..4] for j in 1..5]

``diagonal?(m)`` returns true if the matrix m is square and
diagonal (i.e. all entries of m not on the diagonal are zero) and
false otherwise. ::

  diagonal? matrix [[j**i for i in 0..4] for j in 1..5]

``symmetric?(m)`` returns true if the matrix m is square and
symmetric (i.e. m[i,j] = m[j,i] for all i and j) and false
otherwise. ::

  symmetric? matrix [[j**i for i in 0..4] for j in 1..5]

``antisymmetric?(m)`` returns true if the matrix m is square and
antisymmetric (i.e. m[i,j] = -m[j,i] for all i and j) 
and false otherwise. ::

  antisymmetric? matrix [[j**i for i in 0..4] for j in 1..5]


**Creation**

``zero(m,n)`` returns an m-by-n zero matrix. ::

  z:Matrix(INT):=zero(3,3)

``matrix(l)`` converts the list of lists l to a matrix, where the
list of lists is viewed as a list of the rows of the matrix. ::

  matrix [[1,2,3],[4,5,6],[7,8,9],[1,1,1]]

``scalarMatrix(n,r)`` returns an n-by-n matrix with r's on the
diagonal and zeroes elsewhere. ::

  z:Matrix(INT):=scalarMatrix(3,5)

``diagonalMatrix(l)`` returns a diagonal matrix with the elements
of l on the diagonal. ::

  diagonalMatrix [1,2,3]

``diagonalMatrix([m1,...,mk])`` creates a block diagonal matrix
M with block matrices m1,...,mk down the diagonal,
with 0 block matrices elsewhere.

More precisly: ::
	
 if ri := nrows mi, ci := ncols mi,
 then m is an (r1+..+rk) by (c1+..+ck) - matrix  with entries
   m.i.j = ml.(i-r1-..-r(l-1)).(j-n1-..-n(l-1)), 
 if (r1+..+r(l-1)) < i <= r1+..+rl and
   (c1+..+c(l-1)) < i <= c1+..+cl,
    m.i.j = 0  otherwise.

  diagonalMatrix [matrix [[1,2],[3,4]], matrix [[4,5],[6,7]]]

``coerce(col)`` converts the column col to a column matrix. ::

  coerce([1,2,3])@Matrix(INT)

``transpose(r)`` converts the row r to a row matrix::

  transpose([1,2,3])@Matrix(INT)


**Creation of new matrices from old**

``transpose(m)`` returns the transpose of the matrix m. ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  transpose m

``squareTop(m)`` returns an n-by-n matrix consisting of the first
n rows of the m-by-n matrix m. Error: if m < n. ::

  m:=matrix [[j**i for i in 0..2] for j in 1..5]
  squareTop m

``horizConcat(x,y)`` horizontally concatenates two matrices with
an equal number of rows. The entries of y appear to the right
of the entries of x.  Error: if the matrices
do not have the same number of rows. ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  horizConcat(m,m)

``vertConcat(x,y)`` vertically concatenates two matrices with an
equal number of columns. The entries of y appear below
of the entries of x.  Error: if the matrices
do not have the same number of columns. ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  vertConcat(m,m)

**Part extractions/assignments**

``listOfLists(m)`` returns the rows of the matrix m as a list of lists ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  listOfLists m

``elt(x,rowList,colList)`` returns an m-by-n matrix consisting
of elements of x, where m = # rowList and n = # colList
If rowList = [i<1>,i<2>,...,i<m>] and 
colList = [j<1>,j<2>,...,j<n>], 
then the (k,l)-th entry of elt(x,rowList,colList) is x(i<k>,j<l>). ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  elt(m,3,3)

``setelt(x,rowList,colList,y)`` destructively alters the matrix x.
If y is m-by-n, 
rowList = [i<1>,i<2>,...,i<m>] and
colList = [j<1>,j<2>,...,j<n>], 
then x(i<k>,j<l>)
is set to y(k,l) for k = 1,...,m and l = 1,...,n ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  setelt(m,3,3,10)

``swapRows!(m,i,j)`` interchanges the i-th and j-th
rows of m. This destructively alters the matrix. ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  swapRows!(m,2,4)

``swapColumns!(m,i,j)`` interchanges the i-th and j-th
columns of m. This destructively alters the matrix. ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  swapColumns!(m,2,4)

``subMatrix(x,i1,i2,j1,j2)`` extracts the submatrix [x(i,j)] 
where the index i ranges from i1 to i2
and   the index j ranges from j1 to j2. ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  subMatrix(m,1,3,2,4)

``setsubMatrix(x,i1,j1,y)`` destructively alters the matrix x. 
Here x(i,j) is set to y(i-i1+1,j-j1+1) for
i = i1,...,i1-1+nrows y and j = j1,...,j1-1+ncols y. ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  setsubMatrix!(m,2,2,matrix [[3,3],[3,3]])

----------
Arithmetic
----------

``x + y`` is the sum of the matrices x and y.
It is an error if the dimensions are incompatible. ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  m+m

``x - y`` is the difference of the matrices x and y.
It is an error if the dimensions are incompatible. ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  m-m

``-x`` returns the negative of the matrix x. ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  -m

``x * y`` is the product of the matrices x and y.
It is an error if the dimensions are incompatible. ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  m*m

``r*x`` is the left scalar multiple of the scalar r and the matrix x. ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  1/3*m

``x * r`` is the right scalar multiple of the scalar r and the matrix x. ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  m*1/3

``n * x`` is an integer multiple. ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  3*m

``x * c`` is the product of the matrix x and the column vector c.
It is an **error** if the dimensions are incompatible. ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  c:=coerce([1,2,3,4,5])@Matrix(INT)
  m *c

``r * x`` is the product of the row vector r and the matrix x.
It is an error if the dimensions are incompatible. ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  r:=transpose([1,2,3,4,5])@Matrix(INT)
  r*m

``x ** n`` computes a non-negative integral power of the matrix x.
It is an error if the matrix is not square. ::

  m:=matrix [[j**i for i in 0..4] for j in 1..5]
  m**3

``exquo(m,r)`` computes the exact quotient of the elements
of m by r, returning "failed" if this is not possible. ::

  m:=matrix [[2**i for i in 2..4] for j in 1..5]
  exquo(m,2)

``m/r`` divides the elements of m by r, r must be non-zero. ::

  m:=matrix [[2**i for i in 2..4] for j in 1..5]
  m/4

--------------
Linear algebra
--------------

``rowEchelon(m)`` returns the row echelon form of the matrix m. ::

  rowEchelon matrix [[j**i for i in 0..4] for j in 1..5]

``columnSpace(m)`` returns a sublist of columns of the matrix m ::

  columnSpace matrix [[1,2,3],[4,5,6],[7,8,9],[1,1,1]]

``rank(m)`` returns the rank of the matrix m. ::

  rank matrix [[1,2,3],[4,5,6],[7,8,9]]

``nullity(m)`` returns the nullity of the matrix m. This is
the dimension of the null space of the matrix m. ::

  nullity matrix [[1,2,3],[4,5,6],[7,8,9]]

``nullSpace(m)`` returns a basis for the null space of the matrix m. ::

  nullSpace matrix [[1,2,3],[4,5,6],[7,8,9]]

``determinant(m)`` returns the determinant of the matrix m.
It is an error if the matrix is not square. ::

  determinant matrix [[j**i for i in 0..4] for j in 1..5]

``minordet(m)`` computes the determinant of the matrix m using minors. 
It is an error if the matrix is not square. ::

  minordet matrix [[j**i for i in 0..4] for j in 1..5]

``pfaffian(m)`` returns the Pfaffian of the matrix m.
It is an error if the matrix is not antisymmetric ::

  pfaffian [[0,1,0,0],[-1,0,0,0],[0,0,0,1],[0,0,-1,0]]

``inverse(m)`` returns the inverse of the matrix m.
If the matrix is not invertible, "failed" is returned.
It is an error if the matrix is not square. ::

  inverse matrix [[j**i for i in 0..4] for j in 1..5]

``m**n`` computes an integral power of the matrix m.
It is an error if matrix is not square or 
if the matrix is square but not invertible. ::

  (matrix [[j**i for i in 0..4] for j in 1..5]) ** 2