PARDISO, Sparse RHS

Paweł_J_ · ‎10-08-2015

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.

Gennady_F_Intel · ‎10-08-2015

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

Paweł_J_ · ‎10-08-2015

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