Question

如何调用MKL库

Answer 1

SUBROUTINE pardiso_REDU(JA,A,X,B,IA,N,MAXA,NRHS)

IMPLICIT NONE

C.. Internal solver memory pointer for 64-bit architectures

C.. INTEGER*8 pt(64)

C.. Internal solver memory pointer for 32-bit architectures

C.. INTEGER*4 pt(64)

C.. This is OK in both cases

INTEGER*8 pt(64)

C.. All other variables

INTEGER maxa,maxfct,mnum,mtype,phase,n,nrhs,error,msglvl

INTEGER iparm(64)

INTEGER ia(*)

INTEGER ja(*)

REAL*8 a(*)

REAL*8 b(*)

REAL*8 x(*)

INTEGER i, idum

REAL*8 waltime1, waltime2, ddum

C.. Fill all arrays containing matrix data.

DATA maxfct /1/, mnum /1/

integer omp_get_max_threads

external omp_get_max_threads

C..

C.. Set up PARDISO control parameter

C..

do i=1,n+1

ia(i)=ia(i)+1

enddo

do i = 1, 64

iparm(i) = 0

end do

iparm(1) = 1 ! no solver default

iparm(2) = 2 ! fill-in reordering from METIS

iparm(3) = omp_get_max_threads() ! numbers of processors, value of OM
P_NUM_THREADS

iparm(4) = 61! 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

iparm(13) = 0 ! not in use

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

error = 0 ! initialize error flag

msglvl = 1 ! print statistical information

mtype = 11 ! real unsymmetric

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

write(*,*) 'NTHREADS =', iparm(3)

C.. Reordering and Symbolic Factorization, This step also allocates

C 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,

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

WRITE(*,*) 'Number of nonzeros in factors = ',iparm(18)

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

c do i = 1, n

c b(i) = 1.d0

c end do

CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja,

1 idum, nrhs, iparm, msglvl, b, x, error)

WRITE(*,*) 'Solve completed ... '

WRITE(*,*) 'The solution of the system is '

c DO i = 1, n

c WRITE(*,*) ' x(',i,') = ', x(i)

c END DO

C.. Termination and release of memory

phase = -1 ! release internal memory

CALL pardiso (pt, maxfct, mnum, mtype, phase, n, ddum, idum, idum,

1 idum, nrhs, iparm, msglvl, ddum, ddum, error)

RETURN

END