How to Add Matrices Along Diagonals

[member 428835]

Hi PF!

Let's say I have a 3X3 matrix $K$ and a much larger square matrix, call it $M$. I am trying to add $K$ into $M$ along the diagonals. It's difficult for me to explain (and hence code) so I've attached a picture. Does anyone know a good way of doing this? Notice overlap is only in the corners.

 

[jedishrfu]

Are you trying to do a convolution?

 

[member 428835]

Are you trying to do a convolution?

I am not. I'm ultimately solving a BVP via finite elements. So when constructing a global matrix $M$ I have to superimpose $K$s onto it as described.

 

[jedishrfu]

I guess I’d do this with a couple of for loops. There may be a clever vector scheme but I don’t know Mathematica at all.

For ij=0 to M step K
  For ki=0 to KI
     For kj=0 to KJ
       M(ij+ki, ij+kj) += K(ki,kj)
     End
   End
End

This code might have to be tweaked some like step K-1 instead of step K.

 

[Mark44]

Hi PF!

Let's say I have a 3X3 matrix $K$ and a much larger square matrix, call it $M$. I am trying to add $K$ into $M$ along the diagonals. It's difficult for me to explain (and hence code) so I've attached a picture. Does anyone know a good way of doing this? Notice overlap is only in the corners.

Since you're trying to insert a 3x3 matrix into blocks on a larger matrix, the size of the larger matrix should be 3m x 3m.

Along the lines of what @jedishrfu suggested, you'll need an outer loop that inserts a 3x3 matrix each iteration, and an inner loop that inserts the values for each 3x3 matrix.

The first 3x3 matrix will go in positions (1, 1), (1, 2), (1, 3), (2, 1), (2, 2), (2, 3), (3, 1), (3, 2), (3, 3).

The next 3x3 matrix will go in positions similar to the above, with the upper left index being (4, 4) and the lower right index being (6, 6),

And so on.

The "inner loop" that I mentioned above might actually need another loop nested inside it.

 

[member 428835]

Gotcha! I wasn't sure if Mathematica had some quick internal function that could operate more efficiently than a for loop. So I'm stuck on the syntax; given $K$ nad $M$, in MATLAB I would write

N = 5;
M = zeros(N);
for L = 1:2:N-2
    for i = 1:3
        for j = 1:3
            M(L+i-1,L+j-1) = M(L+i-1,L+j-1) + K(i,j);
        end
    end
end

It should be noted since $M_{33}$ is composed of $K_{11}+K_{33}$, $M$ is actually any matrix with odd dimensions, so not 3m X 3m. I tried showing this in my picture, but I'm not much of an artist. Do either of you know the syntax to compute this triple loop in Mathematica? I've looked it up but the nested loops are confusing for a Mathematic novice like myself.

 

[member 428835]

Perhaps this will clarify my question from above: this is what I have so far in Mathematica

For[L = 1, L <= 1/h - 1, L += 2,
 For[i = 1, i <= 3, i++,
  For[j = 1, j <= 3, j++,
   M[[L + i - 1, L + j - 1]] =
    M[[L + i - 1, L + j - 1]] + K[[i, j]]]]]

and everything works fine so long as I hardcode the dimensions of M. However, what I would like is to have the size of M dependent on a variable h, so something like this
M[h_] := M[h] = Table[0, {i, 1, 1/h + 1}, {j, 1, 1/h + 1}]
The following code executes the previous loop accounting for variable h:

For[L = 1, L <= 1/h - 1, L += 2,
 For[i = 1, i <= 3, i++,
  For[j = 1, j <= 3, j++,
   M[h_][[L + i - 1, L + j - 1]] :=
    M[h][[L + i - 1, L + j - 1]] + K[[i, j]]]]]

But now M is all zeros. Can someone help me here please?

 

[Mark44]

From post #1:

Notice overlap is only in the corners.

This wasn't clear in your first post. So instead of having separate blocks that go down the diagonal, the lower right entry of a block is also the upper left entry of the next block, and this entry is the sum of the corresponding elements from the two blocks. Is that what you're trying to convey?From post #6:

Gotcha! I wasn't sure if Mathematica had some quick internal function that could operate more efficiently than a for loop.

Got me -- I've never written any Mathematica code.

So I'm stuck on the syntax; given $K$ nad $M$, in MATLAB I would write

