Make the meataxe faster by using AddMatrix, AddVector, MultMatrix, MultVector (gap-system#6289)

Consider this setup:

    gap> n:=10;;   # increase this to make the effect even stronger
    gap> G:=GL(56,GF(25));;
    gap> H:=Subgroup(G, Concatenation(GeneratorsOfGroup(G),List([1..n],i->PseudoRandom(G))));;
    gap> dn:=g ->DirectSumGModule(NaturalGModule(g),NaturalGModule(g));;
    gap>
    gap> Reset(GlobalMersenneTwister,1);;Reset(GlobalRandomSource,1);;

Then before this PR:

    gap> MTX.Indecomposition(dn(H));; time;
    9326

After this PR:

    gap> MTX.Indecomposition(dn(H));; time;
    3845

To achieve a substantial speedup beyond this probably requires algorithmic
improvements.

Author: fingolfin

lib/matobj.gi

@@ -1931,6 +1931,18 @@ InstallEarlyMethod( AddMatrix,
     fi;
   end );
 
+InstallMethod( AddMatrix, "for a mutable IsRowListMatrix and a IsRowListMatrix",
+  [ IsRowListMatrix and IsMutable, IsRowListMatrix ],
+  function( dstmat, srcmat )
+    local i;
+    if DimensionsMat(dstmat) <> DimensionsMat(srcmat) then
+      Error("AddMatrix: matrices must have the same dimensions");
+    fi;
+    for i in [1..NrRows(dstmat)] do
+      AddRowVector(dstmat[i], srcmat[i]);
+    od;
+  end );
+
 InstallMethod( AddMatrix, "for a mutable 8bit matrix and an 8bit matrix",
   [ Is8BitMatrixRep and IsMutable, Is8BitMatrixRep ],
   function( dstmat, srcmat )
@@ -1939,7 +1951,7 @@ InstallMethod( AddMatrix, "for a mutable 8bit matrix and an 8bit matrix",
       Error("AddMatrix: matrices must have the same dimensions");
     fi;
     for i in [1..NrRows(dstmat)] do
-      ADD_COEFFS_VEC8BIT_2(dstmat[i], srcmat[i]);
+      ADD_ROWVECTOR_VEC8BITS_2(dstmat[i], srcmat[i]);
     od;
   end );
 
@@ -1977,6 +1989,18 @@ InstallEarlyMethod( AddMatrix,
     fi;
   end );
 
+InstallMethod( AddMatrix, "for a mutable IsRowListMatrix, an IsRowListMatrix, and a scalar",
+  [ IsRowListMatrix and IsMutable, IsRowListMatrix, IsScalar ],
+  function( dstmat, srcmat, scalar )
+    local i;
+    if DimensionsMat(dstmat) <> DimensionsMat(srcmat) then
+      Error("AddMatrix: matrices must have the same dimensions");
+    fi;
+    for i in [1..NrRows(dstmat)] do
+      AddRowVector(dstmat[i], srcmat[i], scalar);
+    od;
+  end );
+
 InstallMethod( AddMatrix, "for a mutable 8bit matrix, an 8bit matrix, and a scalar",
   [ Is8BitMatrixRep and IsMutable, Is8BitMatrixRep, IsFFE ],
   function( dstmat, srcmat, scalar )
@@ -1985,7 +2009,7 @@ InstallMethod( AddMatrix, "for a mutable 8bit matrix, an 8bit matrix, and a scal
       Error("AddMatrix: matrices must have the same dimensions");
     fi;
     for i in [1..NrRows(dstmat)] do
-      ADD_COEFFS_VEC8BIT_3(dstmat[i], srcmat[i], scalar);
+      ADD_ROWVECTOR_VEC8BITS_3(dstmat[i], srcmat[i], scalar);
     od;
   end );
 
@@ -2017,6 +2041,15 @@ InstallEarlyMethod( MultMatrixRight,
     fi;
   end );
 
+InstallMethod( MultMatrixRight, "for a mutable IsRowListMatrix and a scalar",
+  [ IsRowListMatrix and IsMutable, IsScalar ],
+  function( mat, scalar )
+    local i;
+    for i in [1..NrRows(mat)] do
+      MultVectorRight(mat[i], scalar);
+    od;
+  end );
+
 #############################################################################
 ##
 #M  MultMatrixLeft( <mat>, <mult> )
@@ -2045,6 +2078,15 @@ InstallEarlyMethod( MultMatrixLeft,
     fi;
   end );
 
+InstallMethod( MultMatrixLeft, "for a mutable IsRowListMatrix and a scalar",
+  [ IsRowListMatrix and IsMutable, IsScalar ],
+  function( mat, scalar )
+    local i;
+    for i in [1..NrRows(mat)] do
+      MultVectorLeft(mat[i], scalar);
+    od;
+  end );
+
 ############################################################################
 ##  Fallback method for DeterminantMatrix
 

lib/matobj2.gd

@@ -2062,8 +2062,11 @@ DeclareOperationKernel( "SwapMatrixColumns", [ IsMatrixOrMatrixObj and IsMutable
 ##  <Returns>nothing</Returns>
 ##
 ##  <Description>
-##  Computes the calculation <M>M + c \cdot N</M> in-place, storing the result in <A>M</A>.
+##  Computes the calculation <M>M + N \cdot c</M> in-place, storing the result in <A>M</A>.
 ##  If the optional argument <A>c</A> is omitted, then <A>N</A> is added directly.
+##  The matrices must have the same dimensions, otherwise the result is undefined.
+##  Specialized methods may be defined only when <A>M</A> and <A>N</A> have the same
+##  representation.
 ##  If both of the matrices are lists-of-lists, then the operation is delegated
 ##  row by row to <Ref Oper="AddRowVector"/>.
 ##  <Example><![CDATA[

lib/meatauto.gi

@@ -132,10 +132,12 @@ local n, weights, mat, vec, ReduceRow, t,
     if lead > n then
       return rhs;
     fi;
+    lhs := ShallowCopy(lhs);
     for i in [1..Length(weights)] do
       if weights[i] = lead then
         z := lhs[lead];
-        lhs := lhs - z * mat[i]; rhs := rhs - z * vec[i];
+        AddVector(lhs, mat[i], -z);  # lhs is a vector
+        rhs := rhs - z * vec[i];     # rhs is a scalar
         lead := PositionNonZero(lhs, lead);
         if lead > n then
           return rhs;
@@ -217,10 +219,10 @@ local m, n, zero, i, c, j, factor;
     for i in [1..m] do
       c := eqns.weights[i];
       for j in [1..i-1] do
-        if eqns.mat[j][c] <> zero then
+        factor := eqns.mat[j,c];
+        if factor <> zero then
           Info(InfoMtxHom,6,"solveEqns: kill mat[",j,",",c,"]");
-          factor := eqns.mat[j][c];
-          eqns.mat[j] := eqns.mat[j] - factor*eqns.mat[i];
+          AddVector(eqns.mat[j], eqns.mat[i], -factor);
           eqns.vec[j] := eqns.vec[j] - factor*eqns.vec[i];
         fi;
       od;
@@ -306,13 +308,13 @@ local n, coeffs, x, zero, z, i;
     coeffs:=[];
     x:=v;
   else
-    x:=v;
+    x:=ShallowCopy(v);
     zero:=x[1]*0;
     coeffs:=ListWithIdenticalEntries(n, zero);
     for i in [1..n] do
       z:=x[ech[i]];
       if z <> zero then
-        x:=x - z * base[i];
+        AddVector(x, base[i], -z);
         coeffs[i]:=z;
       fi;
     od;
@@ -552,13 +554,13 @@ end);
 # If a linearly dependent set of elements is supplied, this
 # routine will trim it down to a basis.
 BindGlobal("SMTX_EcheloniseMats",function (gens, F)
-local n, m, zero, ech, k, i, j, l;
+local n, m, zero, ech, k, i, j, l, entry;
 
   if Length(gens) = 0 then
     return [ [], [] ];
   fi;
   # copy the list to avoid destroying the original list
-  gens:=ShallowCopy(gens);
+  gens:=List(gens, MutableCopyMatrix);
 
   n:=NrRows(gens[1]);
   m:=NrCols(gens[1]);
@@ -581,12 +583,15 @@ local n, m, zero, ech, k, i, j, l;
       Add(ech, [i,j]);
 
       # First normalise the [i,j] position to 1
-      gens[k]:=gens[k] / gens[k][i,j];
+      entry := gens[k][i,j];
+      if not IsOne(entry) then
+        MultMatrix(gens[k], 1/entry);
+      fi;
 
       # Now zero position [i,j] in all further generators
       for l in [k+1..Length(gens)] do
         if (gens[l][i,j] <> zero) then
-          gens[l]:=gens[l] - gens[k] * gens[l][i,j];
+          AddMatrix(gens[l], gens[k], -gens[l][i,j]);
         fi;
       od;
       k:=k + 1;
@@ -617,7 +622,7 @@ end);
 # The code is heavily commented, and I appreciate suggestions on how to
 # improve it (particularly bits of code).
 BindGlobal("SpinHom",function (V, W)
-local nv, nw, F, zero, zeroW, gV, gW, k, U, echu, r, homs, s, work, ans, v0,
+local nv, nw, F, zero, minusone, zeroW, gV, gW, k, U, echu, r, homs, s, work, ans, v0,
       M, x, pos, z, echm, t, v, echv, a, u, e, start, oldlen, ag, m, uu, ret,
       c, s1, X, mat, uuc, uic, newhoms, hom, Uhom, imv0, imv0c, image, i, j, l;
 
@@ -635,6 +640,7 @@ local nv, nw, F, zero, zeroW, gV, gW, k, U, echu, r, homs, s, work, ans, v0,
   TestModulesFitTogether(V,W);
   F:=V.field;
   zero:=Zero(F);
+  minusone:=-One(F);
 
   zeroW:=ListWithIdenticalEntries(nw,zero);
   zeroW:=ImmutableVector(F,zeroW);
@@ -744,27 +750,26 @@ local nv, nw, F, zero, zeroW, gV, gW, k, U, echu, r, homs, s, work, ans, v0,
           # create new element <x>, with its definition as the
           # difference between <v0^m> and <uu> in <U>.
           x:=v[i] * gV[j];
-          m:=ag;
+          m:=MutableCopyMat(ag);
           ConvertToMatrixRep(m, F);
           uu:=u[i] * gV[j];
 
           ret:=EchResidueCoeffs(U, echu, x,3);
           x:=ret.residue;
-          uu:=uu - ret.projection;
-
+          AddVector(uu, ret.projection, minusone);
           # reduce modulo the new semi-ech basis elements in <v>,
           # storing the coefficients in <c>
           #
           c:=ListWithIdenticalEntries(Length(v),zero);
           for l in [1..Length(v)] do
             z:=x[echv[l]];
             if z <> zero then
-              x:=x - z * v[l];
+              AddVector(x, v[l], -z);
               if Length(m) > 0 then
-                  m:=m - z * a[l];
+                  AddMatrix(m, a[l], -z);
               fi;
               c[l]:=c[l] + z;
-              uu:=uu - z * u[l];
+              AddVector(uu, u[l], -z);
             fi;
           od;
           c:=ImmutableVector(F,c);
@@ -797,14 +802,14 @@ local nv, nw, F, zero, zeroW, gV, gW, k, U, echu, r, homs, s, work, ans, v0,
             for l in [1..Length(v)] do
               if c[l] <> zero then
                 if Length(X) > 0 then
-                  X:=X + c[l] * a[l];
+                  AddMatrix(X, a[l], c[l]);
                 fi;
-                uu:=uu + c[l] * u[l];
+                AddVector(uu, u[l], c[l]);
               fi;
             od;
 
             if Length(X) > 0 then
-              X:=X - ag;
+              AddMatrix(X, ag, minusone);
             fi;
 
             mat:=[];
@@ -814,6 +819,7 @@ local nv, nw, F, zero, zeroW, gV, gW, k, U, echu, r, homs, s, work, ans, v0,
               Add(mat, uuc * homs[l] - uic * homs[l] * gW[j]);
             od;
             Append(mat, X);
+            ConvertToMatrixRep(mat);
             SMTX_AddEqns(e, TransposedMat(mat), zeroW);
           fi;
         od;
@@ -848,26 +854,26 @@ local nv, nw, F, zero, zeroW, gV, gW, k, U, echu, r, homs, s, work, ans, v0,
         ConvertToMatrixRep(Uhom, F);
         for l in [1..s] do
           if ans[i][l] <> zero then
-            Uhom:=Uhom + ans[i][l] * homs[l];
+            AddMatrix(Uhom, homs[l], ans[i][l]);
           fi;
         od;
         for l in [1..r] do
           Add(hom, Uhom[l]);
         od;
       fi;
 
-      imv0:=zeroW * zero;
+      imv0:=ZeroMutable(zeroW);
       for l in [1..t] do
         if ans[i][s+l] <> zero then
-          imv0:=imv0 + ans[i][s+l] * M[l];
+          AddVector(imv0, M[l], ans[i][s+l]);
         fi;
       od;
       imv0c:=EchResidueCoeffs(M, echm, imv0,1);
       for l in [1..Length(v)] do
         if Length(imv0c)=0 then image:=[];
         else image:=imv0c * a[l];fi;
         if r > 0 then
-          image:=image + EchResidueCoeffs(U, echu, u[l],1) * Uhom;
+          AddVector(image, EchResidueCoeffs(U, echu, u[l],1) * Uhom);
         fi;
         Add(hom, image);
       od;
@@ -984,18 +990,21 @@ local proveIndecomposability, addnilpotent, n, F, zero, basis, enddim,
   local i, r, c, k, done, l;
   # NB: <remain> and <nildim> and <cnt> are not local
 
+    a := MutableCopyMatrix(a);
     for i in [1..nildim] do
-      r:=echelon[nilech[i]][1]; c:=echelon[nilech[i]][2];
+      r:=echelon[nilech[i]][1];
+      c:=echelon[nilech[i]][2];
       if a[r,c] <> zero then
-        a:=a - a[r,c] * nilbase[i] / nilbase[i][r,c];
+        AddMatrix(a, nilbase[i], -a[r,c] / nilbase[i][r,c]);
       fi;
     od;
 
     # find which echelon index to remove due to this new element
     k:=1; done:=false;
     while not done and k <= Length(remain) do
       l:=remain[k];
-      r:=echelon[l][1]; c:=echelon[l][2];
+      r:=echelon[l][1];
+      c:=echelon[l][2];
       if a[r,c] <> zero then
         done:=true;
       else
@@ -1377,12 +1386,12 @@ SMTX.IsomorphismModules:=SMTX_IsomorphismModules;
 # running down diagonals below the main diagonal:
 #   [2,1], [3,2], [4,3], ..., [3,1], [4,2], ..., [n-1,1], [n, 2], [n,1]
 BindGlobal("SMTX_EcheloniseNilpotentMatAlg",function (matalg, F)
-local zero, n, base, ech, k, diff, i, j, found, l;
+local zero, n, base, ech, k, diff, i, j, found, l, entry;
 
   zero:=Zero(F);
   n := NrCols(matalg[1]);
 
-  base := ShallowCopy(matalg);
+  base := List(matalg, MutableCopyMatrix);
   ech := [];
   k := 1;
 
@@ -1409,12 +1418,15 @@ local zero, n, base, ech, k, diff, i, j, found, l;
       Add(ech, [i,j]);
 
       # First normalise the [i,j] position to 1
-      base[k] := base[k] / base[k][i,j];
+      entry := base[k][i,j];
+      if not IsOne(entry) then
+        MultMatrix(base[k], 1/entry);
+      fi;
 
       # Now zero position [i,j] in all other basis elements
       for l in [1..Length(base)] do
         if (l <> k) and (base[l][i,j] <> zero) then
-          base[l] := base[l] - base[k] * base[l][i,j];
+          AddMatrix(base[l], base[k], -base[l][i,j]);
         fi;
       od;
       k := k + 1;

tst/testinstall/meataxe.tst

@@ -1,4 +1,4 @@
-#@local G,M,M2,M3,M4,M5,M6,V,bf,bo,cf,homs,m,mat,qf,randM,res,sf,subs,mats,Q,orig,S
+#@local G,M,M2,M3,M4,M5,M6,V,bf,bo,cf,homs,m,mat,qf,randM,res,sf,subs,mats,Q,orig,S,dim,F,i
 gap> START_TEST("meataxe.tst");
 
 #
@@ -71,16 +71,12 @@ Error, generators are different lengths
 #
 #
 #
+gap> G:=SymmetricGroup(3);;
+gap> M:=PermutationGModule(G,GF(2));;
 gap> M2:=TensorProductGModule(M,M);
 rec( IsOverFiniteField := true, dimension := 9, field := GF(2), 
   generators := [ <an immutable 9x9 matrix over GF2>, 
       <an immutable 9x9 matrix over GF2> ], isMTXModule := true )
-gap> M6:=DirectSumGModule(M,M);
-rec( IsOverFiniteField := true, dimension := 6, field := GF(2), 
-  generators := [ <an immutable 6x6 matrix over GF2>, 
-      <an immutable 6x6 matrix over GF2> ], isMTXModule := true )
-gap> M6.generators[1] = DirectSumMat(M.generators[1], M.generators[1]);
-true
 gap> IdGroup(MTX.ModuleAutomorphisms(M2));
 [ 1344, 11301 ]
 gap> cf:=MTX.CompositionFactors(M2);;
@@ -91,6 +87,29 @@ gap> List(Filtered(cf, x -> x.dimension=2), MTX.InvariantQuadraticForm);
 [ <an immutable 2x2 matrix over GF2>, <an immutable 2x2 matrix over GF2>, 
   <an immutable 2x2 matrix over GF2> ]
 
+#
+#
+#
+gap> M6:=DirectSumGModule(M,M);
+rec( IsOverFiniteField := true, dimension := 6, field := GF(2), 
+  generators := [ <an immutable 6x6 matrix over GF2>, 
+      <an immutable 6x6 matrix over GF2> ], isMTXModule := true )
+gap> M6.generators[1] = DirectSumMat(M.generators[1], M.generators[1]);
+true
+
+#
+#
+#
+gap> dim:=10;; F:=GF(25);;
+gap> G:=Sp(dim,F);;
+gap> M:=NaturalGModule(G);;
+gap> for i in [1..3] do
+>      M:=DirectSumGModule(M, NaturalGModule(G^RandomInvertibleMat(dim, F)));
+>    od;
+gap> res:=MTX.Indecomposition(M);;
+gap> Length(res);
+4
+
 #
 #
 #