I have a matrix that I generate with a procedure that looks like this:

blaMX[ptsX_Integer, ptsY_Integer] :=   Table[    With[     {      ix = 1 + Floor[(i - 1)/(ptsY)], jx = 1 + Floor[(j - 1)/(ptsY)],      iy = Mod[i, ptsY, 1], jy = Mod[j, ptsY, 1]      },     If[iy != jy,       0,       k1[ix, jx]       ] +      If[ix != jx,       0,       k2[iy, jy]       ]     ],    {i, ptsX*ptsY},    {j, ptsX*ptsY}    ]; 

where I have two matrices k1 of dimension ptsX and k2 of dimension ptsY

It generates a block-diagonal matrix with a bunch of block-bands. The basic structure looks like:

blaMX[6, 9] /. {_k1 -> 1, _k2 -> 2, c -> 3} // MatrixPlot 

each of those diagonal blocks is k2 duplicated:

blaMX[3, 6][[7 ;; 12, 7 ;; 12]] /. _k1 -> 0 // Grid  k2[1,1] k2[1,2] k2[1,3] k2[1,4] k2[1,5] k2[1,6] k2[2,1] k2[2,2] k2[2,3] k2[2,4] k2[2,5] k2[2,6] k2[3,1] k2[3,2] k2[3,3] k2[3,4] k2[3,5] k2[3,6] k2[4,1] k2[4,2] k2[4,3] k2[4,4] k2[4,5] k2[4,6] k2[5,1] k2[5,2] k2[5,3] k2[5,4] k2[5,5] k2[5,6] k2[6,1] k2[6,2] k2[6,3] k2[6,4] k2[6,5] k2[6,6] 

Each band is made up of blocks, as an unwrapping of k1 with the Band starting at {n+(i-1)*ptsY, n+(j-1)*ptsY} being k1[[i+1, j+1]] (essentially just using Quotient:

With[{n = 1},   Table[    blaMX[6, 9][[       (i - 1)*9 + n,        (j - 1)*9 + n       ]] /. _k2 -> 0,    {i, 6},    {j, 6}    ]   ] // Grid 

Now, my generation procedure is wildly inefficient. I do ptsX^2*ptsY^2 work, when I really only have ptsX*ptsY*(ptsX+ptsY-1) (which is just 2 n^3- n^2 if ptsX==ptsY==n) non-zero elements.

So what is the best way to build a matrix of this form? My current version looks like:

makeKMat[k1_, k2_] :=   With[{     k1SparseRules =      With[{ptsX = Length@k1, ptsY = Length@k2},       Flatten[        Table[         Band[{(i - 1)*ptsY + 1, (j - 1)*ptsY + 1}, {i*ptsY, j*ptsY}] ->          k1[[i, j]],         {i, ptsX},         {j, ptsX}         ],        1        ]       ],     k2SparseRules =      With[{ptsX = Length@k1, ptsY = Length@k2},       {        Band[{1, 1}] ->         ConstantArray[k2, ptsX]        }       ]     },    With[{k3 = SparseArray[k2SparseRules]},     SparseArray[k1SparseRules, Dimensions[k3]] + k3     ]    ]; 

And this works:

k1Pts = 6; k2Pts = 9; k1Test = Array[k1, {k1Pts, k1Pts}]; k2Test = Array[k2, {k2Pts, k2Pts}]; kMatTest = makeKMat[k1Test, k2Test]; kMatTest == blaMX[6, 9]  True 

And is pretty fast, all considered:

k1Test2 = RandomReal[{}, {15, 15}]; k2Test2 = RandomReal[{}, {15, 15}];  {tt, mx} = makeKMat[k1Test2, k2Test2] // RepeatedTiming; tt  0.0071  {tt2, lmx} = lazyMX[k1Test2, k2Test2] // RepeatedTiming; tt2  0.338 

But I can’t help but feel this could be done more efficiently. Can someone propose a better way to build a block-banded SparseArray like this?