Use MeatAxe as an optional back end for dense matrices over `GF(p^n)`, p odd, n>1, `p^n<255`

[SimonKing]

Sage has (or will soon have) fairly good implementations of dense matrices over GF(2), over GF(2^e) (#9562) and over GF(p) (p prime, #4260). However, it uses generic code for dense matrices over GF(p^n), p odd, n>1, p^n<255.

The original suggestion was to use a major modification of MeatAxe Release 2.2.4 instead of the basic implementation. The timings below are with that old version (it is identical with 2.2.3 except for the GPL licence, and 2.2.3 was before 1998).

I now suggest to try and do the same with the latest MeatAxe release 2.4.24, which is from 2011. There also is an experimental 2.5.0 from 2003, but I think we shouldn't rely on that.

Sources

The upstream sources ​http://www.math.rwth-aachen.de/~MTX/meataxe-2.4.24.tar.gz needed to be repackaged, in order to unpack into a single folder. Use attachment:meataxe-2.4.24.tar.gz​

What is done

There is no spkg-check. However, if SAGE_CHECK=yes or of one does sage -i -c meataxe, then a test suite is executed as part of building the package.

It is my experience that the tests pass most of the time. I can not explain why sometimes they don't.

What is missing

Currently, the spkg installs libmtx.a and installs some binaries. However, I also intend to add a Cython wrapper so that one can use MeatAxe matrices in Sage.

Here is the original ticket description:

This is awfully slow:

sage: MS = MatrixSpace(GF(5^3,'y'),2000)
sage: %time A = MS.random_element()
CPU times: user 6.36 s, sys: 0.02 s, total: 6.39 s
Wall time: 6.41 s
sage: type(A)
<type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'>
sage: B = MS.random_element()
sage: %time A*B    # using 6.3% of my computer's memory
CPU times: user 744.20 s, sys: 1.18 s, total: 745.38 s
Wall time: 747.69 s
2000 x 2000 dense matrix over Finite Field in y of size 5^3
sage: %time ~A     # using 10.4% of my computer's memory
CPU times: user 1096.74 s, sys: 1.30 s, total: 1098.05 s
Wall time: 1101.24 s
2000 x 2000 dense matrix over Finite Field in y of size 5^3
sage: %time A.echelon_form() # using 10.4% of my computer's memory
CPU times: user 378.62 s, sys: 0.33 s, total: 378.95 s
Wall time: 380.06 s
2000 x 2000 dense matrix over Finite Field in y of size 5^3

With the optional spkg and the patch, one gets a clear improvement.

sage: MS = MatrixSpace(GF(5^3,'y'),2000)
sage: %time A = MS.random_element()
CPU times: user 0.32 s, sys: 0.00 s, total: 0.32 s
Wall time: 0.33 s
sage: type(A)
<type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
sage: B = MS.random_element()
# The following uses Strassen-Winograd multiplication
sage: %time A*B      # using 3.5% of my computer's memory
CPU times: user 7.68 s, sys: 0.01 s, total: 7.69 s
Wall time: 7.72 s
2000 x 2000 dense matrix over Finite Field in y of size 5^3
# The following is school book multiplication;
# that's more or less the original meataxe speed:
sage: %time A._multiply_classical(B)   # using 3.6% of my computer's memory
CPU times: user 11.68 s, sys: 0.02 s, total: 11.70 s
Wall time: 11.73 s
2000 x 2000 dense matrix over Finite Field in y of size 5^3
# Strassen is not implemented for inversion and echelon form.
sage: %time ~A    # using 3.8% of my computer's memory
CPU times: user 23.55 s, sys: 0.00 s, total: 23.55 s
Wall time: 23.62 s
2000 x 2000 dense matrix over Finite Field in y of size 5^3
sage: %time A.echelon_form() #using 3.9% of my computer's memory
CPU times: user 11.73 s, sys: 0.01 s, total: 11.74 s
Wall time: 11.78 s
2000 x 2000 dense matrix over Finite Field in y of size 5^3

[SimonKing]

Here is more on libmeataxe-1.0.spkg.

I started with Release 2.2.4 of the Aachen C MeatAxe. This is partially because I first came in touch with MeatAxe by old code of David Green, who was using the 2.2.3 release. MeatAxe 2.2.4 was released under GPL2+, which should be fine. There are more recent releases. However, according to Michael Ringe, while the interface changed substantially, the linear algebra is more or less the same.

The spkg contains the unmodified sources in the folder src/, which is not under version control. The file patches/libmeataxe.patch provides my modification. Of course, if one installs the package, the patch will automatically be applied. Here is a summary of my changes:

• Normally, MeatAxe provides some executables that operate on matrices stored in files. I strip it down, such that the executables are not built and only a C library remains (this is ​libmeataxe-1.0.spkg). I wrap the C library using Cython (this is trac12103_meataxe.patch​)
• MeatAxe uses school book multiplication. I added Strassen-Winograd multiplication, using a schedule that minimizes memory consumption. See src/src/window.c.
• The environment variable MTXLIB was supposed to tell MeatAxe where to find multiplication tables. However, this never worked for me. I fix it, so that now multiplication tables in $SAGE_ROOT/local are used, that are created when the package is installed.
• A simpler version of the package has been part of our group cohomology spkg. Bad experiences on Solaris made me change three names: I changed
  • T_STRING into T_STRINGS,
  • matextract into _matextract,
  • matid into MTXmatid.

Please don't ask me why the Solaris linker was mistaking these symbols with symbols in completely different Sage packages, but as a matter of fact it did.

The aim of this ticket is to make the modified MeatAxe an optional back end for dense matrices over GF(p^n), p>2 prime, n>1, p^n<255. Currently, one needs to install an spkg and a patch to the Sage library.

I would prefer to have it all in one, such that sage -i libmeataxe-1.0.spkg would automatically patch the Sage library. Question to Sage experts: Is that possible?

Anyway, I think it can be reviewed...

[SimonKing]

Use MeatAxe as back end for linear algebra over some finite fields.

[SimonKing]

Use MeatAxe as back end for linear algebra over some finite fields. Relative to #11900

[SimonKing]

I have updated both patches, because I wanted to add one method: matrix_from_rows. One could rely on the generic implementation, but in some application I need this to be fast:

            sage: MS = MatrixSpace(GF(5^3,'y'),3,2)
            sage: A = MS(range(6))
            sage: A
            [0 1]
            [2 3]
            [4 0]
            sage: A.matrix_from_rows([2,1])
            [4 0]
            [2 3]
            sage: A.matrix_from_rows([1,1])
            [2 3]
            [2 3]

[SimonKing]

Fix one doctest (coping with a changed class name)

[SimonKing]

Apparently there has been a name change MatrixSpace_generic to MatrixSpace. The result is one doctest failure in my original patch. I am not sure where that name change happened (perhaps it is in one of my patches that aren't in vanilla Sage yet). Therefore I added a second patch, that may or may not be necessary.

Apply trac12103_meataxe_rel11900.patch trac12103_meataxe_docfix.patch

[malb]

The patch applies neither cleanly to 4.8 nor 5.0.beta6. The issue 4.8 is small and easy to fix, 5.0.beta6 are a bit more failures.

[malb]

as of now module_list.py will try to build matrix_modpn_dense regardless of whether MeatAxe? is installed, but it should only build it if MeatAxe? is installed.

[SimonKing]

Replying to malb:

The patch applies neither cleanly to 4.8 nor 5.0.beta6. The issue 4.8 is small and easy to fix, 5.0.beta6 are a bit more failures.

I get

(sage subshell) mpc721:sage king$ hg qpush
applying trac12103_meataxe_rel11900.patch
patching file sage/modules/quotient_module.py
Hunk #12 succeeded at 297 with fuzz 1 (offset 10 lines).
Hunk #13 succeeded at 322 with fuzz 1 (offset 14 lines).
now at: trac12103_meataxe_rel11900.patch
SAGE_ROOT=/home/king/SAGE/sage-5.0.beta6

So, I can not confirm that the first to-be-applied patch is really problematic (ok, fuzz 1 might be fixed at some point, but I don't get actual errors).

How can I test in module_list.py whether or not MeatAxe is installed? I suppose the patch to the sage library should not be applied by the installation script of the spkg.

[SimonKing]

PS: Martin, did you really try to apply trac12103_meataxe_rel11900.patch? Namely, the previous patch (that is not relative to #11900) should be disregarded, as #11900 is already merged.

[malb]

Ah, crap, 11900 was it.

RE: module_list.py I commented on it here: ​https://bitbucket.org/malb/sagemath/pull-request/1/12103-meataxe-wrapper#comment-3510

PS: Did you get my e-mail that I sent to simon.king_uni-jena.de?

[SimonKing]

Replying to malb:

Ah, crap, 11900 was it.

OK, I changed the instructions in the ticket description accordingly.

RE: module_list.py I commented on it here: ​https://bitbucket.org/malb/sagemath/pull-request/1/12103-meataxe-wrapper#comment-3510

PS: Did you get my e-mail that I sent to simon.king_uni-jena.de?

Yes, but I didn't look at bitbucket yet, and I suppose I will not find the time to start and learn the new development before Monday.

[SimonKing]

I think this ticket should be revived. But the aim should be to have a library version of the *latest* MeatAxe upstream sources. And of course it should be a new-style package.

[SimonKing]

The first step will probably be the creation of a new-style package for MeatAxe 2.4.24 (which I take as an exercise for learning new-style packaging).

The README file says:

================================================================================
INSTALLATION (for Unix systems)
================================================================================

1) Copy Makefile.conf.dist to Makefile.conf. Edit Makefile.conf
   and fill in the necessary parameters.

2) To compile all programs an run the test suite, enter

        make check

   You will probably need gmake (GNU make) to build the MeatAxe.

3) If there were no errors, run

        make install

   to install the executables in their final location (as selected
   in the makefile).

I wonder if we should follow that procedure. After all, Makefile.conf.dist just defines a couple of variables that determine the compiler, the flags, the location of the binaries, the location of the library and the implementation (for fields of order up to 255 or for larger fields).

What do people think?

1. a) Should we patch Makefile, by removing the line include Makefile.conf, and define all the variables MTXLIB, MTXBIN etc. from Makefile.conf as environment variables?

or

1. b) Should we keep Makefile, create an empty Makefile.conf, and do the rest with environment variables?

or

2. Should spkg-install first create a fully grown Makefile.conf?

[vbraun]

If all you have to do is touch Makefile.conf then I'd go for that.

[SimonKing]

I am slightly puzzled.

There is a variable ZZZ used in Makefile. But when I do

ZZZ=0
$MAKE

then apparently Makefile doesn't know about ZZZ. What do I need to do to define ZZZ so that Makefile knows about it?

[vbraun]

Unless you export ZZZ it is not passed to subprocesses. Alternatively, you might want to use $MAKE ZZZ=0.

[SimonKing]

Thanks, "export" did the trick.

However, I am afraid that some of the reasons for massively patching MeatAxe 2.2.4 are still present in version 2.4.24:

