hecmw_MKL_wrapper.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
#ifdef HECMW_WITH_MKL
include 'mkl_pardiso.f90'
#endif
 
module m_hecmw_mkl_wrapper
  use hecmw_util
  use m_sparse_matrix
 
#ifdef HECMW_WITH_MKL
  use mkl_pardiso
#endif
 
  implicit none
 
  private                         ! default
  public :: hecmw_mkl_wrapper     ! only entry point of Parallel Direct Solver is public
 
  logical, save :: INITIALIZED = .false.
#ifdef HECMW_WITH_MKL
  type(MKL_PARDISO_HANDLE) :: pt(64)
#endif
  integer maxfct, mnum, mtype, nrhs, msglvl
  integer idum(1), i
  data nrhs /1/, maxfct /1/, mnum /1/
 
  integer, save :: iparm(64)
 
contains
 
  subroutine hecmw_mkl_wrapper(spMAT, phase_start, solx, istat)
    implicit none
    type (sparse_matrix), intent(inout)      :: spmat
    integer(kind=kint), intent(in)           :: phase_start
    integer(kind=kint), intent(out)          :: istat
    real(kind=kreal), pointer, intent(inout) :: solx(:)
 
#ifdef HECMW_WITH_MKL
 
    integer(kind=kint) :: myrank, phase
    real(kind=kreal)   :: t2,t3,t4,t5
 
    myrank=hecmw_comm_get_rank()
 
    if (.not. initialized) then
      do i=1,64
        pt(i)%dummy = 0
      enddo
      iparm(:) = 0
      iparm(1) = 1 ! no solver default
      iparm(2) = 3 ! fill-in reordering from METIS
      iparm(3) = 1
      iparm(10) = 13 ! perturb the pivot elements with 1E-13
      iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
      iparm(13) = 1 ! maximum weighted matching algorithm is switched-off
      iparm(18) = -1
      iparm(19) = -1
      msglvl = 0 ! print statistical information
      if (spmat%loglevel > 1) msglvl = 1
      initialized = .true.
    endif
 
    if ( phase_start == 1 ) then
      phase = -1
      call pardiso(pt, maxfct, mnum, mtype, phase, spmat%N, spmat%A, spmat%IRN, spmat%JCN, idum,  &
        nrhs, iparm, msglvl, spmat%RHS, solx, istat)
      if (istat < 0) then
        write(*,*) 'ERROR: MKL returned with error', istat
        return
      endif
 
      if (spmat%symtype == sparse_matrix_symtype_spd) then
        mtype = 2
      else if (spmat%symtype == sparse_matrix_symtype_sym) then
        mtype = -2
      else
        mtype = 11
      endif
 
      ! ANALYSIS
      t2=hecmw_wtime()
 
      phase = 11
      call pardiso(pt, maxfct, mnum, mtype, phase, spmat%N, spmat%A, spmat%IRN, spmat%JCN, idum,  &
        nrhs, iparm, msglvl, spmat%RHS, solx, istat)
      if (istat < 0) then
        write(*,*) 'ERROR: MKL returned with error', istat
        return
      endif
 
      t3=hecmw_wtime()
      if (myrank==0 .and. spmat%loglevel > 0) &
        write(*,'(A)') ' [Pardiso]: Analysis completed.'
      if (myrank==0 .and. spmat%timelog > 0) &
        write(*,'(A,f10.3)') ' [Pardiso]: Analysis time(sec)=',t3-t2
    endif
 
    ! FACTORIZATION
    if ( phase_start .le. 2 ) then
      t3=hecmw_wtime()
 
      phase = 22
      call pardiso(pt, maxfct, mnum, mtype, phase, spmat%N, spmat%A, spmat%IRN, spmat%JCN, idum,  &
        nrhs, iparm, msglvl, spmat%RHS, solx, istat)
      if (istat < 0) then
        write(*,*) 'ERROR: MKL returned with error', istat
        return
      endif
 
      t4=hecmw_wtime()
      if (myrank==0 .and. spmat%loglevel > 0) &
        write(*,'(A)') ' [Pardiso]: Factorization completed.'
      if (myrank==0 .and. spmat%timelog > 0) &
        write(*,'(A,f10.3)') ' [Pardiso]: Factorization time(sec)=',t4-t3
    endif
 
    ! SOLUTION
    t4=hecmw_wtime()
    phase = 33
    call pardiso(pt, maxfct, mnum, mtype, phase, spmat%N, spmat%A, spmat%IRN, spmat%JCN, idum,  &
      nrhs, iparm, msglvl, spmat%RHS, solx, istat)
    if (istat < 0) then
      write(*,*) 'ERROR: MKL returned with error', istat
      return
    endif
 
    t5=hecmw_wtime()
    if (myrank==0 .and. spmat%loglevel > 0) &
       write(*,'(A)') ' [Pardiso]: Solution completed.'
    if (myrank==0 .and. spmat%timelog > 0) &
       write(*,'(A,f10.3)') ' [Pardiso]: Solution time(sec)=',t5-t4
 
#else
    stop "MKL Pardiso not available"
#endif
  end subroutine hecmw_mkl_wrapper
 
end module m_hecmw_mkl_wrapper