solve_LINEQ_mkl.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 m_solve_lineq_mkl_contact
  use hecmw_util
  use m_sparse_matrix
  use m_sparse_matrix_contact
  use hecmw_matrix_dump
  use m_hecmw_mkl_wrapper
  use m_hecmw_clustermkl_wrapper
  use hecmw_matrix_ass
  use hecmw_matrix_misc
 
  implicit none
 
  private
  public :: solve_lineq_mkl_contact_init
  public :: solve_lineq_mkl_contact
 
  logical, save :: NEED_ANALYSIS = .true.
  type (sparse_matrix), save :: spMAT
 
contains
 
  subroutine solve_lineq_mkl_contact_init(hecMESH,is_sym)
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    logical, intent(in) :: is_sym
 
    integer(kind=kint) :: spmat_type
    integer(kind=kint) :: spmat_symtype
    integer(kind=kint) :: i
 
    call sparse_matrix_finalize(spmat)
 
    if (is_sym) then
      spmat_symtype = sparse_matrix_symtype_sym
    else
      spmat_symtype = sparse_matrix_symtype_asym
    end if
    if(hecmesh%PETOT.GT.1) then
      spmat_type = sparse_matrix_type_coo
    else
      spmat_type = sparse_matrix_type_csr
    endif
    call sparse_matrix_set_type(spmat, spmat_type, spmat_symtype)
 
    need_analysis = .true.
  end subroutine solve_lineq_mkl_contact_init
 
  subroutine solve_lineq_mkl_contact(hecMESH,hecMAT,hecLagMAT,istat,conMAT)
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix    ), intent(inout) :: hecmat
    type (hecmwst_matrix_lagrange), intent(inout) :: heclagmat
    integer(kind=kint), intent(out) :: istat
    type (hecmwst_matrix), intent(in),optional :: conmat
 
    integer(kind=kint)  :: phase_start
    real(kind=kreal)    :: t1,t2
    integer(kind=kint)  :: mpc_method
 
    t1=hecmw_wtime()
 
    mpc_method = hecmw_mat_get_mpc_method(hecmat)
    if (mpc_method < 1 .or. 3 < mpc_method) then
      mpc_method = 1
      call hecmw_mat_set_mpc_method(hecmat,mpc_method)
    endif
    if (mpc_method /= 1) then
      write(*,*) 'ERROR: MPCMETHOD other than penalty is not available for DIRECTmkl solver', &
          ' in contact analysis without elimination'
      stop
    endif
    call hecmw_mat_ass_equation(hecmesh, hecmat)
    call hecmw_mat_ass_equation_rhs(hecmesh, hecmat)
 
    call hecmw_mat_dump(hecmat, hecmesh)
 
    if (need_analysis) then
      !constrtuct new structure
      call sparse_matrix_contact_init_prof(spmat, hecmat, heclagmat, hecmesh)
    endif
 
    !  ----  For Parallel Contact with Multi-Partition Domains
    if(present(conmat)) then
      call sparse_matrix_para_contact_set_vals(spmat, hecmat, heclagmat, conmat)
      call sparse_matrix_para_contact_set_rhs(spmat, hecmat, heclagmat, conmat)
    else
      call sparse_matrix_contact_set_vals(spmat, hecmat, heclagmat)
      !call sparse_matrix_dump(spMAT)
      call sparse_matrix_contact_set_rhs(spmat, hecmat, heclagmat)
    endif
 
    t2=hecmw_wtime()
    if ( hecmw_comm_get_rank() == 0 .and. spmat%timelog > 0) then
      if( hecmesh%PETOT .GT. 1 ) then
         write(*,'(A,f10.3)') ' [Cluster Pardiso]: Setup completed.            time(sec)=',t2-t1
      else
         write(*,'(A,f10.3)') ' [Pardiso]: Setup completed.          time(sec)=',t2-t1
      end if
    endif
 
    phase_start = 2
    if (need_analysis) then
      phase_start = 1
      need_analysis = .false.
    endif
 
    ! SOLVE
    if( hecmesh%PETOT.GT.1 ) then
      call hecmw_clustermkl_wrapper(spmat, phase_start, hecmat%X, istat)
      call sparse_matrix_contact_get_rhs(spmat, hecmat, heclagmat)
    else
      call hecmw_mkl_wrapper(spmat, phase_start, hecmat%X, istat)
      deallocate(spmat%rhs)
    endif
 
    call hecmw_mat_dump_solution(hecmat)
  end subroutine solve_lineq_mkl_contact
 
end module m_solve_lineq_mkl_contact