Hi!

I have the situation that I construct block matrices where I want to avoid to explicitly have them.

For example, with Fourier transforms we can either use slow matrices or fast FFTs:

```
using FFTW
using LinearAlgebra
z = zeros((2,2))
# my true vectors are 1_000_000 elements long
x1 = randn((2,))
x2 = randn((2,))
# FFT plan
p = plan_fft(x1)
p_mat = fft(I(2), (1,))
# slow mat mul
M = [p_mat z; z p_mat]
res1 = M * [x1; x2]
# fast FFT
res2 = [p * x1; p * x2]
####
julia> res1 == res2
true
```

The issue is that my matrices are not necessarily block diagonals, it can be any kind of diagonal matrix. Also my matrix is not quadratic.

And finally I would like to invert the matrix `M`

which obviously does not work straightforwardly since I can’t access matrix elements but only the forward pass `M * x`

.

But maybe we can calculate the pseudoinverse because we should be able to calculate `M' *x`

?

Is there any kind of package allowing to generate the block matrices and efficient calculations?

For example, it probably is perfect if I could write: `(@efficient_block_matrix [p 0; 0 p]) * [x1; x2]`

Best,

Felix

*Existing packages:*

- BlockArrays.jl allowing only dense matrices
- BandedMatrices.jl same here