Re: 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]);

}