I use Pardiso to solve reactive transport problem. I tried to do symbolic factorization at first and then in the following newton iterations, only do numerical factorization and substitution.
For the simple problem, this works fine. But for the complex problem, things are quite different. The results are correct in the beginning, but after many iterations, the result is incorrect with large error. Then I tried to do symbolic factorization every step, it can generate correct results.
The structure of matrix (ia, ja) does not change for the whole simulation, only the coefficient (a) and right hand (b) changes. What's wrong with my setting?
The parameter for pardiso are as follows:
iparm = 0
iparm(1) = 1 ! no solver default
iparm(2) = 3 ! fill-in reordering from METIS ,0-MIN DEGREE, 2-METIS, 3-OPENMP VERSION
iparm(3) = 0 ! numbers of processors. Input the next call mkl_set_dynamic(0), mkl_set_num_threads(n);
iparm(4) = 61 ! 0-no iterative-direct algorithm; 10*L+K, K=1 CGS, K=2 CGS for symmetric, 1.0E-L: stopping criterion
iparm(5) = 0 ! no user fill-in reducing permutation
iparm(6) = 0 ! if == 0, the array of b is replaced with the solution x.
iparm(7) = 0 ! Output, Number of iterative refinement steps performed
iparm(8) = 9 ! numbers of iterative refinement steps, must be 0 if a solution is calculated with separate substitutions (phase = 331, 332, 333)
iparm(9) = 0 ! not in use
iparm(10) = 13 ! Default value 13, perturbe the pivot elements with 1E-13
iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
iparm(12) = 0 ! not in use
iparm(13) = 1 ! maximum weighted matching algorithm is switched-on (default for non-symmetric)
iparm(14) = 0 ! Output: number of perturbed pivots
iparm(15) = 0 ! Output, Peak memory on symbolic factorization.
iparm(16) = 0 ! Output, Permanent memory on symbolic factorization. This value is only computed in phase 1.
iparm(17) = 0 ! Output, Size of factors/Peak memory on numerical factorization and solution.
iparm(18) = 0 ! Input/output. Report the number of non-zero elements in the factors. >= 0 Disable reporting.
iparm(19) = 0 ! Input/output. Report number of floating point operations to factor matrix A. >= 0 Disable reporting.
iparm(20) = 0 ! Output: Numbers of CG Iterations. >0 CGS succeeded, reports the number of completed iterations.
iparm(24) = 1 ! Parallel factorization control, 0: classic algorithm, 1: two-level factorization algorithm, improve scalability on many threads.
iparm(25) = 0 ! Parallel forward/backward solve control. 0: Use parallel algorithm for the solve step; 1: Use the sequential forward/backward solve.
iparm(27) = 0 !check matrix error, 0-without check, 1-check
maxfct = 1
mnum = 1
nrhs = 1
error = 0 ! initialize error flag
msglvl = 0 ! print statistical information
mtype = 11 ! real unsymmetric
Thanks and regards,