mpkrir: Krylov-IR solver

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

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 plain vanilla iterative refinement.

This version is analogous to A\b and combines the factorization and the solve. You start with MPA=MPArray(A) and then pass MPA to mpgeslir and combine the factorization and the solve.

You can also get the multiprecision factorization directly with

MPF=mplu(A)

and then pass MPF to mpgeslir.

I use this to get some timing results and it's also convenient if you want to do factor and solve in one statement.

You can also get this with x = MPA\b.

If you set the kwarg reporting to true you can get the IR residual history. The output of

x = MPA\b

Krylov-IR solver

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) 9.45199e-12

Look at the residual history

julia> solout.rhist 4-element Vector{Float64}: 1.00000e+00 1.27149e-10 9.00036e-12 9.45199e-12

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

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

4-5 Krylovs per iteration.

BiCGSTAB works the same way.