mpkrir: Krylov-IR solver

MultiPrecisionArrays.mpkrirMethod

mpkrir(AF::MPKFact, b; reporting=false, verbose=false, mpdebug=false, termparms=termparms_default)

Krylov-IR solver

I do not export this function. The idea is that you use mpglu and do not touch either the constructor or the solver directly.

Use a multi-precision factorization to solve a linear system with IR-Krylov metods

This is the generic solver used by GMRES-IR and BiCGSTAB-IR. You use the correct MPKFact = Union{MPGFact,MPBFact} structure and mpkrir will do the right thing.

You should not be calling this directly. Use \ to solve linear systems with the multiprecision factorization and to use the optional kwargs.

We overload the backslash operator to call mpkrir for a multiprecision MPGFact factorization. So if MPA is an MPGArray and

AF = mpglu!(MPA)

Then AF\b maps to

mpkrir(AF, b)

which does the GMRES-IR solve. You can also get the multiprecision factoriztion directly with

AF = mpglu(A)

which builds the mutliprcision MPGArray and then factors the low preicsion copy.

Similarly if MPA is an MPBArray. Then

AF = mpblu!(MPA)

Then AF\b maps to

mpkrir(AF, b)

which does the BiCGSTAB-IR solve.

You can also use the \ operator to harvest iteration statistics.

You may not get exactly the same results for this example on different hardware, BLAS, versions of Julia or this package. I am still playing with the termination criteria and the iteration count could grow or shrink as I do that.

Example

```jldoctest julia> using MultiPrecisionArrays.Examples

julia> N=4096; A = I - 800.0 * Gmat(N); b=ones(N);

julia> AF=mpglu(A);

julia> solout=(AF, b; reporting=true);

Correct result?

julia> x=solout.sol; norm(b-A*x,Inf) 8.20877e-12

Look at the residual history

julia> solout.rhist

4-element Vector{Float64}: 1.00000e+00 1.25881e-10 8.58902e-12 8.20877e-12

and the correction norm history. No correction for the final iteration.

julia> solout.dhist 3-element Vector{Float64}: 2.00129e+02 1.05241e-10 2.28629e-11

Stagnation after the 2nd iteration. Now the Krylovs/iteration

julia> solout.khist 3-element Vector{Int64}: 4 5 5

4-5 Krylovs per iteration.

BiCGSTAB works the same way.

source