I can straight forwardly create the matrix C as (where the elements 0, A, B are 2X2 matrices) [ 0 A 0 0 0 0 ... 0 0 B 0 A 0 0 0 ... 0 0 0 B 0 A 0 0 ... 0 0 0 0 B 0 A 0 .... 0 0 .......... 0 0 0 0 0 .... B 0] where 0 is the 2X2 matrix consisting of zeros, A, B are given 2X2 matrices (that are to be placed along the upper and lower diagonals respectively) and the matrix C is of size 2JX2J However, it occurred to me that it is probably more efficient to directly make use of sparse matrix notation and then set the non-zero elements using something like (with J = 5) def TheMatrix(J,alpha): A = [[0.5, alpha],[ alpha,0.5]] B =[[0.5,-alpha],[-alpha,0.5]] upperDiag = A*np.ones(J-1) lowerDiag = B*np.ones(J-1) diagonals = [upperDiag, lowerDiag] # non-zero diagonals ASparse = sps.diags(diagonals,offsets=[1, -1]) ASparse = sps.csc_matrix(ASparse) However, this does not seem to work. Does anyone have any suggestions on how to proceed? In detail, the error reported (for the case J =5) is File "myprg.py", line 16, in TheMatrix(J,alpha) upperDiag = A*np.ones(J-1) ~^~~~~~~~~~~~~ ValueError: operands could not be broadcast together with shapes (2,2) (5,) So, I understand that the need is to somehow replace np.ones(J-1) to something that makes it compatible with A, B being 2X2 matrices and not scalars. Continue reading...
