Hi,
I use matrix multiplication to perform basis transforms in a partial differential equation code. Because the basis transforms are orthogonal for different dimensions, the matrix product can be factorised to only multiply across the dimension being transformed. As a result, the elements used in the matrix multiplication will have a non-unit stride between them. Is it possible to do this in Eigen? I can't really find any information about matrices with strided memory access.

For example, I need to take an array with dimensions 40x30x20 representing my three dimensional grid, and perform matrix multiplications on each dimension. i.e.
B = 40x30x20 'matrix',
A = 30x30 basis transform matrix,
C_{i j k} = sum_{p} A_{j p} B_{i p k}.

This is as opposed to the normal matrix product C_{j k} = sum_{p} A_{j p} B_{p k}

Currently, I implement this using Level 3 BLAS by treating B as a 40x30 matrix with entries stored in memory 20 elements apart and then loop over the offsets for the 20 elements.

Ideally, I would like to be able to calculate the product in one expression and have the whole thing vectorised for me, as this will be more efficient than having a loop, but I suspect I might be dreaming if I want that to happen.

I would appreciate any help about how I might be able to do this within Eigen.

Regards,
Graham

Hi,

since Eigen does not support 3D matrices, the best you can do is to store it to a 2D matrices and perform 20 matrix-matrix product just like you currently do with BLAS. However instead of interleaving the 3th dimension, I would recommend you to assemble your 20 matrices such that each matrix is sequentially stored in memory: this optimize cache and won't prevent vectorization. So for instance, if your matrix is column-major, then put the 20 40x30 matrices into a 40 x (30*20) one and use .block(0,20*k,40,30) to access each sub matrix.

Of course you can add some convenience functions, either using generic global functions or by extending MatrixBase with custom functions (http://eigen.tuxfamily.org/dox/Customiz ... MatrixBase). To give you an idea:

Code: Select all

Block subMatrix(int k, int size) {
  return this->block(0,k*cols()/size,rows(),size);
}


If you really need to have the 3th dimension in the most nested storage order (such that you can see your matrix as a 40x30 matrix of vector of size 20), then you can emulate a "stride" object using Map + block. Map allows you to "reshape" your matrix so that your submatrix can be seen as a block.

Unfortunately I must have my storage in that order as in addition to the 30x30 matrix multiplication, I also need to perform a 40x40 matrix multiplication over the first dimension, and a 20x20 matrix multiplication over the last dimension.

I suppose I'm not looking for direct support for 3D matrices (as actually I need N-D matrices), they are stored in such a way that they can be accessed as 2D matrices, but I was hoping that it might be possible to extend Eigen to support more general matrix products and to automatically vectorise them. This may be a large amount of work though.

I was thinking that the operation I want might be able to be vectorised as the 40x40 matrix multiplication can easily be treated as it is just a (40x40) X (40 x (30 * 20)) matrix multiplication treating the second and third dimensions of the B matrix as a single dimension, and I'm guessing that the 20x20 multiplication will work in Eigen and be vectorised in the same way by treating the first and second dimensions as one dimension, and then using the 'transpose' expression template by doing something like:
C= B * A.transpose();
Actually, Eigen wouldn't need to support vectorising that operation in my case as A doesn't change, so A.transpose() can simply be stored.

'By extension' I was hoping that the case of the 30x30 matrix product might also be able to be handled and vectorised.

Anyway, thanks for your assistance.

I'm not nearly as much of an expert as Gael in this matter, but: perhaps your matrices are already big enough that it's worth extracting them from the 3D matrix into plain matrices (so, no stride), performing the product on these plain matrices, and putting the result back into the 3D matrix with stride?

Of course this is less ideal than working "in place" but might be your only option if you want the matrix product to be vectorized (with this approach at least, the matrix product in question is between plain matrices, so for sure it gets vectorized).

The cost of matrix product is O(n^3) while the cost of extracting the matrices and storing them back into the 3D matrix, is just O(n^2), so, for large enough sizes, the overhead brought by that method should become negligible. However 1) your size 40x30x20 isn't very big, and 2) i don't know how bad the additional memory access traffic is for performance so ... I don't know

Indeed, the size of your matrices is a very important factor. Actually if they are greater than 16x16, Eigen will automatically use is own efficient matrix product implementation which is equivalent to a Level3 BLAS matrix product (it can works with any stride). For such "large" matrix it is always faster to use this implementation rather than a vectorized expression. Internally this matrix-matrix product will do some copy of your data:
- one side is copied into a block format suitable for vectorization (this is done per block, not on the whole matrix at once)
- while the other might be copied into subvectors if the data are not correctly aligned.

So actually you don't have to bother too much, just extract your matrices using a block expression and Eigen will do the rest for you ! Now if you think you are missing something to extract your matrices (like a stride expression) then please tell us precisely what you would like because such an expression is my TODO list and it is always better to a maximal number of use cases and opinions.