[Fortran] 纯文本查看 复制代码
MODULE MMkl_pardiso
USE Constants
IMPLICIT NONE
SAVE
PRIVATE
PUBLIC :: cal_coomatrix_mkl_pardiso ! 非对称阵
CONTAINS
SUBROUTINE cal_coomatrix_mkl_pardiso(nnz, acoo, rowind, colind, m, n, b, x)
IMPLICIT NONE
INTEGER(KINT) :: job(8)
INTEGER(KINT) :: nnz
REAL(KREAL), INTENT(INOUT) :: acoo(:)
INTEGER(KINT), INTENT(INOUT) :: rowind(:)
INTEGER(KINT), INTENT(INOUT) :: colind(:)
INTEGER(KINT), INTENT(IN) :: m
INTEGER(KINT), INTENT(INOUT) :: n
REAL(KREAL), ALLOCATABLE :: acsr(:)
INTEGER(KINT), ALLOCATABLE :: ja(:)
INTEGER(KINT), ALLOCATABLE :: ia(:)
INTEGER(KINT) :: info
INTEGER(KINT) :: acsr_dim
!
REAL(KREAL), INTENT(INOUT) :: b(:)
REAL(KREAL), INTENT(INOUT) :: x(:)
!
!=====================the variables of pardiso====================
integer(kind=4) :: pt(64) ! kind=8:x64; kind=4:win32
integer :: maxfct, mnum, mtype, phase, nrhs, error, msglvl
integer, allocatable :: perm(:)
integer :: iparm(64)
integer, external :: mkl_get_max_threads
!=====================the variables of pardiso====================
job(1) = 1 ! if job(1)=1, the matrix in the coordinate format is converted to the CSR format.
job(2) = 1 ! if job(2)=1, one-based indexing for the matrix in CSR format is used.
job(3) = 1 ! if job(3)=1, one-based indexing for the matrix in coordinate format is used.
job(4) = 0
job(5) = nnz ! nzmax - maximum number of the non-zero elements allowed if job(1)=0.
job(6) = 0 ! all arrays acsr, ja, ia are filled in for the output storage.
acsr_dim = nnz
! n = nnz
ALLOCATE( acsr(acsr_dim) )
ALLOCATE( ja(acsr_dim) )
ALLOCATE( ia(m+1))
ALLOCATE( perm(n))
acsr = REAL_ZERO
ja = INT_ZERO
ia = INT_ZERO
CALL convert_coo_to_csr(job, n, acsr, ja, ia, nnz, acoo, rowind, colind, info)
!===========================solve equations=============================
mtype = 11 ! mtype = -2: real and symmetric indefinite,
! mtype = 2: real and symmetric positive definite
! mtype = 11: real and nonsymmetric
! Initialize Intel MKL PARDISO with default parameters in accordance with the matrix type.
call pardisoinit (pt, mtype, iparm)
! pt: Solver internal data address pointer
maxfct = 1 ! Maximal number of factors in memory
! Generally used value is 1
mnum = 1 ! The number of matrix (from 1 to maxfct) to solve
! Generally used value is 1
! mtype =
! 1 Real and structurally symmetric
! 2 Real and symmetric positive definite
! -2 Real and symmetric indefinite
! 3 Complex and structurally symmetric
! 4 Complex and Hermitian positive definite
! -4 Complex and Hermitian indefinite
! 6 Complex and symmetric matrix
! 11 Real and nonsymmetric matrix
! 13 Complex and nonsymmetric matrix
! phase =
! 11 Analysis
! 12 Analysis, numerical factorization
! 13 Analysis, numerical factorization, solve
! 22 Numerical factorization
! 23 Numerical factorization, solve
! 33 Solve, iterative refinement
! 331 phase=33, but only forward substitution
! 332 phase=33, but only diagonal substitution
! 333 phase=33, but only backward substitution
! 0 Release internal memory for L and U of the matrix number mnum[attach]2427[/attach]
! -1 Release all internal memory for all matrices
! n Number of equations in the sparse linear system A*X = B
! a(*) Contains the non-zero elements of the coefficient matrix A
! ia(n+1) rowIndex array in CSR3 format
! ja(*) columns array in CSR3 format
! perm(n) Holds the permutation vector of size n or specifies elements used for computing a partial solution
! nrhs Number of right-hand sides that need to be solved for
! iparm(64) This array is used to pass various parameters to Intel MKL PARDISO
! and to return some useful information after execution of the solver
! (see pardiso iparm Parameter for more details)
! msglvl 0 Intel MKL PARDISO generates no output
! 1 Intel MKL PARDISO prints statistical information
! b(n*nrhs) Right-hand side vectors
! x(n*nrhs) Solution vectors
! error Error indicator
! call pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, perm, nrhs, iparm, msglvl, b, x, error)
perm = 0
nrhs = 1 ! Number of right-hand sides that need to be solved for.
x = 0.0D0
! iparm(1) = 0 ! iparm use default values
! iparm(3): parallel threads are obtained by MKL. In general, the number of cpus is the same as the number of threads
!call mkl_set_num_threads(4)
iparm(3) = mkl_get_max_threads()
error = 0
msglvl = 0
phase = 12 ! LU decompose, Analysis, numerical factorization
call pardiso( pt, maxfct, mnum, mtype, phase, m, acsr, ia, ja, perm, nrhs, iparm, msglvl, b, x, error )
phase = 33 ! Solve, iterative refinement
call pardiso( pt, maxfct, mnum, mtype, phase, m, acsr, ia, ja, perm, nrhs, iparm, msglvl, b, x, error )
phase = -1 ! Release all internal memory for all matrices
call pardiso (pt, maxfct, mnum, mtype, phase, n, acsr, ia, ja, perm, nrhs, iparm, msglvl, b, x, error)
DEALLOCATE( acsr )
DEALLOCATE( ja )
DEALLOCATE( ia)
DEALLOCATE( perm)
RETURN
END SUBROUTINE cal_coomatrix_mkl_pardiso
SUBROUTINE convert_coo_to_csr(job, n, acsr, ja, ia, nnz, acoo, rowind, colind, info)
IMPLICIT NONE
INTEGER(KINT), INTENT(IN OUT) :: job(8)
INTEGER(KINT), INTENT(IN OUT) :: n
REAL(KREAL) :: acsr(:)
INTEGER(KINT) :: ja(:)
INTEGER(KINT) :: ia(:)
INTEGER(KINT) :: nnz
REAL(KREAL), INTENT(IN OUT) :: acoo(:)
INTEGER(KINT), INTENT(IN OUT) :: rowind(:)
INTEGER(KINT), INTENT(IN OUT) :: colind(:)
INTEGER(KINT) :: info
call mkl_dcsrcoo(job, n, acsr, ja, ia, nnz, acoo, rowind, colind, info)
RETURN
END SUBROUTINE convert_coo_to_csr
END MODULE MMkl_pardiso 
F 币
贡献
权杖
窥视卡
发表于 2019-11-1 08:54:23
微信
收藏
点赞
点踩
提升卡
置顶卡
变色卡
千斤顶
显身卡
发表于 2019-11-1 15:25:13
楼主