4
$\begingroup$

This is sort of a duplicate of Creating overlapped block matrices, but the question wasn't fully answered since the OP never returned to clarify.

I am using Kuba's answer for FEM, but I have been unable to get it to work properly. The problem is that the block matrix isn't being simplified even when the arrays we are working with are known (not arbitrary). Below is sample code. I have used examples for a, b and c, but in my actual code the matrices are larger and much more complicated.

a = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};
b = {{3, 2, 7}, {2, 0, 3}, {4, 1, 2}};
c = {{5, 1, 4}, {4, 5, 8}, {9, 1, 6}};

n = {3, 3};
arrays = Array[#, n] & /@ {a, b, c}; (*arrays to work with*)

app[a1_, a2_, overlap_: 1] := With[{dim = Dimensions@a1},
                                   ArrayPad[a1, {0, n[[1]] - overlap}] + 
                                   ArrayPad[a2, Transpose@{dim - overlap, {0, 0}}]];

Fold[app, First@arrays, Rest@arrays] // MatrixForm

Result

Notice how Mathematica returns Matrix[index] rather than simplifying it. How can I adjust the code and simplify this, and replace e.g. {{1,2,3},{4,5,6},{7,8,9}}[1,1] with 1 since that is row 1, column 1 of the matrix?

$\endgroup$
4
$\begingroup$

Here is a somewhat general routine for building block diagonal matrices with overlaps:

blockOverlap[matList_?ArrayQ, r : (_Integer | {__Integer}) : 1] := 
     With[{spopt = SystemOptions["SparseArrayOptions"]},
          Internal`WithLocalSettings[SetSystemOptions["SparseArrayOptions" ->
                                     {"TreatRepeatedEntries" -> 1}],
          SparseArray[Sort[Flatten[MapThread[
                      ArrayRules[SparseArray[Band[{#1, #1}] -> {#2}]] &,
                      {Accumulate[Prepend[(Length /@ Most[matList]) - r, 1]], 
                       matList}]]]],
          SetSystemOptions[spopt]]]

OP's example:

blockOverlap[{a, b, c}] // MatrixForm

$$\begin{pmatrix} 1 & 2 & 3 & 0 & 0 & 0 & 0 \\ 4 & 5 & 6 & 0 & 0 & 0 & 0 \\ 7 & 8 & 12 & 2 & 7 & 0 & 0 \\ 0 & 0 & 2 & 0 & 3 & 0 & 0 \\ 0 & 0 & 4 & 1 & 7 & 1 & 4 \\ 0 & 0 & 0 & 0 & 4 & 5 & 8 \\ 0 & 0 & 0 & 0 & 9 & 1 & 6 \\ \end{pmatrix}$$

A slightly more general example:

blockOverlap[{a, b, c}, {1, 2}] // MatrixForm

$$\begin{pmatrix} 1 & 2 & 3 & 0 & 0 & 0 \\ 4 & 5 & 6 & 0 & 0 & 0 \\ 7 & 8 & 12 & 2 & 7 & 0 \\ 0 & 0 & 2 & 5 & 4 & 4 \\ 0 & 0 & 4 & 5 & 7 & 8 \\ 0 & 0 & 0 & 9 & 1 & 6 \\ \end{pmatrix}$$

$\endgroup$
  • $\begingroup$ Hi J.M. Thanks a lot for your response! The coupling method you are suggesting requires you to state Band[{k,k}]-> for each array. As I said, I just simplified my FEM matrices. Each one is large, and I can't predict how many matrices there will be, so I need a solution that can apply to multiple matrices of any defined size. This is why Kuba's answer was the closest, for me. The only problem is to get it simplified. Could you please help suggest how to edit the original code, or how to turn yours into a function that admits any number of matrices of any size? $\endgroup$ – Cogicero Mar 9 '16 at 0:13
  • $\begingroup$ Well... 1. are all those matrices on the diagonal? 2. are they all the same size? 3. is the amount of overlap (1 in your OP) the same for all these matrices? $\endgroup$ – J. M. is away Mar 9 '16 at 0:15
  • $\begingroup$ Yes, they are all on the diagonal (creating a tridiagonal matrix, pentadiagonal matrix, etc), all the same size and the overlap is always one (and an addition is always done at the overlap). All I need to supply to the function will be a "list" of the matrices & the size of each matrix, and it will return the sparse block matrix form (using overlaps on the diagonal). $\endgroup$ – Cogicero Mar 9 '16 at 0:16
  • $\begingroup$ If the diagonal blocks are all the same size and the overlap is fixed, it should not be too hard; give me a few... $\endgroup$ – J. M. is away Mar 9 '16 at 0:22
  • 1
    $\begingroup$ Fixed (had used Union[] instead of Sort[]); try it now. $\endgroup$ – J. M. is away Mar 9 '16 at 1:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.