mpgeslir: IR solver

mpgeslir(AF::MPFact, b; reporting=false, verbose=true)

I do not export this function. The idea is that you use mplu 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

Use a multi-precision factorization to solve a linear system with plain vanilla iterative refinement.

MPFact is a union of all the MultiPrecision factorizations in the package. The triangular solver will dispatch on the various types depending on how the interprecision transfers get done.

mpgeslir(MPA::MPArray, b; reporting = false, verbose = true)

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

or

x=MPF\b

is the solition. The output of

mout = \(MPA,b; reporting=true)

or

mout = \(MPF,b; reporting=true)

is a structure. mpout.sol is the solution. mpout.rhist is the residual history. mpout also contains the datatypes TW for high precision and TF for low precision.

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

julia> using MultiPrecisionArrays.Examples

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

julia> MPF=mplu(A);

julia> mout=\(MPF, b; reporting=true);

julia> mout.rhist
6-element Vector{Float64}:
 9.90000e+01
 3.65823e-03
 6.17917e-07
 8.74678e-11
 2.04636e-12
 2.03215e-12

# Stagnation after four IR iterations

julia> [mout.TW mout.TF]
1×2 Matrix{DataType}:
 Float64  Float32