%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%A  matrix.msk                  GAP documentation            Martin Schoenert
%A                                                           Alexander Hulpke
%%
%A  @(#)$Id: matrix.msk,v 1.40.2.5 2006/09/16 19:02:49 jjm Exp $
%%
%Y  (C) 1998 School Math and Comp. Sci., University of St.  Andrews, Scotland
%Y  Copyright (C) 2002 The GAP Group
%%
\Chapter{Matrices}

Matrices are represented in {\GAP} by lists of row vectors (see "Row Vectors").
The vectors must all have the same length, and their elements must lie in
a common ring. However, since checking rectangularness can be expensive
functions and methods of operations for matrices often will not give an error
message for non-rectangular lists of lists -- in such cases the result is
undefined.

Because matrices are just a special case of lists, all operations and
functions for lists are applicable to matrices also (see chapter
"Lists"). This especially includes accessing elements of a matrix (see
"List Elements"), changing elements of a matrix (see "List Assignment"),
and comparing matrices (see "Comparisons of Lists").

Note that, since a matrix is a list of lists, the behaviour of `ShallowCopy'
for matrices is just a special case of `ShallowCopy' for lists
(see~"Duplication of Lists");
called with an immutable matrix <mat>, `ShallowCopy' returns a mutable matrix
whose rows are identical to the rows of <mat>.
In particular the rows are still immutable.
To get a matrix whose rows are mutable,
one can use `List( <mat>, ShallowCopy )'.

\Declaration{InfoMatrix}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{Categories of Matrices}

\Declaration{IsMatrix}
\beginexample
gap> mat:=[[1,2,3],[4,5,6],[7,8,9]];
[ [ 1, 2, 3 ], [ 4, 5, 6 ], [ 7, 8, 9 ] ]
gap> IsMatrix(mat);
true
\endexample

Note also the filter `IsTable' (see section "IsTable")
which  may be more appropriate than `IsMatrix' for some purposes.

Note that the empty list '[ ]' and more complex ``empty'' structures
such as `[[ ]]' are *not* matrices, although special
methods allow them be used in place of matrices in some
situations. See "EmptyMatrix" below.

\beginexample
gap> [[0]]*[[]];
[ [  ] ]
gap> IsMatrix([[]]);
false
\endexample


\Declaration{IsOrdinaryMatrix}
\Declaration{IsLieMatrix}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{Operators for Matrices}

The rules for arithmetic operations involving matrices are in fact
special cases of those for the arithmetic of lists,
given in Section~"Arithmetic for Lists" and the following sections,
here we reiterate that definition, in the language of vectors and matrices.

Note that the additive behaviour sketched below is defined only for lists in
the category `IsGeneralizedRowVector',
and the multiplicative behaviour is defined only for lists in the category
`IsMultiplicativeGeneralizedRowVector'
(see~"Filters Controlling the Arithmetic Behaviour of Lists").


\>`<mat1> + <mat2>'{addition!matrices} O

returns the sum of the two matrices <mat1> and <mat2>,
Probably the most usual situation is that <mat1> and <mat2> have the same
dimensions and are defined over a common field;
in this case the sum is a new matrix over the same field where each entry
is the sum of the corresponding entries of the matrices.

In more general situations, the sum of two matrices need not be a matrix,
for example adding an integer matrix <mat1> and a matrix <mat2> over
a finite field yields the table of pointwise sums,
which will be a mixture of finite field elements and integers if <mat1> has
bigger dimensions than <mat2>.


\>`<scalar> + <mat>'{addition!scalar and matrix} O
\>`<mat> + <scalar>'{addition!matrix and scalar} O

returns the sum of the scalar <scalar> and the matrix <mat>.
Probably the most usual situation is that the entries of <mat> lie in a
common field with <scalar>;
in this case the sum is a new matrix over the same field where each entry
is the sum of the scalar and the corresponding entry of the matrix.

More general situations are for example the sum of an integer scalar and a
matrix over a finite field, or the sum of a finite field element and an
integer matrix.


\>`<mat1> - <mat2>'{subtraction!matrices}
\>`<scalar> - <mat>'{subtraction!scalar and matrix} O
\>`<mat> - <scalar>'{subtraction!matrix and scalar} O

Subtracting a matrix or scalar is defined as adding its additive inverse,
so the statements for the addition hold likewise.


\>`<scalar> * <mat>'{multiplication!scalar and matrix} O
\>`<mat> * <scalar>'{multiplication!matrix and scalar} O

returns the product of the scalar <scalar> and the matrix <mat>.
Probably the most usual situation is that the elements of <mat> lie in a
common field with <scalar>;
in this case the product is a new matrix over the same field where each
entry is the product of the scalar and the corresponding entry of the matrix.

More general situations are for example the product of an integer scalar and
a matrix over a finite field,
or the product of a finite field element and an integer matrix.


\>`<vec> * <mat>'{multiplication!vector and matrix} O

returns the product of the row vector <vec> and the matrix <mat>.
Probably the most usual situation is that <vec> and <mat> have the same
lengths and are defined over a common field,
and that all rows of <mat> have the same length $m$, say;
in this case the product is a new row vector of length $m$ over the same
field which is the sum of the scalar multiples of the rows of <mat> with the
corresponding entries of <vec>.

More general situations are for example the product of an integer vector and
a matrix over a finite field,
or the product of a vector over a finite field and an integer matrix.


\>`<mat> * <vec>'{multiplication!matrix and vector} O

returns the product of the matrix <mat> and the row vector <vec>.
(This is the standard product of a matrix with a *column* vector.)
Probably the most usual situation is that the length of <vec> and of all rows
of <mat> are equal, and that the elements of <mat> and <vec> lie in a common
field;
in this case the product is a new row vector of the same length as <mat> and
over the same field which is the sum of the scalar multiples of the columns
of <mat> with the corresponding entries of <vec>.

More general situations are for example the product of an integer matrix and
a vector over a finite field,
or the product of a matrix over a finite field and an integer vector.


\>`<mat1> * <mat2>'{multiplication!matrices} O

This form evaluates to the (Cauchy) product of the two matrices <mat1> and
<mat2>.
Probably the most usual situation is that the number of columns of <mat1>
equals the number of rows of <mat2>,
and that the elements of <mat> and <vec> lie in a common field;
if <mat1> is a matrix with $m$ rows and $n$ columns, say, and <mat2> is a
matrix with $n$ rows and $o$ columns, the result is a new matrix with $m$
rows and $o$ columns.
The element in row $i$ at position $j$ of the product is the sum of
$<mat1>[i][l] \* <mat2>[l][j]$, with $l$ running from $1$ to $n$.


\>`Inverse( <mat> )'{inverse!matrix} O

returns the inverse of the matrix <mat>,
which must be an invertible square matrix.
If <mat> is not invertible then `fail' is returned.


\>`<mat1> / <mat2>'{quotient!matrices} O
\>`<scalar> / <mat>'{quotient!scalar and matrix} O
\>`<mat> / <scalar>'{quotient!matrix and scalar} O
\>`<vec> / <mat>'{quotient!vector and matrix} O

In general, `<left> / <right>' is defined as `<left> * <right>^-1'.
Thus in the above forms the right operand must always be invertible.


\>`<mat> ^ <int>'{power!matrix} O
\>`<mat1> ^ <mat2>'{conjugate!matrix} O
\>`<vec> ^ <mat>'{image!vector under matrix} O

Powering a square matrix <mat> by an integer <int> yields the <int>-th power
of <mat>; if <int> is negative then <mat> must be invertible,
if <int> is `0' then the result is the identity matrix `One( <mat> )',
even if <mat> is not invertible.

Powering a square matrix <mat1> by an invertible square matrix <mat2> of the
same dimensions yields the conjugate of <mat1> by <mat2>, i.e.,
the matrix `<mat2>^-1 * <mat1> * <mat2>'.

Powering a row vector <vec> by a matrix <mat> is in every respect equivalent
to `<vec> * <mat>'.
This operations reflects the fact that matrices act naturally on row vectors
by multiplication from the right, and that the powering operator is {\GAP}'s
standard for group actions.


\>`Comm( <mat1>, <mat2> )'{matrices!commutator} O

returns the commutator of the square invertible matrices <mat1> and <mat2>
of the same dimensions and over a common field,
which is the matrix `<mat1>^-1 * <mat2>^-1 * <mat1> * <mat2>'.


The following cases are still special cases of the general list arithmetic
defined in~"Arithmetic for Lists".


\>`<scalar> + <matlist>'{addition!scalar and matrix list} O
\>`<matlist> + <scalar>'{addition!scalar and matrix list} O
\>`<scalar> - <matlist>'{subtraction!scalar and matrix list} O
\>`<matlist> - <scalar>'{subtraction!scalar and matrix list} O
\>`<scalar> * <matlist>'{multiplication!scalar and matrix list} O
\>`<matlist> * <scalar>'{multiplication!scalar and matrix list} O
\>`<matlist> / <scalar>'{quotient!scalar and matrix list} O

A scalar <scalar> may also be added, subtracted, multiplied with, or
divided into a list <matlist> of matrices. The result is a new list
of matrices where each matrix is the result of performing the operation
with the corresponding matrix in <matlist>.


\>`<mat> * <matlist>'{multiplication!matrix and matrix list} O
\>`<matlist> * <mat>'{multiplication!matrix and matrix list} O

A matrix <mat> may also be multiplied with a list <matlist> of matrices.
The result is a new list of matrices, where each entry is the product of
<mat> and the corresponding entry in <matlist>.


\>`<matlist> / <mat>'{quotient!matrix and matrix list} O

Dividing a list <matlist> of matrices by an invertible matrix <mat>
evaluates to `<matlist> * <mat>^-1'.


\>`<vec> * <matlist>'{multiplication!vector and matrix list} O

returns the product of the vector <vec> and the list of matrices <mat>.
The lengths <l> of <vec> and <matlist> must be equal.
All matrices in <matlist> must have the same dimensions. The elements of
<vec> and the elements of the matrices in <matlist> must lie in a common
ring. The product is the sum over `<vec>[<i>] * <matlist>[<i>]' with
<i> running from 1 to <l>.


For the mutability of results of arithmetic operations,
see~"Mutability and Copyability".


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{Properties and Attributes of Matrices}

\Declaration{DimensionsMat}
\beginexample
gap> DimensionsMat([[1,2,3],[4,5,6]]);
[ 2, 3 ]
\endexample
\Declaration{DefaultFieldOfMatrix}
\beginexample
gap> DefaultFieldOfMatrix([[Z(4),Z(8)]]);
GF(2^6)
\endexample

\indextt{Trace!of a matrix}
\Declaration{TraceMat}
\beginexample
gap> TraceMat([[1,2,3],[4,5,6],[7,8,9]]);
15
\endexample
\Declaration{DeterminantMat}
\Declaration{DeterminantMatDestructive}
\beginexample
gap> DeterminantMat([[1,2],[2,1]]);
-3
gap> mm:= [[1,2],[2,1]];;
gap> DeterminantMatDestructive( mm );
-3
gap> mm;
[ [ 1, 2 ], [ 0, -3 ] ]
\endexample

\Declaration{DeterminantMatDivFree}

\Declaration{IsMonomialMatrix}
\beginexample
gap> IsMonomialMatrix([[0,1],[1,0]]);
true
\endexample

\Declaration{IsDiagonalMat}
\Declaration{IsUpperTriangularMat}
\Declaration{IsLowerTriangularMat}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{Matrix Constructions}

\Declaration{IdentityMat}
\Declaration{NullMat}
\beginexample
gap> IdentityMat(3,1);
[ [ 1, 0, 0 ], [ 0, 1, 0 ], [ 0, 0, 1 ] ]
gap> NullMat(3,2,Z(3));
[ [ 0*Z(3), 0*Z(3) ], [ 0*Z(3), 0*Z(3) ], [ 0*Z(3), 0*Z(3) ] ]
\endexample
\Declaration{EmptyMatrix}
\beginexample
gap> EmptyMatrix(5);
EmptyMatrix( 5 )
gap> AsList(last);
[  ]
\endexample

\Declaration{DiagonalMat}
\beginexample
gap> DiagonalMat([1,2,3]);
[ [ 1, 0, 0 ], [ 0, 2, 0 ], [ 0, 0, 3 ] ]
\endexample
\Declaration{PermutationMat}
\beginexample
gap> PermutationMat((1,2,3),4,1);
[ [ 0, 1, 0, 0 ], [ 0, 0, 1, 0 ], [ 1, 0, 0, 0 ], [ 0, 0, 0, 1 ] ]
\endexample

\Declaration{TransposedMatImmutable}
\Declaration{TransposedMatDestructive}
\beginexample
gap> TransposedMat([[1,2,3],[4,5,6],[7,8,9]]);
[ [ 1, 4, 7 ], [ 2, 5, 8 ], [ 3, 6, 9 ] ]
gap> mm:= [[1,2,3],[4,5,6],[7,8,9]];;
gap> TransposedMatDestructive( mm );
[ [ 1, 4, 7 ], [ 2, 5, 8 ], [ 3, 6, 9 ] ]
gap> mm;
[ [ 1, 4, 7 ], [ 2, 5, 8 ], [ 3, 6, 9 ] ]
\endexample

\Declaration{KroneckerProduct}
\beginexample
gap> KroneckerProduct([[1,2]],[[5,7],[9,2]]);
[ [ 5, 7, 10, 14 ], [ 9, 2, 18, 4 ] ]
\endexample


\Declaration{ReflectionMat}

\Declaration{PrintArray}

\Declaration{MutableIdentityMat}

\Declaration{MutableNullMat}

\Declaration{MutableCopyMat}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{Random Matrices}

\Declaration{RandomMat}
\Declaration{RandomInvertibleMat}
\Declaration{RandomUnimodularMat}
\beginexample
gap> RandomMat(2,3,GF(3));
[ [ Z(3), Z(3), 0*Z(3) ], [ Z(3), Z(3)^0, Z(3) ] ]
gap> RandomInvertibleMat(4);
[ [ 1, -2, -1, 0 ], [ 1, 0, 1, -1 ], [ 0, 2, 0, 4 ], [ -1, -3, 1, -4 ] ]
\endexample

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{Matrices Representing Linear Equations and the Gaussian Algorithm}

\atindex{Gaussian algorithm}{@Gaussian algorithm}
\Declaration{RankMat}
\beginexample
gap> mat:=[[1,2,3],[4,5,6],[7,8,9]];;
gap> RankMat(mat);
2
\endexample
\Declaration{TriangulizeMat}
\beginexample
gap> m:=TransposedMatMutable(mat);
[ [ 1, 4, 7 ], [ 2, 5, 8 ], [ 3, 6, 9 ] ]
gap> TriangulizeMat(m);m;
[ [ 1, 0, -1 ], [ 0, 1, 2 ], [ 0, 0, 0 ] ]
\endexample
\Declaration{NullspaceMat}
\Declaration{NullspaceMatDestructive}
\beginexample
gap> mat:=[[1,2,3],[4,5,6],[7,8,9]];;
gap> NullspaceMat(mat);
[ [ 1, -2, 1 ] ]
gap> mm:=[[1,2,3],[4,5,6],[7,8,9]];;
gap> NullspaceMatDestructive( mm );
[ [ 1, -2, 1 ] ]
gap> mm;
[ [ 1, 2, 3 ], [ 0, -3, -6 ], [ 0, 0, 0 ] ]
\endexample
\Declaration{SolutionMat}
\Declaration{SolutionMatDestructive}
\beginexample
gap> mat:=[[1,2,3],[4,5,6],[7,8,9]];;
gap> SolutionMat(mat,[3,5,7]);
[ 5/3, 1/3, 0 ]
gap> mm:=[[1,2,3],[4,5,6],[7,8,9]];;
gap> SolutionMatDestructive( mm, [3,5,7] );
[ 5/3, 1/3, 0 ]
gap> mm;
[ [ 1, 2, 3 ], [ 0, -3, -6 ], [ 0, 0, 0 ] ]
\endexample
\Declaration{BaseFixedSpace}
\beginexample
gap> BaseFixedSpace([[[1,2],[0,1]]]);
[ [ 0, 1 ] ]
\endexample

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{Eigenvectors and eigenvalues}

\Declaration{GeneralisedEigenvalues}
\Declaration{GeneralisedEigenspaces}
\Declaration{Eigenvalues}
\Declaration{Eigenspaces}
\Declaration{Eigenvectors}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{Elementary Divisors}

See also chapter "Integral Matrices and Lattices".

\Declaration{ElementaryDivisorsMat}
\beginexample
gap> mat:=[[1,2,3],[4,5,6],[7,8,9]];;
gap> ElementaryDivisorsMat(mat);
[ 1, 3, 0 ]
gap> x:=X(Rationals,"x");;
gap> mat:=mat*One(x)-x*mat^0;       
[ [ -x+1, 2, 3 ], [ 4, -x+5, 6 ], [ 7, 8, -x+9 ] ]
gap> ElementaryDivisorsMat(PolynomialRing(Rationals,1),mat);
[ 1, 1, x^3-15*x^2-18*x ]
gap> mat:=KroneckerProduct(CompanionMat((x-1)^2),CompanionMat((x^3-1)*(x-1)));;
gap> mat:=mat*One(x)-x*mat^0;                                                 
[ [ -x, 0, 0, 0, 0, 0, 0, 1 ], [ 0, -x, 0, 0, -1, 0, 0, -1 ], 
  [ 0, 0, -x, 0, 0, -1, 0, 0 ], [ 0, 0, 0, -x, 0, 0, -1, -1 ], 
  [ 0, 0, 0, -1, -x, 0, 0, -2 ], [ 1, 0, 0, 1, 2, -x, 0, 2 ], 
  [ 0, 1, 0, 0, 0, 2, -x, 0 ], [ 0, 0, 1, 1, 0, 0, 2, -x+2 ] ]
gap> ElementaryDivisorsMat(PolynomialRing(Rationals,1),mat);
[ 1, 1, 1, 1, 1, 1, x-1, x^7-x^6-2*x^4+2*x^3+x-1 ]
\endexample
\Declaration{DiagonalizeMat}
\beginexample
gap> m:=[[1,2],[2,1]];;
gap> DiagonalizeMat(Integers,m);m;
[ [ 1, 0 ], [ 0, 3 ] ]
\endexample

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{Echelonized Matrices}

\Declaration{SemiEchelonMat}
\Declaration{SemiEchelonMatDestructive}
\beginexample
gap> mm:=[[1,2,3],[4,5,6],[7,8,9]];;
gap> SemiEchelonMatDestructive( mm );
rec( heads := [ 1, 2, 0 ], vectors := [ [ 1, 2, 3 ], [ 0, 1, 2 ] ] )
gap> mm;
[ [ 1, 2, 3 ], [ 0, 1, 2 ], [ 0, 0, 0 ] ]
\endexample

\Declaration{SemiEchelonMatTransformation}
\beginexample
gap> SemiEchelonMatTransformation([[1,2,3],[0,0,1]]);
rec( heads := [ 1, 0, 2 ], vectors := [ [ 1, 2, 3 ], [ 0, 0, 1 ] ], 
  coeffs := [ [ 1, 0 ], [ 0, 1 ] ], relations := [  ] )
\endexample
\Declaration{SemiEchelonMats}
\Declaration{SemiEchelonMatsDestructive}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{Matrices as Basis of a Row Space}

\Declaration{BaseMat}
\Declaration{BaseMatDestructive}
\beginexample
gap> mat:=[[1,2,3],[4,5,6],[7,8,9]];;
gap> BaseMat(mat);
[ [ 1, 2, 3 ], [ 0, 1, 2 ] ]
gap> mm:= [[1,2,3],[4,5,6],[5,7,9]];;
gap> BaseMatDestructive( mm );
[ [ 1, 2, 3 ], [ 0, 1, 2 ] ]
gap> mm;
[ [ 1, 2, 3 ], [ 0, 1, 2 ], [ 0, 0, 0 ] ]
\endexample
\Declaration{BaseOrthogonalSpaceMat}

\Declaration{SumIntersectionMat}
\beginexample
gap> SumIntersectionMat(mat,[[2,7,6],[5,9,4]]);
[ [ [ 1, 2, 3 ], [ 0, 1, 2 ], [ 0, 0, 1 ] ], [ [ 1, -3/4, -5/2 ] ] ]
\endexample

\Declaration{BaseSteinitzVectors}
\beginexample
gap> BaseSteinitzVectors(IdentityMat(3,1),[[11,13,15]]);
rec( factorspace := [ [ 0, 1, 15/13 ], [ 0, 0, 1 ] ], 
  factorzero := [ 0, 0, 0 ], subspace := [ [ 1, 13/11, 15/11 ] ], 
  heads := [ -1, 1, 2 ] )
\endexample

See also chapter "Integral Matrices and Lattices"


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{Triangular Matrices}

\Declaration{DiagonalOfMat}
\beginexample
gap> DiagonalOfMat([[1,2],[3,4]]);
[ 1, 4 ]
\endexample
\Declaration{UpperSubdiagonal}
\beginexample
gap> UpperSubdiagonal(mat,1);
[ 2, 6 ]
\endexample
\Declaration{DepthOfUpperTriangularMatrix}
\beginexample
gap> DepthOfUpperTriangularMatrix([[0,1,2],[0,0,1],[0,0,0]]);
1
gap> DepthOfUpperTriangularMatrix([[0,0,2],[0,0,0],[0,0,0]]);
2
\endexample

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{Matrices as Linear Mappings}

\Declaration{CharacteristicPolynomial}
\beginexample
gap> CharacteristicPolynomial( [ [ 1, 1 ], [ 0, 1 ] ] );
x^2-2*x+1
gap> mat := [[0,1],[E(4)-1,E(4)]];;
gap> CharacteristicPolynomial( mat );
x^2+(-E(4))*x+(1-E(4))
gap> CharacteristicPolynomial( Rationals, CF(4), mat );
x^4+3*x^2+2*x+2
gap> mat:= [ [ E(4), 1 ], [ 0, -E(4) ] ];;
gap> CharacteristicPolynomial( mat );
x^2+1
gap> CharacteristicPolynomial( Rationals, CF(4), mat );
x^4+2*x^2+1
\endexample

\Declaration{JordanDecomposition}
\beginexample
gap> mat:=[[1,2,3],[4,5,6],[7,8,9]];;
gap> JordanDecomposition(mat);
[ [ [ 1, 2, 3 ], [ 4, 5, 6 ], [ 7, 8, 9 ] ], 
  [ [ 0, 0, 0 ], [ 0, 0, 0 ], [ 0, 0, 0 ] ] ]
\endexample

\Declaration{BlownUpMat}
\Declaration{BlownUpVector}
\beginexample
gap> B:= Basis( CF(4), [ 1, E(4) ] );;
gap> mat:= [ [ 1, E(4) ], [ 0, 1 ] ];;  vec:= [ 1, E(4) ];;
gap> bmat:= BlownUpMat( B, mat );;  bvec:= BlownUpVector( B, vec );;
gap> Display( bmat );  bvec;
[ [   1,   0,   0,   1 ],
  [   0,   1,  -1,   0 ],
  [   0,   0,   1,   0 ],
  [   0,   0,   0,   1 ] ]
[ 1, 0, 0, 1 ]
gap> bvec * bmat = BlownUpVector( B, vec * mat );
true
\endexample

\Declaration{CompanionMat}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{Matrices over Finite Fields}

Just as for row vectors, (see section "Row Vectors over Finite
Fields"), {\GAP} has a special representation for matrices over small
finite fields.

To be eligible to be represented in this way, each row of a matrix
must be able to be represented as a compact row vector of the same
length  over *the same* finite field.

\beginexample
gap> v := Z(2)*[1,0,0,1,1];
[ Z(2)^0, 0*Z(2), 0*Z(2), Z(2)^0, Z(2)^0 ]
gap> ConvertToVectorRep(v,2);
2
gap> v;
<a GF2 vector of length 5>
gap> m := [v];; ConvertToMatrixRep(m,GF(2));; m;
<a 1x5 matrix over GF2>
gap> m := [v,v];; ConvertToMatrixRep(m,GF(2));; m;
<a 2x5 matrix over GF2>
gap> m := [v,v,v];; ConvertToMatrixRep(m,GF(2));; m;
<a 3x5 matrix over GF2>
gap> v := Z(3)*[1..8];
[ Z(3), Z(3)^0, 0*Z(3), Z(3), Z(3)^0, 0*Z(3), Z(3), Z(3)^0 ]
gap> ConvertToVectorRep(v);
3
gap> m := [v];; ConvertToMatrixRep(m,GF(3));; m;
[ [ Z(3), Z(3)^0, 0*Z(3), Z(3), Z(3)^0, 0*Z(3), Z(3), Z(3)^0 ] ]
gap> RepresentationsOfObject(m);
[ "IsPositionalObjectRep", "Is8BitMatrixRep" ]
gap> m := [v,v,v,v];; ConvertToMatrixRep(m,GF(3));; m;
< mutable compressed matrix 4x8 over GF(3) >
\endexample

All compressed matrices over GF(2) are viewed as `\<a <n>x<m> matrix
over GF2>', while over fields GF(q) for q between 3 and 256, matrices
with 25 or more entries are viewed in this way, and smaller ones as
lists of lists.

Matrices can be converted to this special representation via
the following functions.

\Declaration{ImmutableMatrix}

\Declaration{ConvertToMatrixRep}

Note that the main advantage of this special representation of
matrices is in low dimensions, where various overheads can be
reduced. In higher dimensions, a list of compressed vectors will be
almost as fast. Note also that list access and assignment will be
somewhat slower for compressed matrices than for plain lists.

In order to form a row of a compressed matrix a vector must accept
certain restrictions. Specifically, it cannot change its length or
change the field over which it is compressed. The main consequences of
this are: that only elements of the appropriate field can be assigned
to entries of the vector, and only to positions between 1 and the
original length; that the vector cannot be shared between two matrices
compressed over different fields.

This is enforced by the filter `IsLockedRepresentationVector'.  When a
vector becomes part of a compressed matrix, this filter is set for it.
Assignment, `Unbind',  `ConvertToVectorRep' and `ConvertToMatrixRep'
are all prevented from altering a vector with this filter.

\beginexample
gap> v := [Z(2),Z(2)];; ConvertToVectorRep(v,GF(2));; v;
<a GF2 vector of length 2>
gap> m := [v,v]; 
[ <a GF2 vector of length 2>, <a GF2 vector of length 2> ]
gap> ConvertToMatrixRep(m,GF(2)); 
2
gap> m2 := [m[1], [Z(4),Z(4)]]; # now try and mix in some GF(4)
[ <a GF2 vector of length 2>, [ Z(2^2), Z(2^2) ] ]
gap> ConvertToMatrixRep(m2); # but m2[1] is locked
#I  ConvertToVectorRep: locked vector not converted to different field
fail
gap> m2 := [ShallowCopy(m[1]), [Z(4),Z(4)]]; # a fresh copy of row 1
[ <a GF2 vector of length 2>, [ Z(2^2), Z(2^2) ] ]
gap> ConvertToMatrixRep(m2); # now it works
4
gap> m2;
[ [ Z(2)^0, Z(2)^0 ], [ Z(2^2), Z(2^2) ] ]
gap> RepresentationsOfObject(m2);
[ "IsPositionalObjectRep", "Is8BitMatrixRep" ]
\endexample

Arithmetic operations (see~"Arithmetic for Lists" and the following
sections) preserve the compression status of matrices in the sense that
if all arguments are compressed matrices written over the same field and
the result is a matrix then also the result is a compressed matrix
written over this field.

There are also two operations that are only available for matrices
written over finite fields.

\Declaration{ProjectiveOrder}
\beginexample
gap> ProjectiveOrder([[1,4],[5,2]]*Z(11)^0);
[ 5, Z(11)^5 ]
\endexample

\Declaration{SimultaneousEigenvalues}



Finally, there are two operations that deal with matrices over a ring,
but only care about the residues of their entries modulo some ring element.
In the case of the integers and a prime number $p$, say,
this is effectively computation in a matrix over the prime field
in characteristic $p$.

\Declaration{InverseMatMod}
\beginexample
gap> mat:= [ [ 1, 2 ], [ 3, 4 ] ];;  inv:= InverseMatMod( mat, 5 );
[ [ 3, 1 ], [ 4, 2 ] ]
gap> mat * inv;
[ [ 11, 5 ], [ 25, 11 ] ]
\endexample

\Declaration{NullspaceModQ}
\beginexample
gap> mat:= [ [ 1, 3 ], [ 1, 2 ], [ 1, 1 ] ];;  NullspaceModQ( mat, 5 );
[ [ 0, 0, 0 ], [ 1, 3, 1 ], [ 2, 1, 2 ], [ 4, 2, 4 ], [ 3, 4, 3 ] ]
\endexample

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{Special Multiplication Algorithms for Matrices over GF(2)}

When multiplying two compressed matrices $M$ and $N$ over GF(2) of dimensions
$a\times b$ and $b\times c$, say, where $a$, $b$ and $c$ are all
greater than or equal to 128, \GAP\ by default uses a more
sophisticated matrix multiplication algorithm, in which linear
combinations of groups of 8 rows of $M$ are remembered and re-used in
constructing various rows of the product. This is called level 8
grease. To optimise memory access patterns, these combinations are
stored for $(b+255)/256$ sets of 8 rows at once. This number is called
the blocking level.

These levels of grease and blocking are found experimentally to give
good performance across a range of processors and matrix sizes, but
other levels may do even better in some cases. You can control the
levels exactly using the functions below:

\>PROD_GF2MAT_GF2MAT_SIMPLE( <m1>, <m2> ) F

This function performs the standard unblocked and ungreased matrix
multiplication for matrices of any size.

\>PROD_GF2MAT_GF2MAT_ADVANCED( <m1>, <m2>, <g>, <b> ) F

This function computes the product of <m1> and <m2>, which must be
compressed matrices over GF(2) of compatible dimensions, using level <g>
grease and level <b> blocking.

We plan to include greased blocked matrix multiplication for other
finite fields, and greased blocked algorithms for inversion and other
matrix operations in a future release.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{Block Matrices}

\FileHeader{matblock}[1]

\Declaration{AsBlockMatrix}
\Declaration{BlockMatrix}
\beginexample
gap> M := BlockMatrix([[1,1,[[1, 2],[ 3, 4]]],
>                      [1,2,[[9,10],[11,12]]],
>                      [2,2,[[5, 6],[ 7, 8]]]],2,2);
<block matrix of dimensions (2*2)x(2*2)>
gap> Display(M);
[ [   1,   2,   9,  10 ],
  [   3,   4,  11,  12 ],
  [   0,   0,   5,   6 ],
  [   0,   0,   7,   8 ] ]
\endexample

\Declaration{MatrixByBlockMatrix}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%E