|
# SPDX-License-Identifier: GPL-2.0-or-later
# Gauss: Extended Gauss functionality for GAP
#
# Implementations
#
## Implementation stuff for Gauss with sparse matrices.
##
DeclareRepresentation( "IsSparseMatrixRep",
IsSparseMatrix, [ "nrows", "ncols", "indices", "entries", "ring" ] );
##
DeclareRepresentation( "IsSparseMatrixGF2Rep",
IsSparseMatrix, [ "nrows", "ncols", "indices", "ring" ] );
##
BindGlobal( "TheFamilyOfSparseMatrices",
NewFamily( "TheFamilyOfSparseMatrices" ) );
##
BindGlobal( "TheTypeSparseMatrix",
NewType( TheFamilyOfSparseMatrices, IsSparseMatrixRep ) );
##
BindGlobal( "TheTypeSparseMatrixGF2",
NewType( TheFamilyOfSparseMatrices, IsSparseMatrixGF2Rep ) );
## <#GAPDoc Label="SparseMatrix">
## <ManSection >
## <Func Arg="mat[, R]" Name="SparseMatrix" Label="constructor using gap matrices" />
## <Returns>a sparse matrix over the ring <A>R</A></Returns>
## <Func Arg="nrows, ncols, indices" Name="SparseMatrix" Label="constructor using indices" />
## <Returns>a sparse matrix over GF(2)</Returns>
## <Func Arg="nrows, ncols, indices, entries[, R]" Name="SparseMatrix" Label="constructor using indices and entries" />
## <Returns>a sparse matrix over the ring <A>R</A></Returns>
## <Description>
## The sparse matrix constructor.
## In the one-argument form the SparseMatrix constructor converts a &GAP; matrix
## to a sparse matrix. If not provided the base ring <A>R</A> is found automatically.
## For the multi-argument form <A>nrows</A> and <A>ncols</A> are the dimensions
## of the matrix. <A>indices</A> must be a list of length <A>nrows</A> containing
## lists of the column indices of the matrix in ascending order.
##
## <Example><![CDATA[
## gap> M := [ [ 0 , 1 ], [ 3, 0 ] ] * One( GF(2) );
## [ [ 0*Z(2), Z(2)^0 ], [ Z(2)^0, 0*Z(2) ] ]
## gap> SM := SparseMatrix( M );
## <a 2 x 2 sparse matrix over GF(2)>
## gap> IsSparseMatrix( SM );
## true
## gap> Display( SM );
## . 1
## 1 .
## gap> SN := SparseMatrix( 2, 2, [ [ 2 ], [ 1 ] ] );
## <a 2 x 2 sparse matrix over GF(2)>
## gap> SN = SM;
## true
## gap> SN := SparseMatrix( 2, 3,
## > [ [ 2 ], [ 1, 3 ] ],
## > [ [ 1 ], [ 3, 2 ] ] * One( GF(5) ) );
## <a 2 x 3 sparse matrix over GF(5)>
## gap> Display( SN );
## . 1 .
## 3 . 2
## ]]></Example>
## </Description>
## </ManSection>
## <#/GAPDoc>
##
InstallGlobalFunction( SparseMatrix,
function( arg )
local nargs, M, nrows, ncols, i, j, indices, entries, ring;
nargs := Length( arg );
if IsRing( arg[ nargs ] ) then
ring := arg[ nargs ];
nargs := nargs - 1;
fi;
if nargs >= 3 then
nrows := arg[1];
ncols := arg[2];
indices := arg[3];
if not IsList(indices) or Length(indices) <> nrows then
Error("<indices> must be a list of length <nrows>");
fi;
if not ForAll(indices, l -> IsSet(l) and IsPositionsList(l)) then
Error("<indices> must consist of strictly ascending lists of positive integers");
fi;
fi;
if nargs = 3 then
return Objectify( TheTypeSparseMatrixGF2, rec( nrows := nrows, ncols := ncols, indices := indices, ring := GF(2) ) );
elif nargs = 4 then
if not IsBound( ring ) then
ring := FindRing( arg[4] );
fi;
return Objectify( TheTypeSparseMatrix, rec( nrows := nrows, ncols := ncols, indices := indices, entries := arg[4], ring := ring ) );
elif nargs = 5 then
return Objectify( TheTypeSparseMatrix, rec( nrows := nrows, ncols := ncols, indices := indices, entries := arg[4], ring := arg[5] ) );
fi;
if nargs > 1 then
Error( "wrong number of arguments! SparseMatrix expects nrows, ncols, indices, [entries], [ring]!" );
elif not IsList( arg[ 1 ] ) then
Error( "SparseMatrix constructor with 1 argument expects a (dense) Matrix or List!" );
fi;
M := arg[1];
nrows := Length( M );
if nrows = 0 then
return SparseMatrix( 0, 0, [], [], "unknown" );
fi;
ncols := Length( M[1] );
if ncols = 0 then
return SparseMatrix( nrows, 0, List( [ 1 .. nrows ], i -> [] ), List( [ 1 .. nrows ], i -> [] ), "unknown" );
fi;
if not IsBound( ring ) then
ring := FindRing( List( M, i -> Filtered( i, j -> not IsZero(j) ) ) );
fi;
if ring = GF(2) then
indices := List( [ 1 .. nrows ], i -> Filtered( [ 1 .. ncols ], j -> IsOne( M[i][j] ) ) );
return SparseMatrix( nrows, ncols, indices );
fi;
indices := [];
entries := [];
for i in [ 1..nrows ] do
indices[i] := [];
entries[i] := [];
for j in [ 1..ncols ] do
if not IsZero( M[i][j] ) then
Add( indices[i], j );
Add( entries[i], M[i][j] );
fi;
od;
od;
return SparseMatrix( nrows, ncols, indices, entries, ring );
end
);
## <#GAPDoc Label="ConvertSparseMatrixToMatrix">
## <ManSection >
## <Meth Arg="sm" Name="ConvertSparseMatrixToMatrix" />
## <Returns>a &GAP; matrix, [], or a list of empty lists</Returns>
## <Description>
## This function converts the sparse matrix <A>sm</A> into a &GAP; matrix.
## In case of <C>nrows(sm)=0</C> or <C>ncols(sm)=0</C> the return value is the
## empty list or a list of empty lists, respectively.
## <Example><![CDATA[
## gap> M := [ [ 0 , 1 ], [ 3, 0 ] ] * One( GF(3) );
## [ [ 0*Z(3), Z(3)^0 ], [ 0*Z(3), 0*Z(3) ] ]
## gap> SM := SparseMatrix( M );
## <a 2 x 2 sparse matrix over GF(3)>
## gap> N := ConvertSparseMatrixToMatrix( SM );
## [ [ 0*Z(3), Z(3)^0 ], [ 0*Z(3), 0*Z(3) ] ]
## gap> M = N;
## true
## ]]></Example>
## </Description>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( ConvertSparseMatrixToMatrix,
[ IsSparseMatrix ],
function( SM )
local indices, entries, i, j, ring, M;
if SM!.nrows = 0 then
return [ ];
elif SM!.ncols = 0 then
return List( [ 1 .. SM!.nrows ], i -> [] );
fi;
ring := SM!.ring;
if ring = "unknown" then
ring := Rationals;
fi;
indices := SM!.indices;
entries := SM!.entries;
M := NullMat( SM!.nrows, SM!.ncols, ring );
for i in [ 1 .. SM!.nrows ] do
for j in [ 1 .. Length( indices[i] ) ] do
M[ i ][ indices[i][j] ] := entries[i][j];
od;
od;
return M;
end
);
## <#GAPDoc Label="CopyMat">
## <ManSection >
## <Meth Arg="sm" Name="CopyMat" />
## <Returns>a shallow copy of the sparse matrix <A>sm</A></Returns>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( CopyMat,
[ IsSparseMatrix ],
function( M )
local indices, entries, i;
indices := [];
entries := [];
for i in [ 1 .. M!.nrows ] do
indices[i] := ShallowCopy( M!.indices[i] );
entries[i] := ShallowCopy( M!.entries[i] );
od;
return SparseMatrix( M!.nrows, M!.ncols, indices, entries, M!.ring );
end
);
## <#GAPDoc Label="GetEntry">
## <ManSection >
## <Meth Arg="sm, i, j" Name="GetEntry" />
## <Returns>a ring element.</Returns>
## <Description>
## This returns the entry <C>sm[i,j]</C> of the sparse matrix <A>sm</A>
## </Description>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( GetEntry,
[ IsSparseMatrix, IsInt, IsInt ],
function( M, i, j )
local p;
p := PositionSet( M!.indices[i], j );
if p = fail then
return Zero( M!.ring );
else
return M!.entries[i][p];
fi;
end
);
## <#GAPDoc Label="SetEntry">
## <ManSection >
## <Meth Arg="sm, i, j, elm" Name="SetEntry" />
## <Returns>nothing.</Returns>
## <Description>
## This sets the entry <C>sm[i,j]</C> of the sparse matrix <A>sm</A> to <A>elm</A>.
## </Description>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( SetEntry,
[ IsSparseMatrix, IsInt, IsInt, IsRingElement ],
function( M, i, j, e )
local ring, pos, res;
ring := M!.ring;
if not e in ring then
Error( "the element has to be in ", ring, "!" );
fi;
pos := PositionSorted( M!.indices[i], j );
if IsBound( M!.indices[i][pos] ) and M!.indices[i][pos] = j then
if e = Zero( ring ) then
Remove( M!.indices[i], pos );
Remove( M!.entries[i], pos );
else
M!.entries[i][pos] := e;
fi;
else
if e <> Zero( ring ) then
Add( M!.indices[i], j, pos );
Add( M!.entries[i], e, pos );
fi;
fi;
end
);
## <#GAPDoc Label="AddToEntry">
## <ManSection >
## <Meth Arg="sm, i, j, elm" Name="AddToEntry" />
## <Returns><K>true</K> or a ring element</Returns>
## <Description>
## AddToEntry adds the element <A>elm</A> to the sparse matrix <A>sm</A> at the
## <A>(i,j)</A>-th position. This is a Method with a side effect which
## returns true if you tried to add zero or the sum of <C>sm[i,j]</C> and
## <A>elm</A> otherwise. Please use this method whenever possible.
## </Description>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( AddToEntry,
[ IsSparseMatrix, IsInt, IsInt, IsRingElement ],
function( M, i, j, e )
local ring, pos, res;
ring := M!.ring;
if not e in ring then
Error( "the element has to be in ", ring, "!" );
fi;
if e = Zero( ring ) then
return true;
fi;
pos := PositionSorted( M!.indices[i], j );
if IsBound( M!.indices[i][pos] ) and M!.indices[i][pos] = j then
res := M!.entries[i][pos] + e;
if res = Zero( ring ) then
Remove( M!.indices[i], pos );
Remove( M!.entries[i], pos );
else
M!.entries[i][pos] := res;
fi;
return res;
else
Add( M!.indices[i], j, pos );
Add( M!.entries[i], e, pos );
return e;
fi;
end
);
###############################
## View and Display methods: ##
###############################
##
InstallMethod( ViewObj,
[ IsSparseMatrix ],
function( M )
Print( "<a ", M!.nrows, " x ", M!.ncols, " sparse matrix over ", M!.ring, ">" );
end
);
##
InstallMethod( Display,
[ IsSparseMatrix ],
function( M )
local str, ws, i, last, j;
if M!.nrows = 0 or M!.ncols = 0 then
Print( "(an empty ", M!.nrows, " x ", M!.ncols, " matrix)" );
elif Characteristic( M!.ring ) = 0 or ( HasDegreeOverPrimeField( M!.ring ) and DegreeOverPrimeField( M!.ring ) > 1 ) then
Display( ConvertSparseMatrixToMatrix( M ) );
else
str := "";
ws := ListWithIdenticalEntries( Length( String( Int( - One( M!.ring ) ) ) ), ' ' );
for i in [ 1 .. M!.nrows ] do
last := 0;
for j in [ 1 .. Length( M!.indices[i] ) ] do
str := Concatenation( str, Concatenation( ListWithIdenticalEntries( M!.indices[i][j] - 1 - last, Concatenation( ws, "." ) ) ), ws{ [ 1 .. Length( ws ) + 1 - Length( String( Int( M!.entries[i][j] ) ) ) ] }, String( Int( M!.entries[i][j] ) ) );
last := M!.indices[i][j];
od;
str := Concatenation( str, Concatenation( ListWithIdenticalEntries( M!.ncols - last, Concatenation( ws, "." ) ) ), "\n" );
od;
Print( str );
fi;
end
);
##
InstallMethod( PrintObj,
[ IsSparseMatrix ],
function( M )
Print( "SparseMatrix( ", M!.nrows, ", ", M!.ncols, ", ", M!.indices, ", ", M!.entries, ", ", M!.ring, " )" );
end
);
###############################
##
InstallMethod( FindRing,
[ IsList ],
function( entries )
local found_ring, i, ring;
found_ring := false;
i := 1;
while found_ring = false and i <= Length( entries ) do
if IsBound( entries[i][1] ) then
ring := DefaultRing( entries[i][1] );
found_ring := true;
fi;
i := i + 1;
od;
if found_ring = true then
return ring;
else
return "unknown";
fi;
end
);
## <#GAPDoc Label="SparseZeroMatrix">
## <ManSection >
## <Func Arg="nrows[, ring]" Name="SparseZeroMatrix" Label="constructor using number of rows" />
## <Returns>a sparse <<A>nrows</A> x <A>nrows</A>> zero matrix over the ring <A>ring</A></Returns>
## <Func Arg="nrows, ncols[, ring]" Name="SparseZeroMatrix" Label="constructor using number of rows and columns"/>
## <Returns>a sparse <<A>nrows</A> x <A>ncols</A>> zero matrix over the ring <A>ring</A></Returns>
## </ManSection>
## <#/GAPDoc>
##
InstallGlobalFunction( SparseZeroMatrix,
function( arg )
local nargs, ring;
nargs := Length( arg );
if IsRing( arg[ nargs ] ) or arg[ nargs ] = "unknown" then
ring := arg[ Length( arg ) ];
nargs := nargs - 1;
else
ring := "unknown";
fi;
if ring = GF(2) then
if nargs = 1 then
return SparseMatrix( arg[1], arg[1], List( [ 1 .. arg[1] ], i -> [] ) );
elif nargs = 2 then
return SparseMatrix( arg[1], arg[2], List( [ 1 .. arg[1] ], i -> [] ) );
fi;
elif nargs = 1 then
return SparseMatrix( arg[1], arg[1], List( [ 1 .. arg[1] ], i -> [] ), List( [ 1 .. arg[1] ], i -> [] ), ring );
elif nargs = 2 then
return SparseMatrix( arg[1], arg[2], List( [ 1 .. arg[1] ], i -> [] ), List( [ 1 .. arg[1] ], i -> [] ), ring );
fi;
Error( "wrong number of arguments in SparseZeroMatrix!" );
end
);
## <#GAPDoc Label="SparseIdentityMatrix">
## <ManSection >
## <Func Arg="dim[, ring]" Name="SparseIdentityMatrix" />
## <Returns>a sparse <<A>dim</A> x <A>dim</A>> identity matrix over the ring <A>ring</A>.
## If no ring is specified (one should try to avoid this if possible)
## the Rationals are the default ring.</Returns>
## </ManSection>
## <#/GAPDoc>
##
InstallGlobalFunction( SparseIdentityMatrix,
function( arg )
local nargs, ring, indices, entries;
nargs := Length( arg );
if IsRing( arg[ nargs ] ) then
ring := arg[ Length( arg ) ];
nargs := nargs - 1;
else
ring := Rationals;
fi;
if nargs = 1 then
indices := List( [ 1 .. arg[1] ], i -> [i] );
if ring = GF(2) then
return SparseMatrix( arg[1], arg[1], indices );
else
entries := List( [ 1 .. arg[1] ], i -> [ One( ring ) ] );
return SparseMatrix( arg[1], arg[1], indices, entries, ring );
fi;
else
Error( "wrong number of arguments in SparseIdentityMatrix!" );;
fi;
end
);
##
InstallMethod( \=,
[ IsSparseMatrix, IsSparseMatrix ],
function( A, B )
return A!.nrows = B!.nrows and
A!.ncols = B!.ncols and
A!.ring = B!.ring and
A!.indices = B!.indices and
A!.entries = B!.entries;
end
);
## <#GAPDoc Label="TransposedSparseMat">
## <ManSection >
## <Meth Arg="sm" Name="TransposedSparseMat" />
## <Returns>the transposed matrix of the sparse matrix <A>sm</A></Returns>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( TransposedSparseMat,
[ IsSparseMatrix ],
function( M )
local T, i, j;
T := SparseZeroMatrix( M!.ncols, M!.nrows, M!.ring );
for i in [ 1 .. M!.nrows ] do
for j in [ 1 .. Length( M!.indices[i] ) ] do
Add( T!.indices[ M!.indices[i][j] ], i );
Add( T!.entries[ M!.indices[i][j] ], M!.entries[i][j] );
od;
od;
return T;
end
);
## <#GAPDoc Label="CertainRows">
## <ManSection >
## <Meth Arg="sm, list" Name="CertainRows" />
## <Returns>the submatrix <C>sm{[list]}</C> as a sparse matrix</Returns>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( CertainRows,
[ IsSparseMatrix, IsList ],
function( M, L )
return SparseMatrix( Length( L ), M!.ncols, M!.indices{ L }, M!.entries{ L }, M!.ring );
end
);
## <#GAPDoc Label="CertainColumns">
## <ManSection >
## <Meth Arg="sm, list" Name="CertainColumns" />
## <Returns>the submatrix <C>sm{[1..nrows(sm)]}{[list]}</C> as a sparse matrix</Returns>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( CertainColumns,
[ IsSparseMatrix, IsList ],
function( M, L )
local indices, entries, list, i, j, column, p;
indices := List( [ 1 .. M!.nrows ], i -> [] );
entries := List( [ 1 .. M!.nrows ], i -> [] );
for i in [ 1 .. M!.nrows ] do
for j in [ 1 .. Length( L ) ] do
column := L[j];
p := PositionSet( M!.indices[i], column);
if p <> fail then
Add( indices[i], j );
Add( entries[i], M!.entries[i][p] );
fi;
od;
od;
return SparseMatrix( M!.nrows, Length( L ), indices, entries, M!.ring );
end
);
## <#GAPDoc Label="SparseUnionOfRows">
## <ManSection>
## <Func Arg="L" Name="SparseUnionOfRows" Label="for a list of sparse matrices"/>
## <Returns>a sparse matrix</Returns>
## <Description>
## Stack the sparse matrices in the non-empty list <A>L</A>.
## </Description>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( SparseUnionOfRows,
"of a list of sparse matrices",
[ IsList ],
function( L )
local R, nr_cols, indices;
R := L[1]!.ring;
nr_cols := L[1]!.ncols;
indices := Concatenation( List( L, x -> x!.indices ) );
if IsSparseMatrixGF2Rep( L[1] ) then
return SparseMatrix( Sum( L, x -> x!.nrows ), nr_cols, indices, R );
else
return SparseMatrix( Sum( L, x -> x!.nrows ), nr_cols, indices, Concatenation( List( L, x -> x!.entries ) ), R );
fi;
end
);
## <#GAPDoc Label="SparseUnionOfColumns">
## <ManSection>
## <Func Arg="L" Name="SparseUnionOfColumns" Label="for a list of sparse matrices"/>
## <Returns>a sparse matrix</Returns>
## <Description>
## Augment the sparse matrices in the non-empty list <A>L</A>.
## </Description>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( SparseUnionOfColumns,
"of a list of sparse matrices",
[ IsList ],
function( L )
local R, nr_rows, indices, c, l, i;
R := L[1]!.ring;
nr_rows := L[1]!.nrows;
indices := List( [ 1 .. nr_rows ], x -> [] );
c := 0;
for l in L do
for i in [ 1 .. nr_rows ] do
indices[i] := Concatenation( indices[i], l!.indices[i] + c );
od;
c := c + l!.ncols;
od;
if IsSparseMatrixGF2Rep( L[1] ) then
return SparseMatrix( nr_rows, Sum( L, x -> x!.ncols ), indices, R );
else
return SparseMatrix( nr_rows, Sum( L, x -> x!.ncols ), indices, List( [ 1 .. nr_rows ], i -> Concatenation( List( L, x -> x!.entries[i] ) ) ), R );
fi;
end
);
## <#GAPDoc Label="SparseDiagMat">
## <ManSection >
## <Func Arg="list" Name="SparseDiagMat" />
## <Returns>the block diagonal matrix composed of the
## sparse matrices in <A>list</A></Returns>
## </ManSection>
## <#/GAPDoc>
##
InstallGlobalFunction( SparseDiagMat,
function( e )
if Length( e ) = 1 then
return e[1];
elif Length( e ) > 2 then
return SparseDiagMat( Concatenation( [ SparseDiagMat( e{[1,2]} ) ], e{[ 3 .. Length( e ) ]} ) );
else
if IsSparseMatrixGF2Rep( e[1] ) then
return SparseMatrix( e[1]!.nrows + e[2]!.nrows, e[1]!.ncols + e[2]!.ncols, Concatenation( e[1]!.indices, e[2]!.indices + e[1]!.ncols ) );
else
return SparseMatrix( e[1]!.nrows + e[2]!.nrows, e[1]!.ncols + e[2]!.ncols, Concatenation( e[1]!.indices, e[2]!.indices + e[1]!.ncols ), Concatenation( e[1]!.entries, e[2]!.entries ), e[1]!.ring );
fi;
fi;
end
);
##
InstallMethod( \*,
[ IsSparseMatrix, IsRingElement ],
function( A, a )
local i, m;
if IsZero( a ) then
return SparseZeroMatrix( A!.nrows, A!.ncols, A!.ring );
else
return SparseMatrix( A!.nrows, A!.ncols, A!.indices, A!.entries * a, A!.ring );
fi;
end
);
##
InstallMethod( \*,
[ IsRingElement, IsSparseMatrix ],
function( a, A )
local i, m;
if IsZero( a ) then
return SparseZeroMatrix( A!.nrows, A!.ncols, A!.ring );
else
return SparseMatrix( A!.nrows, A!.ncols, A!.indices, A!.entries * a, A!.ring );
fi;
end
);
##
InstallMethod( \*,
[ IsSparseMatrix, IsSparseMatrix ],
function( A, B )
local C, i, j, rownr, m;
if A!.ncols <> B!.nrows or A!.ring <> B!.ring then
return fail;
fi;
C := SparseZeroMatrix( A!.nrows, B!.ncols, A!.ring );
for i in [ 1 .. C!.nrows ] do
for j in [ 1 .. Length( A!.indices[i] ) ] do
rownr := A!.indices[i][j];
m := MultRow( B!.indices[ rownr ], B!.entries[ rownr ], A!.entries[i][j] );
AddRow( m.indices, m.entries, C!.indices, C!.entries, i );
od;
od;
return C;
end
);
##
InstallMethod( \+,
[ IsSparseMatrix, IsSparseMatrix ],
function( A, B )
local C, i;
C := CopyMat( A );
for i in [ 1 .. C!.nrows ] do
AddRow( B!.indices[i], B!.entries[i], C!.indices, C!.entries, i );
od;
return C;
end
);
##
InstallMethod( \-,
[ IsSparseMatrix, IsSparseMatrix ],
function( A, B )
return A + ( - One( B!.ring ) ) * B;
end
);
InstallOtherMethod( AdditiveInverseMutable,# strangely, doesn't have the filters for arith.op
[IsSparseMatrix],
function(A)
return SparseMatrix( A!.nrows, A!.ncols, A!.indices, -A!.entries, A!.ring );
end);
InstallOtherMethod(ComplexConjugate,
[IsSparseMatrix],
function(A)
return SparseMatrix(A!.nrows,A!.ncols,A!.indices,ComplexConjugate(A!.entries),A!.ring);
end);
# no need to give it a new name -- here make the usual one available
InstallOtherMethod(TransposedMat, [IsSparseMatrix], TransposedSparseMat);
## <#GAPDoc Label="Nrows">
## <ManSection >
## <Meth Arg="sm" Name="Nrows" />
## <Returns>the number of rows of the sparse matrix <A>sm</A>.
## This should be preferred to the equivalent <C>sm!.nrows</C>.</Returns>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( Nrows,
[ IsSparseMatrix ],
function( M )
return M!.nrows;
end
);
## <#GAPDoc Label="Ncols">
## <ManSection >
## <Meth Arg="sm" Name="Ncols" />
## <Returns>the number of columns of the sparse matrix <A>sm</A>.
## This should be preferred to the equivalent <C>sm!.ncols</C>.</Returns>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( Ncols,
[ IsSparseMatrix ],
function( M )
return M!.ncols;
end
);
## <#GAPDoc Label="IndicesOfSparseMatrix">
## <ManSection >
## <Meth Arg="sm" Name="IndicesOfSparseMatrix" />
## <Returns>the indices of the sparse matrix <A>sm</A> as a ListList.
## This should be preferred to the equivalent <C>sm!.indices</C>.</Returns>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( IndicesOfSparseMatrix,
[ IsSparseMatrix ],
function( M )
return M!.indices;
end
);
## <#GAPDoc Label="EntriesOfSparseMatrix">
## <ManSection >
## <Meth Arg="sm" Name="EntriesOfSparseMatrix" />
## <Returns>the entries of the sparse matrix <A>sm</A> as a ListList of ring elements.
## This should be preferred to the equivalent <C>sm!.entries</C> and has the additional
## advantage of working for sparse matrices over GF(2) which do not have any entries.</Returns>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( EntriesOfSparseMatrix,
[ IsSparseMatrix ],
function( M )
if IsSparseMatrixGF2Rep( M ) then
return List( M!.indices, row -> ListWithIdenticalEntries( Length( row ), One( GF(2) ) ) );
else
return M!.entries;
fi;
end
);
## <#GAPDoc Label="RingOfDefinition">
## <ManSection >
## <Meth Arg="sm" Name="RingOfDefinition" />
## <Returns>the base ring of the sparse matrix <A>sm</A>.
## This should be preferred to the equivalent <C>sm!.ring</C>.</Returns>
## </ManSection>
## <#/GAPDoc>
##
InstallMethod( RingOfDefinition,
[ IsSparseMatrix ],
function( M )
return M!.ring;
end
);
##
InstallMethod( IsSparseZeroMatrix,
[ IsSparseMatrix ],
function( M )
return List( [ 1 .. Length( M!.indices ) ], i -> Length( M!.indices[i] ) ) = List( [ 1 .. Length( M!.indices ) ], i -> 0 );
end
);
##
InstallMethod( IsSparseIdentityMatrix,
[ IsSparseMatrix ],
function( M )
local one, i;
one := One( M!.ring );
for i in [ 1 .. M!.nrows ] do
if M!.indices[i] <> [i] or M!.entries[i] <> [ one ] then
return false;
fi;
od;
return true;
end
);
##
InstallMethod( IsSparseDiagonalMatrix,
[ IsSparseMatrix ],
function( M )
local i;
for i in [ 1 .. M!.nrows ] do
if M!.indices[i] <> [i] then
return false;
fi;
od;
return true;
end
);
##
InstallMethod( SparseKroneckerProduct,
[ IsSparseMatrix, IsSparseMatrix ],
function( A, B )
local indices, entries, i1, i2, rowindex, j1, j2, prod;
indices := [];
entries := [];
for i1 in [ 1 .. A!.nrows ] do
for i2 in [ 1 .. B!.nrows ] do
rowindex := ( i1 - 1 ) * B!.nrows + i2;
indices[ rowindex ] := [];
entries[ rowindex ] := [];
for j1 in [ 1 .. Length( A!.indices[i1] ) ] do
for j2 in [ 1.. Length( B!.indices[i2] ) ] do
prod := A!.entries[i1][j1] * B!.entries[i2][j2];
if not IsZero( prod ) then
Add( indices[ rowindex ], ( A!.indices[i1][j1] - 1 ) * B!.ncols + B!.indices[i2][j2] );
Add( entries[ rowindex ], prod );
fi;
od;
od;
od;
od;
return SparseMatrix( A!.nrows * B!.nrows, A!.ncols * B!.ncols, indices, entries, A!.ring );
end
);
##
InstallMethod( SparseDualKroneckerProduct,
[ IsSparseMatrix, IsSparseMatrix ],
function( B, A )
local indices, entries, i1, i2, rowindex, j1, j2, prod;
indices := [];
entries := [];
for i1 in [ 1 .. A!.nrows ] do
for i2 in [ 1 .. B!.nrows ] do
rowindex := ( i1 - 1 ) * B!.nrows + i2;
indices[ rowindex ] := [];
entries[ rowindex ] := [];
for j1 in [ 1 .. Length( A!.indices[i1] ) ] do
for j2 in [ 1.. Length( B!.indices[i2] ) ] do
prod := B!.entries[i2][j2] * A!.entries[i1][j1];
if not IsZero( prod ) then
Add( indices[ rowindex ], ( A!.indices[i1][j1] - 1 ) * B!.ncols + B!.indices[i2][j2] );
Add( entries[ rowindex ], prod );
fi;
od;
od;
od;
od;
return SparseMatrix( A!.nrows * B!.nrows, A!.ncols * B!.ncols, indices, entries, A!.ring );
end
);
##
InstallMethod( SparseZeroRows,
[ IsSparseMatrix ],
function( M )
return Filtered( [ 1 .. M!.nrows ], i -> Length( M!.indices[i] ) = 0 );
end
);
##
InstallMethod( SparseZeroColumns,
[ IsSparseMatrix ],
function( M )
return Difference( [ 1 .. M!.ncols ], Union( M!.indices ) );
end
);
##
InstallMethod( MultRow, #no side effect!
[ IsList, IsList, IsRingElement ],
function( indices, entries, x )
local prod, list;
if IsUnit( x ) then
return rec( indices := indices, entries := entries * x );
else
prod := entries * x;
list := Filtered( [ 1 .. Length( prod ) ], i -> not IsZero( prod[i] ) );
return rec( indices := indices{ list }, entries := prod{ list } );
fi;
end
);
InstallOtherMethod(AddRow,[IsList,IsList,IsList,IsList,IsInt],
function (row1_indices, row1_entries, p_sum_indices, p_sum_entries, i)
# (row1_indices, row1_entries) is row1
# (p_sum_indices[i], p_sum_entries[i]) is row2
# modifies p_sum_indices[i] and p_sum_entriesi[] to make them equal row1 + row2
local pos1, pos2, len1, len2, row2_indices, row2_entries, sum_indices, sum_entries, index1, index2, s;
len1 := Length(row1_indices);
if len1 = 0 then return; fi;
row2_indices := p_sum_indices[i]; row2_entries := p_sum_entries[i];
len2 := Length(row2_indices);
pos1 := 1;
pos2 := 1;
p_sum_indices[i] := []; sum_indices := p_sum_indices[i];
p_sum_entries[i] := []; sum_entries := p_sum_entries[i];
while true do
if pos1 > len1 then
Append(sum_indices, row2_indices{[pos2..len2]});
Append(sum_entries, row2_entries{[pos2..len2]});
return;
fi;
if pos2 > len2 then
Append(sum_indices, row1_indices{[pos1..len1]});
Append(sum_entries, row1_entries{[pos1..len1]});
return;
fi;
index1 := row1_indices[pos1];
index2 := row2_indices[pos2];
if index1 > index2 then
Add(sum_indices, index2);
Add(sum_entries, row2_entries[pos2]);
pos2 := pos2 + 1;
elif index1 < index2 then
Add(sum_indices, index1);
Add(sum_entries, row1_entries[pos1]);
pos1 := pos1 + 1;
else
s := row1_entries[pos1] + row2_entries[pos2];
if not IsZero(s) then
Add(sum_indices, index1);
Add(sum_entries, s);
fi;
pos1 := pos1 + 1;
pos2 := pos2 + 1;
fi;
od;
end);
##
InstallMethod( AddRow, #with desired side effect!
[ IsList, IsList, IsList, IsList ],
function( row1_indices, row1_entries, row2_indices, row2_entries )
local m, zero, j, i, index1, index2;
m := Length( row1_indices );
if m = 0 then
return rec( indices := row2_indices, entries := row2_entries );
fi;
zero := Zero( row1_entries[1] );
i := 1;
j := 1;
while i <= m do
if j > Length( row2_indices ) then
Append( row2_indices, row1_indices{[ i .. m ]} );
Append( row2_entries, row1_entries{[ i .. m ]} );
break;
fi;
index1 := row1_indices[i];
index2 := row2_indices[j];
if index1 > index2 then
j := j + 1;
elif index1 < index2 then
Add( row2_indices, index1, j );
Add( row2_entries, row1_entries[i], j );
i := i + 1;
elif index1 = index2 then
row2_entries[j] := row2_entries[j] + row1_entries[i];
if row2_entries[j] = zero then
Remove( row2_entries, j );
Remove( row2_indices, j );
fi;
i := i + 1;
fi;
od;
return rec( indices := row2_indices, entries := row2_entries );
end
);
InstallMethod( IsRowOfSparseMatrix,
[ IsSparseMatrix, IsSparseMatrix ],
function(mat, row)
local pos;
pos := Position(mat!.indices, row!.indices[1]);
if pos = fail then
return false;
fi;
if mat!.entries[pos] = row!.entries[1] then
return true;
else
return false;
fi;
end );
[ Dauer der Verarbeitung: 0.35 Sekunden
(vorverarbeitet)
]
|
