ksfengjia 发表于 2019-11-1 08:54:23

intel openmp和intel mkl冲突

环境:vs2015+ivf2017
背景:程序采用mkl库pardiso并行进行矩阵计算,准备采用openmp对其他代码进行局部优化。
环境变量:NUMBER_OF_PROCESSORS=4,OMP_STACKSIZE=1000000
问题:单打开/Qmkl:parallel,程序正常运行
          同时打开/Qopenmp_stubs和/Qmkl:parallel编译无问题,正常运行没有问题
          同时打开/Qopenmp和/Qmkl:parallel编译无问题,在pardiso计算后,报stack损坏
现状:已搜索查询很多网站,还是没有发现如何规避此问题,求助大神解答。

报错

代码:
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
      !         -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   
   
    SUBROUTINEconvert_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   


li913 发表于 2019-11-1 15:25:13

本帖最后由 li913 于 2019-11-1 15:30 编辑

你给的代码不全,没法执行。试试关闭嵌套并行。

ksfengjia 发表于 2019-11-19 13:50:08

integer(kind=4)         :: pt(64)                                       ! kind=8:x64; kind=4:win32
最后发现并非mkl和openmp冲突,是由于64位编译下pt需要用kind=8,问题解决
页: [1]
查看完整版本: intel openmp和intel mkl冲突