OpenMP threaded nonlinear solver for complex numbers

Question

Problem: I have translated Jacobian-Free Newton-Krylov solver written by C. T. Kelley to Fortran and now want to parallelize it on a shared-memory system with OpenMP. In addition, I want to precondition the system with ILU0 or ILUT preconditioners.

Considered solution: I want to use FGMRES from Intel MKL library - I expect that it is highly optimized and threaded. Since it does not support complex numbers, I will follow Intel's solution to split the problem:

$$ \left(\array{A_r & -A_i \\A_i &A_r}\right)\left(\array{x_r\\x_i}\right)=\left(\array{f_r\\f_i}\right) $$

For ILU0 precondition, I have decided to use a piece of code from Yousef Saad's book Iterative Methods for Sparse Linear Systems Second Edition and parallelize it.

Then I will use ?getri from MKL which "Computes the inverse of an LU-factored general matrix".

Question: Does it sounds reasonable? Any other solver for shared-memory systems?

Answer 1

What other preconditioners have you considered? In my experience with ILU, it's more difficult to parallelize effectively than other methods like alternating Schwarz or algebraic multigrid. As always, your mileage may vary: experiment and profile before trusting nincompoops on the internet.

Using alternating Schwarz breaks up the problem into several sub-problems that are solved serially, for which ILU + a Krylov subspace solver or a direct LU factorization can be very effective depending on size.

As regards OpenMP, the best-of-breed libraries (PETSc, Trilinos) for this sort of thing use message passing via MPI before shared-memory parallelism. The aforementioned libraries can link to Intel MKL;if you choose to use them, you can still take full advantage of all the optimized linear algebra kernels. However, you'll have to weigh the time sunk into learning how to use the library vs. how long it would take you to write it yourself.

Finally, shared memory parallelism lends itself like gangbusters to bugs that are way harder to diagnose than message-passing. What happens when library A uses OpenMP and library B uses OpenMP and they short-circuit each other? That only becomes a consideration when you write big applications, but you never know when some script that you wrote because you couldn't sleep on a long flight metamorphoses into like... half of your Ph.D. thesis*.

*Speaking from experience.