topic Re: Help! How to use Pardiso in this way? in Intel® oneAPI Math Kernel Library

Help! How to use Pardiso in this way?

zlh007 — Sat, 24 Oct 2009 06:06:29 GMT

Hello,
Mycode needs to repeatly solve three problems: A(1)*x=b(1),A(2)*x=b(2),A(3)*x=b(3), where A(1),A(2) and A(3) are the coefficient matrices, b(1),b(2) and b(3)are the right hand sides of equations.And matrix A(i)(i=1,2,3) are real unsymmetric. At each time calling Pardiso,the coefficient matrix A(i) is the same (although A(1), A(2) and A(3) are not the same), only the right sides bi(i=1,2,3) change.Buteverytime solving A(i)x=b(i)using the pardiso, the solving proceduresincludethree steps, i.e.phase=11, phase=22and phase=33, although the factorized coefficient matrices have absolutely the same structure and values. Hence,it takes too much time to factorizing the same coefficient matrix at eachcall of Pardiso. Is thereany methods to factorize the coefficient matrix only once, andsubsequently call pardiso onlyto solve the equation withdifferent b(i) without facotrization steps ?Or, how can I obtain the L and U matrices from pardiso because the L and U matrices are saved in the internal memory before I use phase=-1 to clear them? Because the coefficient matrix is too large and sparse, the Lapack LU decompositionroutine does not work.

Re: Help! How to use Pardiso in this way?

zlh007 — Sat, 24 Oct 2009 06:15:02 GMT

To clear the problem, the code now using is:

type(mkl_pardiso_handle) :: pt(64)
INTEGER maxfct, mnum, mtype, phase, neq, nrhs, error, msglvl
parameter (N=3, neq=nx*ny,nrhs=1,maxfct=1,mnum=1)
INTEGER iparm(64)
INTEGER ia(neq+1,N)
INTEGER ja((nx-2)*(ny-2)*5+2*(nx+ny-2),N)
REAL*8 a((nx-2)*(ny-2)*5+2*(nx+ny-2),N)
REAL*8 b(neq,N)
REAL*8 x(neq,N)
INTEGER idum(neq)
REAL*8 waltime1, waltime2, ddum(neq)
maxfct=1
do i = 1, 64
iparm(i) = 0
end do
iparm(2) = 2 ! fill-in reordering from METIS
iparm(3) = 2 ! numbers of processors
iparm(4) = 0 ! no iterative-direct algorithm
iparm(5) = 0 ! no user fill-in reducing permutation
iparm(6) = 0 ! =0 solution on the first n compoments of x
iparm(7) = 0 ! not in use
iparm(8) = 9 ! numbers of iterative refinement steps
iparm(9) = 0 ! not in use
iparm(10) = 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) = 0 ! not in use
iparm(14) = 0 ! Output: number of perturbed pivots
iparm(15) = 0 ! not in use
iparm(16) = 0 ! not in use
iparm(17) = 0 ! not in use
iparm(18) = -1 ! Output: number of nonzeros in the factor LU
iparm(19) = -1 ! Output: Mflops for LU factorization
iparm(20) = 0 ! Output: Numbers of CG Iterations
iparm(25) = 1 !Parallel Forward/Backward Solve.
iparm(28)=1 !Parallel Reordering for METIS
error = 0 ! initialize error flag
msglvl = 0 ! print statistical information
mtype = 11 ! real unsymmetric
do it=1,10000
call Bdef(N,nx,ny,0,delta,r,b) !bchangesat each 'it' step, but a, ia and jaremainthe same values
do k=1,N
pt(:)=mkl_pardiso_handle(0)

phase = 13 ! only reordering and symbolic factorization
call pardiso (pt, maxfct, mnum, mtype, phase, neq, a(:,k), ia(:,k), ja(:,k),&
idum, nrhs, iparm, msglvl, b(:,k), x(:,k), error)

phase = -1 ! release internal memory
call pardiso (pt, maxfct, mnum, mtype, phase, neq, a(:,k), ia(:,k), ja(:,k),&
idum, nrhs, iparm, msglvl, b(:,k), x(:,k), error)
enddo
enddo

Re: Help! How to use Pardiso in this way?

Sergey_P_Intel2 — Tue, 27 Oct 2009 06:34:30 GMT

Hello!

There are two possible ways of solving several matrices by PARDISO at the same time.

First way is for matrices with the same sparsity structure. In this case PARDISO can store in the same handle pt the sparcity structure of LU factorizations (it is the same for all LU factors) and values of LU factorization for different input matrices. For doing so PARDISO provides maxfct and mnum parameters in its interface.

In the case of different sparsity structure the LU factorization of matrices can be kept in memory with different memory address pointers pt. Please see MKL manuals for details.

For your three matrices with identical sparsity structure both variants are applicable but first one is better from the memory utilization and performance points of view. For example, you can set maxfct=3 in your application, call reordering step in PARDISO for matrix A(1), then call factorization step for A(1) /with mnum = 1/, A(2) /with mnum = 2/ and A(3) /with mnum = 3/.
As the result you can call solving PARDISO step for each of the matrices A(1), A(2) or A(3) by using mnum = 1, 2 or 3, respectively.

