Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
6956 Discussions

Help! How to use Pardiso in this way?

zlh007
Beginner
381 Views
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.
0 Kudos
7 Replies
zlh007
Beginner
381 Views
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
0 Kudos
Sergey_P_Intel2
Employee
381 Views
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
0 Kudos
zlh007
Beginner
381 Views
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

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
0 Kudos
zlh007
Beginner
381 Views
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.
0 Kudos
zlh007
Beginner
381 Views
Does anyone know this problem?
0 Kudos
Alexander_K_Intel2
381 Views
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
0 Kudos
Gennady_F_Intel
Moderator
381 Views

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

0 Kudos
Reply