Pardiso for complex symmetrc matrix giving few NaN elements

Ines_F_ · ‎12-05-2016

Dear all,

I use the following code to solve complex symmetric matrix by Pardiso and get result with few NaN with an output error = 0.

I'm using intel composer xe 13.1.3 on linux.

could anyone please help me to find out the problem ? Thanks !

INCLUDE 'mkl_pardiso.f90'

SUBROUTINE pardiso_solver(n,nrhs,nnz,ia,ja,a,b,x)

  USE mkl_pardiso
  IMPLICIT NONE
  
  !IN/OUT
  INTEGER, INTENT(IN) :: n,nnz,nrhs
  INTEGER, Dimension(n+1), INTENT(IN) :: ia
  INTEGER, Dimension(nnz), INTENT(IN) :: ja  
  COMPLEX, Dimension(nnz), INTENT(IN) :: a
  COMPLEX, Dimension(n,nrhs), INTENT(INOUT) :: b
  COMPLEX, Dimension(n,nrhs), INTENT(INOUT) :: x    
  
  !.. Internal solver memory pointer 
  TYPE(MKL_PARDISO_HANDLE), ALLOCATABLE  :: pt(:)
  !.. All other variables
  INTEGER maxfct, mnum, mtype, phase, error, msglvl
  INTEGER error1
  INTEGER, ALLOCATABLE :: iparm( : )
  
  INTEGER i,j,c, idum(1)
  COMPLEX ddum(1)
  
  !..
  !.. Set up PARDISO control parameter
  !..
  ALLOCATE( iparm ( 64 ) )
  
  do i = 1, 64
     iparm(i) = 0
  end do 
  
  iparm(1) = 1 ! no solver default
  iparm(2) = 3 ! fill-in reordering from METIS
  iparm(3) = 0 ! Not in use
  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) = 2 ! numbers of iterative refinement steps
  iparm(10) = 13 ! perturbe the pivot elements with 1E-13
  iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
  iparm(13) = 1 ! maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm(13) = 1 in case of inappropriate accuracy
  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(21) = 1;  !Apply 1x1 and 2x2 Bunch and Kaufman pivoting during the factorization process
  iparm(24) = 1;  !PARDISO uses new two - level factorization algorithm
  
  maxfct = 1 !Maximum number of numerical factorizations
  mnum = 1  !Which factorization to use
  
  error  = 0 ! initialize error flag
  msglvl = 1 ! print statistical information
  mtype  = 6 !complex symmetric matrix
  !mtype  = 13      ! complex unsymmetric matrix
  
  !.. Initiliaze the internal solver memory pointer. This is only
  ! necessary for the FIRST call of the PARDISO solver.
  
  ALLOCATE ( pt ( 64 ) )
  do i = 1, 64
     pt( i )%DUMMY =  0 
  end do
  
  !.. Reordering and Symbolic Factorization, This step also allocates
  ! all memory that is necessary for the factorization
  
  phase = 11 ! only reordering and symbolic factorization
  
  CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
  idum, nrhs, iparm, msglvl, ddum, ddum, error)
     
  WRITE(*,*) 'Reordering completed ... '
  IF (error /= 0) THEN
     WRITE(*,*) 'The following ERROR was detected: ', error
     GOTO 1000
  END IF
  WRITE(*,*) 'Number of nonzeros in factors = ',iparm(18)
  WRITE(*,*) 'Number of factorization MFLOPS = ',iparm(19)
  
  !.. Factorization.
  phase = 22 ! only factorization
  CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
  idum, nrhs, iparm, msglvl, ddum, ddum, error)
  WRITE(*,*) 'Factorization completed ... '
  IF (error /= 0) THEN
     WRITE(*,*) 'The following ERROR was detected: ', error
     GOTO 1000
  ENDIF
  
  !.. Back substitution and iterative refinement
  phase = 33 ! only factorization
  
  CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
  idum, nrhs, iparm, msglvl, b, x, error)
  WRITE(*,*) 'Solve completed ... '
  IF (error /= 0) THEN
     WRITE(*,*) 'The following ERROR was detected: ', error
     GOTO 1000
  ENDIF
  WRITE(*,*) 'The solution of the system is '
  DO i = 1, n
    Do c=1,nrhs
     WRITE(*,*) ' x(',i,',',c,') = ', x(i,c)
    End DO
  END DO
        
  1000 CONTINUE
  !.. Termination and release of memory
  phase = -1 ! release internal memory
  CALL pardiso (pt, maxfct, mnum, mtype, phase, n, ddum, idum, idum, &
  idum, nrhs, iparm, msglvl, ddum, ddum, error1)
  
  IF ( ALLOCATED( iparm ) )   DEALLOCATE( iparm )
  
  IF (error1 /= 0) THEN
     WRITE(*,*) 'The following ERROR on release stage was detected: ', error1
  ENDIF
  
  IF ( error /= 0 ) STOP 1
  STOP 0

END SUBROUTINE

mecej4 · ‎12-05-2016

With the current version of MKL (MKL 2017), the document says that iparm(28) = 1 needs to be specified if the matrices contain single-precision REAL or COMPLEX elements. You are setting iparm(28)=0 in the loop soon after allocating the array iparm. That is for double-precision.