#12103 closed defect (fixed)
Use MeatAxe as an optional back end for dense matrices over `GF(p^n)`, p odd, n>1, `p^n<255`
Reported by:  SimonKing  Owned by:  jason, was 

Priority:  major  Milestone:  sage7.0 
Component:  packages: optional  Keywords:  linear algebra, MeatAxe 
Cc:  malb, jdemeyer, vbraun  Merged in:  
Authors:  Simon King  Reviewers:  Jeroen Demeyer, Travis Scrimshaw 
Report Upstream:  Reported upstream. No feedback yet.  Work issues:  
Branch:  74edf19 (Commits, GitHub, GitLab)  Commit:  
Dependencies:  #19240  Stopgaps: 
Description (last modified by )
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.rwthaachen.de/~MTX/meataxe2.4.24.tar.gz needed to be repackaged, in order to unpack into a single folder. Use attachment:meataxe2.4.24.tar.gz
What is done
There is no spkgcheck. 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 StrassenWinograd 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
Attachments (4)
Change History (177)
comment:1 Changed 9 years ago by
 Description modified (diff)
 Status changed from new to needs_review
Changed 9 years ago by
Use MeatAxe
as back end for linear algebra over some finite fields. Relative to #11900
comment:2 Changed 9 years ago by
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]
comment:3 Changed 9 years ago by
 Description modified (diff)
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
comment:4 followup: ↓ 6 Changed 9 years ago by
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.
comment:5 Changed 9 years ago by
comment:6 in reply to: ↑ 4 Changed 9 years ago by
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/sage5.0.beta6
So, I can not confirm that the first tobeapplied 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.
comment:7 Changed 9 years ago by
comment:8 followup: ↓ 9 Changed 9 years ago by
Ah, crap, 11900 was it.
RE: module_list.py I commented on it here: https://bitbucket.org/malb/sagemath/pullrequest/1/12103meataxewrapper#comment3510
PS: Did you get my email that I sent to simon.king_unijena.de?
comment:9 in reply to: ↑ 8 Changed 9 years ago by
 Description modified (diff)
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/pullrequest/1/12103meataxewrapper#comment3510
PS: Did you get my email that I sent to simon.king_unijena.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.
comment:10 Changed 8 years ago by
 Milestone changed from sage5.11 to sage5.12
comment:11 Changed 7 years ago by
 Milestone changed from sage6.1 to sage6.2
comment:12 Changed 7 years ago by
 Status changed from needs_review to needs_work
comment:13 Changed 7 years ago by
 Milestone changed from sage6.2 to sage6.3
comment:14 Changed 7 years ago by
 Milestone changed from sage6.3 to sage6.4
comment:15 Changed 6 years ago by
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 newstyle package.
comment:16 Changed 6 years ago by
 Description modified (diff)
comment:17 Changed 6 years ago by
 Cc jdemeyer vbraun added; mringe@… david.green@… removed
The first step will probably be the creation of a newstyle package for MeatAxe 2.4.24 (which I take as an exercise for learning newstyle 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?
 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
 b) Should we keep Makefile, create an empty Makefile.conf, and do the rest with environment variables?
or
 Should spkginstall first create a fully grown Makefile.conf?
comment:18 Changed 6 years ago by
If all you have to do is touch Makefile.conf
then I'd go for that.
comment:19 Changed 6 years ago by
 Component changed from linear algebra to packages: experimental