Regards,
Sergey

Re: Help! How to use Pardiso in this way?

zlh007 — Wed, 28 Oct 2009 18:39:14 GMT

Quoting - Sergey Pudov (Intel)

Thanks for your replying! I tried your method. There still exist some problems confusing me.
Firstly, I set maxfct=3, and follow the method as your suggestion, i.e. reordering once, factorizing 3 times and then solving3 times,the obtained results are wrong. The code is:
pt(:)=mkl_pardiso_handle(0)
phase = 11 ! only reordering and symbolic factorization
mnum=1
call pardiso (pt, maxfct, mnum, mtype, phase, neq, a(:,1), ia(:,1), ja(:,1),&
idum, nrhs, iparm, msglvl, ddum, ddum,error)
do k=1,3
phase = 22 ! factorization
mnum=k
call pardiso (pt, maxfct, mnum, mtype, phase, neq, a(:,k), ia(:,k), ja(:,k),&
idum, nrhs, iparm, msglvl, ddum, ddum, error)
enddo
do k=1,3
phase = 33 !solve
mnum=k
call pardiso (pt, maxfct, mnum, mtype, phase, neq, a(:,k), ia(:,k), ja(:,k),&
idum, nrhs, iparm, msglvl, b(:,k), x(:,k), error)
enddo
But when I change the calling sequence as: reordering once for A(1), factorizingA(1) and solving A(1),factorizingA(2) and solving A(2),factorizingA(3) and solving A(3), the results iscorrect. Butthe procedurehasadditional factorization steps.The code is:
pt(:)=mkl_pardiso_handle(0)
phase = 11 ! only reordering and symbolic factorization
mnum=1
call pardiso (pt, maxfct, mnum, mtype, phase, neq, a(:,1), ia(:,1), ja(:,1),&
idum, nrhs, iparm, msglvl, ddum, ddum,error)
do k=1,3
mnum=k
phase = 22 !factorization
call pardiso (pt, maxfct, mnum, mtype, phase, neq, a(:,k), ia(:,k), ja(:,k),&
idum, nrhs, iparm, msglvl, ddum, ddum, error)
phase = 33 !solve
call pardiso (pt, maxfct, mnum, mtype, phase, neq, a(:,k), ia(:,k), ja(:,k),&
idum, nrhs, iparm, msglvl, b(:,k), x(:,k), error)
enddo
I think the code will run much more quickly without repeatly calling the factorization steps, i.e. only the first three calling of factorization steps for A(1),A(2) and A(3) is needed.Because theequations with the three coefficient matrices will besolved many times (such as Ntime>1000) only withdifferent right hand sides.Is thereany method toimplement only three solving steps everytime in the Ntime loop?
Regards,
Linhao

Re: Help! How to use Pardiso in this way?

zlh007 — Wed, 28 Oct 2009 19:25:53 GMT

In fact, my problem is quite similar to that in
http://software.intel.com/en-us/forums/showthread.php?t=65254
Using one step 1, three step 2, and then three step 3 can'tproduced the correct result, althoughsome results can be obtained without any error warning.

Re: Help! How to use Pardiso in this way?

zlh007 — Thu, 29 Oct 2009 17:59:09 GMT

Does anyone know this problem?

Re: Help! How to use Pardiso in this way?

Alexander_K_Intel2 — Tue, 17 Nov 2009 06:07:34 GMT

Quoting - zlh007

Does anyone know this problem?

Hi!

Try to use nextconstruction, solving step you could use several times:

Pt1(:,1)=mkl_pardiso_handle1(0)

Pt2(:,2)=mkl_pardiso_handle2(0)

Pt3(:,3)=mkl_pardiso_handle3(0)

mnum=1
do k=1,3
phase = 11 ! only reordering and symbolic factorization
call pardiso (pt(:,k), maxfct, mnum, mtype, phase, neq, a(:,k), ia(:,k), ja(:,k),&
idum, nrhs, iparm, msglvl, ddum, ddum, error)
enddo

do k=1,3
phase = 22 ! factorization
call pardiso (pt(:,k), maxfct, mnum, mtype, phase, neq, a(:,k), ia(:,k), ja(:,k),&
idum, nrhs, iparm, msglvl, ddum, ddum, error)
enddo

do k=1,3
phase = 33 !solve
call pardiso (pt(:,k), maxfct, mnum, mtype, phase, neq, a(:,k), ia(:,k), ja(:,k),&
idum, nrhs, iparm, msglvl, b(:,k), x(:,k), error)
enddo

With best regards,

Alexander Kalinkin

Re: Help! How to use Pardiso in this way?

Gennady_F_Intel — Wed, 24 Feb 2010 12:51:00 GMT

Hi zlh007,

The problem you reported has been fixed into the latest MKL 10.2 Update 4 version released recently.

You can Download this version from Intel Registration Center: <>. See announcement about that here. Please let us know if the problem is still exist.

--Gennady