JuliaSparse / Pardiso.jl

Calling the PARDISO library from Julia
Other
100 stars 27 forks source link

Should iterative refinement be opt-in? #72

Open ViralBShah opened 3 years ago

ViralBShah commented 3 years ago

It seems like iterative refinement is the default. https://github.com/JuliaSparse/Pardiso.jl/blob/81c8391e86cebaf71e936363664d878b22b1fd9a/src/Pardiso.jl#L242

Discussion on iterative refinement for Julia sparse solvers: https://github.com/JuliaLang/julia/issues/31105#issuecomment-465030232

cc @andreasnoack

ViralBShah commented 3 years ago

cc @ranjanan

ranjanan commented 3 years ago

According to https://software.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-fortran/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-interface/onemkl-pardiso-parameters-in-tabular-form.html, it looks like 23 corresponds to just factorization + solve? Only 33 is iterative refinement.

ViralBShah commented 3 years ago

It seems that the intel table that @ranjanan links above is not consistent with pardiso params.

According to the pardiso params, parameter 13 is Analysis, numerical factorization, solve, iterative refinement. In fact, I can't quite tell if there is an easy way to solve without refinement.

j-fu commented 3 years ago

Just my 2ct: yes, iterative refinement should be opt-in. May be we want to do e.g. iterative refinement with IterativeSolvers, or the solution is part of a bigger block preconditioner...

With pardiso, setting iparm[8]=0 should prevent iterative refinement unless pivot perturbation was activated.