hecmw_solver_direct_MUMPS.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
module hecmw_solver_direct_mumps
  use hecmw_util
  use m_sparse_matrix
  use m_sparse_matrix_hec
  use m_hecmw_mumps_wrapper
  use hecmw_matrix_ass
  use hecmw_matrix_dump
 
  private
  public :: hecmw_solve_direct_mumps
 
  logical, save :: INITIALIZED = .false.
  type (sparse_matrix), save :: spMAT
 
contains
 
  subroutine hecmw_solve_direct_mumps(hecMESH,hecMAT)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix    ), intent(inout) :: hecmat
    integer(kind=kint) :: spmat_type
    integer(kind=kint) :: spmat_symtype
    integer(kind=kint) :: mumps_job
    integer(kind=kint) :: istat,myrank
    real(kind=kreal) :: t1,t2,t3
 
    call hecmw_mat_dump(hecmat, hecmesh)
 
    t1=hecmw_wtime()
    myrank=hecmw_comm_get_rank()
 
    if (initialized .and. hecmat%Iarray(98) .eq. 1) then
      mumps_job=-2
      call hecmw_mumps_wrapper(spmat, mumps_job, istat)
      if (istat < 0) then
        write(*,*) 'ERROR: MUMPS returned with error', istat
        stop
      endif
      call sparse_matrix_finalize(spmat)
      initialized = .false.
    endif
 
    if (.not. initialized) then
      spmat_type = sparse_matrix_type_coo
      if (hecmat%symmetric) then
        spmat_symtype = sparse_matrix_symtype_spd
      else
        spmat_symtype = sparse_matrix_symtype_asym
      end if
      call sparse_matrix_set_type(spmat, spmat_type, spmat_symtype)
      mumps_job = -1
      call hecmw_mumps_wrapper(spmat, mumps_job, istat)
      if (istat < 0) then
        write(*,*) 'ERROR: MUMPS returned with error', istat
        stop
      endif
      initialized = .true.
      hecmat%Iarray(98) = 1
    endif
 
    !* Flag to activate symbolic factorization: 1(yes) 0(no)  hecMESH%Iarray(98)
    !* Flag to activate numeric  factorization: 1(yes) 0(no)  hecMESH%Iarray(97)
 
    if (hecmat%Iarray(98) .eq. 1) then
      ! ANALYSIS and FACTORIZATION
      call sparse_matrix_hec_init_prof(spmat, hecmat, hecmesh)
      call sparse_matrix_hec_set_vals(spmat, hecmat)
      !call sparse_matrix_dump(spMAT)
      mumps_job=4
      call hecmw_mumps_wrapper(spmat, mumps_job, istat)
      if (istat < 0) then
        write(*,*) 'ERROR: MUMPS returned with error', istat
        stop
      endif
      if (myrank==0 .and. (spmat%iterlog > 0 .or. spmat%timelog > 0)) &
           write(*,*) ' [MUMPS]: Analysis and Factorization completed.'
      hecmat%Iarray(98) = 0
      hecmat%Iarray(97) = 0
    endif
    if (hecmat%Iarray(97) .eq. 1) then
      ! FACTORIZATION
      call sparse_matrix_hec_set_vals(spmat, hecmat)
      !call sparse_matrix_dump(spMAT)
      mumps_job=2
      call hecmw_mumps_wrapper(spmat, mumps_job, istat)
      if (istat < 0) then
        write(*,*) 'ERROR: MUMPS returned with error', istat
        stop
      endif
      if (myrank==0 .and. (spmat%iterlog > 0 .or. spmat%timelog > 0)) &
           write(*,*) ' [MUMPS]: Factorization completed.'
      hecmat%Iarray(97) = 0
    endif
 
    t2=hecmw_wtime()
 
    ! SOLUTION
    call sparse_matrix_hec_set_rhs(spmat, hecmat)
    mumps_job=3
    call hecmw_mumps_wrapper(spmat, mumps_job, istat)
    if (istat < 0) then
      write(*,*) 'ERROR: MUMPS returned with error', istat
      stop
    endif
    call sparse_matrix_hec_get_rhs(spmat, hecmat)
    if (myrank==0 .and. (spmat%iterlog > 0 .or. spmat%timelog > 0)) &
         write(*,*) ' [MUMPS]: Solution completed.'
 
    t3=hecmw_wtime()
    if (myrank==0 .and. spmat%timelog > 0) then
      write(*,*) 'setup time : ',t2-t1
      write(*,*) 'solve time : ',t3-t2
    endif
 
    call hecmw_mat_dump_solution(hecmat)
 
    !call sparse_matrix_finalize(spMAT)
  end subroutine hecmw_solve_direct_mumps
 
end module hecmw_solver_direct_mumps