Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
Announcements
The Intel sign-in experience is changing in February to support enhanced security controls. If you sign in, click here for more information.
6680 Discussions

## PARDISO, Sparse RHS

Beginner
225 Views

Hi,

I have sparse system of equations around 800 000 unknowns. After factorization I want to solve it with multiple sparse RHS. I set array parm, as described in manual, and I set iparm(31). But it doesnt matter if iparm(31) =1 or iparm(31)=2 I always get the same solving times as I was using iparm(31)=0. I set parm array before first call to pardiso. Do I miss something?

Best regards,
Pawel J.

2 Replies
Moderator
225 Views

Pawel, am I understand you right, that when you make the call of Pardiso on the 33 stage, you set the proper number   right-hand side, like nrhs = 100 ( as an example)?

if yes, then, if you want to see the benefits of using iparm(31) [ quote from FR: "which can be used if only a few components of the solution vectors are needed or if you want to reduce the computation cost at the solve step by utilizing the sparsity of the right-hand sides" ]., please tell us what is the sparsity of your right-hands?

From our expectation, in the case of big number of zeroes, the solving step would be accelerate significantly.

it also would be much better if the reproducer the problem will be provided for checking the problem on our side.

thanks

Beginner
225 Views

Thank you for response!

The sparsity of RHS is around 0.006%. It will be hard to provide you the exact code, because it's a part of whole complicated system, but I extracted some main parts that call PARDISO:

iparm(1) = 1

! iparm(2) = 0
! iparm(2) = 2
iparm(2) = 3
iparm(3) = 0

iparm(4) = 0

iparm(5) = 0
! iparm(5) = 1

iparm(6) = 0

iparm(7) = 0

iparm(8) = 0

iparm(9) = 0

iparm(10) = 8

iparm(11) = 1
iparm(12) = 0

iparm(13) = 1
iparm(14) = 0
iparm(15) = 0
iparm(16) = 0
iparm(17) = 0
iparm(18) = -1
iparm(19) = 0

iparm(20) = 0
iparm(21) = 1

iparm(22) = 0
iparm(23) = 0

iparm(24) = 0

iparm(25) = 0

iparm(27) = 0
iparm(28) = 0

iparm(30) = 0

iparm(31) = 1

iparm(35) = 0   ! Fortran style indexing.

iparm(60) = 0

!multiple sparse rhs
nrhs = 12

allocate(manyx(neq * nrhs))
allocate(manyx1(neq * nrhs))
allocate(permSparse(neq * nrhs))

do i = 0, neq*nrhs
manyx1(i) = 0.d0
permSparse(i) = 0
enddo
call createManyFullRHS2(manyx, permSparse, neq, nrhs)

!end of multiple sparse rhs

phase = 11 ! only reordering and symbolic factorization
CALL pardiso (pt, maxfct, mnum, mtype, phase, neq,   ! (PT, MAXFCT, MNUM, MTYPE, PHASE, N,
1           a,ia,ja),        !  A, IA, JA,
2           permSparse, nrhs, iparm, msglvl, ddum, ddum, mklerr)   !  PERM, NRHS, IPARM, MSGLVL, B, X, ERROR)

phase = 22 ! numerical factorization

CALL pardiso (pt, maxfct, mnum, mtype, phase, neq,   ! (PT, MAXFCT, MNUM, MTYPE, PHASE, N,
1           a,ia,ja),        !  A, IA, JA,
2           permSparse, nrhs, iparm, msglvl, ddum, ddum, mklerr)   !  PERM, NRHS, IPARM, MSGLVL, B, X, ERROR)

phase = 33 ! only solve

CALL pardiso (pt, maxfct, mnum, mtype, phase, neq,
1           a,ia,ja),        !  A, IA, JA,
2      parmSparse, nrhs, iparm, msglvl, manyx, manyx1, mklerr)

subroutine createManyFullRHS2(x1,perm,neq,no)

implicit none
integer :: neq, i,no,k
real*8    :: x1(neq*no)
integer    :: perm(neq*no)
integer :: howMany

howMany = INT(neq * no * 0.00006)

do k=1, no
do i = 1, neq
if(Mod(k+i,neq*no/howMany).eq.0) then
x1(i+(k-1)*no) = -1.0d0/(3.0d0*k)
perm(i+(k-1)*no) = 1
else
x1(i+(k-1)*no) = 0.0d0
endif
enddo
enddo

end