GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
############################################################################
##
## Sparse.gi Gauss package Simon Goertzen
##
## Copyright 2007-2008 Lehrstuhl B f��r Mathematik, RWTH Aachen
##
## Implementation stuff for performing algorithms on sparse matrices.
## Also includes the documentation of these algorithms.
##
#############################################################################
## <#GAPDoc Label="EchelonMat">
## <ManSection Label="echelonmat">
## <Meth Arg="mat" Name="EchelonMat"/>
## <Returns>a record that contains information about an echelonized
## form of the matrix <A>mat</A>.<P/>
## The components of this record are<P/>
## `vectors'<P/>
## the reduced row echelon / hermite form of the matrix <A>mat</A>
## without zero rows.<P/>
## `heads'<P/>
## list that contains at position <i>, if nonzero, the
## number of the row for that the pivot element is in column <i>.
## </Returns>
## <Description>
## computes the reduced row echelon form RREF of a dense or sparse matrix
## <A>mat</A> over a field,
## or the hermite form of a sparse matrix <A>mat</A> over <M>&ZZ; / < p^n ></M>.
## <Example><![CDATA[
## gap> M := [[0,0,0,1,0],[0,1,1,1,1],[1,1,1,1,0]] * One( GF(2) );;
## gap> Display(M);
## . . . 1 .
## . 1 1 1 1
## 1 1 1 1 .
## gap> EchelonMat(M);
## rec( heads := [ 1, 2, 0, 3, 0 ],
## vectors := [ <an immutable GF2 vector of length 5>,
## <an immutable GF2 vector of length 5>,
## <an immutable GF2 vector of length 5> ] )
## gap> Display( last.vectors );
## 1 . . . 1
## . 1 1 . 1
## . . . 1 .
## gap> SM := SparseMatrix( M );
## <a 3 x 5 sparse matrix over GF(2)>
## gap> EchelonMat( SM );
## rec( heads := [ 1, 2, 0, 3, 0 ], vectors := <a 3 x 5 sparse matrix over GF(
## 2)> )
## gap> Display(last.vectors);
## 1 . . . 1
## . 1 1 . 1
## . . . 1 .
## gap> SM := SparseMatrix( [[7,4,5],[0,0,6],[0,4,4]] * One( Integers mod 8 ) );
## <a 3 x 3 sparse matrix over (Integers mod 8)>
## gap> Display( SM );
## 7 4 5
## . . 6
## . 4 4
## gap> EchelonMat( SM );
## rec( heads := [ 1, 2, 3 ],
## vectors := <a 3 x 3 sparse matrix over (Integers mod 8)> )
## gap> Display( last.vectors );
## 1 . 1
## . 4 .
## . . 2
## ]]></Example>
## </Description>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( EchelonMat,
"method for sparse matrices",
[ IsSparseMatrix ],
function( mat )
if IsField( mat!.ring ) then
return EchelonMatDestructive( CopyMat( mat ) );
else
return HermiteMatDestructive( CopyMat( mat ) );
fi;
end
);
## <#GAPDoc Label="EchelonMatTransformation">
## <ManSection>
## <Meth Arg="mat" Name="EchelonMatTransformation" />
## <Returns>a record that contains information about an echelonized
## form of the matrix <A>mat</A>.<P/>
## The components of this record are<P/>
## `vectors'<P/>
## the reduced row echelon / hermite form of the matrix
## <A>mat</A> without zero rows.<P/>
## `heads'<P/>
## list that contains at position <i>, if nonzero, the
## number of the row for that the pivot element is in column <i>.<P/>
## `coeffs'<P/>
## the transformation matrix needed to obtain the RREF from <A>mat</A>.<P/>
## `relations'<P/>
## the kernel of the matrix <A>mat</A> if RingOfDefinition(<A>mat</A>) is a field.
## Otherwise these are only the obvious row relations of <A>mat</A>,
## there might be more kernel vectors - &see; <Ref Func="KernelMat" Style="Number"/>.
## </Returns>
## <Description>
## computes the reduced row echelon form RREF of a dense or sparse matrix
## <A>mat</A> over a field,
## or the hermite form of a sparse matrix <A>mat</A> over <M>&ZZ; / < p^n ></M>.
## In either case, the transformation matrix <M>T</M> is calculated
## as the row union of `coeffs' and `relations'.
## <Example><![CDATA[
## gap> M := [[1,0,1],[1,1,0],[1,0,1],[1,1,0],[1,1,1]] * One( GF(2) );;
## gap> EchelonMatTransformation( M );
## rec(
## coeffs := [ <an immutable GF2 vector of length 5>,
## <an immutable GF2 vector of length 5>,
## <an immutable GF2 vector of length 5> ], heads := [ 1, 2, 3 ],
## relations :=
## [ <an immutable GF2 vector of length 5>,
## <an immutable GF2 vector of length 5> ],
## vectors := [ <an immutable GF2 vector of length 3>,
## <an immutable GF2 vector of length 3>,
## <an immutable GF2 vector of length 3> ] )
## gap> Display(last.vectors);
## 1 . .
## . 1 .
## . . 1
## gap> Display(last.coeffs);
## 1 1 . . 1
## 1 . . . 1
## . 1 . . 1
## gap> Display(last.relations);
## 1 . 1 . .
## . 1 . 1 .
## gap> Display( Concatenation( last.coeffs, last.relations ) * M );
## 1 . .
## . 1 .
## . . 1
## . . .
## . . .
## gap> SM := SparseMatrix( M );
## <a 5 x 3 sparse matrix over GF(2)>
## gap> EchelonMatTransformation( SM );
## rec( coeffs := <a 3 x 5 sparse matrix over GF(2)>,
## heads := [ 1, 2, 3 ],
## relations := <a 2 x 5 sparse matrix over GF(2)>,
## vectors := <a 3 x 3 sparse matrix over GF(2)> )
## gap> Display(last.vectors);
## 1 . .
## . 1 .
## . . 1
## gap> Display(last.coeffs);
## 1 1 . . 1
## 1 . . . 1
## . 1 . . 1
## gap> Display(last.relations);
## 1 . 1 . .
## . 1 . 1 .
## gap> Display( UnionOfRows( last.coeffs, last.relations ) * SM );
## 1 . .
## . 1 .
## . . 1
## . . .
## . . .
## ]]></Example>
## </Description>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( EchelonMatTransformation,
"method for sparse matrices",
[ IsSparseMatrix ],
function( mat )
if IsField( mat!.ring ) then
return EchelonMatTransformationDestructive( CopyMat( mat ) );
else
return HermiteMatTransformationDestructive( CopyMat( mat ) );
fi;
end
);
## <#GAPDoc Label="ReduceMat">
## <ManSection>
## <Meth Arg="A, B" Name="ReduceMat" />
## <Returns>a record with a single component `reduced_matrix' := M.
## M is created by reducing <A>A</A> with <A>B</A>, where <A>B</A> must
## be in Echelon/Hermite form. M will have the same dimensions as <A>A</A>.</Returns>
## <Description>
## <Example><![CDATA[
## gap> M := [[0,0,0,1,0],[0,1,1,1,1],[1,1,1,1,0]] * One( GF(2) );;
## gap> Display(M);
## . . . 1 .
## . 1 1 1 1
## 1 1 1 1 .
## gap> N := [[1,1,0,0,0],[0,0,1,0,1]] * One( GF(2) );;
## gap> Display(N);
## 1 1 . . .
## . . 1 . 1
## gap> ReduceMat(M,N);
## rec(
## reduced_matrix := [ <a GF2 vector of length 5>, <a GF2 vector of length 5>,
## <a GF2 vector of length 5> ] )
## gap> Display(last.reduced_matrix);
## . . . 1 .
## . 1 . 1 .
## . . . 1 1
## gap> SM := SparseMatrix(M); SN := SparseMatrix(N);
## <a 3 x 5 sparse matrix over GF(2)>
## <a 2 x 5 sparse matrix over GF(2)>
## gap> ReduceMat(SM,SN);
## rec( reduced_matrix := <a 3 x 5 sparse matrix over GF(2)> )
## gap> Display(last.reduced_matrix);
## . . . 1 .
## . 1 . 1 .
## . . . 1 1
## ]]></Example>
## </Description>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( ReduceMat,
"for sparse matrices over a ring, second argument must be in REF",
[ IsSparseMatrix, IsSparseMatrix ],
function( mat, N )
if IsField( mat!.ring ) then
return ReduceMatWithEchelonMat( mat, N );
else
return ReduceMatWithHermiteMat( mat, N );
fi;
end
);
## <#GAPDoc Label="ReduceMatTransformation">
## <ManSection>
## <Meth Arg="A, B" Name="ReduceMatTransformation" />
## <Returns>a record with a component `reduced_matrix' := M.
## M is created by reducing <A>A</A> with <A>B</A>, where <A>B</A> must
## be in Echelon/Hermite form. M will have the same dimensions as <A>A</A>.
## In addition to the (identical) output as ReduceMat this record also
## includes the component `transformation', which stores the row operations
## that were needed to reduce <A>A</A> with <A>B</A>. This differs from "normal"
## transformation matrices because only rows of <A>B</A> had to be moved.
## Therefore, the transformation matrix solves M = A + T * B.</Returns>
## <Description>
## <Example><![CDATA[
## gap> M := [[0,0,0,1,0],[0,1,1,1,1],[1,1,1,1,0]] * One( GF(2) );;
## gap> Display(M);
## . . . 1 .
## . 1 1 1 1
## 1 1 1 1 .
## gap> N := [[1,1,0,0,0],[0,0,1,0,1]] * One( GF(2) );;
## gap> Display(N);
## 1 1 . . .
## . . 1 . 1
## gap> ReduceMatTransformation(M,N);
## rec(
## reduced_matrix :=
## [ <a GF2 vector of length 5>, <a GF2 vector of length 5>,
## <a GF2 vector of length 5> ],
## transformation := [ <a GF2 vector of length 2>,
## <a GF2 vector of length 2>, <a GF2 vector of length 2> ] )
## gap> Display(last.reduced_matrix);
## . . . 1 .
## . 1 . 1 .
## . . . 1 1
## gap> Display(last.transformation);
## . .
## . 1
## 1 1
## gap> Display( M + last.transformation * N );
## . . . 1 .
## . 1 . 1 .
## . . . 1 1
## gap> SM := SparseMatrix(M); SN := SparseMatrix(N);
## <a 3 x 5 sparse matrix over GF(2)>
## <a 2 x 5 sparse matrix over GF(2)>
## gap> ReduceMatTransformation(SM,SN);
## rec( reduced_matrix := <a 3 x 5 sparse matrix over GF(2)>,
## transformation := <a 3 x 2 sparse matrix over GF(2)> )
## gap> Display(last.reduced_matrix);
## . . . 1 .
## . 1 . 1 .
## . . . 1 1
## gap> Display(last.transformation);
## . .
## . 1
## 1 1
## gap> Display( SM + last.transformation * SN );
## . . . 1 .
## . 1 . 1 .
## . . . 1 1
## ]]></Example>
## </Description>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( ReduceMatTransformation,
"for sparse matrices over a ring, second argument must be in REF",
[ IsSparseMatrix, IsSparseMatrix ],
function( mat, N )
if IsField( mat!.ring ) then
return ReduceMatWithEchelonMatTransformation( mat, N );
else
return ReduceMatWithHermiteMatTransformation( mat, N );
fi;
end
);
## <#GAPDoc Label="KernelMat">
## <ManSection Label="kernelmat">
## <Func Arg="M" Name="KernelMat"/>
## <Returns>a record with a single component `relations'.</Returns>
## <Description>
## If <A>M</A> is a matrix over a field this is the same output as
## <Ref Meth="EchelonMatTransformation"/> provides in the
## `relations' component, but with less memory and CPU usage.
## If the base ring of <A>M</A> is a non-field, the Kernel might
## have additional generators, which are added to the output.
## <Example><![CDATA[
## gap> M := [[2,1],[0,2]];
## [ [ 2, 1 ], [ 0, 2 ] ]
## gap> SM := SparseMatrix( M * One( GF(3) ) );
## <a 2 x 2 sparse matrix over GF(3)>
## gap> KernelMat(SM);
## rec( relations := <a 0 x 2 sparse matrix over GF(3)> )
## gap> SN := SparseMatrix( M * One( Integers mod 4 ) );
## <a 2 x 2 sparse matrix over (Integers mod 4)>
## gap> KernelMat(SN);
## rec( relations := <a 1 x 2 sparse matrix over (Integers mod 4)> )
## gap> Display(last.relations);
## 2 1
## ]]></Example>
## </Description>
## </ManSection>
## <#/GAPDoc>
##
InstallGlobalFunction( KernelMatSparse,
function( arg )
local M;
M := CopyMat( arg[1] );
if IsField( M!.ring ) and Length( arg ) = 1 then
return KernelEchelonMatDestructive( M, [ 1 .. M!.nrows ] );
elif IsField( M!.ring ) and Length( arg ) > 1 then
return KernelEchelonMatDestructive( M, arg[2] );
elif Length( arg ) = 1 then
return KernelHermiteMatDestructive( M, [ 1 .. M!.nrows ] );
elif Length( arg ) > 1 then
return KernelHermiteMatDestructive( M, arg[2] );
fi;
end
);
## <#GAPDoc Label="Rank">
## <ManSection>
## <Meth Arg="sm[, boundary]" Name="Rank" />
## <Returns>the rank of the sparse matrix <A>sm</A>. Only works for fields.</Returns>
## <Description>
## Computes the rank of a sparse matrix. If the optional argument
## <A>boundary</A> is provided, some algorithms take into account
## the fact that Rank(<A>sm</A>) <= <A>boundary</A>, thus
## possibly terminating earlier.
## <Example><![CDATA[
## gap> M := SparseDiagMat( ListWithIdenticalEntries( 10,
## > SparseMatrix( [[1,1],[1,1]] * One( GF(5) ) ) ) );
## <a 20 x 20 sparse matrix over GF(5)>
## gap> Display(M);
## 1 1 . . . . . . . . . . . . . . . . . .
## 1 1 . . . . . . . . . . . . . . . . . .
## . . 1 1 . . . . . . . . . . . . . . . .
## . . 1 1 . . . . . . . . . . . . . . . .
## . . . . 1 1 . . . . . . . . . . . . . .
## . . . . 1 1 . . . . . . . . . . . . . .
## . . . . . . 1 1 . . . . . . . . . . . .
## . . . . . . 1 1 . . . . . . . . . . . .
## . . . . . . . . 1 1 . . . . . . . . . .
## . . . . . . . . 1 1 . . . . . . . . . .
## . . . . . . . . . . 1 1 . . . . . . . .
## . . . . . . . . . . 1 1 . . . . . . . .
## . . . . . . . . . . . . 1 1 . . . . . .
## . . . . . . . . . . . . 1 1 . . . . . .
## . . . . . . . . . . . . . . 1 1 . . . .
## . . . . . . . . . . . . . . 1 1 . . . .
## . . . . . . . . . . . . . . . . 1 1 . .
## . . . . . . . . . . . . . . . . 1 1 . .
## . . . . . . . . . . . . . . . . . . 1 1
## . . . . . . . . . . . . . . . . . . 1 1
## gap> Rank(M);
## 10
## ]]></Example>
## </Description>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( Rank,
"method for sparse matrices",
[ IsSparseMatrix ],
function( mat )
if mat!.ring = GF(2) then
return Length( RankOfIndicesListList( mat!.indices ).vectors );
elif IsField( mat!.ring ) then
return RankDestructive( CopyMat( mat ), mat!.ncols );
else
Error( "no Rank method for matrices over ", mat!.ring, "!" );
fi;
end
);
##
InstallOtherMethod( Rank,
"method for sparse matrices with upper boundary",
[ IsSparseMatrix, IsInt ],
function( mat, upper_boundary )
if IsField( mat!.ring ) then
return RankDestructive( CopyMat( mat ), Minimum( upper_boundary, mat!.ncols ) );
else
Error( "no Rank method for matrices over ", mat!.ring, "!" );
fi;
end
);