# Is there an iterative solver that is able to operate on block-matrices (non-concatenated)

**URL:** <https://discourse.julialang.org/t/is-there-an-iterative-solver-that-is-able-to-operate-on-block-matrices-non-concatenated/1516>\
**Category:** Numerics\
**Created:** [January 16, 2017, 3:15pm UTC](https://discourse.julialang.org/t/is-there-an-iterative-solver-that-is-able-to-operate-on-block-matrices-non-concatenated/1516 "2017-01-16T15:15:26Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![rleegates](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rleegates/32/1029_2.png) [@rleegates](https://discourse.julialang.org/u/rleegates)\
**Post date:** [January 16, 2017, 3:15pm UTC](https://discourse.julialang.org/t/is-there-an-iterative-solver-that-is-able-to-operate-on-block-matrices-non-concatenated/1516/1 "2017-01-16T15:15:26Z")

</div>

Dear users,

I haven’t been able to find an implementation of an iterative solver (invertible matrix, no further specialization) general enough to work on matrices of matrices, i.e. A = [B C D…; E F G…; H I J…; …], b = [K; L; M; …], x = [N; O; P; …], A is block sparse, b, x is block dense. I would like to avoid concatenation in my use-case.

Shouldn’t this be generally possible using GMRES? Perhaps there’s an implementation out there that I’ve missed, or a way to use IterativeSolvers.jl/KrylovMethods.jl more appropriately than I have tried?

Thanks, Robert

---

<div class="post-metadata">

**Author:** ![cortner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cortner/32/204_2.png) [@cortner](https://discourse.julialang.org/u/cortner)\
**Post date:** [January 16, 2017, 8:40pm UTC](https://discourse.julialang.org/t/is-there-an-iterative-solver-that-is-able-to-operate-on-block-matrices-non-concatenated/1516/2 "2017-01-16T20:40:56Z")

</div>

I’ve never tried it but have been planning to do so with each block being a `StaticArray`. At least `IterativeSolvers` has no restrictions on the type of the matrix elements so it should _in principle_ work.

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [January 17, 2017, 8:35am UTC](https://discourse.julialang.org/t/is-there-an-iterative-solver-that-is-able-to-operate-on-block-matrices-non-concatenated/1516/3 "2017-01-17T08:35:26Z")

</div>

If you wrap at least one of the individual blocks in a linear operator, LinearOperators.jl will let you build such a block operator using the standard bracket notation. You can then compute products with the result as if it were a usual matrix. Krylov.jl uses linear operators extensively, but it doesn’t have GMRES.

---

<div class="post-metadata">

**Author:** ![rleegates](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rleegates/32/1029_2.png) [@rleegates](https://discourse.julialang.org/u/rleegates)\
**Post date:** [January 17, 2017, 10:29am UTC](https://discourse.julialang.org/t/is-there-an-iterative-solver-that-is-able-to-operate-on-block-matrices-non-concatenated/1516/4 "2017-01-17T10:29:47Z")

</div>

Hi Dominique,  
so we’re talking about a matrix of LinearOperators? A matrix of matrices can also be multiplied by a vector of matrices, however the solvers appear to fail since they don’t know how to initialize the elements to zero (or one). Do you perhaps have a small example using LinearOperators?

---

<div class="post-metadata">

**Author:** ![cortner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cortner/32/204_2.png) [@cortner](https://discourse.julialang.org/u/cortner)\
**Post date:** [January 18, 2017, 1:29am UTC](https://discourse.julialang.org/t/is-there-an-iterative-solver-that-is-able-to-operate-on-block-matrices-non-concatenated/1516/5 "2017-01-18T01:29:33Z")

</div>

you need to define `zero(::T)` and `one(::T)` where `T` is the type of your blocks.

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [January 18, 2017, 6:27am UTC](https://discourse.julialang.org/t/is-there-an-iterative-solver-that-is-able-to-operate-on-block-matrices-non-concatenated/1516/6 "2017-01-18T06:27:26Z")

</div>

A “matrix of linear operators” will avoid concatenating your matrices. I didn’t realize you have multiple right-hand sides. Neither Krylov.jl nor KrylovMethods.jl accommodate that. At first glance, IterativeSolvers.jl might but I’m not sure it’s a good idea in general for iterative methods. In Krylov.jl, the rhs should also be contiguous because we use BLAS calls. Blas.dot is quite a bit faster than dot.

---

<div class="post-metadata">

**Author:** ![rleegates](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rleegates/32/1029_2.png) [@rleegates](https://discourse.julialang.org/u/rleegates)\
**Post date:** [January 18, 2017, 8:58am UTC](https://discourse.julialang.org/t/is-there-an-iterative-solver-that-is-able-to-operate-on-block-matrices-non-concatenated/1516/7 "2017-01-18T08:58:59Z")

</div>

@dpo : Actually, that was supposed to symbolize a vector of matrices, i.e. a block partitioned vector, where each block has `size(b,2) >= 1`. I looked in literature on the solution of such systems and didn’t find all too much except DUNE-ISTL. For now, I abandoned the partitioning of the RHS and solution vector, i.e. I now have a sparse matrix of matrices `A`, as well as contiguous `x` and `b`. I overloaded `*(::SparseMatrixCSC{Matrix{T},Int}, ::Vector{T})`, `size(::SparseMatrixCSC{Matrix{T},Int}, ::Int)` and `eltype(::SparseMatrixCSC{Matrix{T},Int})` in line with the `IterativeSolvers` interface requirements. That did the trick.

@cortner : redefining `zero` and `one` doesn’t work conceptually as the size of the blocks is not constant.