Matrix multiplication alglib

Question

How do I multiply two matrices with AlgLib

Answer 1

Disclaimer: I haven't used AlgLib; I'm just going by what the documentation seems to say. I'd be happy to be corrected by someone more expert.

Anyway, I'm afraid the answer seems to be that you need to use cmatrixgemm or rmatrixgemm (which one depends on whether your matrices are real or complex), like this:

rmatrixgemm(m,n,k, 1, A,0,0,0, B,0,0,0, 0, C,0,0);

where:

• m,n,k are the sizes of the matrices (A is m by k, B is k by n, C is m by n)
• the 1 is what to multiply the product by (if you happen to want, say, 3AB instead of AB, put 3 there instead)
• the A,0,0,0 and B,0,0,0 groups are: matrix, row offset, column offset, operation type
  • the operation type is 0 to use A or B as it is, 1 to use the transpose, and 2 to use the conjugate transpose (of course you can't use 2 for rmatrixgemm)
• the next 0 says to add 0*C to the result (if you put 0 here then the initial values in C are completely ignored)
• the two 0s after C are a row and column offset

You might perhaps think that that level of generality is overkill, and that there should be a simpler function that provides those defaults. I wouldn't disagree with that, but so far as I can see there's no such simpler function in AlgLib. You might want to write your own (which would just call rmatrixgemm or cmatrixgemm).

(Why so much generality? Because doing an efficient matrix multiplication requires quite sophisticated code, and it's essentially the same quite sophisticated code as you need to do the more general C=a.f(A).g(B)+b.C operation that *matrixgemm does, and sometimes that more general operation is useful.)

EDITED to add a few more remarks that might be useful.

• The offsets are so that you can do things with submatrices. Being able to do this is useful in some numerical algorithms. I assume m,n,k are the sizes of the submatrices you're using; in the common case, they'll be the same as the dimensions of your arrays, and the offsets will be zero.
• The arrays themselves need to exist and be of appropriate sizes before you call rmatrixgemm or cmatrixgemm. At least, A and B certainly do; C is passed as a ref so it's possible that these functions will create it if it's null on entry.
• You might think from the signature of rmatrixgemm or cmatrixgemm that A and B get copied whereas C is passed by reference, but if I'm not totally confused about C#'s semantics they're all effectively passed by (object) reference.

Answer 2

Just to confirm what Garech wrote:

double[,] a = new double[,] {
    {1,2,3},
    {4,5,6}
};
double[,] b = new double[,] {
    {7,8,9,10},
    {11,12,13,14},
    {15,16,17,18}
};
int m = a.GetLength(0);
int n = b.GetLength(1);
int k = a.GetLength(1);
double[,] c = new double[m,n];
alglib.rmatrixgemm(m, n, k, 1, a, 0,0,0, b,0,0,0, 0, ref c, 0,0);
//c = {{74, 80, 86, 92}, {173, 188, 203, 218}}