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

PARDISO, Sparse RHS

Paweł_J_
Beginner
386 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.

0 Kudos
2 Replies
Gennady_F_Intel
Moderator
386 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

 

0 Kudos
Paweł_J_
Beginner
386 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 

0 Kudos
Reply