Fortran Coder

标题: intel openmp和intel mkl冲突 [打印本页]

作者: ksfengjia    时间: 2019-11-1 08:54
标题: 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损坏
现状:已搜索查询很多网站,还是没有发现如何规避此问题,求助大神解答。

报错

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




作者: li913    时间: 2019-11-1 15:25
本帖最后由 li913 于 2019-11-1 15:30 编辑

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

作者: ksfengjia    时间: 2019-11-19 13:50
integer(kind=4)           :: pt(64)                                         ! kind=8:x64; kind=4:win32
最后发现并非mkl和openmp冲突,是由于64位编译下pt需要用kind=8,问题解决





欢迎光临 Fortran Coder (http://bbs.fcode.cn/) Powered by Discuz! X3.2