comment:20 Changed 6 years ago by
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?
comment:21 Changed 6 years ago by
Unless you export ZZZ
it is not passed to subprocesses. Alternatively, you might want to use $MAKE ZZZ=0
.
comment:22 followups: ↓ 23 ↓ 25 ↓ 26 Changed 6 years ago by
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 spkginstall 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 callsmake clean
therebut 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 spkginstall ALWAYS do make check
? Or should that be solely the job of spkgcheck (with the drawback that then spkgcheck would rebuild meataxe)?
comment:23 in reply to: ↑ 22 Changed 6 years ago by
Replying to SimonKing:
Question on
make check
: Would it be OK to let spkginstall ALWAYS domake check
? Or should that be solely the job of spkgcheck (with the drawback that then spkgcheck would rebuild meataxe)?
Or should spkginstall test whether SAGE_CHECK=yes, and choose between make
and make check
accordingly?
comment:24 Changed 6 years ago by
Do I understand correctly that for automatically downloading the sources from upstream, I should create spkgsrc similarly to what I find in build/pkg/singular/spkgsrc?
comment:25 in reply to: ↑ 22 ; followup: ↓ 32 Changed 6 years ago by
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 csources so that MeatAxe would be able to find the tables in a fixed directory. Now I cannot see any of these tables.
comment:26 in reply to: ↑ 22 Changed 6 years ago by
Replying to SimonKing:
Question on
make check
: Would it be OK to let spkginstall ALWAYS domake check
? Or should that be solely the job of spkgcheck (with the drawback that then spkgcheck 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 spkgcheck
.
comment:27 Changed 6 years ago by
You don't need spkgsrc, its just a bandaid to modify tarballs (don't do it if you can avoid it)
comment:28 Changed 6 years ago by
 Branch set to u/SimonKing/meataxe
comment:29 Changed 6 years ago by
 Commit set to e732f0ee8fa0dc631f9cffd62b4dd1f6f2c10c2d
 Description modified (diff)
New commits:
e732f0e  An optional MeatAxe package

comment:30 Changed 6 years ago by
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?
comment:31 Changed 6 years ago by
PS: Note that MeatAxe has a nice documentation, and as usual it is installed if SAGE_SPKG_INSTALL_DOCS=yes.
comment:32 in reply to: ↑ 25 ; followup: ↓ 33 Changed 6 years ago by
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 csources 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/?
comment:33 in reply to: ↑ 32 Changed 6 years ago by
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?
comment:34 Changed 6 years ago by
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.
comment:35 Changed 6 years ago by
PS: It says "Otherwise, the value of the environment variable MTXLIB is used." That definitely does not work.
comment:36 Changed 6 years ago by
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.
comment:37 Changed 6 years ago by
Internally, MeatAxe stores MTXLIB in MtxLibDir?. But directly defining MtxLibDir? at run time didn't help. Sigh.
comment:38 Changed 6 years ago by
 Commit changed from e732f0ee8fa0dc631f9cffd62b4dd1f6f2c10c2d to 6c324700abeb37c3e32ebf36c7cf9027c8e70588
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

comment:39 followup: ↓ 41 Changed 6 years ago by
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/meataxe2.4.24/
.
comment:40 Changed 6 years ago by
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 IOfunctions that accept an argument saying that the file is a library filebut in the IOfunctions, 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 IOargument for "library files" is used wrongly. So, after fixing the IOfunction, 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.
comment:41 in reply to: ↑ 39 ; followup: ↓ 43 Changed 6 years ago by
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 ofDOT_SAGE
, probablyDOT_SAGE/meataxe2.4.24/
.
Right, that could be a problem. The tables should indeed be permanent. So, DOT_SAGE/meataxe
(without version number) it is.
comment:42 Changed 6 years ago by
 Commit changed from 6c324700abeb37c3e32ebf36c7cf9027c8e70588 to d9675fbca6cf8647548f6f37522d2c388812669d
Branch pushed to git repo; I updated commit sha1. New commits:
d9675fb  Use DOT_SAGE, not SAGE_LOCAL, for storing multiplication tables

comment:43 in reply to: ↑ 41 ; followup: ↓ 44 Changed 6 years ago by
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 ofDOT_SAGE
, probablyDOT_SAGE/meataxe2.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.)
comment:44 in reply to: ↑ 43 ; followup: ↓ 45 Changed 6 years ago by
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
andmatplotlib
, which is why their directories inDOT_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.
comment:45 in reply to: ↑ 44 Changed 6 years ago by
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.
comment:46 Changed 6 years ago by
 Description modified (diff)
comment:47 Changed 6 years ago by
 Commit changed from d9675fbca6cf8647548f6f37522d2c388812669d to 20a16e18db945d6e123f457a44ee5835427a1351
