I struggle a bit to understand (circumvent if possible) the pardiso performance drop measured as "seconds per right hand side" when decreasing the number of right hand sides.
The sparse matrix feed into pardiso has 2,562,698 rows/columns and the factorization yields 168,330,156 non-zero elements (iparm(18)).
After calling pardiso with phase=12, it is called with phase=33. The iparm vector was left with default values. When using 10 r-h-s, the solution was obtained after 2.21 seconds (phase 33 only), when using 750 r-h-s, the solution was obtained after 16.7 seconds (phase 33 only). Thus, in the first case pardiso needed 0.221 seconds per r-h-s, in the second case 0.022 seconds per r-h-s, which is a decrease of almost factor 10.
The timing was obtained on a Intel(R) Xeon(R) CPU E5-2697 v4 @ 2.30GHz with 36 real cores and 256GB ram. The number of MKL cores was 36.
The iparm vector for the 10 r-h-s case was:
for the 750 r-h-s case:
Any suggestion about how I can achieve a time per r-h-s for the 10 r-h-s case similar to that for the 750 r-h-s case are highly appreciated.
I've got an idea (unfortunately, not yet the out-of-box pardiso improvement). Can you try to set iparm = 2 and tell us if you see any improvement? This option would allow pardiso to parallelize over matrix, not over the righthand sides.
I'd expect this should help improve the performance.