MKL Pardiso Factorize once for many solves

[greg_rowe]

Can I use Pardiso to factorize the matrix once and then use that factorization to solve successive rhs vectors?

[ArturGuzik]

Quoting - greg.rowe

Can I use Pardiso to factorize the matrix once and then use that factorization to solve successive rhs vectors?

Hi Greg,

you can. You control this (as well as other stages of computations) by setting phase code. From manual:

Phase 1: Fill-reduction analysis and symbolic factorization
Phase 2: Numerical factorization
Phase 3: Forward and Backward solve including iterative refinements
This phase can be divided into two or three separate substitutions: forward, backward, and diagonal (see above).
Termination and Memory Release Phase (phase 0)

phase Solver Execution Steps
11 Analysis
12 Analysis, numerical factorization
13 Analysis, numerical factorization, solve, iterative refinement
22 Numerical factorization
23 Numerical factorization, solve, iterative refinement
33 Solve, iterative refinement
331 like phase=33, but only forward substitution
332 like phase=33, but only diagonal substitution
333 like phase=33, but only backward substitution
0 Release internal memory for L and U matrix number mnum
-1 Release all internal memory for all matrices

So, for factorization only you'll have something along the lines below:

! Factorization.
phase = 22 ! only factorization
CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, 1 idum, nrhs, iparm, msglvl, ddum, ddum, error)

See also Examples for Sparse Symmetric Linear Systems in the manual.

As any solver it handles multiple RHS vectors at one call as well (see nrhs parameter).

A.

[greg_rowe]

Thanks for the info. Do you know of any clear example where the factorization is done first, followed by multiple solves. Multiple RHS solves using matrix of vectors is not exactly what I want because each subsequent RHS vector in my problem depends on the previous solve. Thanks again, Greg.

(To clear this up, we want to call pardiso once with a phase of 12, then call it multiple times with a phase of 33 using a RHS vector which depends on the result of the previous pardiso call)

[ArturGuzik]

Quoting - greg.rowe

Thanks for the info. Do you know of any clear example where the factorization is done first, followed by multiple solves. Multiple RHS solves using matrix of vectors is not exactly what I want because each subsequent RHS vector in my problem depends on the previous solve. Thanks again, Greg.

(To clear this up, we want to call pardiso once with a phase of 12, then call it multiple times with a phase of 33 using a RHS vector which depends on the result of the previous pardiso call)

Greg,

in the mklexamplessolversource you'll find pardiso_sym_f.f file. Just make this working and you can easily try to factorize first and then call again many times phase = 33 with other B (rhs vector(s)).

Say...

!.. Factorization.
phase = 22 ! only factorization
CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, idum, nrhs, iparm, msglvl, ddum, ddum, error)
WRITE(*,*) 'Factorization completed ... '
IF (error .NE. 0) THEN
WRITE(*,*) 'The following ERROR was detected: ', error
STOP
ENDIF
!.. Back substitution and iterative refinement
iparm(8) = 2 ! max numbers of iterative refinement steps
phase = 33 ! Solve, iterative refinement
do i = 1, n
b(i) = 1.d0
end do

CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, idum, nrhs, iparm, msglvl, b, x, error)
WRITE(*,*) 'Solve completed ... '
WRITE(*,*) 'The solution of the system is '

DO i = 1, n
WRITE(*,*) ' x(',i,') = ', x(i)
END DO


! modify RHS
b =2.d0
x =0.d0

CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, idum, nrhs, iparm, msglvl, b, x, error)
WRITE(*,*) 'Solve completed ... '
WRITE(*,*) 'The solution of the system is '

DO i = 1, n
WRITE(*,*) ' x(',i,') = ', x(i)
END DO

A.

[Sergey_K_Intel1]

Quoting - greg.rowe

Can I use Pardiso to factorize the matrix once and then use that factorization to solve successive rhs vectors?

PARDISO can process several matrices with identical matrix sparsity pattern and is able to store the factors of these matrices at the same time. Please look at the desciption of inputparameters maxfct and mnum (second and third parameters in PARDISO calling sequence).

The factors of matrices with different sparsity structures can be kept in memory with different memory address pointers pt.

All the best
Sergey