Overhaul IsPlistMatrixRep and IsPlistVectorRep

[fingolfin]

Progress towards issue #2148 : This PR is an incomplete proof of concept (see #2776 for its genesis), I am putting it out here to make sure people are aware of this possible direction, and in the vague hope that somebody else might help out with it.

The idea is to tune IsPlistMatrixRep as much as possible so that it becomes a viable replacement for lists-of-lists (= lol) "matrices". Indeed, with this hackish POC, I already benchmarked some cases where the IsPlistMatrixRep code is faster than lists-of-lists, probably in parts because of avoided method selection overhead; another part is that the optimized code benefits from kernel optimizations for multiplication of lol matrices

To give an idea about the performance of IsPlistMatrixRep before and after this PR, I run some mini benchmarks, using the benchmarking code in this gist, plus the following code:

M1:=IdentityMatrix(GF(4),10);
M2:=IdentityMat(10)*Z(2); # "decompress" the matrix
M3:=IdentityMatrix(IsPlistMatrixRep, GF(4),10);
I1:=IdentityMat(10);
I2:=IdentityMatrix(IsPlistMatrixRep,Integers,10);

inputs := [
   [ M1, GF(4), Is8BitMatrixRep ],
   [ M2, GF(4), IsPlistRep ],
   [ M3, GF(4), IsPlistMatrixRep ],
   [ I1, Integers, IsPlistRep ],
   [ I2, Integers, IsPlistMatrixRep ],
];

f := function(mat) local i, x; for i in [1..10000] do x:=mat*mat; od; end;;
g := function(mat) local i, x; for i in [1..10000] do x:=mat^5+mat; od; end;;
opts := rec(mintime:=1000);

for test in [ f, g ] do
    Print("-----\n");
    Print("Testing function ", NameFunction(test), "\n");
    Print("-----\n");

    for input in inputs do
        Print("... on ", NameFunction(input[3]), " matrix over ", input[2], "...\n");
        Assert(0, input[3](input[1])); # verify filter matches
        Benchmark(function() test(input[1]); end, opts);;
        Print("\n");
    od;

    Print("\n");
od;

The results of this are as follows with master (using commit 5b7263cc69a93f0dc9a067d5df90c4e756d46bbb from 2019-07-11), on a Ubuntu server in Gießen:

-----
Testing function f
-----
... on Is8BitMatrixRep matrix over GF(2^2)...
Performed 128 iterations, taking 1008. milliseconds.
average: 7.87 +/- 0.02 (+/- 0.%)
median: 7.870849609375

... on IsPlistRep matrix over GF(2^2)...
Performed 5 iterations, taking 1517. milliseconds.
average: 303.45 +/- 0.12 (+/- 0.%)
median: 303.48095703125

... on IsPlistMatrixRep matrix over GF(2^2)...
Performed 5 iterations, taking 5801. milliseconds.
average: 1160.23 +/- 0.89 (+/- 0.%)
median: 1159.994140625

... on IsPlistRep matrix over Integers...
Performed 18 iterations, taking 1001. milliseconds.
average: 55.62 +/- 0.07000000000000001 (+/- 0.%)
median: 55.6290283203125

... on IsPlistMatrixRep matrix over Integers...
Performed 5 iterations, taking 1900. milliseconds.
average: 380.09 +/- 0.33 (+/- 0.%)
median: 380.011962890625


-----
Testing function g
-----
... on Is8BitMatrixRep matrix over GF(2^2)...
Performed 29 iterations, taking 1010. milliseconds.
average: 34.84 +/- 0.02 (+/- 0.%)
median: 34.8369140625

... on IsPlistRep matrix over GF(2^2)...
Performed 5 iterations, taking 4797. milliseconds.
average: 959.4400000000001 +/- 0.21 (+/- 0.%)
median: 959.5048828125

... on IsPlistMatrixRep matrix over GF(2^2)...
Performed 5 iterations, taking 15796. milliseconds.
average: 3159.28 +/- 1.58 (+/- 0.%)
median: 3159.17724609375

... on IsPlistRep matrix over Integers...
Performed 5 iterations, taking 1014. milliseconds.
average: 202.89 +/- 0.09 (+/- 0.%)
median: 202.889892578125

... on IsPlistMatrixRep matrix over Integers...
Performed 5 iterations, taking 6131. milliseconds.
average: 1226.19 +/- 0.49 (+/- 0.%)
median: 1226.10791015625

With this PR (commit cb79cdf58b85ae83b3fbf218332e36a16970cdfb), the numbers for Is8BitMatrixRep and IsPlistRep stay essentially the same, but for IsPlistMatrixRep improve to more or less match those for IsPlistRep (they are a bit slower, but not that much):

-----
Testing function f
-----
... on Is8BitMatrixRep matrix over GF(2^2)...
Performed 127 iterations, taking 1001. milliseconds.
average: 7.88 +/- 0.06 (+/- 1.%)
median: 7.870849609375

... on IsPlistRep matrix over GF(2^2)...
Performed 5 iterations, taking 1513. milliseconds.
average: 302.69 +/- 0.09 (+/- 0.%)
median: 302.649169921875

... on IsPlistMatrixRep matrix over GF(2^2)...
Performed 5 iterations, taking 1582. milliseconds.
average: 316.49 +/- 0.06 (+/- 0.%)
median: 316.454833984375

... on IsPlistRep matrix over Integers...
Performed 18 iterations, taking 1004. milliseconds.
average: 55.77 +/- 0.06 (+/- 0.%)
median: 55.7733154296875

