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

## Pardiso for complex symmetrc matrix giving few NaN elements Beginner
159 Views

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.

```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
```  