Shifted inverse iteration example

[gflegar]

As promised, here is an example which uses the newly added linear combination support (#131). It implements the shifted inverse iteration method which needs to solve a system with (A - zI) in each iteration. For the stopping criterion, we need multiplication with A, so both matrices are needed. To avoid storing the matrix twice (and enable the solution of a linear system no matter in what format A is stored), we use gko::Combination to represent (A - zI).

TODO:

• merge Linear combination of linear operators #131, as this PR also includes commits from there

[hartwiganzt]

very nice job!

[hartwiganzt]

Not sure where, but I feel this would also be a nice general routine part of the library. Obviously, it does not really fit the "linear operator" scheme...

[gflegar]

Not sure where, but I feel this would also be a nice general routine part of the library. Obviously, it does not really fit the "linear operator" scheme...

It may fit the factory scheme. If you look at the eigenvalue decomposition A V = V L, then A = V L V^-1, which means you could treat it as any other factorization algorithm, which will be represented as factories that produce a combined operator containing the factors. That's easy for normal matrices (A^*A = AA^*), where V^--1 = V^*, so there's no additional computation. If V is not full, but still normal, we could also do A -> V L V^*, where the resulting operator represents a low-rank approximation of A.

The only problem is what to do if the matrix is not normal, so V does not have to be orthogonal. Then, one could do a low-rank approximation as A ~ V L V^+, but computing the pseudo-inverse might be expensive. We'll have to think about it more when we decide to support such matrices.

Concerning inverse iterations specifically, not sure we want to have it - it's not a very robust algorithm. It can only find 1 eigenvalue, and only if it's unique. Do people actually need this in practice, or would you always do something a bit more advanced? The natural next step would be to iterate a multidimensional subspace (leading to QR iterations if adding the full domain, not sure if there's an equivalent for sparse).

I'll create an issue later today to have a place to discuss how we may introduce eigensolvers to Ginkgo.

[gflegar]

I was afraid to test it, but everything seems to work like a charm in complex arithmetic (even though the unit tests only check double, something that should be on our TODO list).
So I updated it to use complex. At least then we then have an example of how to use complex with Ginkgo.