Reshaping block matrix in block diagonal form

[Daniel_Paul]

Hello,

I have a block matrix of the form A = {A1, A2, ..., AN} with symmetric matrices A1, .., AN each of the same order 9*9 (N could be as large as 400). I want to shape the matrix in the block-diagonal form B = {C 0 0, 0 D 0, ..., 0 0 E} with matrices of the form C = {2A1 A1 A1, A1 2A1 A1, A1 A1 2A1} and D and E of the same form with A2 vs. AN instead of A1. If I allow loops, I can figure out a solution, but I try to avoid loops. Any suggestions?

Bye, Daniel 

[Rick_SAS · Accepted solution]

Use Kronecker products to construct the C, D, and E matrices from A1, A2, and A3.

Use the BLOCK function to get the block diagonal matrix of C, D, and E.

proc iml;
N = 3;
N2 = N*(N+1)/2;
A1 = sqrsym(1:N2);
A2 = sqrsym(N2:1);

M = j(N, N, 1) + I(N);
C = M @ A1;
D = M @ A2;
print C, D;

B = block(C, D);
print B;

[Daniel_Paul]

I see the principal idee, but N can be large, so I have to loop in the block part. I try to avoid that loop.

proc iml;
p = 3;
N = 2;

start TBlock(AN, p); 
M = j(p,p,1) + I(p);
C = M @ AN;
return (C);
finish TBlock;

start AD(N, A, p);
t = 1;
do rs=1 to N;
	B = A[t:(t+2),];
	BB = TBlock(B, p);
	AF = block(AF, BB);
	t = t + p;
end;
return (AF);
finish AD;

A1 = {1 2 3, 
	  2 4 5,
      3 5 7};
A2 = {1 7 9,
	  7 8 1,
	  9 1 1};

AC = A1//A2;

AF = AD(N, AC, p);

print AF;
quit;

[IanWakeling]

It looks too complicated to be able to avoid the loop.  What might be slowing things down, is the fact that you are 'growing' the matrix AF by steps within the loop.  After each iteration a new version of AF is declared, and ever more data needs to be moved around in memory.  It would be more efficient to declare AF as a zero matrix, at it's final size, before the loop, and then write the block diagonal sub-matrices using row and column indexing.