hello
As you can see in the picture, there is an error when I try to use pardiso to solve eqs. I don't know why it is. What should I do ?
SUBROUTINE SOL_PARDISO(n,lent,a,ia,ja,b,mtype)
C
C PARDISO : PARALLEL SPARSE DIRECT LINEAR SOLVER
C Based on '/sw/sdev/cmkl/10.0.2.18/examples/solver/source/
C pardiso_sym_f.f'
C
use omp_lib
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
INTEGER*8 PT(64)
INTEGER, PARAMETER :: dp = KIND(1.0D0)
C.. All other variables
INTEGER MAXFCT,MNUM,MTYPE,PHASE,N,NRHS,ERROR,error1,MSGLVL,lent
INTEGER iparm(64),idum(1)
REAL(KIND=DP) ddum(1)
dimension a(lent),ia(n+1),ja(lent),b(n),x(n)
DATA nrhs /1/, maxfct /1/, mnum /1/
C.. Set up PARDISO control parameter
do i = 1, 64
iparm(i) = 0
end do
iparm(1) = 3 ! no solver default
iparm(2) = 2 ! fill-in reordering from METIS
iparm(3) = 1 ! numbers of processors
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) = 9 ! numbers of iterative refinement steps
iparm(9) = 0 ! not in use
iparm(10) = 13 ! perturbe the pivot elements with 1E-13
iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
iparm(12) = 0 ! not in use
! if(mtype .eq. -2) then !Improved accuracy using (non-) symmetric weighted matching
iparm(13) = 0 ! Disable matching. Default for symmetric indefinite matrices
if(mtype .eq. 11) then
iparm(13) = 1 ! Enable matching. Default for nonsymmetric matrices.
endif
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(24) =1
iparm(27) = 1 !turn on matrix checker on your matrix
error = 0 ! initialize error flag
msglvl = 0 !0/1 ! print statistical information
c mtype = -2 ! symmetric, indefinite
!! mtype = 11 ! real and nonsymmetric
C.. Initiliaze the internal solver memory pointer. This is only
C necessary for the FIRST call of the PARDISO solver.
do i = 1, 64
pt(i) = 0
end do
call mkl_set_dynamic(0)
call mkl_set_num_threads(4)
call omp_set_num_threads(4)
c phase = -1
phase = 11 ! only reordering and symbolic factorization
CALL pardiso(pt,maxfct,mnum,mtype,phase,n,a,ia,ja,
1 idum,nrhs,iparm,msglvl,ddum,ddum,error)
! WRITE(*,*) 'Reordering completed ... '
IF (error .NE. 0) THEN
WRITE(*,*) 'The following ERROR was detected: ', error
STOP
END IF
c WRITE(*,*) 'Number of nonzeros in factors = ',iparm(18)
c WRITE(*,*) 'Number of factorization MFLOPS = ',iparm(19)
C.. Factorization.
phase=22 ! only factorization
CALL pardiso(pt,maxfct,mnum,mtype,phase,n,a,ia,ja,
1 idum,nrhs,iparm,msglvl,ddum,ddum,error)
! WRITE(*,*) 'Factorization completed ... '
IF (error .NE. 0)THEN
WRITE(*,*)'The following ERROR was detected: ',error
STOP
ENDIF
C.. Back substitution and iterative refinement
iparm(8) = 2 ! max numbers of iterative refinement steps
phase = 33 ! only factorization
do i=1,n
b(i)=1.d0
end do
CALL pardiso(pt, maxfct, mnum, mtype, phase, n, a, ia, ja,
1 idum, nrhs,iparm, msglvl, b, x, error)
! WRITE(*,*) 'Solve completed ... '
b=x
C.. Termiation and release of memory
c phase = -1 ! release internal memory
phase=11
CALL pardiso(pt,maxfct,mnum,mtype,phase,n,ddum,idum,idum,
1 idum,nrhs,iparm,msglvl,ddum,ddum,error)
IF (error .NE. 0)THEN
WRITE(*,*)'The following ERROR on release stage was detected: ',error
STOP
ENDIF
return
end
C
Thank you for all your reply.
Link Copied
