How to tune numerical D? (matrix multiplication is faster in g++ vs gdc)

[J]

On Monday, 4 March 2013 at 08:02:46 UTC, J wrote:
> That's a really good point. I wonder if there is a canonical 
> matrix that would be preferred?

I'm not sure if they are the recommended/best practice for matrix 
handling in D at the moment (please advise if they are not), but 
with a little searching, I found that the SciD library has nice 
matrices and MattrixViews (column-major storage, LAPACK 
compatible).

Now I like MatrixViews because they let me beat the original 
(clearly non-optimal) C matrix multiplication by a couple 
seconds, and the D code with operator overloading in place makes 
matrix multiplication look elegant.

One shootout benchmark down, eleven to go. :-)

- J

p.s. Am I right in concluding there are no multimethods (multiple 
dispatch) in D... it seemed a little awkward to have to wrap the 
MatrixView in a new struct solely in order to overload 
multiplication. Is there a better way that I've missed?

#beating the C/malloc/shootout based code that ran in 1m32sec:

$ time ./scid_matrix_test
-1015380632 859379360 -367726792 -1548829944

real    1m27.967s
user    1m27.930s
sys 0m0.030s
$ time ./scid_matrix_test
-1015380632 859379360 -367726792 -1548829944

real    1m28.259s
user    1m28.230s
sys 0m0.020s
$

module scid_matrix_test;

// compile: gdmd -O -release -inline -noboundscheck 
scid_matrix_test.d -L-L/home/you/pkg/scid/scid/ -L-lscid 
-L-L/usr/lib/x86_64-linux-gnu/ -L-lgfortran -L-lblas -L-llapack 
-L-L.

import scid.matrix;

import std.stdio, std.string, std.array, std.conv;
const int SIZE = 2000;

import std.stdio;


// Doesn't seem to have multiple dispatch / multimethods, so
// I guess we have to wrap MatrixView to implement 
opBinary!"*"(lhs,rhs)

struct Multipliable (T) {
   MatrixView!(T) _m;
   alias _m this;
   this(MatrixView!(T) src) {
     _m = src;
   }

   Multipliable!(T) opBinary(string op : "*")(ref Multipliable 
rhs) if (op == "*") {
     auto r = Multipliable!int(matrix!int(_m.rows,rhs.cols));
     return mmult(this, rhs, r);
   }

}


Multipliable!(int) mkmatrix(int rows, int cols)
{
     auto m = Multipliable!int();
     m._m = matrix!int(rows,cols);

     int count = 1;
     for(int i = 0; i < rows; ++i)
     {
         for(int j = 0; j < cols; ++j)
         {
             m[i,j] = count++;
         }
     }
     return(m);
}



Multipliable!(int) mmult(ref Multipliable!(int) m1, ref 
Multipliable!(int) m2, ref Multipliable!(int) m3) {
     int i, j, k, val;
     ulong rows = m1.rows;
     ulong cols = m1.cols;
     for (i=0; i<rows; i++) {
     for (j=0; j<cols; j++) {
         val = 0;
         for (k=0; k<cols; k++) {
         val += m1[i,k] * m2[k,j];
         }
         m3[i,j] = val;
     }
     }
     return(m3);
}


void main()
{
     auto m1 = mkmatrix(SIZE,SIZE);
     auto m2 = mkmatrix(SIZE,SIZE);
     auto mm = mkmatrix(SIZE,SIZE);

     mm = m1 * m2;

     writefln("%d %d %d %d",mm[0,0],mm[2,3],mm[3,2],mm[4,4]);

}