[数学库] 求教Pardiso 6.0的用法

subroutine RealPardiso(matrix,b,x,N,error)
    use libpardiso600
    implicit none
    integer                   :: n  !// depend on your equation
    integer                   :: i, j, mm, tmp, nn, fileid, first, num, ncolumn, count
    real(kind=8)              :: matrix(n,n), b(n), x(n),dparm(64)
    real(kind=8), allocatable :: aa(:)
    integer, allocatable      :: ia(:), ja(:)
  
    !// =====================the variables of pardiso====================
    integer(kind=8)           :: pt(64)   !// kind=8:x64; kind=4:win32
    integer                   :: maxfct, mtype, phase, nrhs, error, msglvl, mnum,  solver
    integer, allocatable      :: perm(:)
    integer                   :: iparm(64)
    integer, external         :: mkl_get_max_threads
    !// =====================the variables of pardiso====================
  
    
  
    !// input data
    allocate( perm( n ) )
 
    !// =========================store data by CSR format=======================

    num=0
    do I=1,N
        do J=1,N
            if ( abs(matrix(i,j)) /= 0.0d0 ) then
                num=num+1
            endif
        enddo
    enddo

    allocate( aa(num), ja(num), ia(n+1) )  !// store the no-zero element

    ia(1)=1
    count=1
    do I=1,N
        ncolumn=0
        do J=1,N
            if ( abs(matrix(I,J)) /= 0.0d0 ) then
                ncolumn=ncolumn+1
                ja(count)=J
                aa(count)=abs(matrix(I,J))
                count=count+1
            endif
        enddo
        ia(I+1)=count
    enddo

    !// ===========================solve equations=============================
    pt = 0  !// pointer initialization
    maxfct = 1; mnum = 1; mtype = 11  !// mtype = -2: symmetric nonpositive definite matrix, mtype = 2: symmetric positive definite matrix
    perm = 0; nrhs = 1
    x = 0.d0
    iparm(1) = 0  !// iparm use default values
    iparm(3) = mkl_get_max_threads()  !// iparm(3): parallel threads are obtained by MKL. In general, the number of cpus is the same as the number of threads
    error = 0  
    msglvl = 1
    solver=0
    
!     iparm(32)=1
!     iparm(33)=1
    phase = 12  !// LU decompose
    call pardiso( pt, maxfct, mnum, mtype, phase, n, aa, ia, ja, perm, nrhs, iparm, msglvl, b, x, error,dparm  )
!     write(*,*)dparm(32),dparm(33),iparm(33)
    
    
    iparm(8)=100
    iparm(32)=1
!     iparm(33)=1
    solver=0
    phase = 33  !// solve equations
    call pardiso( pt, maxfct, mnum, mtype, phase, n, aa, ia, ja, perm, nrhs, iparm, msglvl, b, x, error,dparm )
    write(*,*)iparm(7),iparm(20),dparm(35)
    do I=1,64
        write(*,*)dparm(I)
    enddo

    phase = -1
    call pardiso (pt, maxfct, mnum, mtype, phase, n, aa, ia, ja, perm, nrhs, iparm, msglvl, b, x, error,dparm)
    
    deallocate( aa, ja, ia, perm )
    !// ===========================stop solving equations=============================
  
end