• The environment variables MTXLIB and MTXBIN are ignored during build. I.e., the library is put into ./tmp/, the binaries are put into ./bin/

Suggested solution: Let spkg-install move the created binaries and libmtx.a to the desired location

• make clean doesn't work, as it also goes into the directory test/ and calls make clean there---but there is no test/Makefile.

Not relevant here, as the build directory is deleted anyway.

• make check rebuilds everything in ./tmp/ and ./bin/ before testing. Hence, the test suite actually does NOT test what is in $SAGE_ROOT/local/lib and $SAGE_ROOT/local/bin.

Question on make check: Would it be OK to let spkg-install ALWAYS do make check? Or should that be solely the job of spkg-check (with the drawback that then spkg-check would rebuild meataxe)?

[SimonKing]

Replying to SimonKing:

Question on make check: Would it be OK to let spkg-install ALWAYS do make check? Or should that be solely the job of spkg-check (with the drawback that then spkg-check would rebuild meataxe)?

Or should spkg-install test whether SAGE_CHECK=yes, and choose between make and make check accordingly?

[SimonKing]

Do I understand correctly that for automatically downloading the sources from upstream, I should create spkg-src similarly to what I find in build/pkg/singular/spkg-src?

[SimonKing]

Replying to SimonKing:

However, I am afraid that some of the reasons for massively patching MeatAxe 2.2.4 are still present in version 2.4.24:

But one aspect has improved: MeatAxe used to rely on certain tables that were put (and when needed also created) in the current directory. And I had to patch the c-sources so that MeatAxe would be able to find the tables in a fixed directory. Now I cannot see any of these tables.

[jhpalmieri]

Replying to SimonKing:

Question on make check: Would it be OK to let spkg-install ALWAYS do make check? Or should that be solely the job of spkg-check (with the drawback that then spkg-check would rebuild meataxe)?

In my opinion, if make check doesn't add too much time (a minute?) to the build, then you can always do it, but if it's more than that, it should be done only in spkg-check.

[vbraun]

You don't need spkg-src, its just a bandaid to modify tarballs (don't do it if you can avoid it)

[SimonKing]

New commits:

​e732f0e

An optional MeatAxe package

[SimonKing]

Question on how to continue:

My aim is to provide a Cython wrapper, so, I need to copy meataxe.h to some include location.

But how to deal with the Cython wrapper? Apparently it would be in the Sage library, but apparently it will not compile if meataxe is not installed.

So, should the package install the cython code into the Sage library? Or is it possible to have code permanently in the Sage library that is only compiled if a certain package is installed, and otherwise is left untouched?

[SimonKing]

PS: Note that MeatAxe has a nice documentation, and as usual it is installed if SAGE_SPKG_INSTALL_DOCS=yes.

[SimonKing]

Replying to SimonKing:

But one aspect has improved: MeatAxe used to rely on certain tables that were put (and when needed also created) in the current directory. And I had to patch the c-sources so that MeatAxe would be able to find the tables in a fixed directory. Now I cannot see any of these tables.

I am afraid I stand corrected. Meanwhile I have a Cython wrapper on my laptop, and when I create a matrix over GF(2), then a file p002.zzz (defining multiplication table) is created in the current directory.

That's a very nasty property of MeatAxe. What shall we do about it?

Shall we leave it like that?

Shall I try to port my patch from 2.2.4 to 2.4.24 and ensure that all multiplication tables should be created in SAGE_ROOT/local/lib/meataxe/?

[SimonKing]

Replying to SimonKing:

Shall I try to port my patch from 2.2.4 to 2.4.24 and ensure that all multiplication tables should be created in SAGE_ROOT/local/lib/meataxe/?

If "yes": Where shall the tables be placed? In SAGE_LOCAL/share/meataxe/ perhaps?

[SimonKing]

From the MeatAxe sources:

/**
 ** MeatAxe Library Directory.
 ** This variable contains the name of the MeatAxe library directory.
 ** Arithmetic table files are searched in this directory. The value of 
 ** MtxLibDir can be set on the command line with the "-L" option. Otherwise, 
 ** the value of the environment variable MTXLIB is used. If neither "-L" nor 
 ** MTXLIB are defined, the default directory, which was selected when 
 ** building the MeatAxe, is used.
 ** @see  MtxBinDir 
 **/

So, apparently the multiplication tables are supposed to be in MTXLIB. But that has NEVER worked for me.

Or rather, it could be that MTXLIB has to be an environment variable at running time. Would it be acceptable that a module sage.matrix.meataxe sets an environment variable for MeatAxe? I'll test if that works.

[SimonKing]

PS: It says "Otherwise, the value of the environment variable MTXLIB is used." That definitely does not work.

[SimonKing]

Setting the environment variable at running time doesn't work either. That's a bug that has always been present in MeatAxe. So, I guess I need to patch the sources.

[SimonKing]

Internally, MeatAxe stores MTXLIB in MtxLibDir?. But directly defining MtxLibDir? at run time didn't help. Sigh.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​71c9b94	Install meataxe header; no separate meataxe folder in SAGE_LOCAL/lib/
​0320de9	Declare that meataxe is an optional package
​e0ca6d0	Use -fPIC for creating libmtx.a
​6e6cdd8	Choose a location for multiplication tables
​6c32470	Patch MeatAxe so that MTXLIB is used to store library files

[jhpalmieri]

You shouldn't assume that the Sage directory is writeable. So the multiplication tables should be stored in SAGE_TMP or (if you want to keep them) in some other subdirectory of DOT_SAGE, probably DOT_SAGE/meataxe-2.4.24/.

[SimonKing]

Recent changes:

• Since we want to wrap the code, it is needed to make the meataxe header available.
• I found that packages are supposed to give a file called "type". I declare there that the package is supposed to be optional.
• Since we move libmtx.a from a temporary location to somewhere in SAGE_LOCAL, it is needed to use "-fPIC".
• The multiplication tables shouldn't pollute the user's current directory. I propose to store them in SAGE_SHARE/meataxe/.
• MeatAxe has IO-functions that accept an argument saying that the file is a library file---but in the IO-functions, it is always first looked into the current directory, and the library location is only considered if the current directory didn't work.
• In at least one function, the IO-argument for "library files" is used wrongly. So, after fixing the IO-function, the wrong argument needed to be fixed, too.

In meataxe/patches is a patch that fixes these two upstream bugs.

Current status:

• The package builds fine.
• When running the testsuite, I sometimes see a wrong checksum, but when I repeat the testsuite, everything is fine. I guess upstream's checksum function is fragile. Not sure if I'm going to fix that.
• With the Cython wrapper on my laptop, I can create MeatAxe matrices in a Sage session, and the current directory is not polluted by multiplication tables.

TODO:

• Finish the Cython wrapper.

[SimonKing]

Replying to jhpalmieri:

You shouldn't assume that the Sage directory is writeable. So the multiplication tables should be stored in SAGE_TMP or (if you want to keep them) in some other subdirectory of DOT_SAGE, probably DOT_SAGE/meataxe-2.4.24/.

Right, that could be a problem. The tables should indeed be permanent. So, DOT_SAGE/meataxe (without version number) it is.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​d9675fb

Use DOT_SAGE, not SAGE_LOCAL, for storing multiplication tables

[jhpalmieri]

Replying to SimonKing:

Replying to jhpalmieri:

You shouldn't assume that the Sage directory is writeable. So the multiplication tables should be stored in SAGE_TMP or (if you want to keep them) in some other subdirectory of DOT_SAGE, probably DOT_SAGE/meataxe-2.4.24/.

Right, that could be a problem. The tables should indeed be permanent. So, DOT_SAGE/meataxe (without version number) it is.

As long as you're pretty confident that the tables will be compatible with future versions. (We ran into this problem with IPython and matplotlib, which is why their directories in DOT_SAGE include the version number.)

[SimonKing]

Replying to jhpalmieri:

As long as you're pretty confident that the tables will be compatible with future versions. (We ran into this problem with IPython and matplotlib, which is why their directories in DOT_SAGE include the version number.)

The format hasn't changed since version 2.2.3, and the current version already is 4 years old. So, I don't think there will be changes any time soon.

Moreover, if we really need to change it, we can create all multiplication tables (we have bounded field size) in a fraction of a second during package installation.

[SimonKing]

Replying to SimonKing:

Moreover, if we really need to change it, we can create all multiplication tables (we have bounded field size) in a fraction of a second during package installation.

Aaahm, if we use SAGE_SHARE.

Anyway, if things really go bad in future meataxe versions, it should be possible to patch it so that invalid old multiplication tables are automatically removed when accidentally opening them.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​20a16e1

Implement and use Strassen-Winograd matrix multiplication in MeatAxe

[SimonKing]

In the recent commit, patches StrassenWinogradImplementation.patch and StrassenWinogradUsage.patch are provided.

There, I implement the Winograd-Strassen multiplication algorithm in MeatAxe. I optimized the cutoff parameter using a Cython wrapper that I didn't post yet (needs doctests and needs to wrap echelonisation).

I tested on matrices of different sizes and different fields. First of all, the school book multiplication in MeatAxe is vastly faster than the current matrix multiplication in Sage for all non-prime finite fields of odd characteristic (of course with field order < 256).

Moreover, it turns out that the Winograd-Strassen implementation is not only asymptotically fast, but even on small examples (200x200) is as fast as MeatAxe's school book multiplication. So, I think it makes sense to use it not only in my Cython wrapper but also in the MeatAxe functions/executables, which is done in StrassenWinogradUsage.patch.

The fact that MeatAxe's test suite still works is an indicative that my Winograd-Strassen implementation ok.

Of course it will be even clearer in the Cython wrapper that will be subject of the next commit. Then, I'll show timings.

[jdemeyer]

The tarball is badly made: there should be top-level directory like meataxe-2.4.24 with everything contained in that directory.