Branch pushed to git repo; I updated commit sha1. New commits:
20a16e1  Implement and use StrassenWinograd matrix multiplication in MeatAxe

comment:48 Changed 6 years ago by
In the recent commit, patches StrassenWinogradImplementation.patch
and StrassenWinogradUsage.patch
are provided.
There, I implement the WinogradStrassen 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 nonprime finite fields of odd characteristic (of course with field order < 256).
Moreover, it turns out that the WinogradStrassen 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 WinogradStrassen 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.
comment:49 followup: ↓ 51 Changed 6 years ago by
The tarball is badly made: there should be toplevel directory like meataxe2.4.24
with everything contained in that directory.
Also, we no longer list maintainers in SPKG.txt
(#19238)
comment:50 followup: ↓ 52 Changed 6 years ago by
Patches need to be documented (preferably in the patch itself). In particular: where do the patches come from?
comment:51 in reply to: ↑ 49 ; followups: ↓ 53 ↓ 146 Changed 6 years ago by
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.
comment:52 in reply to: ↑ 50 ; followup: ↓ 56 Changed 6 years ago by
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 IOfixes, but failed for the other two.
In particular: where do the patches come from?
My name suffices?
comment:53 in reply to: ↑ 51 ; followup: ↓ 55 Changed 6 years ago by
comment:54 followup: ↓ 57 Changed 6 years ago by
I guess nobody except myself uses the branch. Hence, I will rebase and forcepush the changes requested by Jeroen, unless someone tells me not to.
comment:55 in reply to: ↑ 53 Changed 6 years ago by
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.
comment:56 in reply to: ↑ 52 ; followup: ↓ 64 Changed 6 years ago by
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!
comment:57 in reply to: ↑ 54 ; followup: ↓ 62 Changed 6 years ago by
Replying to SimonKing:
Hence, I will rebase and forcepush 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
)
comment:58 Changed 6 years ago by
After replacing the old tarball by a new one called meataxe2.4.24.p0.tar.gz
, I did ./sage fixpkgchecksums
. However, it did not produce a new checksum.ini
in build/pkgs/meataxe/
.
What can I do now?
comment:59 Changed 6 years ago by
Aha. Apparently it was needed to call it meataxe2.4.24.tar.gz
without .p0.
comment:60 Changed 6 years ago by
 Commit changed from 20a16e18db945d6e123f457a44ee5835427a1351 to a722d4e9b8208f900c7850e8ddf8adda2e297550
comment:61 Changed 6 years ago by
 Description modified (diff)
