Add blockdiagind() indexing function

[imciner2]

The scalar elements on the diagonals of arrays can be accessed using diagind (https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.diagind) to get the AbstractRange of an entire diagonal and operate on it (e.g. doing D[diagind( D, -1 )] .= 1 to set all values on the diagonal 1 below the main to 1). That also works on a BlockArray for accessing the scalar elements on a diagonal. It would be nice to have a similar function for accessing the block diagonal elements (e.g. access entire blocks on a diagonal) instead of just the scalar elements.

Example of desired behavior:

julia> D = BlockArray( randn(4,4), [2, 2], [2, 2] )                                                                                                                                                                                                       
2×2-blocked 4×4 BlockArray{Float64,2}:
  0.456552   0.508835  │  -1.36404    0.0981592
  0.579405   1.1095    │  -0.171994   1.02463  
 ──────────────────────┼───────────────────────
  0.189181  -0.593364  │   0.931672  -0.648214 
 -0.414329   1.02497   │  -0.990573   0.154431 

julia> D[diagind(D,-1)]                                                                                                                                                                                                                                   
3-element Array{Float64,1}:
  0.5794046395673269
 -0.5933636927794533
 -0.9905727125760019

julia> D[blockdiagind(D,-1)]              
2×2 Array{Float64,2}:
  0.189181  -0.593364
 -0.414329   1.02497

Having this available will simplify some operations such as:

 # Create a NxN block matrix with identity matrices of size nx x nx on the diagonal
otherDiag = kron( I(N), to_matrix( mt, I, nx ) )

D = zeros( mt, N, N )
D[diagind( D, -1 )] .= 1    

# Populate the block diagonal 1 below the main block diagonal with the nx x nx matrix sys.A
otherDiag = otherDiag + kron( D, -sys.A )

could become

 # Create a NxN block matrix with identity matrices of size nx x nx on the diagonal
otherDiag = kron( I(N), to_matrix( mt, I, nx ) )

# Populate the block diagonal 1 below the main block diagonal with the nx x nx matrix sys.A
otherDiag[blockdiagind( otherDiag, -1 )] .= -sys.A

[dlfivefifty]

Will merge a PR