How To Use Block_diag Repeatedly

I have rather simple question but still couldn´t make it work. I want a block diagonal n^2*n^2 matrix. The blocks are sparse n*n matrices with just the diagonal, first off diagon

Solution 1:

A simple way of constructing a long list of DD2 matrices, is with a list comprehension:

In [128]: sparse.block_diag([DD2 for _ in range(20)]).A
Out[128]: 
array([[-2,  1,  0, ...,  0,  0,  0],
       [ 1, -2,  1, ...,  0,  0,  0],
       [ 0,  1, -2, ...,  0,  0,  0],
       ..., 
       [ 0,  0,  0, ..., -2,  1,  0],
       [ 0,  0,  0, ...,  1, -2,  1],
       [ 0,  0,  0, ...,  0,  1, -2]])

In [129]: _.shape
Out[129]: (80, 80)

At least in my version, block_diag wants a list of arrays, not *args:

In [133]: sparse.block_diag(DD2,DD2,DD2,DD2)
...
TypeError: block_diag() takes at most 3 arguments (4 given)

In [134]: sparse.block_diag([DD2,DD2,DD2,DD2])
Out[134]: 
<16x16 sparse matrix of type '<type 'numpy.int32'>'
    with 40 stored elements in COOrdinate format>

This probably isn't the fastest way to construct such a block diagonal array, but it's a start.

================

Looking at the code for sparse.block_mat I deduce that it does:

In [145]: rows=[]
In [146]: for i inrange(4):
    arow=[None]*4
    arow[i]=DD2
    rows.append(arow)
   .....:     

In [147]: rows
Out[147]: 
[[<4x4 sparse matrix of type'<type 'numpy.int32'>'with10 stored elements (5 diagonals) in DIAgonal format>,
  None,
  None,
  None],
 [None,
  <4x4 sparse matrix of type'<type 'numpy.int32'>'
  ...
  None,
  <4x4 sparse matrix of type'<type 'numpy.int32'>'with10 stored elements (5 diagonals) in DIAgonal format>]]

In other words, rows is a 'matrix' of None with DD2 along the diagonals. It then passes these to sparse.bmat.

In [148]: sparse.bmat(rows)
Out[148]: 
<16x16 sparse matrix of type'<type 'numpy.int32'>'with40 stored elements in COOrdinate format>

bmat in turn collects the data,rows,cols from the coo format of all the input matricies, joins them into master arrays, and builds a new coo matrix from them.

So an alternative is to construct those 3 arrays directly.