Question: Efficiently adding submatrices?

This question is related to the Question How to change only the diagonal values in matrix?

Hello, Maple wizards.

Can we talk about how to perform arithmetic efficiently on submatrices? To illustrate, let matrices A and B be represented by a block structure using submatrices, so that A and B are

A11  A12
A21  A22

and

B11 B12
B21 B22

respectively. I want to efficiently replace the lower submatrix A21 A22  with A21+B21   A22+B22.

Here's the catch (to me, anyway): the indices representing the partition between upper and lower submatrices is changing during execution, as are the contents of B21 and B22.I considered trying

  • a LinearAlgebra[MatrixAdd] with the submatrices as operands but that didn't fly
  • zeroing out or restoring the upper block of B before doing the inplace MatrixAdd, but prep will overcome efficiency of the MatrixAdd
  • element-wise addition on the lower block of the two matrices -- works, but doesn't seem to get me any parallelism from the dgemm in the MKL runtime (I'm on a 6-way AMD processor in Maple 15 - haven't tried on Maple16 yet)
  • creating a matrix alias (assuming they are stored row-major) for the addends and doing in place MatrixAdd -- but creating a new handle for the alias every time can't be cheap,

They aren't huge, maybe 100K x 2, but I'm going to be doing a bunch of them.

A recent answer that acer posted to serena88 made me consider maintaining an Array where the arithmetic could occur, and then use ArrayTools[Copy] to update the lower block of A, but think I need to treat the Array as a Vector to get fast pipelined arithmetic (is that true?), but the varying offset and size gets me back to the aliasing problem. The matrices are stored as rtables, nothing tricky, and except for the zero block, they are definitely not sparse.

Acer brought it last February with a nice inplace matrix multiply, if I need to resort to something like that, I will try. But elegant in-the-Maple-syntax would be nice.

The right approach is probably obvious to any Maple wizard, but from my earlier posts, we know that isn't me.

So, I would welcome your suggestion, and hope this is not an RT*M question.

Looking forward to your thoughts,

- jimmyinhmb

Please Wait...