I would like to use ISS from MKL for solving linear systems of equations, originating from discretization of partial differential equations used in CFD.
I found out that CG and FGMRES are available, but not everything is clear from the documentation:
- Threading: Can ISS be threaded? I am using Pardiso and it can be made parallel easily.
- Preconditioning: In case of CG, RCI_request might equal 3, which indicates to "apply the preconditioner C inverse to tmp(1:n,3) and put the result in tmp(1:n,4)". This is marked an OPTIONAL in the scheme of using the RCI CG Routines, however, it is stated that "Indicates that you MUST apply the preconditioner to tmp(:, 3), put the result in the tmp(:, 4), and return control back to the routine dcg" in the description of dcg. Finally, using preconditioners is not recommended for CG as it generally produces non-symmetric matrices.
By the way, can Pardiso as a direct solver work also in an iterative "mode", meaning that it would not calculate the solution, but just a user-defined number of iterations?
I'll try to answer your questions.
1) When you're using Iterative Sparse Solvers in Intel MKL, they work through RCI which means that the user should supply building blocks for the main operations involved in solving Ax=b iteratively, which are computing the matrix-vector product A*v, computing action of the preconditoner P * v (optional) and also checking the stopping criteria.
You can see an example of doing that in examples/solverc/fgmres_full_funct_c.c, for instance. There, for RCI_request = 1 just mkl_sparse_d_mv is used, for RCI_request=3 a custom preconditoner action is hardcoded. The user is free to choose any method for the computational routines. To know what should be the input and output of the user-provided functions, one can look in the documentation, e.g. https://software.intel.com/en-us/mkl-developer-reference-c-dfgmres for dfgmres. So, the user can provide a parallel routine to perform each of the RCI_request (in the considered example, mkl_sparse_d_mv in RCI_request = 1 will be parallel if a correct threading is used when linking mkl).
2) As for the preconditioner in CG method (RCI_request=3), parameter ipar controls whether preconditioner will be used or not, see https://software.intel.com/en-us/mkl-developer-reference-c-cg-interface-description#RCI_ID_COMMONPARAMETERS. If you set ipar=0, then CG will never call RCI_request=3 and thus you don't need ti supply any routine for it.
3) Last, but not least: PARDISO is a direct solver but it supports using preconditioned CGS/CG in the mode when you want to solve multiple systems. Once a matrix is factored, you can use the factorization as a preconditioner when solving next systems (for example, if the matrix changes slightly and you expect that the previous matrix is a good preconditoner), see description of iparm at https://software.intel.com/en-us/mkl-developer-reference-c-pardiso-iparm-parameter. I believe this is not what you need though.