Also, we no longer list maintainers in SPKG.txt (#19238)

[jdemeyer]

Patches need to be documented (preferably in the patch itself). In particular: where do the patches come from?

[SimonKing]

Replying to jdemeyer:

The tarball is badly made:

Complain to upstream. The tarball is directly from their website. I'm not going to fix that.

[SimonKing]

Replying to jdemeyer:

Patches need to be documented (preferably in the patch itself).

Right, I did so (if I recall correctly) for the first patch concerning IO-fixes, but failed for the other two.

In particular: where do the patches come from?

My name suffices?

[SimonKing]

Replying to SimonKing:

Replying to jdemeyer:

The tarball is badly made:

Complain to upstream. The tarball is directly from their website. I'm not going to fix that.

Or would it be ok to repack the upstream tarball? I thought the purpose of the new-style packages is to use clean upstream sources?

[SimonKing]

I guess nobody except myself uses the branch. Hence, I will rebase and force-push the changes requested by Jeroen, unless someone tells me not to.

[jdemeyer]

Replying to SimonKing:

Or would it be ok to repack the upstream tarball?

Yes. In this case, it is clearly needed.

I am going to complain to upstream by the way.

[jdemeyer]

Replying to SimonKing:

My name suffices?

Yes, it should somehow be clear that you wrote the patch, that it's not an official patch from upstream. Have the patches been submitted upstream? They should be!

[jdemeyer]

Replying to SimonKing:

Hence, I will rebase and force-push the changes requested by Jeroen, unless someone tells me not to.

Go ahead. If you're rebasing, you might as well squash everything in one commit (with the f option in git rebase -i)

[SimonKing]

After replacing the old tarball by a new one called meataxe-2.4.24.p0.tar.gz, I did ./sage -fix-pkg-checksums. However, it did not produce a new checksum.ini in build/pkgs/meataxe/.

What can I do now?

[SimonKing]

Aha. Apparently it was needed to call it meataxe-2.4.24.tar.gz without .p0.

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​dea8e06	An optional MeatAxe package
​a722d4e	Implement and use Strassen-Winograd matrix multiplication in MeatAxe

[SimonKing]

Repackaged upstream sources

[SimonKing]

Replying to jdemeyer:

Go ahead. If you're rebasing, you might as well squash everything in one commit (with the f option in git rebase -i)

Done. I kept two commits of largely different purpose. The first makes MeatAxe build and install its tables into DOT_SAGE. The second commit implements and uses Winograd-Strassen multiplication in MeatAxe.

The self-tests (by sage -f -c meataxe) should pass. I don't know why they are a little flaky. Perhaps a race condition? Anyway, I think that's an upstream problem. Most of the time, the self-tests just pass fine.

The third commit will be cython wrapper for MeatAxe with which I can demonstrate that things work correct and very fast. And then I have to deal with doctests.

[SimonKing]

Replying to jdemeyer:

Replying to SimonKing:

My name suffices?

Yes, it should somehow be clear that you wrote the patch, that it's not an official patch from upstream. Have the patches been submitted upstream? They should be!

Long time ago, I have submitted a Winograd-Strassen patch for version 2.2.4 to upstream. But too much has changed since version 2.2.4, the old patches are useless.

But for sure it makes sense to submit the new three patches to upstream. Even the first one is useful. I wonder why upstream keeps telling that multiplication tables are stored in MTX_LIB when in the last 10 years they have always been installed in the current directory.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​539b134

A very basic MeatAxe Cython wrapper

[jdemeyer]

Replying to SimonKing:

Long time ago, I have submitted a Winograd-Strassen patch for version 2.2.4 to upstream. But too much has changed since version 2.2.4, the old patches are useless.

I don't understand. In this ticket, you are adding patches to a package which you claim yourself are useless???

[SimonKing]

The last commit really isn't finished yet (most tests and much documentation is missing). But for now it is enough to demonstrate the progress compared with the current Sage beta.

Please make sure that you use attachment:meataxe-2.4.24.tar.gz​, not the source tarball from the upstream download page.

In the latest development version, we have

sage: MS = MatrixSpace(GF(5^3,'y'),2000)
sage: %time A = MS.random_element()
CPU times: user 6.71 s, sys: 20 ms, total: 6.73 s
Wall time: 6.79 s
sage: type(A)
<type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'>
sage: B = MS.random_element()
sage: %time A*B
CPU times: user 23min 49s, sys: 490 ms, total: 23min 50s
Wall time: 23min 55s # using 7.5% of my computer's memory
2000 x 2000 dense matrix over Finite Field in y of size 5^3 (use the '.str()' method to see the entries)
sage: %time ~A
CPU times: user 23min 42s, sys: 842 ms, total: 23min 43s
Wall time: 23min 47s # using 11.3% of my computer's memory
2000 x 2000 dense matrix over Finite Field in y of size 5^3 (use the '.str()' method to see the entries)

23 minutes for a single multiplication or inversion of a not so big matrix over a small field!!

You see why I see the need of a different backend?

Here is the same example when the optional meataxe package (and its wrapper) is used:

sage: MS = MatrixSpace(GF(5^3,'y'),2000)
sage: %time A = MS.random_element()
CPU times: user 320 ms, sys: 5 ms, total: 325 ms
Wall time: 351 ms
sage: type(A)
<type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
sage: B = MS.random_element()
sage: %time A*B  # using 4.8% of my computer's memory
CPU times: user 14.4 s, sys: 1 ms, total: 14.4 s
Wall time: 14.4 s
2000 x 2000 dense matrix over Finite Field in y of size 5^3 (use the '.str()' method to see the entries)
sage: %time ~A   # using 5.9% of my computer's memory
CPU times: user 32.7 s, sys: 69 ms, total: 32.8 s
Wall time: 33.2 s
2000 x 2000 dense matrix over Finite Field in y of size 5^3 (use the '.str()' method to see the entries)
sage: %time _*A == A*_ == 1
CPU times: user 28.7 s, sys: 60 ms, total: 28.8 s
Wall time: 29 s
True

The last example indicates that the (Strassen) multiplication is vastly faster than what we have now in Sage.

Admittedly, the patched MeatAxe is still a lot slower than linbox for prime fields or M4RIE for characteristic 2. That's why I only suggest to use it for non-prime fields of odd characteristic.

By the way, MeatAxe is supposed to have a different kernel for fields of size > 255. But strangely, the respective source file is not in the upstream source tarball.

[jdemeyer]

The meataxe part of src/sage/matrix/matrix_modpn_dense.pxd should probably be moved to src/sage/libs/meataxe.pxd or src/sage/libs/meataxe/multiple_files_here.pxd

Also, can you use the standard heading template from ​http://doc.sagemath.org/html/en/developer/coding_basics.html#headings-of-sage-library-code-files

[SimonKing]

Replying to jdemeyer:

Replying to SimonKing:

Long time ago, I have submitted a Winograd-Strassen patch for version 2.2.4 to upstream. But too much has changed since version 2.2.4, the old patches are useless.

I don't understand. In this ticket, you are adding patches to a package which you claim yourself are useless???

The patches FOR THE OLD VERSION 2.2.4 that I send to upstream a couple of years ago are useless FOR THE CURRENT VERSION 2.4.24. All function and struct names in MeatAxe have changed between version 2.2.4 and 2.4.24.

[jdemeyer]

Replying to SimonKing:

The patches FOR THE OLD VERSION 2.2.4 that I send to upstream a couple of years ago are useless FOR THE CURRENT VERSION 2.4.24. All function and struct names in MeatAxe have changed between version 2.2.4 and 2.4.24.

OK, but does it make sense to send the new patches to upstream again?

[SimonKing]

Replying to jdemeyer:

src/sage/matrix/matrix_modpn_dense.pxd should probably be moved to src/sage/libs/meataxe.pxd or src/sage/libs/meataxe/multiple_files_here.pxd

I don't think that meataxe.pxd is a good name. It is a dense implementation for matrices over GF(p^n), and thus it is matrix_modpn_dense.pxd. This is totally analogous to matrix_mod2e_dense.pxd; or do you suggest to rename the latter by m4rie.pxd?

Also, can you use the standard heading template from ​http://doc.sagemath.org/html/en/developer/coding_basics.html#headings-of-sage-library-code-files

I think I have copied the heading from some existing pxd files. But I can change it.

[jdemeyer]

If I understand things correctly, the modpn in the filename refers to a finite field of size p^n, right? That's a very confusing name, can you use matrix_gfpn_dense.pyx or something?

[SimonKing]

Replying to jdemeyer:

Replying to SimonKing:

The patches FOR THE OLD VERSION 2.2.4 that I send to upstream a couple of years ago are useless FOR THE CURRENT VERSION 2.4.24. All function and struct names in MeatAxe have changed between version 2.2.4 and 2.4.24.

OK, but does it make sense to send the new patches to upstream again?

As I said in comment:64: "But for sure it makes sense to submit the new three patches to upstream."

[jdemeyer]

Replying to SimonKing:

This is totally analogous to matrix_mod2e_dense.pxd; or do you suggest to rename the latter by m4rie.pxd?

That file actually exists: src/sage/libs/m4rie.pxd

It's cimported by matrix_mod2e_dense.pxd, which is what you should do for meataxe also: separate the Sage matrix stuff from the meataxe library stuff.

[SimonKing]

Replying to jdemeyer:

If I understand things correctly, the modpn in the filename refers to a finite field of size p^n, right? That's a very confusing name, can you use matrix_gfpn_dense.pyx or something?

Do you suggest to rename Martin Albrecht's matrix_mod2e_dense.pyx to matrix_gf2e_dense.pyx? By the way, IIRC that's where I got the heading from.

I think we should be consistent with the existing naming scheme. We currently have matrix_modn_dense, matrix_mod2_dense, matrix_mod2e_dense, and the obvious continuation of that scheme is matrix_modpn_dense.

[SimonKing]

Replying to jdemeyer:

Replying to SimonKing:

This is totally analogous to matrix_mod2e_dense.pxd; or do you suggest to rename the latter by m4rie.pxd?

That file actually exists: src/sage/libs/m4rie.pxd

Aha! Now I understand.

OK, that makes more sense. Have a src/sage/libs/meataxe.pxd without a corresponding .pyx file, and do the rest in src/sage/matrix/matrix_modpn_dense.p??, which imports from sage.libs.meataxe.

[jdemeyer]

Replying to SimonKing:

Do you suggest to rename Martin Albrecht's matrix_mod2e_dense.pyx to matrix_gf2e_dense.pyx?

Yes!! It's a very confusing name.

[jdemeyer]

Replying to SimonKing:

Have a src/sage/libs/meataxe.pxd without a corresponding .pyx file, and do the rest in src/sage/matrix/matrix_modpn_dense.p??, which imports from sage.libs.meataxe.

Yes, exactly. It's good practice to separate the library (which has nothing to do with Sage) from the Sage interface.

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​a8f1847

A very basic MeatAxe Cython wrapper

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​da665d9

A very basic MeatAxe Cython wrapper

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​106aeef

A very basic MeatAxe Cython wrapper

[SimonKing]

Replying to jdemeyer:

Replying to SimonKing:

Do you suggest to rename Martin Albrecht's matrix_mod2e_dense.pyx to matrix_gf2e_dense.pyx?

Yes!! It's a very confusing name.

OK, if you are so determined, then I better change my module and class name, which is done in the previous (forced) commit.

I think I will now let it rest for a few hours, and then I'll tell how I found the default cutoff for Winograd-Strassen multiplication.

[SimonKing]

I tested that classical multiplication yields the same result as Strassen multiplication, for various fields (GF(q) for q=4,8,9,25,27,32,125,243). No problem. Then I tested GF(2)---and got different results. Sigh. I don't know what's wrong with 2, but I need to debug it.

[jdemeyer]

In meataxe.pxd, I assume that ctypedef struct FILE refers to the <stdio.h> type FILE? In that case, it's better to use

from libc.stdio cimport FILE

(or even remove that line altogether, since you don't actually use FILE)

[jdemeyer]

In spkg-install, why this?

PATCHES=../patches
for patch in "$PATCHES"/*.patch; do

It's not wrong, but usually one simply does

for patch in ../patches/*.patch; do

[SimonKing]

Replying to jdemeyer:

In meataxe.pxd, I assume that ctypedef struct FILE refers to the <stdio.h> type FILE? In that case, it's better to use

from libc.stdio cimport FILE

(or even remove that line altogether, since you don't actually use FILE)

I refer to whatever MeatAxe refers to. But you are right, I don't use it. Similarly I don't use Ulong, Ushort etc. I left them in for completeness.

[SimonKing]

Replying to jdemeyer:

In spkg-install, why this?

PATCHES=../patches
for patch in "$PATCHES"/*.patch; do

It's not wrong, but usually one simply does

for patch in ../patches/*.patch; do

I guess I copied it from somewhere. But "somewhere" probably means a package in which the "patches" folder is referred to repeatedly.

[SimonKing]

Meanwhile I know that the wrong results of Strassen multiplication in field order 2 come from the function FfAddMapRowWindow, which in fact has a special case for the field of order 2. This function is a modification of the MeatAxe function FfMapRow. So, I have to find out where my modifications went wrong.

[jdemeyer]

Replying to SimonKing:

I refer to whatever MeatAxe refers to.

I don't understand what you mean with "refers to".

I hope that meataxe does not redefine FILE (otherwise it's more broken than I thought).

Similarly I don't use Ulong, Ushort etc.

This is different, since these are meataxe-specific.

[SimonKing]

Replying to jdemeyer:

Replying to SimonKing:

I refer to whatever MeatAxe refers to.

I don't understand what you mean with "refers to".

Hm. Now I am puzzled myself. I thought for a few minutes that I found it in meataxe.h of version 2.2.4, and when I ported my wrapper to version 2.4.24 I didn't check if it is still there. But now I checked: FILE is *not* defined in meataxe.h, neither in the old nor the new version.

I hope that meataxe does not redefine FILE (otherwise it's more broken than I thought).

It doesn't. It does #include <stdio.h> both in the old and the new version.

[jdemeyer]

Replying to SimonKing:

It doesn't. It does #include <stdio.h> both in the old and the new version.

Excellent! Then there is certainly no reason to define FILE the way you did in meataxe.pxd.

[SimonKing]

Replying to SimonKing:

Meanwhile I know that the wrong results of Strassen multiplication in field order 2 come from the function FfAddMapRowWindow, which in fact has a special case for the field of order 2. This function is a modification of the MeatAxe function FfMapRow.

Got it! There is a pointer that is moved through the matrix by which we are multiplying. In case of order two, if an addition of a row happens then the pointer is moved to the next item after the current row of the matrix.

In FfMapRow, this position automatically is the first item of the next row of the matrix. In FfAddMapRowWindow, the matrix is in fact just a window of a larger matrix, and this position generally is outside of the window. Hence, I have to manually move the pointer to the first position of the next line of the window.

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​de61563	An optional MeatAxe package
​f4873f7	Implement and use Strassen-Winograd matrix multiplication in MeatAxe
​c1d5026	A very basic MeatAxe Cython wrapper

[SimonKing]

I force-pushed the new version.

Changes: The bug is fixed, and FILE and PATCHES are removed.

Now for timings and the quest for a good (default) cutoff.

When using Winograd-Strassen, one needs to define a matrix size such that divide-and-conquer is applied if and only if all resulting blocks have at least that size.

Meataxe packs matrix coefficients densely. I.e, if you have a field of size 3, then it will pack 5 matrix coefficients into one byte (2^5<256).

Each meataxe matrix is stored in a block of memory, row by row. Internally, meataxe uses "long" operations in some functions. That's why each row uses a byte size that is a multiple of sizeof(long), using up to sizeof(long)-1 trailing zero bytes ("padding bytes").

Winograd-Strassen operates on matrix windows. It is rather obvious that the window shouldn't end in the middle of a byte, i.e., the number of columns of a window should be a multiple of the "number of coefficients per byte".

In my old patches for version 2.2.4, in fact I used "byte" as unit. But in the new patches for version 2.4.24, I adopt meataxe's convention that row length should be divisible by sizeof(long). So, in my Winograd-Strassen implementation, the number of columns of matrix windows is a multiple of "sizeof(long)*(number of coefficients)".

Consequently, the cutoff parameter in meataxe is given in the unit "sizeof(long)" (in my current cython wrapper the unit is byte, but that's likely to change).

With the following function used in a Sage session, I have determined the optimal cutoff parameter for different field sizes and different matrix dimensions. Also I have tested for what matrix dimensions the school book algorithm is faster than the Winograd-Strassen algorithm, and I checked that both algorithms yield the same result.

from sage.matrix.matrix_gfpn_dense import Matrix_gfpn_dense
def find_cutoff(q, D):
    print "GF({})".format(q)
    for n in [100,500,1000,2000]:
        print "  n = {}:".format(n),
        MS = MatrixSpace(GF(q,'x'), n)
        MS._MatrixSpace__matrix_class = Matrix_gfpn_dense
        M = MS.random_element()
        t = cputime()
        N = M._multiply_classical(M)
        D[q,n] = set()
        D[q,n].add((cputime(t),0))
        for c in range(1,9):
            print c*8,
            sys.stdout.flush()
            t = cputime()
            N1 = M._multiply_strassen(M, c*8)
            D[q,n].add((cputime(t), c))
            assert N1==N
        print
        print "->",sorted(list((D[q,n])))[:5]

The result was

sage: D = {}
sage: for q in [2,3,8,9,25,27,32,49,81,125,243]: find_cutoff(q, D)
GF(2)
  n = 100: 8 16 24 32 40 48 56 64
-> [(0.0, 0), (0.0, 1), (0.0, 2), (0.0, 3), (0.0, 4)]
  n = 500: 8 16 24 32 40 48 56 64
-> [(0.002999999999985903, 5), (0.002999999999985903, 6), (0.002999999999985903, 7), (0.0030000000000427463, 4), (0.0030000000000427463, 8)]
  n = 1000: 8 16 24 32 40 48 56 64
-> [(0.019000000000005457, 0), (0.01999999999998181, 8), (0.02199999999999136, 4), (0.02199999999999136, 5), (0.02199999999999136, 7)]
  n = 2000: 8 16 24 32 40 48 56 64
-> [(0.12000000000000455, 0), (0.14599999999995816, 8), (0.1560000000000059, 4), (0.15699999999998226, 5), (0.15699999999998226, 6)]
GF(3)
  n = 100: 8 16 24 32 40 48 56 64
-> [(0.0, 0), (0.0, 4), (0.0, 7), (0.0009999999999763531, 2), (0.0009999999999763531, 3)]
  n = 500: 8 16 24 32 40 48 56 64
-> [(0.04899999999997817, 0), (0.05000000000001137, 7), (0.05099999999998772, 8), (0.051999999999964075, 5), (0.05299999999999727, 4)]
  n = 1000: 8 16 24 32 40 48 56 64
-> [(0.3509999999999991, 4), (0.3509999999999991, 5), (0.3509999999999991, 6), (0.36199999999996635, 8), (0.36299999999999955, 7)]
  n = 2000: 8 16 24 32 40 48 56 64
-> [(2.4659999999999513, 6), (2.4669999999999845, 4), (2.4669999999999845, 5), (2.5469999999999686, 7), (2.548000000000002, 8)]
GF(8)
  n = 100: 8 16 24 32 40 48 56 64
-> [(0.0009999999999763531, 3), (0.0009999999999763531, 6), (0.0009999999999763531, 8), (0.0010000000000331966, 1), (0.0019999999999527063, 0)]
  n = 500: 8 16 24 32 40 48 56 64
-> [(0.1169999999999618, 3), (0.117999999999995, 2), (0.1209999999999809, 5), (0.12100000000003774, 4), (0.12100000000003774, 6)]
  n = 1000: 8 16 24 32 40 48 56 64
-> [(0.9129999999999541, 7), (0.9130000000000109, 4), (0.9130000000000109, 6), (0.9139999999999873, 5), (0.9519999999999982, 3)]
  n = 2000: 8 16 24 32 40 48 56 64
-> [(6.291999999999973, 7), (6.29200000000003, 4), (6.293000000000006, 5), (6.293000000000006, 6), (6.569000000000017, 2)]
GF(9)
  n = 100: 8 16 24 32 40 48 56 64
-> [(0.0009999999999763531, 7), (0.0010000000000331966, 2), (0.0010000000000331966, 4), (0.0019999999999527063, 0), (0.0019999999999527063, 3)]
  n = 500: 8 16 24 32 40 48 56 64
-> [(0.12299999999999045, 3), (0.1239999999999668, 7), (0.12400000000002365, 2), (0.125, 4), (0.125, 5)]
  n = 1000: 8 16 24 32 40 48 56 64
-> [(0.9429999999999836, 4), (0.9429999999999836, 6), (0.9430000000000405, 7), (0.9440000000000168, 5), (0.9830000000000041, 8)]
  n = 2000: 8 16 24 32 40 48 56 64
-> [(6.517999999999972, 6), (6.522999999999968, 7), (6.524000000000001, 5), (6.5330000000000155, 4), (6.793999999999983, 8)]
GF(25)
  n = 100: 8 16 24 32 40 48 56 64
-> [(0.002999999999985903, 0), (0.002999999999985903, 2), (0.002999999999985903, 3), (0.002999999999985903, 4), (0.002999999999985903, 6)]
  n = 500: 8 16 24 32 40 48 56 64
-> [(0.28300000000001546, 4), (0.2839999999999918, 6), (0.2839999999999918, 7), (0.2840000000001055, 5), (0.2949999999999591, 8)]
  n = 1000: 8 16 24 32 40 48 56 64
-> [(1.9640000000000555, 7), (1.9650000000000318, 6), (1.9669999999999845, 4), (1.969000000000051, 5), (2.044999999999959, 8)]
  n = 2000: 8 16 24 32 40 48 56 64
-> [(13.799999999999955, 4), (13.805000000000064, 6), (13.807000000000016, 7), (13.81499999999994, 5), (14.360000000000014, 8)]
GF(27)
  n = 100: 8 16 24 32 40 48 56 64
-> [(0.0029999999999290594, 0), (0.0029999999999290594, 4), (0.0029999999999290594, 6), (0.0030000000000427463, 2), (0.0030000000000427463, 3)]
  n = 500: 8 16 24 32 40 48 56 64
-> [(0.28099999999994907, 5), (0.28099999999994907, 7), (0.28100000000006276, 4), (0.28100000000006276, 6), (0.29300000000000637, 8)]
  n = 1000: 8 16 24 32 40 48 56 64
-> [(1.9509999999999081, 4), (1.9510000000000218, 6), (1.9510000000000218, 7), (1.9520000000001119, 5), (2.033999999999878, 8)]
  n = 2000: 8 16 24 32 40 48 56 64
-> [(13.709999999999923, 6), (13.718000000000075, 7), (13.731999999999971, 5), (13.764999999999873, 4), (14.337999999999965, 8)]
GF(32)
  n = 100: 8 16 24 32 40 48 56 64
-> [(0.0029999999999290594, 2), (0.0029999999999290594, 5), (0.0029999999999290594, 7), (0.0030000000000427463, 0), (0.0030000000000427463, 3)]
  n = 500: 8 16 24 32 40 48 56 64
-> [(0.2700000000000955, 7), (0.27099999999995816, 4), (0.27099999999995816, 6), (0.2720000000000482, 5), (0.2799999999999727, 2)]
  n = 1000: 8 16 24 32 40 48 56 64
-> [(1.8730000000000473, 4), (1.893000000000029, 5), (1.8960000000000719, 7), (1.8980000000000246, 6), (1.9220000000000255, 2)]
  n = 2000: 8 16 24 32 40 48 56 64
-> [(13.370999999999981, 4), (13.513000000000034, 5), (13.535999999999945, 2), (13.625, 3), (13.687000000000012, 7)]
GF(49)
  n = 100: 8 16 24 32 40 48 56 64
-> [(0.0029999999999290594, 4), (0.0029999999999290594, 7), (0.0030000000000427463, 0), (0.0030000000000427463, 2), (0.0030000000000427463, 3)]
  n = 500: 8 16 24 32 40 48 56 64
-> [(0.28300000000001546, 4), (0.2860000000000582, 6), (0.2869999999999209, 5), (0.28700000000003456, 7), (0.2939999999999827, 8)]
  n = 1000: 8 16 24 32 40 48 56 64
-> [(1.9269999999999072, 6), (1.9300000000000637, 7), (1.9329999999999927, 5), (1.9470000000000027, 4), (2.038000000000011, 8)]
  n = 2000: 8 16 24 32 40 48 56 64
-> [(13.607000000000198, 6), (13.632000000000062, 7), (13.685999999999922, 5), (13.788999999999987, 4), (14.29099999999994, 2)]
GF(81)
  n = 100: 8 16 24 32 40 48 56 64
-> [(0.0029999999999290594, 2), (0.0029999999999290594, 4), (0.0029999999999290594, 5), (0.0029999999999290594, 7), (0.0029999999999290594, 8)]
  n = 500: 8 16 24 32 40 48 56 64
-> [(0.2819999999999254, 4), (0.2819999999999254, 7), (0.2840000000001055, 6), (0.28899999999998727, 5), (0.29500000000007276, 2)]
  n = 1000: 8 16 24 32 40 48 56 64
-> [(1.9349999999999454, 6), (1.9409999999998035, 7), (1.9610000000000127, 4), (1.9629999999999654, 5), (2.0399999999999636, 8)]
  n = 2000: 8 16 24 32 40 48 56 64
-> [(13.505999999999858, 7), (13.623000000000047, 5), (13.644000000000233, 6), (13.719000000000278, 4), (14.055000000000064, 8)]
GF(125)
  n = 100: 8 16 24 32 40 48 56 64
-> [(0.0029999999999290594, 0), (0.0029999999999290594, 2), (0.0029999999999290594, 4), (0.0029999999999290594, 5), (0.0029999999999290594, 7)]
  n = 500: 8 16 24 32 40 48 56 64
-> [(0.29599999999982174, 5), (0.2960000000000491, 4), (0.2960000000000491, 6), (0.2960000000000491, 7), (0.3009999999999309, 8)]
  n = 1000: 8 16 24 32 40 48 56 64
-> [(2.0429999999998927, 5), (2.04300000000012, 7), (2.0449999999998454, 4), (2.0450000000000728, 6), (2.0960000000000036, 8)]
  n = 2000: 8 16 24 32 40 48 56 64
-> [(14.353999999999814, 7), (14.366000000000213, 5), (14.383000000000266, 4), (14.414999999999964, 6), (14.596000000000004, 8)]
GF(243)
  n = 100: 8 16 24 32 40 48 56 64
-> [(0.0029999999999290594, 4), (0.0029999999999290594, 6), (0.003000000000156433, 0), (0.003000000000156433, 2), (0.003000000000156433, 8)]
  n = 500: 8 16 24 32 40 48 56 64
-> [(0.3369999999999891, 5), (0.3369999999999891, 7), (0.33799999999996544, 4), (0.3380000000001928, 6), (0.34799999999995634, 8)]
  n = 1000: 8 16 24 32 40 48 56 64
-> [(2.336999999999989, 4), (2.336999999999989, 6), (2.336999999999989, 7), (2.341999999999871, 5), (2.4020000000000437, 8)]
  n = 2000: 8 16 24 32 40 48 56 64
-> [(17.700999999999794, 4), (17.958000000000084, 5), (17.998000000000047, 7), (18.02800000000002, 2), (18.03599999999983, 6)]

Evaluation of these timings

• My Winograd-Strassen implementation gives the same result as meataxe's school book implementation. --> Hooray
• There is not a single case (not even for small matrices) where meataxe's school book implementation is significantly faster than Winograd-Strassen --> Winograd-Strassen should be the default algorithm in matrix_gfpn_dense
• It seems that a cutoff between 4 and 7 in the unit "sizeof(long)" are about the same, and are generally better than other cutoffs.

The last statement is backed up by another study, where I consider 3601x3601-matrices over different fields, varying the cutoff parameter. I got these timings for multiplication:

cutoff       1       2        3       4        5        6        7         8
GF(4):     15.5 s   12.2 s  12.3 s   11 s     11 s     11 s     11 s      11.5 s
GF(9):     41.3 s   33.7 s  33.5 s   31.6 s   31.5 s   31.5 s   31.5 s    34 s
GF(16):    26.4 s   21.1 s   21 s    19.9 s   19.9 s   19.9 s   19.9 s    21.4 s
GF(25):   1min 34s 1min 12s 1min 12s 1min 8s  1min 8s  1min 8s  1min 8s  1min 12s
GF(32):   1min 27s 1min 6s  1min 6s  1min 4s  1min 4s  1min 4s  1min 4s  1min 10s
GF(243):  1min 55s 1min 27s 1min 27s 1min 18s 1min 18s 1min 18s 1min 19s 1min 24s

Now, the real question is: What is the meaning of "4"? My guess is that "4" is "sizeof(long)/2", so that the default cutoff parameter now is "sizeof(long)/2". But I could be mistaken. Perhaps on a 32-bit machine, cutoff 4 is better than cutoff "2==sizeof(long)/2", and on a 128-bit machine, cutoff 4 is better than cutoff "8=sizeof(long)/2"? Or perhaps it all depends on L1/L2 cache size?

Request to the reviewer

• Do you have access to machines of different architecture (32-bit, 128-bit, different cache sizes)? Could you try to use the function above to find the optimal cutoff parameters?
• Meataxe uses an environment variable ASM_MMX in the file kernel-0.c. It means that in some functions it uses assembler code in characteristic two. By copying I use it as well, in FfAddMapRowWindow mentioned above. In some comment, it says that "this assumes Intel with 4 bytes per long, but MMX implies Intel anyway."

My questions:

• Is ASM_MMX some kind of standard variable?
• Can you test the package on a machine where ASM_MMX is used?

Note that I do have Intel, but my longs are 8 bytes, not 4.

[SimonKing]

Just to see if I can improve the code a bit, I did

sage: M = MatrixSpace(GF(125,'x'),2000).random_element()
sage: %crun N = M*M
PROFILE: interrupts/evictions/bytes = 1466/1183/718104
Using local file /home/king/Sage/git/sage/local/bin/python.
Using local file /home/king/.sage/temp/linux-va3e.site/18759/tmp_efM3tZ.perf.
Total: 1466 samples
       0   0.0%   0.0%     1466 100.0% MatMulStrassen
       0   0.0%   0.0%     1466 100.0% PyEval_EvalCode
       0   0.0%   0.0%     1466 100.0% PyEval_EvalCodeEx
       0   0.0%   0.0%     1466 100.0% PyEval_EvalFrameEx
       0   0.0%   0.0%     1466 100.0% PyNumber_Multiply
       0   0.0%   0.0%     1466 100.0% PyObject_Call
       0   0.0%   0.0%     1466 100.0% PyRun_FileExFlags
       0   0.0%   0.0%     1466 100.0% PyRun_SimpleFileExFlags
       0   0.0%   0.0%     1466 100.0% PyRun_StringFlags
       0   0.0%   0.0%     1466 100.0% Py_Main
       0   0.0%   0.0%     1466 100.0% StrassenStep
       0   0.0%   0.0%     1466 100.0% __libc_start_main
       0   0.0%   0.0%     1466 100.0% __pyx_f_4sage_6matrix_17matrix_gfpn_dense_17Matrix_gfpn_dense__matrix_times_matrix_
       0   0.0%   0.0%     1466 100.0% __pyx_f_4sage_6matrix_17matrix_gfpn_dense_17Matrix_gfpn_dense__multiply_strassen
       0   0.0%   0.0%     1466 100.0% __pyx_pf_4sage_9structure_7element_6Matrix_2__mul__ (inline)
       0   0.0%   0.0%     1466 100.0% __pyx_pw_4sage_9structure_7element_6Matrix_3__mul__
       0   0.0%   0.0%     1466 100.0% _start
       0   0.0%   0.0%     1466 100.0% binary_op1 (inline)
       0   0.0%   0.0%     1466 100.0% call_function (inline)
       0   0.0%   0.0%     1466 100.0% exec_statement (inline)
       0   0.0%   0.0%     1466 100.0% ext_do_call (inline)
       0   0.0%   0.0%     1466 100.0% fast_function (inline)
       0   0.0%   0.0%     1466 100.0% function_call
       0   0.0%   0.0%     1466 100.0% run_mod (inline)
       2   0.1%   0.1%     1382  94.3% WindowAddMul
    1368  93.3%  93.5%     1368  93.3% FfAddMapRowWindow
       0   0.0%  93.5%       46   3.1% WindowDif
       2   0.1%  93.6%       31   2.1% WindowSum
      30   2.0%  95.6%       30   2.0% FfSetNoc
      23   1.6%  97.2%       23   1.6% FfAddRowPartial
      23   1.6%  98.8%       23   1.6% FfSubRowPartial
       8   0.5%  99.3%        8   0.5% FfSubRowPartialReverse
       0   0.0%  99.3%        4   0.3% WindowClear
       4   0.3%  99.6%        4   0.3% __GI_memset
       0   0.0%  99.6%        3   0.2% FfAlloc
       0   0.0%  99.6%        3   0.2% MatAlloc
       0   0.0%  99.6%        3   0.2% WindowAlloc
       2   0.1%  99.7%        2   0.1% FfStepPtr
       0   0.0%  99.7%        2   0.1% SysMalloc
       1   0.1%  99.8%        2   0.1% __GI___libc_malloc
       0   0.0%  99.8%        1   0.1% FfAddRow
       0   0.0%  99.8%        1   0.1% FfGetPtr
       1   0.1%  99.9%        1   0.1% FfMulRow
       1   0.1%  99.9%        1   0.1% __memcpy_sse2_unaligned
       0   0.0%  99.9%        1   0.1% _int_malloc
       1   0.1% 100.0%        1   0.1% malloc_consolidate
sage: %crun N = M._multiply_classical(M)
PROFILE: interrupts/evictions/bytes = 2242/24/17464
Using local file /home/king/Sage/git/sage/local/bin/python.
Using local file /home/king/.sage/temp/linux-va3e.site/18759/tmp_R9Qzph.perf.
Total: 2242 samples
       0   0.0%   0.0%     2242 100.0% PyEval_EvalCode
       0   0.0%   0.0%     2242 100.0% PyEval_EvalCodeEx
       0   0.0%   0.0%     2242 100.0% PyEval_EvalFrameEx
       0   0.0%   0.0%     2242 100.0% PyObject_Call
       0   0.0%   0.0%     2242 100.0% PyRun_FileExFlags
       0   0.0%   0.0%     2242 100.0% PyRun_SimpleFileExFlags
       0   0.0%   0.0%     2242 100.0% PyRun_StringFlags
       0   0.0%   0.0%     2242 100.0% Py_Main
       0   0.0%   0.0%     2242 100.0% __libc_start_main
       0   0.0%   0.0%     2242 100.0% __pyx_f_4sage_6matrix_17matrix_gfpn_dense_17Matrix_gfpn_dense__multiply_classical
       0   0.0%   0.0%     2242 100.0% __pyx_pf_4sage_6matrix_17matrix_gfpn_dense_17Matrix_gfpn_dense_32_multiply_classical (inline)
       0   0.0%   0.0%     2242 100.0% __pyx_pw_4sage_6matrix_17matrix_gfpn_dense_17Matrix_gfpn_dense_33_multiply_classical
       0   0.0%   0.0%     2242 100.0% _start
       0   0.0%   0.0%     2242 100.0% call_function (inline)
       0   0.0%   0.0%     2242 100.0% exec_statement (inline)
       0   0.0%   0.0%     2242 100.0% ext_do_call (inline)
       0   0.0%   0.0%     2242 100.0% fast_function (inline)
       0   0.0%   0.0%     2242 100.0% function_call
       0   0.0%   0.0%     2242 100.0% run_mod (inline)
    2241 100.0% 100.0%     2241 100.0% FfMapRow
       0   0.0% 100.0%     2241 100.0% MatMul
       0   0.0% 100.0%        1   0.0% FfAlloc
       1   0.0% 100.0%        1   0.0% FfMulRow
       0   0.0% 100.0%        1   0.0% MatAlloc
       0   0.0% 100.0%        1   0.0% MatDup

So, if I want to optimise something then I should have a look at FfMapRow/FfAddMapRowWindow.

[SimonKing]

I tried some tweaks, but it didn't result in a clear improvement. So, I'll now focus on wrapping echelon computation and then on the doctests.

[jdemeyer]

It would be convenient if you would base this ticket on #19240 (and review that), because there will be conflicts.

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​1e6264d	Rename matrix_mod2e_dense to matrix_gf2e_dense
​26f359e	Fix unpickling old Matrix_mod2e_dense instances
​32b2cb1	Fix docstring by making it a raw string
​1e47929	Change old into new style raise statement
​8e72fb5	An optional MeatAxe package
​c7d75fb	Implement and use Strassen-Winograd matrix multiplication in MeatAxe
​9e2a8c6	A very basic MeatAxe Cython wrapper

[SimonKing]

I don't know why "git trac push 12103 --force" does not return, although apparently it has succeeded and all commits have arrived on trac. Anyway, I have rebased the branch on top of #19240.

[SimonKing]

By the way, I asked upstream about the "big" kernel (for field orders larger than 255) and got the reply that they had it in version 2.2.4, but removed in the most recent versions because apparently nobody has used it. So, if I want to use it, then I have to port the "big" kernel from 2.2.4 to 2.4.24 myself.

Very probably I will not do it. For my applications, field orders up to 255 are enough.

[SimonKing]

In order to have a fully-fledged matrix implementation, I will have a look at sage.matrix.matrix_generic_dense to see what methods need to be implemented.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​4bdd285

A full wrapper for MeatAxe matrices

[SimonKing]

The implementation of _matrix_times_matrix_ is now inherited from sage.matrix.matrix0, i.e., instead of overriding it I created a method _strassen_default_cutoff that returns 0 (which means that Strassen-Winograd is used for all matrices).

Also, I have wrapped the computation of echelon form. Meataxe only computes a semi-echelon form. Thus, I need to work a little harder to get the reduced row echelon form. Again, it is a lot faster than what currently is in Sage for nonprime finite fields of odd order.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​b237c5f

Improve echelon computation in MeatAxe

[SimonKing]

In Gaussian elimination, one adds a scalar multiple of one row to another row. In fact, one knows that pivot of the row to add, i.e., in fact one only needs to add a *part* of the row (namely the part after the pivot). However, MeatAxe did always add the complete rows, including the leading zeroes.

That's obviously a waste of time, and thus I added another patch to MeatAxe in my recent commit.

With the patch:

sage: M = MatrixSpace(GF(25,'x'),5000).random_element(density=0.001)
sage: %time N = M.echelon_form()
CPU times: user 1min 9s, sys: 156 ms, total: 1min 9s
Wall time: 1min 9s

Without the patch, it is 1min 31s. I guess that's something.

I am unsure if I will implement echelon computation à la Strassen. If I do, it is essential to be able to add *parts* of a row, because one operates on matrix windows and not on matrices.

[jdemeyer]

Following up on your sage-devel post, you should add a file build/pkgs/meataxe/dependencies containing

# no dependencies

In that case, you need to create $SAGE_LOCAL/include to avoid

cp: cannot create regular file '/usr/local/src/sage-config/local/include/': Not a directory
Error copying MeatAxe header.

[SimonKing]

Replying to jdemeyer:

Following up on your sage-devel post, you should add a file build/pkgs/meataxe/dependencies containing

# no dependencies

Thank you! Doing it now.

In that case, you need to create $SAGE_LOCAL/include to avoid

cp: cannot create regular file '/usr/local/src/sage-config/local/include/': Not a directory
Error copying MeatAxe header.

Pardon? is $SAGE_LOCAL/include not a standard Sage folder? Or do you mean that I need to do mkdir -p "$SAGE_LOCAL/include" because the local/include folder is only created during installation of Sage, making sage a dependency?

[SimonKing]

Replying to SimonKing:

Replying to jdemeyer:

Following up on your sage-devel post, you should add a file build/pkgs/meataxe/dependencies containing

# no dependencies

Thank you! Doing it now.

Didn't help (edit: ... against the three occurrences of make -j2 configure. It did help against sage -b).

[jdemeyer]

Replying to SimonKing:

Or do you mean that I need to do mkdir -p "$SAGE_LOCAL/include" because the local/include folder is only created during installation of Sage

Yes, exactly.

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​db6b7fe

Improve echelon computation in MeatAxe, and fix some compiler warnings

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​191477e

Improve echelon computation in MeatAxe, and fix some compiler warnings

[SimonKing]

Sorry for force-pushing again...

In the modified version of the previous commit, I do the "mkdir -p" thingy for all folders that I use ($SAGE_LOCAL/include and $SAGE_LOCAL/bin).

Also, I remove all but one compiler warnings for meataxe. The only remaining warning can not be removed. I need that the function MatrxiToWindow has a const argument, because it is used in a multiplication function which has const arguments. However, in MatrixToWindow, I assign the const argument to a pointer target type.

I think it is fine to ignore the warning. The matrix definitely is not altered inside of MatrixToWindow (but may of course be altered by a function that *later* modifies the window).

[SimonKing]

I think there is a very severe problem: Error handling.

I thought that meataxe would set an error and then return an error value, for example in matinv.c in the function zmatinv:

        if (f1 == FF_ZERO) 
        {
            MTX_ERROR1("%E",MTX_ERR_DIV0);
            return -1;
        }

which should propage to the function MatInverse:

    if (zmatinv(tmp,dest->Data) != 0) 
    {
        MatFree(dest);
        return NULL;
    }
    return dest;

And in my wrapper, I also test for the error return value:

        sig_on()
        OUT.Data = MatInverse(self.Data)
        sig_off()
        if OUT.Data != NULL:
            return OUT
        raise ArithmeticError("This matrix is not invertible")

However, it doesn't work. I get an immediate crash of Sage without any traceback:

sage: M = MatrixSpace(GF(27,'x'),5).random_element(density=0.4)
sage: M.rank()
4
sage: ~M
matinv.c(50):Division by zero

Hence, I need to study how meataxe manages to crash everything, even though it pretends to handle errors.

[SimonKing]

I think it is the default error handler doing the following:

static void DefaultHandler(const MtxErrorRecord_t *e)
{
    static int count = 0;
    if (LogFile == NULL)
        LogFile = stderr;

    if (e->FileInfo != NULL)
        fprintf(LogFile,"%s(%d):",e->FileInfo->BaseName,e->LineNo);
    fprintf(LogFile,"%s\n",e->Text);
    if (--count <= 0)
        exit(255);
}

Apparently Sage's sig_on()/sig_off() cannot deal with the line exit(255);, right?

So, should I patch the default error handler, making it signal something that sig_on()/off() can deal with? What would that be?

Or perhaps I know something better: I think the usual MeatAxe binaries could very well keep using the current error handler. However, in my wrapper, I could instead define a new error handler and set it with MtxSetErrorHandler.

[jdemeyer]

The fact that the error handler is called DefaultHandler seems to imply that it's possible to change the error handler...

For a simple example of a custom error handler for a C library, see the error_handler function in src/sage/libs/gap/util.pyx

[SimonKing]

Replying to jdemeyer:

The fact that the error handler is called DefaultHandler seems to imply that it's possible to change the error handler...

Sure. I already mentioned MtxSetErrorHandler.

For a simple example of a custom error handler for a C library, see the error_handler function in src/sage/libs/gap/util.pyx

Thank you, I'll have a look!

[jdemeyer]

In case you're wondering: the function PyErr_SetObject() is a low-level version of raise: it "raises" an exception but execution continues to the sig_error() which is handled by sig_on().

[SimonKing]

I am not sure if the error handler should do sig_error().

There are some places in meataxe where a function A calls a function B, and if B returns an error value then A deallocates temporary data before returning an error value as well.

Wouldn't we prevent meataxe from these clean-up operations (causing a memory leak on error) if we would raise sig_error()?

[jdemeyer]

Replying to SimonKing:

Wouldn't we prevent meataxe from these clean-up operations (causing a memory leak on error) if we would raise sig_error()?

Yes.

But since I don't know the specifics of the meataxe error handling system, I cannot comment further.

[SimonKing]

Replying to jdemeyer:

Replying to SimonKing:

Wouldn't we prevent meataxe from these clean-up operations (causing a memory leak on error) if we would raise sig_error()?

Yes.

The question is if we care. But perhaps we should, since catching errors and coping with it is typical for Python, and that has caused problems in the past (I recall some problem with the Gap interface, because a tight loop occasionally catching errors eventually exceeded the allowed number of errors in Gap).

Does meataxe care itself? That's unclear. I have seen instances of the cleaning-up-on-error, but I have seen other cases in which it does not clean up.

But since I don't know the specifics of the meataxe error handling system, I cannot comment further.

I am not sure if one can call it a "system". Anyway. There is a function that defines the error handler to-be-used. The function MtxError creates an error record (containing information on file name and line number and error code) and calls the error handler with it.

In sage.libs.gap.util, the error handler sets an exception and signals sig_error(). An alternative approach would be an error handler that just sets an exception. If the meataxe functions are wrapped with an error value (such as cdef Matrix_t *MatInverse(Matrix_t M) except NULL) then this error would be raised, wouldn't it?

To me, it seems like a cleaner approach, as it allows meataxe to cope with the exception. Plus, we would see if meataxe really respects its own error values (otherwise we would see crashes in the wrapper) and could fix that. Meanwhile I have THREE patches that should be send upstream, and I wouldn't mind a fourth patch.

[jdemeyer]

Replying to SimonKing:

In sage.libs.gap.util, the error handler sets an exception and signals sig_error(). An alternative approach would be an error handler that just sets an exception. If the meataxe functions are wrapped with an error value (such as cdef Matrix_t *MatInverse(Matrix_t M) except NULL) then this error would be raised, wouldn't it?

At least in theory, that would work. You do need to be sure that the result is NULL if and only if the error handler was called. I agree that it would be cleaner if you can make it work because that would completely bypass the sig_on() system.

[SimonKing]

Replying to jdemeyer:

At least in theory, that would work. You do need to be sure that the result is NULL if and only if the error handler was called.

At least the meataxe documentation occasionally mentions that a certain value is returned on error.

And to be on the safe side, one could perhaps do cdef Matrix_t *MatInverse(Matrix_t M) except? NULL. Then, Cython would check if there is an exception set, and otherwise it would just continue.

[jdemeyer]

Replying to SimonKing:

otherwise it would just continue.

Would that make sense? If you're doing something with the result, you need that it's not NULL, so the except NULL can actually be an additional sanity check (IIRC, Cython raises a SystemError when there is no exception set if there should be one. That's safer that a segfault because you forgot to check for a NULL).

[SimonKing]

Replying to jdemeyer:

Replying to SimonKing:

otherwise it would just continue.

Would that make sense? If you're doing something with the result, you need that it's not NULL, so the except NULL can actually be an additional sanity check (IIRC, Cython raises a SystemError when there is no exception set if there should be one. That's safer that a segfault because you forgot to check for a NULL).

Correct. The edits I am currently doing generally assume an error if an error value occurs (accepting Cython's SystemError as an answer). There is one function where the docs say that the value -1 means that PERHAPS there was an exception, and I'll use except? only in those cases.

[SimonKing]

I was mistaken in one point: MtxErrorRecord_t does not contain the error code. It only contains a message that was created from the error record. The error messages are taken from a "static struct" that is part of the file message.c

Since the list of these error messages is finite and is in one-to-one correspondence to the error codes, I will translate each mesage back into the error code and use the appropriate error (ZeroDivisionError, MemoryError, ArithmeticError etc).

If you think that it isn't worth it, I could instead create a new error type MtxError, that is raised whenever an error occurs in meataxe. But I will go for the "long" solution, if you don't stop me...

[SimonKing]

It seems to work. With the commit that I am preparing now and that also contains some doctests, I obtain:

sage: M = MatrixSpace(GF(27,'x'),5).random_element(density=0.4)
sage: M.rank()
5
sage: ~M
[            1 x^2 + 2*x + 2             0             0             0]
[    x^2 + 2*x     2*x^2 + 2           x^2             0   x^2 + x + 2]
[      x^2 + 2           x^2         2*x^2     2*x^2 + 1       x^2 + 2]
[            1       x^2 + x           2*x             0             0]
[    2*x^2 + x             x             0             0             0]
sage: M = MatrixSpace(GF(27,'x'),5).random_element(density=0.4)
sage: M.rank()
4
sage: ~M
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-6-1888fd9efc66> in <module>()
----> 1 ~M

/home/king/Sage/git/sage/src/sage/matrix/matrix_gfpn_dense.pyx in sage.matrix.matrix_gfpn_dense.Matrix_gfpn_dense.__invert__ (build/cythonized/sage/matrix/matrix_gfpn_dense.c:12186)()
   1032         OUT._cache = {}
   1033         sig_on()
-> 1034         OUT.Data = MatInverse(self.Data)
   1035         sig_off()
   1036         if OUT.Data != NULL:

ZeroDivisionError: Division by zero in file matinv.c (line 50)

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​d889e0b	Improve echelon computation in MeatAxe, and fix some compiler warnings
​2e64257	Doctests and error handling for MeatAxe

[SimonKing]

With the new commits, there is full doctest coverage for the MeatAxe wrapper. Both

sage -t --optional=sage,meataxe src/sage/matrix/matrix_gfpn_dense.pyx

and

sage -t src/sage/matrix/matrix_gfpn_dense.pyx

pass.

Moreover, as I have announced above, I can now handle errors in MeatAxe.

Another change: I have removed the __pow__ method: The default __pow__ method automatically uses Strassen multiplication and is thus faster than the previous method, which called MeatAxe's MatPower (which doesn't use Strassen).

Probably I have to work on several doctests in other files: All tests depending on a random matrix may fail with meataxe, thus, for these cases I need to force usage of Matrix_generic_dense.

[SimonKing]

Problem:

sage: K.<a> = GF(25)
sage: MS = MatrixSpace(K, 2, 4)
sage: M = MS([4, 4, 1, 0, 0, 2*a+1, a+2, 1])
sage: M.echelon_form()
[      4       4       1       0]
[      0 2*a + 1   a + 2       1]

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​7c38969	Fix computation of row-reduced echelon form
​55a278d	Fix doctests when meataxe is installed

[SimonKing]

I have fixed the computation of row-reduced echelon form.

• With the current branch, meataxe in Sage is substantially improved compared with upstream (fixes for multiplication table IO, some speed-ups for echelon computation, Strassen-Winograd, error handling).
• Matrix arithmetic over small non-prime fields of odd order will be vastly faster with meataxe than with the previously used matrix implementation in Sage.
• Doctest coverage of the Cython wrapper is 100%.
• make test should now pass (I need to re-run "make test" after fixing the errors, though).

Hence: Needs review!

[SimonKing]

PS: I suggest that meataxe will be optional, not experimental.

[tscrim]

+1 for making it optional.

[SimonKing]

Sigh. I shot into my own foot: After installing this package, I couldn't install my cohomology package (because the latter relies on meataxe-2.2.4).

[SimonKing]

Motivated by the attempt to rebase my cohomology package on the new meataxe version, I try to be a bit more careful concerning errors in meataxe.

I made a thorough list of all functions in meataxe that use MTX_ERROR (or some variation thereof). I found that not all of them have error values set. Some of them have a documentation defining an error value, but the function in some cases does not return them. And if a function returning an error value is called by a second function, then the error value is not always dealt with.

It would be a very messy work to fix error handling thoroughly. I thus try to go a middle way:

• I will ignore MeatAxe's stand-alone programs, because they are *supposed* to immediately exit on error. In particular, I will also ignore functions checking command line arguments.
• I will try to only deal with functions that are relevant for matrix arithmetic, because that's the only part of meataxe that I am wrapping here.

Hence, if someone in future wants to wrap the more interesting parts of MeatAxe (all the representation theory stuff), then he or she has to deal with it.

[SimonKing]

PS: Also ignore the test suite, as that is supposed to have a non-default error handler anyway. And as long as the test suite passes, I could care less about its error handling.

[SimonKing]

To summarise, it seems to me that everything concerning

• stand-alone programs
• test suite
• "bit strings"
• polynomials and their factorisation (I currently do not wrap characteristic polynomials!)
• "greased" matrices (but see below!)
• integer matrices (I am only interested in finite fields)
• representations
• lattices
• matrix sets
• sets
• permutations
• operating system stuff. I am not going to fix setting a timeout in meataxe, because we can set a timeout in Sage, with alarm(...), provided that we wrap meataxe calls in sig_on()/sig_off().
• word generators

is not relevant at that point. I'll focus on what remains.

Greased matrices

Something that I have ignored but perhaps should consider in future: "Greased matrices". From the documentation, it seems that greased matrices use more memory, but have faster row operations achieved by precomputing linear combinations of blocks of rows. So, if meataxe becomes an optional package with Matrix_gfpn_dense wrapping meataxe matrices, then a further development would be a wrapper for greased matrices.

[SimonKing]

Also I will not do error propagation for MatFree.

[SimonKing]

I have a patch that "should work" (TM). But the test suite consistently fails in a part that I didn't change (concerning the stand-alone program zmw, which is supposed to "Make word"). The program relies on computation of null space which (I thinK) I haven't changed, except for error handling. Sigh.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​f733377

Use and propagate specific return values on error in matrix-related MeatAxe functions.

[SimonKing]

The problems with error propagation are solved, tests with the current patch should still pass, and I think it is good that with the fifth patch MeatAxe becomes a lot more consequent in using specific return values indicating an error. Thus, it an error handler that doesn't exit on error became even more useful.

For now, I have no idea for further improvements (but perhaps you have ideas?). So, I will continue to try and rebase David Green's programs for my group cohomology package.

[SimonKing]

Replying to SimonKing:

For now, I have no idea for further improvements ...

One detail comes to my mind: MeatAxe does not only have MTX_ERROR, but also MTX_ASSERT and MTX_VERIFY. I made sure that after all instances of MTX_ERROR outside of stand-alone executables an error return value is returned. This is yet to be done for MTX_ASSERT and MTX_VERIFY.

[SimonKing]

OTOH, I feel that I am overdoing it. So, I will *not* deal with MTX_ASSERT and MTX_VERIFY, unless requested by the referee.

[jdemeyer]

Replying to SimonKing:

Complain to upstream. The tarball is directly from their website. I'm not going to fix that.

Upstream's answer:

I am currently working on a 2.4 release that will solve some known 
problems and hopefully not break anthing else (it's been a long
time since I last touched the code...). The next TAR release will
contain a properly named root directory as you suggested. 

I am not sure, however, when the release will be ready.

[SimonKing]

Replying to jdemeyer:

Replying to SimonKing:

Complain to upstream. The tarball is directly from their website. I'm not going to fix that.

Upstream's answer:

I am currently working on a 2.4 release that will solve some known 
problems and hopefully not break anthing else (it's been a long
time since I last touched the code...). The next TAR release will
contain a properly named root directory as you suggested. 

I am not sure, however, when the release will be ready.

OMG. I see coming that my patches won't apply on the new release.

[SimonKing]

I have sent my patches upstream a minute ago.

[jdemeyer]

Any answer from upstream?

Do you have an idea who could review those patches, because I seems to me that it should only be reviewed by somebody who knows meataxe.

One detail I spotted in the code: replace except: by except BaseException: (I assume that people write except: by mistake, so it's better to be explicit that you really want to catch everything).

[SimonKing]

Replying to jdemeyer:

Any answer from upstream?

Unfortunately not.

Do you have an idea who could review those patches, because I seems to me that it should only be reviewed by somebody who knows meataxe.

This might be difficult. But the following should be possible even without deeper knowledge of meataxe:

• Does the package keeps the promise that multiplication tables (name pXYZ.zzz, where XYZ is the order of the field) will not pollute the current directory?
• Are there memory leaks? Combine various random matrices by stacking, slicing, multiplication, copying etc, delete everything, repeat in a loop, and see if any memory is lost.
• Manually create the same large random matrices once in meataxe and once in the current sage backends. That includes the case of fields of characteristic 2 and prime fields, where meataxe will not be used by default. Then do arithmetic, including echolonization, and compare (using M.list()) if the results coincide in all implementations. There are examples of that type in the docs.
• In addition to the previous class of tests, compare the timings. You will hopefully find that for non-prime fields of odd characteristic meataxe is clearly superior over the current Sage implementation.
• Concerning Strassen-Winograd (which is not in upstream meataxe): Compare the results of multiplication using school book and Strassen-Winograd algorithm, with different cutoff and different matrix dimensions. Here, both the results (via M.list()) and the performance (memory consumption) should be compared. Also it can be tested if the default cutoff I chose is optimal even on 32bit machines (I only tested on 64 bit).

I think a review based on such test would be valid. I doubt that for our current standard backends the review was more thorough than what I suggest here for an optional backend.

One detail I spotted in the code: replace except: by except BaseException: (I assume that people write except: by mistake, so it's better to be explicit that you really want to catch everything).

OK, can do.

[jdemeyer]

Replying to SimonKing:

I doubt that for our current standard backends the review was more thorough than what I suggest here for an optional backend.

A fair comment, but those standard backends are not heavily patched with patches which are not accepted by upstream.

Upstream code does not need a review, but these patches do need a review.

[jdemeyer]

I think you are still mis-using some of the error handling.

For example, you have

int MatEchelonize(Matrix_t *mat) except -1

but you still write

        if MatEchelonize(self.Data) == -1:
            raise ArithmeticError("Error echelonizing this matrix")

The latter should just be

MatEchelonize(self.Data)

[jdemeyer]

You can remove these parts from spkg-install:

+# Just to be sure, we also create other folders, although
+# they are standard SageMath folders
+
+mkdir -p $MTXBIN
+
+if [ $? -ne 0 ]; then
+    echo >&2 "Error creating directory for meataxe binaries."
+    exit 1
+fi
+
+mkdir -p "$SAGE_LOCAL/include"
+
+if [ $? -ne 0 ]; then
+    echo >&2 "Error creating SageMath's include directory."
+    exit 1
+fi
+
+mkdir -p "$SAGE_LOCAL/lib"
+
+if [ $? -ne 0 ]; then
+    echo >&2 "Error creating SageMath's lib folder."
+    exit 1
+fi

(I know I asked you to add these, but in the mean time the build system was improved)

[SimonKing]

Replying to jdemeyer:

I think you are still mis-using some of the error handling.

For example, you have

int MatEchelonize(Matrix_t *mat) except -1

but you still write

        if MatEchelonize(self.Data) == -1:
            raise ArithmeticError("Error echelonizing this matrix")

The latter should just be

MatEchelonize(self.Data)

Right. I guess I wrote it before I changed to using error return values. How important is it to change that?

[jdemeyer]

Replying to SimonKing:

How important is it to change that?

That's a very philosophical question. It's not important. This ticket is not important. Even Sage itself is not really that important.

Seriously: it's a small thing which makes your code more readable, so why wouldn't you fix it?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​710668d

Remove overcautious commands in spkg-install; rely on default error return values in matrix_gfpn_dense

[SimonKing]

In fact there were several places in which I explicitly tested for error return values, where I should better let Cython do the job. Fixed now. Of course, I can't easily solve the problem that upstream didn't comment on (or just acknowledges receive of) my patches.

[SimonKing]

I thought that I did test pickling. But now I obtain crashes when unpickling a pickle.

sage: M = MatrixSpace(GF(9,'x'),10,10)(0)
sage: M[2,3] = M[3,7] = 1
sage: f = tmp_filename()
sage: save(M, f)
sage: load(f+".sobj")
------------------------------------------------------------------------
/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/ext/interrupt/interrupt.so(+0x4125)[0x7f5a7dd08125]
/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/ext/interrupt/interrupt.so(+0x4177)[0x7f5a7dd08177]
/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/ext/interrupt/interrupt.so(+0x60c9)[0x7f5a7dd0a0c9]
/lib64/libpthread.so.0(+0xf890)[0x7f5a85922890]
/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/matrix/matrix_gfpn_dense.so(+0xd9a0)[0x7f5a46fad9a0]
/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/matrix/matrix_dense.so(+0x629a)[0x7f5a6290029a]
/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/matrix/matrix0.so(+0x39496)[0x7f5a62d99496]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyObject_Call+0x43)[0x7f5a85b84a23]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_CallObjectWithKeywords+0x47)[0x7f5a85c35317]
/home/king/Sage/git/sage/local/lib/python2.7/lib-dynload/cPickle.so(+0x4f4e)[0x7f5a80250f4e]
/home/king/Sage/git/sage/local/lib/python2.7/lib-dynload/cPickle.so(+0xa669)[0x7f5a80256669]
/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/structure/sage_object.so(+0x21fa2)[0x7f5a7e83bfa2]
/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/structure/sage_object.so(+0x77df)[0x7f5a7e8217df]
/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/structure/sage_object.so(+0x25e10)[0x7f5a7e83fe10]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x48a9)[0x7f5a85c3a1a9]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x80d)[0x7f5a85c3b9cd]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCode+0x32)[0x7f5a85c3bb02]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x4dae)[0x7f5a85c3a6ae]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x80d)[0x7f5a85c3b9cd]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x43a1)[0x7f5a85c39ca1]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x80d)[0x7f5a85c3b9cd]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x43a1)[0x7f5a85c39ca1]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x80d)[0x7f5a85c3b9cd]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x43a1)[0x7f5a85c39ca1]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x80d)[0x7f5a85c3b9cd]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x43a1)[0x7f5a85c39ca1]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x80d)[0x7f5a85c3b9cd]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x43a1)[0x7f5a85c39ca1]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x80d)[0x7f5a85c3b9cd]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x43a1)[0x7f5a85c39ca1]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x80d)[0x7f5a85c3b9cd]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCode+0x32)[0x7f5a85c3bb02]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyRun_FileExFlags+0x92)[0x7f5a85c666d2]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyRun_SimpleFileExFlags+0xd9)[0x7f5a85c67c09]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(Py_Main+0xc4d)[0x7f5a85c7d99d]
/lib64/libc.so.6(__libc_start_main+0xf5)[0x7f5a84e85b05]
python[0x4007be]
------------------------------------------------------------------------
Attaching gdb to process id 4561.
...

[SimonKing]

What can I change in my local git configuration so that I can push to trac if port 22 is blocked?

[tscrim]

You could try the solution posted here: ​http://stackoverflow.com/questions/4891527/git-protocol-blocked-by-company-how-can-i-get-around-that, but I'm not sure that will work because it is a different port (and IDK about how the trac server is setup and which ports it is listening to).

[SimonKing]

Replying to tscrim:

You could try the solution posted here: ​http://stackoverflow.com/questions/4891527/git-protocol-blocked-by-company-how-can-i-get-around-that, but I'm not sure that will work because it is a different port (and IDK about how the trac server is setup and which ports it is listening to).

It didn't work; in particular, after the change, the error messages still mentioned port 22, which is not available to me.

[vbraun]

not access to port 22 == no access to the internet. Complain to whoever is blocking you that you can't do your work and their service is useless to you

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​c67cb42

Fix pickling of meataxe matrices

[SimonKing]

The problem was: The default implementation of unpickling relies on set_unsafe(i,j,val). In my implementation, set_unsafe relies on a converter from finite field elements to the internal representation---and the converter is not created during the unpickling.

Hence, I created a custom __reduce__ relying on some unpickling function.

There, I used the occasion to return to what I did in the old group cohomology spkg: Pickling is *not* by conversion of the matrix to a list of entries, but it dumps the densely packed internal representation. Advantage: It is a lot faster!

For the group cohomology spkg, it was checked that the internal representation does not depend on details such as big versus little endian machines: A pickle created on a big-endian machine can correctly be unpickled with a little-endian machine.

Another advantage of the new approach will be apparent in #18514: Since the new pickling/unpickling scheme for meataxe matrices is a lot closer to the scheme used in the old group cohomology spkg, the *new* group cohomology package should relatively easily be able to read data from existing cohomology data bases.

[tscrim]

I get 2 doctest failures with MeatAxe not installed. Both I think are because they should be marked optional:

File "../../matrix/matrix_gfpn_dense.pyx", line 213, in sage.matrix.matrix_gfpn_dense.PrimeFieldConverter_class.int_to_field
Failed example:
    C.int_to_field(int(2))
Exception raised:
    Traceback (most recent call last):
      File "/home/travis/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 496, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/travis/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 858, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.matrix.matrix_gfpn_dense.PrimeFieldConverter_class.int_to_field[1]>", line 1, in <module>
        C.int_to_field(int(Integer(2)))
    NameError: name 'C' is not defined

**********************************************************************

File "../../matrix/matrix_space.py", line 992, in sage.matrix.matrix_space.MatrixSpace._get_matrix_class
Failed example:
    type(matrix(GF(125,'z'), 2, range(4)))
Expected:
    <type 'sage.matrix.matrix_gfpn_dense.Matrix_gfpn_dense'>
Got:
    <type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'>

[tscrim]

However, all tests pass in src/matrix with it installed.

[git]

Branch pushed to git repo; I updated commit sha1. Last 10 new commits:

​c2e6fe5	A full wrapper for MeatAxe matrices
​395aa9a	Improve echelon computation in MeatAxe, and fix some compiler warnings
​c5f328c	Doctests and error handling for MeatAxe
​c29d1f8	Fix computation of row-reduced echelon form
​f816e41	Fix doctests when meataxe is installed
​80af75e	Use and propagate specific return values on error in matrix-related MeatAxe functions.
​6649c82	Remove overcautious commands in spkg-install; rely on default error return values in matrix_gfpn_dense
​d24260c	Fix pickling of meataxe matrices
​bddc09d	Merge branch 'u/SimonKing/meataxe' of trac.sagemath.org:sage into t/12103/meataxe
​74edf19	Make two doctests optional

[SimonKing]

I made the two failing tests optional and at that occasion merged the latest develop.

Is it fine now?

[tscrim]

Looks good to me now. Thanks.

[dimpase]

there is no tarball on the mirrors, apparently

[vbraun]

Fixed