N = 5;
M = zeros(N);
for L = 1:2:N-2
    for i = 1:3
        for j = 1:3
            M(L+i-1,L+j-1) = M(L+i-1,L+j-1) + K(i,j);
        end
    end
end

It should be noted since $M_{33}$ is composed of $K_{11}+K_{33}$, $M$ is actually any matrix with odd dimensions, so not 3m X 3m.

No, I think you'll need to be more specific about how the blocks (submatrices) are laid out, including the sizes of each block and where they go.
Otherwise this doesn't make much sense. Is M a 5x5 matrix? And is K a 3x3 matrix? For the moment, let's assume that M is initialized to all 0s and let's assume that J is a 3x3 matrix initialized to all 1's and K is 3x3 initialized to all 2's.
If I'm understanding what you're trying to say, then $M_{3,3} = J_{3, 3} + K_{1, 1} = 1 + 2 = 3$.
As I see it, the complete matrix would look like this:
$\begin{bmatrix} 1&1&1&0&0\\ 1&1&1&0&0\\ 1&1&3 & 2&2\\ 0&0&2&2&2\\ 0&0&2&2&2 \end{bmatrix}$
The upper left 3x3 block is J and the lower right 3x3 block is K. Where they intersect the entry is the sum of $J_{3, 3}$ and $K_{1, 1}$

The MATLAB code you have above doesn't appear to be doing this.

I tried showing this in my picture, but I'm not much of an artist. Do either of you know the syntax to compute this triple loop in Mathematica? I've looked it up but the nested loops are confusing for a Mathematic novice like myself.

 

[Mark44]

Perhaps this will clarify my question from above:
this is what I have so far in Mathematica

For[L = 1, L <= 1/h - 1, L += 2,

What is the value of h?

 

[member 428835]

From post #1: This wasn't clear in your first post. So instead of having separate blocks that go down the diagonal, the lower right entry of a block is also the upper left entry of the next block, and this entry is the sum of the corresponding elements from the two blocks. Is that what you're trying to convey?

Yes! Thanks for phrasing it so concisely.It should work. Try this:

h = 0.25;
N = 1/h + 1;                
M = zeros(N);
K = ones
for L = 1:2:N-2
    for i = 1:3
        for j = 1:3
            M(L+i-1,L+j-1) = M(L+i-1,L+j-1) + K(i,j);
        end
    end
end

Matrix should look like this
$\begin{bmatrix} 1&1&1&0&0\\ 1&1&1&0&0\\ 1&1&2 & 1&1\\ 0&0&1&1&1\\ 0&0&1&1&1 \end{bmatrix}$

The issue I face now is solely Mathematica syntax. Specifically, how to allow the size of M to be a function of h in the loop?

 

[Mark44]

Here's a first cut with some working C++ code that produces the matrix I showed in post #3. I'm sure there's a lot of room for improving my code, to make it more general. I'm not familiar with Mathematica, but possibly you can translate this into Mathematica code.

// Test.cpp

#include <iostream>

const int N = 5; // Matrix M is N x N

const int cornerIdx = 2;
void printArray(int arr[][N], int);

int J[][3] = { {1, 1, 1}, {1, 1, 1}, {1, 1, 1} };
int K[][3] = { {2, 2, 2}, {2, 2, 2}, {2, 2, 2} };
int M[N][N] = { 0 }; // Initializes all N x N elements to 0.int main()
{
    int i, j;

    // Insert matrix J as the upper left 3x3 block
    for (i = 0; i <= cornerIdx; i++)
       for (j = 0; j <= cornerIdx; j++) 
          M[i][j] = J[i][j];

    // Insert matrix K as the lower right 3x3 block
    for (i = cornerIdx ; i < N; i++)
       for (j = cornerIdx; j < N; j++)
          M[i][j] = K[i - cornerIdx][j - cornerIdx];

    // Add the lower right corner of J to the upper left corner of K,
    // and insert into the overlap spot of M
    M[cornerIdx][cornerIdx] = J[cornerIdx][cornerIdx] + K[1][1];

    // Print the array
    printArray(M, N);
    return 0;
}

void printArray(int arr[][N], int numRows)
{
    for (int i = 0; i < numRows; i++)
    {
       for (int j = 0; j < N; j++)
          std::cout << arr[i][j] << " ";
       std::cout << std::endl;
    }
    std::cout << std::endl;
}

Here's the output:

1 1 1 0 0
1 1 1 0 0
1 1 3 2 2
0 0 2 2 2
0 0 2 2 2