... on IsPlistMatrixRep matrix over Integers...
Performed 15 iterations, taking 1008. milliseconds.
average: 67.19 +/- 0.16 (+/- 0.%)
median: 67.14599609375


-----
Testing function g
-----
... on Is8BitMatrixRep matrix over GF(2^2)...
Performed 29 iterations, taking 1003. milliseconds.
average: 34.6 +/- 0.03 (+/- 0.%)
median: 34.587158203125

... on IsPlistRep matrix over GF(2^2)...
Performed 5 iterations, taking 4775. milliseconds.
average: 955. +/- 0.5600000000000001 (+/- 0.%)
median: 955.198974609375

... on IsPlistMatrixRep matrix over GF(2^2)...
Performed 5 iterations, taking 4858. milliseconds.
average: 971.5700000000001 +/- 0.7000000000000001 (+/- 0.%)
median: 971.367919921875

... on IsPlistRep matrix over Integers...
Performed 5 iterations, taking 1015. milliseconds.
average: 202.91 +/- 0.36 (+/- 0.%)
median: 202.7880859375

... on IsPlistMatrixRep matrix over Integers...
Performed 5 iterations, taking 1140. milliseconds.
average: 228.01 +/- 0.09 (+/- 0.%)
median: 227.99169921875

Of course proper benchmarks would include a range of matrix sizes and ground fields; more interesting / realistic operations; and graph all of that nicely...

[codecov[bot]]

Codecov Report

Merging #2973 into master will decrease coverage by 12.35%. The diff coverage is 48.76%.

@@             Coverage Diff             @@
##           master    #2973       +/-   ##
===========================================
- Coverage   84.65%   72.29%   -12.36%     
===========================================
  Files         698      635       -63     
  Lines      345626   308665    -36961     
===========================================
- Hits       292582   223150    -69432     
- Misses      53044    85515    +32471

Impacted Files	Coverage Δ
lib/grpmat.gi	66.56% <ø> (-7.44%)	:arrow_down:
lib/matobj2.gd	100% <ø> (ø)	:arrow_up:
lib/matobjplist.gd	100% <100%> (ø)	:arrow_up:
lib/matrix.gi	74.76% <100%> (-12.98%)	:arrow_down:
lib/vecmat.gi	83.2% <100%> (-3.56%)	:arrow_down:
lib/vecmat.gd	100% <100%> (ø)	:arrow_up:
lib/grpffmat.gi	91.88% <100%> (-4.04%)	:arrow_down:
lib/matobjplist.gi	56.17% <35.63%> (+5.69%)	:arrow_up:
lib/matobj.gi	65.66% <53.84%> (-0.54%)	:arrow_down:
lib/teachm2.g	15.38% <0%> (-84.62%)	:arrow_down:
... and 395 more

[coveralls]

Coverage increased (+0.03%) to 84.699% when pulling 0e7599b3cd18a3317be9dbe0ba29f33dba26f1ea on fingolfin:mh/IsPlistMatrixRep into d2db1bf548a0efaf27cfa497f53bae10465d11f0 on gap-system:master.

[fingolfin]

I think this is getting close to being mergeable. I am sure there are bugs in my changes, but to flush those out, we should write more and better systematic tests for the interfaces.

[fingolfin]

@ssiccha @ThomasBreuer thank you for your feedback

@ThomasBreuer you wrote:

The documentation of IsPlistMatrixRep states that a plain list of IsPlistVectorRep objects is stored, should this be changed?

Yes, thank you. In fact, there is more in it that needs to be changed: it implies that IsPlistMatrixRep is a row list matrix, which was the case before, but now isn't / shouldn't be... I'll try to explain my train of thought here (and this definitely should be explained in the commit message(s), too): for optimal performance (and I run a lot of benchmarks back then), we really don't want to store a IsPlistMatrixRep as a list of IsPlistVectorRep; this just adds a ton of overhead. It's much better to simply wrap an old-style IsMatrix and IsPlistRep matrix, which this PR does. Also, for testing purposes, it is good to have an IsMatrixObj implementation which is not an IsRowListMatrix (@danielrademacher asked for this, too). So that's a useful application for IsPlistMatrixRep.

HOWEVER another useful application for IsPlistMatrixRep would be to use it as drop-in replacement for any old-style IsMatrix and IsPlistRep; i.e., one may want to change code producing these to instead produce IsPlistMatrixRep (not all at once, but one at a time, if/when it seems useful and safe). But for this application, one wants those matrices to work as closely as possible to the old-style matrices; so in particular, here we want an IsRowListMatrix.

So here we have a clash of goal. So my thought here now is to introduce IsRowPlistMatrixRep (or something like that) which specializes IsPlistMatrixRep and is also in IsRowListMatrix. I think all that is needed to make that happen is to install suitable [] and []:= methods for IsRowPlistMatrixRep (and perhaps a few more, dunno). Then mat[i] could simply return an IsVectorObj which wraps the list mat![ROWSPOS][i] (and not a shallow copy; it really must point to exactly that list). This way, mat[i][j] := 1 will still work. Of course this incurs a performance penalty; I think that at least the [] method should contain a Info(InfoPerformance, ...) message. If we want to optimize it a bit, we could add a WeakPointerObj as cache for these row wrapper objects (some care needs to be done: replacing a row via []:= should invalidate the respective cache entry).

[fingolfin]

Err, I've started to add IsRowPlistMatrixRep, but completing it will be tedious, as I need to once again implement NewZeroMatrix, NewIdentityMatrix, ... -- we ought to resolve #4366