Add InverseMap for lazy inverse of a linear map.

[fredrikekre]

Idea from the inverse_map class in deal.ii: https://www.dealii.org/current/doxygen/deal.II/group__LAOperators.html#ga87e38fbde431397c069a88692bd24ae7 with example usage in https://www.dealii.org/current/doxygen/deal.II/step_20.html#SolvingusingtheSchurcomplement

cc @kimauth

[codecov]

Codecov Report

Merging #184 (9d67671) into master (99bbc0e) will increase coverage by 0.00%.
The diff coverage is 100.00%.

❗ Current head 9d67671 differs from pull request most recent head 63f290f. Consider uploading reports for the commit 63f290f to get more accurate results

@@           Coverage Diff           @@
##           master     #184   +/-   ##
=======================================
  Coverage   98.60%   98.61%           
=======================================
  Files          16       17    +1     
  Lines        1365     1373    +8     
=======================================
+ Hits         1346     1354    +8     
  Misses         19       19           

Impacted Files	Coverage Δ
src/LinearMaps.jl	100.00% <ø> (ø)
src/inversemap.jl	100.00% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 99bbc0e...63f290f. Read the comment docs.

[dkarrasch]

Thanks for that very nice contribution! I have a few requests:

• Could you please add a short feature announcement in the history.md, perhaps together with the two verbal examples from the docstring: factorization and ldiv!, or LinearMap and iterative solver.
• Speaking of these two use cases, it would be nice to add two such examples to the title page of the documentation, where we currently briefly demonstrate how to use iterative eigensolvers.
• This PR currently adds very few methods for InverseMap, and therefore relies on many generic fallbacks: issymmetric, ishermitian, isposdef, transpose and adjoint. I wonder if we can do a little better. The first three should inherit that property from A; transpose and adjoint can be thrown on A because inversion and transposition commute. It would be good then to test that these work properly under multiplication.

[fredrikekre]

Thanks for the review, I think I updated everything (also took the liberty to bump the version number).

[dkarrasch]

Should we add a little convenience with

InverseMap(A::AbstractMatrix) = InverseMap(factorize(A))

? I know it may not work in all cases, but as a default maybe not too bad?

[fredrikekre]

I thought about it too. However, ldiv! itself doesn't support it deliberately and forces you to factorize yourself (I can't find a reference to this now, but I believe that method even threw a descriptive error message about this in some version of Julia). Maybe the same argumentation doesn't apply here though.

Edit:
And also, that would remove some control from the user, and wouldn't let you do e.g.

iA = InverseMap(A::AbstractMatrix; solver=IterativeSolvers.cg!)

since that would try to factorize A.

Edit 2: I guess you could do

iA = InverseMap(LinearMap(A::AbstractMatrix); solver=IterativeSolvers.cg!)

to avoid it, but I feel like it is better to ask the user to provide an appropriate factorization.

[dkarrasch]

That's all true. Also, implicit conversion/modification is not done anywhere in this package (only wrapping), and factorize is even an expensive operation, so that was a pretty stupid idea of mine. 🤦

[fredrikekre]

I added some tests for 3- and 5 arg mul! too.

[dkarrasch]

Thanks Fredrik once again. The tests are awesome, especially the last test example which showcases the seamless integration with IterativeSolvers.jl nicely. I wonder if we should make that one more visible? Otherwise this is good to go.

[fredrikekre]

It is basically the same as the example I included on the index page. I can expand on it I suppose.

[dkarrasch]

I liked the extra preconditioner bit and the "context", i.e., the comparison to the original blocked problem. But I guess you're right. That relates more to IterativeSolvers.jl, and interested users will very likely know that context anyway.

[fredrikekre]

I'm happy to write a more comprehensive example, but doesn't feel like it fits on the index page where you just quickly want to showcase the package.

[dkarrasch]

If you don't mind the effort, then you could perhaps add a new page to the documentation on which we demonstrate in some more detail "advanced applications" (or something like this). That would be cool. The "original" blocked operator, for instance, could be built both as a matrix or as a LinearMap (the latter may be advantageous when the building matrices are of very inhomogeneous type, sparse and dense etc.). There are a few other nice applications that we don't demonstrate, that users need to see themselves by looking at the bare offered functionality/docstrings.

[dkarrasch]

@fredrikekre I think I was asking for too much. Let's keep the PR as is. We should pursuit the tutorial idea further, but perhaps first collect ideas on how to do it and what to present. One option could be to keep the presentation short and reference downstream packages. Another option could be to present some nice examples explicitly and in detail without necessarily being relevant in concrete use cases. In any case, I think some more thought is required on how to make it useful.

[fredrikekre]

I started on it already, but I will make a follow up PR instead then :)