- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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

Link Copied

1 Reply

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page