comment:62 in reply to: ↑ 57 Changed 6 years ago by
Replying to jdemeyer:
Go ahead. If you're rebasing, you might as well squash everything in one commit (with the
f
option ingit 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 WinogradStrassen multiplication in MeatAxe.
The selftests (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 selftests 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.
comment:63 Changed 6 years ago by
 Dependencies #9562 #4260 deleted
comment:64 in reply to: ↑ 56 ; followup: ↓ 66 Changed 6 years ago by
 Dependencies set to #9562 #4260
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 WinogradStrassen 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.
comment:65 Changed 6 years ago by
 Commit changed from a722d4e9b8208f900c7850e8ddf8adda2e297550 to 539b1349359004b78712a67d2739cce9ce0227da
Branch pushed to git repo; I updated commit sha1. New commits:
539b134  A very basic MeatAxe Cython wrapper

comment:66 in reply to: ↑ 64 ; followup: ↓ 69 Changed 6 years ago by
Replying to SimonKing:
Long time ago, I have submitted a WinogradStrassen 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???
comment:67 Changed 6 years ago by
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:meataxe2.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 nonprime 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.
comment:68 followup: ↓ 71 Changed 6 years ago by
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#headingsofsagelibrarycodefiles
comment:69 in reply to: ↑ 66 ; followup: ↓ 70 Changed 6 years ago by
Replying to jdemeyer:
Replying to SimonKing:
Long time ago, I have submitted a WinogradStrassen 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.
comment:70 in reply to: ↑ 69 ; followup: ↓ 73 Changed 6 years ago by
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?
comment:71 in reply to: ↑ 68 ; followup: ↓ 74 Changed 6 years ago by
Replying to jdemeyer:
src/sage/matrix/matrix_modpn_dense.pxd
should probably be moved tosrc/sage/libs/meataxe.pxd
orsrc/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#headingsofsagelibrarycodefiles
I think I have copied the heading from some existing pxd files. But I can change it.
comment:72 followup: ↓ 76 Changed 6 years ago by
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?
comment:73 in reply to: ↑ 70 Changed 6 years ago by
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."
comment:74 in reply to: ↑ 71 ; followup: ↓ 77 Changed 6 years ago by
Replying to SimonKing:
This is totally analogous to
matrix_mod2e_dense.pxd
; or do you suggest to rename the latter bym4rie.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.
comment:75 Changed 6 years ago by
 Dependencies #9562 #4260 deleted
comment:76 in reply to: ↑ 72 ; followup: ↓ 78 Changed 6 years ago by
 Dependencies set to #9562 #4260
Replying to jdemeyer:
If I understand things correctly, the
modpn
in the filename refers to a finite field of sizep^n
, right? That's a very confusing name, can you usematrix_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
.
comment:77 in reply to: ↑ 74 ; followup: ↓ 79 Changed 6 years ago by
Replying to jdemeyer:
Replying to SimonKing:
This is totally analogous to
matrix_mod2e_dense.pxd
; or do you suggest to rename the latter bym4rie.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
.
comment:78 in reply to: ↑ 76 ; followup: ↓ 83 Changed 6 years ago by
Replying to SimonKing:
Do you suggest to rename Martin Albrecht's
matrix_mod2e_dense.pyx
tomatrix_gf2e_dense.pyx
?
Yes!! It's a very confusing name.
comment:79 in reply to: ↑ 77 Changed 6 years ago by
 Dependencies #9562 #4260 deleted
Replying to SimonKing:
Have a
src/sage/libs/meataxe.pxd
without a corresponding.pyx
file, and do the rest insrc/sage/matrix/matrix_modpn_dense.p??
, which imports fromsage.libs.meataxe
.
Yes, exactly. It's good practice to separate the library (which has nothing to do with Sage) from the Sage interface.
comment:80 Changed 6 years ago by
 Commit changed from 539b1349359004b78712a67d2739cce9ce0227da to a8f18474d1b95afb024a53c36ca59b311519e88d
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
a8f1847  A very basic MeatAxe Cython wrapper

comment:81 Changed 6 years ago by
 Commit changed from a8f18474d1b95afb024a53c36ca59b311519e88d to da665d9f683dc4ac77df6270e18f8f494d9118ac
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
da665d9  A very basic MeatAxe Cython wrapper

comment:82 Changed 6 years ago by
 Commit changed from da665d9f683dc4ac77df6270e18f8f494d9118ac to 106aeef1f449f7db3a617e986d2e697865ee34ac
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
106aeef  A very basic MeatAxe Cython wrapper

comment:83 in reply to: ↑ 78 Changed 6 years ago by
Replying to jdemeyer:
Replying to SimonKing:
Do you suggest to rename Martin Albrecht's
matrix_mod2e_dense.pyx
tomatrix_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 WinogradStrassen multiplication.
comment:84 Changed 6 years ago by
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.
comment:85 followup: ↓ 87 Changed 6 years ago by
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
)
comment:86 followup: ↓ 88 Changed 6 years ago by
In spkginstall
, 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
comment:87 in reply to: ↑ 85 ; followup: ↓ 90 Changed 6 years ago by
Replying to jdemeyer:
In
meataxe.pxd
, I assume thatctypedef struct FILE
refers to the<stdio.h>
typeFILE
? In that case, it's better to usefrom 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.
comment:88 in reply to: ↑ 86 Changed 6 years ago by
Replying to jdemeyer:
In
spkginstall
, why this?PATCHES=../patches for patch in "$PATCHES"/*.patch; doIt'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.
comment:89 followup: ↓ 93 Changed 6 years ago by
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.
comment:90 in reply to: ↑ 87 ; followup: ↓ 91 Changed 6 years ago by
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 meataxespecific.
comment:91 in reply to: ↑ 90 ; followup: ↓ 92 Changed 6 years ago by
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.
comment:92 in reply to: ↑ 91 Changed 6 years ago by
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
.
comment:93 in reply to: ↑ 89 Changed 6 years ago by
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 functionFfMapRow
.
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.
comment:94 Changed 6 years ago by
 Commit changed from 106aeef1f449f7db3a617e986d2e697865ee34ac to c1d5026ade5e1bb454401a34a734c91c15cbf75f
comment:95 Changed 6 years ago by
I forcepushed 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 WinogradStrassen, one needs to define a matrix size such that divideandconquer 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").
WinogradStrassen 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 WinogradStrassen 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 WinogradStrassen 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 WinogradStrassen 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 WinogradStrassen > WinogradStrassen 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 3601x3601matrices 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 32bit machine, cutoff 4 is better than cutoff "2==sizeof(long)/2", and on a 128bit 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 (32bit, 128bit, 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 kernel0.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.
comment:96 Changed 6 years ago by
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/linuxva3e.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/linuxva3e.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
.
comment:97 Changed 6 years ago by
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.
comment:98 Changed 6 years ago by
It would be convenient if you would base this ticket on #19240 (and review that), because there will be conflicts.
comment:99 Changed 6 years ago by
 Dependencies set to #19240
comment:100 Changed 6 years ago by
 Commit changed from c1d5026ade5e1bb454401a34a734c91c15cbf75f to 9e2a8c6027a1034ec21ea52dd2881965c05292d3
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 StrassenWinograd matrix multiplication in MeatAxe

9e2a8c6  A very basic MeatAxe Cython wrapper

comment:101 Changed 6 years ago by
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.
comment:102 Changed 6 years ago by
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.
comment:103 Changed 6 years ago by
In order to have a fullyfledged matrix implementation, I will have a look at sage.matrix.matrix_generic_dense to see what methods need to be implemented.
comment:104 Changed 6 years ago by
 Commit changed from 9e2a8c6027a1034ec21ea52dd2881965c05292d3 to 4bdd285cd6320851aa4ed7aee8fb820960005692
Branch pushed to git repo; I updated commit sha1. New commits:
4bdd285  A full wrapper for MeatAxe matrices

comment:105 Changed 6 years ago by
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 StrassenWinograd is used for all matrices).
Also, I have wrapped the computation of echelon form. Meataxe only computes a semiechelon 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.
comment:106 Changed 6 years ago by
 Commit changed from 4bdd285cd6320851aa4ed7aee8fb820960005692 to b237c5f4a500f313e573728d4a9b55cca902a255
Branch pushed to git repo; I updated commit sha1. New commits:
b237c5f  Improve echelon computation in MeatAxe

comment:107 Changed 6 years ago by
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.
comment:108 followup: ↓ 109 Changed 6 years ago by
Following up on your sagedevel 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/sageconfig/local/include/': Not a directory Error copying MeatAxe header.
comment:109 in reply to: ↑ 108 ; followups: ↓ 110 ↓ 111 Changed 6 years ago by
Replying to jdemeyer:
Following up on your sagedevel 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 avoidcp: cannot create regular file '/usr/local/src/sageconfig/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?
comment:110 in reply to: ↑ 109 Changed 6 years ago by
Replying to SimonKing:
Replying to jdemeyer:
Following up on your sagedevel post, you should add a file
build/pkgs/meataxe/dependencies
containing# no dependenciesThank you! Doing it now.
Didn't help.
comment:111 in reply to: ↑ 109 Changed 6 years ago by
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.
comment:112 Changed 6 years ago by
 Commit changed from b237c5f4a500f313e573728d4a9b55cca902a255 to db6b7fe1f49c5412a5331b7eaa3c9584fcf0d636
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

comment:113 Changed 6 years ago by
 Commit changed from db6b7fe1f49c5412a5331b7eaa3c9584fcf0d636 to 191477e697d5fb02c0e6bf7f8b80850e1092d4f6
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

comment:114 Changed 6 years ago by
Sorry for forcepushing 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).
comment:115 Changed 6 years ago by
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.
comment:116 Changed 6 years ago by
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
.
comment:117 followup: ↓ 118 Changed 6 years ago by
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
comment:118 in reply to: ↑ 117 Changed 6 years ago by
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 insrc/sage/libs/gap/util.pyx
Thank you, I'll have a look!
comment:119 Changed 6 years ago by
In case you're wondering: the function PyErr_SetObject()
is a lowlevel version of raise
: it "raises" an exception but execution continues to the sig_error()
which is handled by sig_on()
.
comment:120 followup: ↓ 121 Changed 6 years ago by
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 cleanup operations (causing a memory leak on error) if we would raise sig_error()?
comment:121 in reply to: ↑ 120 ; followup: ↓ 122 Changed 6 years ago by
Replying to SimonKing:
Wouldn't we prevent meataxe from these cleanup 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.
comment:122 in reply to: ↑ 121 ; followup: ↓ 123 Changed 6 years ago by
Replying to jdemeyer:
Replying to SimonKing:
Wouldn't we prevent meataxe from these cleanup 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 cleaninguponerror, 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 tobeused. 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.
comment:123 in reply to: ↑ 122 ; followup: ↓ 124 Changed 6 years ago by
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.
comment:124 in reply to: ↑ 123 ; followup: ↓ 125 Changed 6 years ago by
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.
comment:125 in reply to: ↑ 124 ; followup: ↓ 126 Changed 6 years ago by
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
).
comment:126 in reply to: ↑ 125 Changed 6 years ago by
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 theexcept NULL
can actually be an additional sanity check (IIRC, Cython raises aSystemError
when there is no exception set if there should be one. That's safer that a segfault because you forgot to check for aNULL
).
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.
comment:127 Changed 6 years ago by
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 onetoone 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...
comment:128 Changed 6 years ago by
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) <ipythoninput61888fd9efc66> 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)
comment:129 Changed 6 years ago by
 Commit changed from 191477e697d5fb02c0e6bf7f8b80850e1092d4f6 to 2e6425793607152a39296a95c5c50f63cf796dea
comment:130 Changed 6 years ago by
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
.
comment:131 Changed 6 years ago by
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]
comment:132 Changed 6 years ago by
 Commit changed from 2e6425793607152a39296a95c5c50f63cf796dea to 55a278da06ba77fdfde839aa2e45d43a6806f2fb
comment:133 Changed 6 years ago by
 Component changed from packages: experimental to packages: optional
 Status changed from needs_work to needs_review
I have fixed the computation of rowreduced echelon form.
 With the current branch, meataxe in Sage is substantially improved compared with upstream (fixes for multiplication table IO, some speedups for echelon computation, StrassenWinograd, error handling).
 Matrix arithmetic over small nonprime 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 rerun "make test" after fixing the errors, though).
Hence: Needs review!
comment:134 Changed 6 years ago by
PS: I suggest that meataxe will be optional, not experimental.
comment:135 Changed 6 years ago by
 Milestone changed from sage6.4 to sage6.9
+1 for making it optional.
comment:136 Changed 6 years ago by
Sigh. I shot into my own foot: After installing this package, I couldn't install my cohomology package (because the latter relies on meataxe2.2.4).
comment:137 Changed 6 years ago by
 Status changed from needs_review to needs_work
 Work issues set to improve error handling in meataxe (within reason)
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 standalone 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.
comment:138 Changed 6 years ago by
PS: Also ignore the test suite, as that is supposed to have a nondefault error handler anyway. And as long as the test suite passes, I could care less about its error handling.
comment:139 Changed 6 years ago by
To summarise, it seems to me that everything concerning
 standalone 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.
comment:140 Changed 6 years ago by
Also I will not do error propagation for MatFree
.
comment:141 Changed 6 years ago by
I have a patch that "should work" (TM). But the test suite consistently fails in a part that I didn't change (concerning the standalone 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.
comment:142 Changed 6 years ago by
 Commit changed from 55a278da06ba77fdfde839aa2e45d43a6806f2fb to f73337711df114571d1be0fa61140a41628703b7
Branch pushed to git repo; I updated commit sha1. New commits:
f733377  Use and propagate specific return values on error in matrixrelated MeatAxe functions.

comment:143 followup: ↓ 144 Changed 6 years ago by
 Status changed from needs_work to needs_review
 Work issues improve error handling in meataxe (within reason) deleted
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.
comment:144 in reply to: ↑ 143 Changed 6 years ago by
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 standalone executables an error return value is returned. This is yet to be done for MTX_ASSERT and MTX_VERIFY.
comment:145 Changed 6 years ago by
OTOH, I feel that I am overdoing it. So, I will *not* deal with MTX_ASSERT and MTX_VERIFY, unless requested by the referee.
comment:146 in reply to: ↑ 51 ; followup: ↓ 147 Changed 6 years ago by
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.
comment:147 in reply to: ↑ 146 Changed 6 years ago by
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.
comment:148 Changed 6 years ago by
I have sent my patches upstream a minute ago.
comment:149 followup: ↓ 150 Changed 6 years ago by
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).
comment:150 in reply to: ↑ 149 ; followup: ↓ 151 Changed 6 years ago by
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 nonprime fields of odd characteristic meataxe is clearly superior over the current Sage implementation.
 Concerning StrassenWinograd (which is not in upstream meataxe): Compare the results of multiplication using school book and StrassenWinograd 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:
byexcept BaseException:
(I assume that people writeexcept:
by mistake, so it's better to be explicit that you really want to catch everything).
OK, can do.
comment:151 in reply to: ↑ 150 Changed 6 years ago by
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.
comment:152 followup: ↓ 154 Changed 6 years ago by
I think you are still misusing 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)
comment:153 Changed 6 years ago by
You can remove these parts from spkginstall
:
+# 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)
comment:154 in reply to: ↑ 152 ; followup: ↓ 155 Changed 6 years ago by
Replying to jdemeyer:
I think you are still misusing some of the error handling.
For example, you have
int MatEchelonize(Matrix_t *mat) except 1but 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?
comment:155 in reply to: ↑ 154 Changed 6 years ago by
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?
comment:156 Changed 6 years ago by
 Commit changed from f73337711df114571d1be0fa61140a41628703b7 to 710668d45f5f1e75672482db62a00b57c29cd16a
Branch pushed to git repo; I updated commit sha1. New commits:
710668d  Remove overcautious commands in spkginstall; rely on default error return values in matrix_gfpn_dense

comment:157 Changed 6 years ago by
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.
comment:158 Changed 6 years ago by
 Report Upstream changed from None of the above  read trac for reasoning. to Reported upstream. No feedback yet.
comment:159 Changed 6 years ago by
 Status changed from needs_review to needs_work
 Work issues set to Pickling
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/sitepackages/sage/ext/interrupt/interrupt.so(+0x4125)[0x7f5a7dd08125] /home/king/Sage/git/sage/local/lib/python2.7/sitepackages/sage/ext/interrupt/interrupt.so(+0x4177)[0x7f5a7dd08177] /home/king/Sage/git/sage/local/lib/python2.7/sitepackages/sage/ext/interrupt/interrupt.so(+0x60c9)[0x7f5a7dd0a0c9] /lib64/libpthread.so.0(+0xf890)[0x7f5a85922890] /home/king/Sage/git/sage/local/lib/python2.7/sitepackages/sage/matrix/matrix_gfpn_dense.so(+0xd9a0)[0x7f5a46fad9a0] /home/king/Sage/git/sage/local/lib/python2.7/sitepackages/sage/matrix/matrix_dense.so(+0x629a)[0x7f5a6290029a] /home/king/Sage/git/sage/local/lib/python2.7/sitepackages/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/libdynload/cPickle.so(+0x4f4e)[0x7f5a80250f4e] /home/king/Sage/git/sage/local/lib/python2.7/libdynload/cPickle.so(+0xa669)[0x7f5a80256669] /home/king/Sage/git/sage/local/lib/python2.7/sitepackages/sage/structure/sage_object.so(+0x21fa2)[0x7f5a7e83bfa2] /home/king/Sage/git/sage/local/lib/python2.7/sitepackages/sage/structure/sage_object.so(+0x77df)[0x7f5a7e8217df] /home/king/Sage/git/sage/local/lib/python2.7/sitepackages/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. ...
comment:160 Changed 6 years ago by
What can I change in my local git configuration so that I can push to trac if port 22 is blocked?
comment:161 followup: ↓ 162 Changed 6 years ago by
You could try the solution posted here: http://stackoverflow.com/questions/4891527/gitprotocolblockedbycompanyhowcanigetaroundthat, 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).
comment:162 in reply to: ↑ 161 Changed 6 years ago by
Replying to tscrim:
You could try the solution posted here: http://stackoverflow.com/questions/4891527/gitprotocolblockedbycompanyhowcanigetaroundthat, 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.
comment:163 Changed 6 years ago by
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
comment:164 Changed 6 years ago by
 Commit changed from 710668d45f5f1e75672482db62a00b57c29cd16a to c67cb42042b25d17aeb17999ac7430aafef0f142
Branch pushed to git repo; I updated commit sha1. New commits:
c67cb42  Fix pickling of meataxe matrices

comment:165 Changed 6 years ago by
 Status changed from needs_work to needs_review
 Work issues Pickling deleted
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 representationand 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 bigendian machine can correctly be unpickled with a littleendian 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.
comment:166 Changed 5 years ago by
 Milestone changed from sage6.9 to sage6.11
 Status changed from needs_review to needs_work
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/sitepackages/sage/doctest/forker.py", line 496, in _run self.compile_and_execute(example, compiler, test.globs) File "/home/travis/sage/local/lib/python2.7/sitepackages/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'>
comment:167 Changed 5 years ago by
However, all tests pass in src/matrix
with it installed.
comment:168 Changed 5 years ago by
 Commit changed from c67cb42042b25d17aeb17999ac7430aafef0f142 to 74edf19ac9217428c482cef93e77226cca84aab3
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 rowreduced echelon form

f816e41  Fix doctests when meataxe is installed

80af75e  Use and propagate specific return values on error in matrixrelated MeatAxe functions.

6649c82  Remove overcautious commands in spkginstall; 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

comment:169 Changed 5 years ago by
 Status changed from needs_work to needs_review
I made the two failing tests optional and at that occasion merged the latest develop.
Is it fine now?
comment:170 Changed 5 years ago by
 Reviewers set to Jeroen Demeyer, Travis Scrimshaw
 Status changed from needs_review to positive_review
Looks good to me now. Thanks.
comment:171 Changed 5 years ago by
 Branch changed from u/SimonKing/meataxe to 74edf19ac9217428c482cef93e77226cca84aab3
 Resolution set to fixed
 Status changed from positive_review to closed
comment:172 Changed 5 years ago by
 Commit 74edf19ac9217428c482cef93e77226cca84aab3 deleted
there is no tarball on the mirrors, apparently
comment:173 Changed 5 years ago by
Fixed
Here is more on
libmeataxe1.0.spkg
.I started with Release 2.2.4 of the Aachen
C MeatAxe
. This is partially because I first came in touch withMeatAxe
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 filepatches/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: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 libmeataxe1.0.spkg). I wrap the C library using Cython (this is trac12103_meataxe.patch)MeatAxe
uses school book multiplication. I added StrassenWinograd multiplication, using a schedule that minimizes memory consumption. Seesrc/src/window.c
.MTXLIB
was supposed to tellMeatAxe
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.T_STRING
intoT_STRINGS
,MTXmatid
.The aim of this ticket is to make the modified
MeatAxe
an optional back end for dense matrices overGF(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 libmeataxe1.0.spkg
would automatically patch the Sage library. Question to Sage experts: Is that possible?Anyway, I think it can be reviewed...