Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Highlighted
Beginner
12 Views

PARDISO, Sparse RHS

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
Highlighted
Moderator
12 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
Highlighted
Beginner
12 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