hecmw_matrix_misc.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_matrix_misc
  use hecmw_util
  use m_hecmw_comm_f
  implicit none
 
  private
  public :: hecmw_mat_clear
  public :: hecmw_mat_clear_b
  public :: hecmw_mat_init
  public :: hecmw_mat_finalize
  public :: hecmw_mat_copy_profile
  public :: hecmw_mat_copy_val
 
  public :: hecmw_mat_set_iter
  public :: hecmw_mat_get_iter
  public :: hecmw_mat_set_method
  public :: hecmw_mat_get_method
  public :: hecmw_mat_set_precond
  public :: hecmw_mat_get_precond
  public :: hecmw_mat_set_nset
  public :: hecmw_mat_get_nset
  public :: hecmw_mat_set_iterpremax
  public :: hecmw_mat_get_iterpremax
  public :: hecmw_mat_set_nrest
  public :: hecmw_mat_get_nrest
  public :: hecmw_mat_set_nbfgs
  public :: hecmw_mat_get_nbfgs
  public :: hecmw_mat_set_scaling
  public :: hecmw_mat_get_scaling
  public :: hecmw_mat_set_penalized
  public :: hecmw_mat_get_penalized
  public :: hecmw_mat_set_penalized_b
  public :: hecmw_mat_get_penalized_b
  public :: hecmw_mat_set_mpc_method
  public :: hecmw_mat_get_mpc_method
  public :: hecmw_mat_set_estcond
  public :: hecmw_mat_get_estcond
  public :: hecmw_mat_set_recompute_residual
  public :: hecmw_mat_get_recompute_residual
  public :: hecmw_mat_set_contact_elim
  public :: hecmw_mat_get_contact_elim
  public :: hecmw_mat_set_iterlog
  public :: hecmw_mat_get_iterlog
  public :: hecmw_mat_set_timelog
  public :: hecmw_mat_get_timelog
  public :: hecmw_mat_set_loglevel
  public :: hecmw_mat_get_loglevel
  public :: hecmw_mat_set_dump
  public :: hecmw_mat_get_dump
  public :: hecmw_mat_set_dump_exit
  public :: hecmw_mat_get_dump_exit
  public :: hecmw_mat_set_usejad
  public :: hecmw_mat_get_usejad
  public :: hecmw_mat_set_matvec_impl
  public :: hecmw_mat_get_matvec_impl
  public :: hecmw_mat_set_precond_impl
  public :: hecmw_mat_get_precond_impl
  public :: hecmw_mat_set_ncolor_in
  public :: hecmw_mat_get_ncolor_in
  public :: hecmw_mat_set_maxrecycle_precond
  public :: hecmw_mat_get_maxrecycle_precond
  public :: hecmw_mat_get_nrecycle_precond
  public :: hecmw_mat_reset_nrecycle_precond
  public :: hecmw_mat_incr_nrecycle_precond
  public :: hecmw_mat_set_flag_numfact
  public :: hecmw_mat_get_flag_numfact
  public :: hecmw_mat_set_flag_symbfact
  public :: hecmw_mat_get_flag_symbfact
  public :: hecmw_mat_clear_flag_symbfact
  public :: hecmw_mat_set_solver_type
  public :: hecmw_mat_get_solver_type
 
  public :: hecmw_mat_set_method2
  public :: hecmw_mat_get_method2
  public :: hecmw_mat_set_flag_converged
  public :: hecmw_mat_get_flag_converged
  public :: hecmw_mat_set_flag_diverged
  public :: hecmw_mat_get_flag_diverged
 
  public :: hecmw_mat_set_solver_opt
  public :: hecmw_mat_get_solver_opt
  public :: hecmw_mat_set_solver_ropt
  public :: hecmw_mat_get_solver_ropt
 
  public :: hecmw_mat_set_resid
  public :: hecmw_mat_get_resid
  public :: hecmw_mat_set_sigma_diag
  public :: hecmw_mat_get_sigma_diag
  public :: hecmw_mat_set_omega
  public :: hecmw_mat_get_omega
  public :: hecmw_mat_set_thresh
  public :: hecmw_mat_get_thresh
  public :: hecmw_mat_set_filter
  public :: hecmw_mat_get_filter
  public :: hecmw_mat_set_penalty
  public :: hecmw_mat_get_penalty
  public :: hecmw_mat_set_penalty_alpha
  public :: hecmw_mat_get_penalty_alpha
 
  public :: hecmw_matvec_impl_bsr
  public :: hecmw_matvec_impl_csr
  public :: hecmw_matvec_impl_sblas
  public :: hecmw_precond_impl_bsr
  public :: hecmw_precond_impl_csr
  public :: hecmw_precond_impl_sblas
  public :: hecmw_matvec_impl_default
  public :: hecmw_precond_impl_default
 
  public :: hecmw_mat_diag_max
  public :: hecmw_mat_diag
  public :: hecmw_mat_recycle_precond_setting
  public :: hecmw_mat_substitute
  public :: hecmw_mat_integrate
 
  integer, parameter :: IDX_I_ITER               = 1
  integer, parameter :: IDX_I_METHOD             = 2
  integer, parameter :: IDX_I_PRECOND            = 3
  integer, parameter :: IDX_I_NSET               = 4
  integer, parameter :: IDX_I_ITERPREMAX         = 5
  integer, parameter :: IDX_I_NREST              = 6
  integer, parameter :: IDX_I_NBFGS              = 60
  integer, parameter :: IDX_I_SCALING            = 7
  integer, parameter :: IDX_I_PENALIZED          = 11
  integer, parameter :: IDX_I_PENALIZED_B        = 12
  integer, parameter :: IDX_I_MPC_METHOD         = 13
  integer, parameter :: IDX_I_ESTCOND            = 14
  integer, parameter :: IDX_I_CONTACT_ELIM       = 15
  integer, parameter :: IDX_I_RECOMPUTE_RESIDUAL = 16
  integer, parameter :: IDX_I_ITERLOG            = 21
  integer, parameter :: IDX_I_TIMELOG            = 22
  integer, parameter :: IDX_I_LOGLEVEL           = 24   ! 23 is steplog (svIarray)
  integer, parameter :: IDX_I_DUMP               = 31
  integer, parameter :: IDX_I_DUMP_EXIT          = 32
  integer, parameter :: IDX_I_USEJAD             = 33
  integer, parameter :: IDX_I_NCOLOR_IN          = 34
  integer, parameter :: IDX_I_MAXRECYCLE_PRECOND = 35
  integer, parameter :: IDX_I_MATVEC_IMPL        = 36
  integer, parameter :: IDX_I_PRECOND_IMPL       = 37
  integer, parameter :: IDX_I_NRECYCLE_PRECOND   = 96
  integer, parameter :: IDX_I_FLAG_NUMFACT       = 97
  integer, parameter :: IDX_I_FLAG_SYMBFACT      = 98
  integer, parameter :: IDX_I_SOLVER_TYPE        = 99
 
  integer, parameter :: IDX_I_METHOD2            = 8
  integer, parameter :: IDX_I_FLAG_CONVERGED     = 81
  integer, parameter :: IDX_I_FLAG_DIVERGED      = 82
 
  integer, parameter :: IDX_I_SOLVER_OPT_S       = 41
  integer, parameter :: IDX_I_SOLVER_OPT_E       = 50
 
  integer, parameter :: IDX_R_RESID         = 1
  integer, parameter :: IDX_R_SIGMA_DIAG    = 2
  integer, parameter :: IDX_R_OMEGA         = 3
  integer, parameter :: IDX_R_THRESH        = 4
  integer, parameter :: IDX_R_FILTER        = 5
  integer, parameter :: IDX_R_PENALTY       = 11
  integer, parameter :: IDX_R_PENALTY_ALPHA = 12
  ! real-valued solver options, mirroring the integer SOLVER_OPT block (41:50)
  integer, parameter :: IDX_R_SOLVER_OPT_S  = 41
  integer, parameter :: IDX_R_SOLVER_OPT_E  = 50
 
  ! storage format the matvec and the preconditioner work in, selected by
  ! IDX_I_MATVEC_IMPL / IDX_I_PRECOND_IMPL.  BSR is hecmwST_matrix's own format and is
  ! the only one every build carries; the others keep a private copy in their own layout
  ! and are built only for the architecture cmake -DARCH= names, so the dispatchers fall
  ! back to BSR for a format this build does not have.  The names match the MATRIXFORMAT
  ! values accepted by the cnt file.
  integer(kind=kint), parameter :: hecmw_matvec_impl_bsr    = 0
  integer(kind=kint), parameter :: hecmw_matvec_impl_csr    = 1
  integer(kind=kint), parameter :: hecmw_matvec_impl_sblas  = 2
  integer(kind=kint), parameter :: hecmw_precond_impl_bsr   = 0
  integer(kind=kint), parameter :: hecmw_precond_impl_csr   = 1
  integer(kind=kint), parameter :: hecmw_precond_impl_sblas = 2
 
  ! cmake -DARCH= names the constants above through the format it defaults to; without
  ! it the build defaults to BSR.  Both initializers of hecmwST_matrix read these, so
  ! the fallback stays in one place.
#ifndef HECMW_ARCH_DEFAULT_MATVEC_IMPL
#define HECMW_ARCH_DEFAULT_MATVEC_IMPL HECMW_MATVEC_IMPL_BSR
#endif
#ifndef HECMW_ARCH_DEFAULT_PRECOND_IMPL
#define HECMW_ARCH_DEFAULT_PRECOND_IMPL HECMW_PRECOND_IMPL_BSR
#endif
  integer(kind=kint), parameter :: hecmw_matvec_impl_default  = hecmw_arch_default_matvec_impl
  integer(kind=kint), parameter :: hecmw_precond_impl_default = hecmw_arch_default_precond_impl
 
contains
 
  subroutine hecmw_mat_clear( hecMAT )
    type(hecmwst_matrix) :: hecmat
 
    hecmat%D = 0.0d0
    hecmat%AL = 0.0d0
    hecmat%AU = 0.0d0
    call hecmw_mat_set_penalized( hecmat, 0 )
    call hecmw_mat_set_penalty_alpha( hecmat, 0.d0 )
  end subroutine hecmw_mat_clear
 
  subroutine hecmw_mat_clear_b( hecMAT )
    type(hecmwst_matrix) :: hecmat
 
    hecmat%B = 0.0d0
    call hecmw_mat_set_penalized_b( hecmat, 0 )
  end subroutine hecmw_mat_clear_b
 
  subroutine hecmw_mat_init( hecMAT )
    type(hecmwst_matrix) :: hecmat
 
    call hecmw_nullify_matrix( hecmat )
 
    hecmat%Iarray = 0
    hecmat%Rarray = 0.d0
 
    call hecmw_mat_set_iter( hecmat, 100 )
    call hecmw_mat_set_method( hecmat, 1 )
    call hecmw_mat_set_precond( hecmat, 1 )
    call hecmw_mat_set_nset( hecmat, 0 )
    call hecmw_mat_set_iterpremax( hecmat, 1 )
    call hecmw_mat_set_nrest( hecmat, 10 )
    call hecmw_mat_set_nbfgs( hecmat, 0 )
    call hecmw_mat_set_scaling( hecmat, 0 )
    call hecmw_mat_set_iterlog( hecmat, 0 )
    call hecmw_mat_set_timelog( hecmat, 0 )
    call hecmw_mat_set_loglevel( hecmat, -1 )   ! -1 = unset: consumers fall back to their default
    call hecmw_mat_set_dump( hecmat, 0 )
    call hecmw_mat_set_dump_exit( hecmat, 0 )
    call hecmw_mat_set_usejad( hecmat, 0 )
    call hecmw_mat_set_matvec_impl( hecmat, hecmw_matvec_impl_default )
    call hecmw_mat_set_precond_impl( hecmat, hecmw_precond_impl_default )
    call hecmw_mat_set_ncolor_in( hecmat, 10 )
    call hecmw_mat_set_estcond( hecmat, 0 )
    call hecmw_mat_set_recompute_residual( hecmat, 0 )   ! 0 = unset: each iterative solver falls back to its own period
    call hecmw_mat_set_maxrecycle_precond( hecmat, 3 )
 
    call hecmw_mat_set_resid( hecmat, 1.d-8 )
    call hecmw_mat_set_sigma_diag( hecmat, 1.d0 )
    call hecmw_mat_set_omega( hecmat, 1.d0 )
    call hecmw_mat_set_thresh( hecmat, 0.10d0 )
    call hecmw_mat_set_filter( hecmat, 0.10d0 )
 
    call hecmw_mat_set_penalized( hecmat, 0 )
    call hecmw_mat_set_penalty( hecmat, 1.d+4 )
    call hecmw_mat_set_penalty_alpha( hecmat, 0.d0 )
    call hecmw_mat_set_mpc_method( hecmat, 0 )
 
    call hecmw_mat_reset_nrecycle_precond( hecmat )
    call hecmw_mat_set_flag_numfact( hecmat, 1 )
    call hecmw_mat_set_flag_symbfact( hecmat, 1 )
    call hecmw_mat_set_solver_type( hecmat, 1 )
  end subroutine hecmw_mat_init
 
  subroutine hecmw_mat_finalize( hecMAT )
    type(hecmwst_matrix) :: hecmat
    if (associated(hecmat%D)) deallocate(hecmat%D)
    if (associated(hecmat%B)) deallocate(hecmat%B)
    if (associated(hecmat%X)) deallocate(hecmat%X)
    if (associated(hecmat%AL)) deallocate(hecmat%AL)
    if (associated(hecmat%AU)) deallocate(hecmat%AU)
    if (associated(hecmat%indexL)) deallocate(hecmat%indexL)
    if (associated(hecmat%indexU)) deallocate(hecmat%indexU)
    if (associated(hecmat%itemL)) deallocate(hecmat%itemL)
    if (associated(hecmat%itemU)) deallocate(hecmat%itemU)
#ifdef _OPENACC
    if (associated(hecmat%A)) deallocate(hecmat%A)
    if (associated(hecmat%indexA)) deallocate(hecmat%indexA)
    if (associated(hecmat%itemA)) deallocate(hecmat%itemA)
#endif
  end subroutine hecmw_mat_finalize
 
  subroutine hecmw_mat_copy_profile( hecMATorg, hecMAT )
    type(hecmwst_matrix), intent(in) :: hecmatorg
    type(hecmwst_matrix), intent(inout) :: hecmat
    hecmat%N    = hecmatorg%N
    hecmat%NP   = hecmatorg%NP
    hecmat%NDOF = hecmatorg%NDOF
    hecmat%NPL  = hecmatorg%NPL
    hecmat%NPU  = hecmatorg%NPU
    allocate(hecmat%indexL(0:size(hecmatorg%indexL)-1))
    allocate(hecmat%indexU(0:size(hecmatorg%indexU)-1))
    allocate(hecmat%itemL (size(hecmatorg%itemL )))
    allocate(hecmat%itemU (size(hecmatorg%itemU )))
    allocate(hecmat%D (size(hecmatorg%D )))
    allocate(hecmat%AL(size(hecmatorg%AL)))
    allocate(hecmat%AU(size(hecmatorg%AU)))
    allocate(hecmat%B (size(hecmatorg%B )))
    allocate(hecmat%X (size(hecmatorg%X )))
    hecmat%indexL = hecmatorg%indexL
    hecmat%indexU = hecmatorg%indexU
    hecmat%itemL  = hecmatorg%itemL
    hecmat%itemU  = hecmatorg%itemU
    hecmat%D  = 0.d0
    hecmat%AL = 0.d0
    hecmat%AU = 0.d0
    hecmat%B  = 0.d0
    hecmat%X  = 0.d0
  end subroutine hecmw_mat_copy_profile
 
  subroutine hecmw_mat_copy_val( hecMATorg, hecMAT )
    type(hecmwst_matrix), intent(in) :: hecmatorg
    type(hecmwst_matrix), intent(inout) :: hecmat
    integer(kind=kint) :: ierr
    integer(kind=kint) :: i
    ierr = 0
    if (hecmat%N    /= hecmatorg%N) ierr = 1
    if (hecmat%NP   /= hecmatorg%NP) ierr = 1
    if (hecmat%NDOF /= hecmatorg%NDOF) ierr = 1
    if (hecmat%NPL  /= hecmatorg%NPL) ierr = 1
    if (hecmat%NPU  /= hecmatorg%NPU) ierr = 1
    if (ierr /= 0) then
      write(0,*) 'ERROR: hecmw_mat_copy_val: different profile'
      stop
    endif
    do i = 1, size(hecmat%D)
      hecmat%D(i)  = hecmatorg%D(i)
    enddo
    do i = 1, size(hecmat%AL)
      hecmat%AL(i) = hecmatorg%AL(i)
    enddo
    do i = 1, size(hecmat%AU)
      hecmat%AU(i) = hecmatorg%AU(i)
    enddo
  end subroutine hecmw_mat_copy_val
 
  subroutine hecmw_mat_set_iter( hecMAT, iter )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: iter
 
    hecmat%Iarray(idx_i_iter) = iter
  end subroutine hecmw_mat_set_iter
 
  function hecmw_mat_get_iter( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_iter
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_iter = hecmat%Iarray(idx_i_iter)
  end function hecmw_mat_get_iter
 
  subroutine hecmw_mat_set_method( hecMAT, method )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: method
 
    hecmat%Iarray(idx_i_method) = method
  end subroutine hecmw_mat_set_method
 
  function hecmw_mat_get_method( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_method
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_method = hecmat%Iarray(idx_i_method)
  end function hecmw_mat_get_method
 
  subroutine hecmw_mat_set_method2( hecMAT, method2 )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: method2
 
    hecmat%Iarray(idx_i_method2) = method2
  end subroutine hecmw_mat_set_method2
 
  function hecmw_mat_get_method2( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_method2
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_method2 = hecmat%Iarray(idx_i_method2)
  end function hecmw_mat_get_method2
 
  subroutine hecmw_mat_set_precond( hecMAT, precond )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: precond
 
    hecmat%Iarray(idx_i_precond) = precond
  end subroutine hecmw_mat_set_precond
 
  function hecmw_mat_get_precond( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_precond
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_precond = hecmat%Iarray(idx_i_precond)
  end function hecmw_mat_get_precond
 
  subroutine hecmw_mat_set_nset( hecMAT, nset )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: nset
 
    hecmat%Iarray(idx_i_nset) = nset
  end subroutine hecmw_mat_set_nset
 
  function hecmw_mat_get_nset( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_nset
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_nset = hecmat%Iarray(idx_i_nset)
  end function hecmw_mat_get_nset
 
  subroutine hecmw_mat_set_iterpremax( hecMAT, iterpremax )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: iterpremax
 
    if (iterpremax.lt.0) iterpremax= 0
    if (iterpremax.gt.4) iterpremax= 4
 
    hecmat%Iarray(idx_i_iterpremax) = iterpremax
  end subroutine hecmw_mat_set_iterpremax
 
  function hecmw_mat_get_iterpremax( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_iterpremax
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_iterpremax = hecmat%Iarray(idx_i_iterpremax)
  end function hecmw_mat_get_iterpremax
 
  subroutine hecmw_mat_set_nrest( hecMAT, nrest )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: nrest
 
    hecmat%Iarray(idx_i_nrest) = nrest
  end subroutine hecmw_mat_set_nrest
 
  function hecmw_mat_get_nrest( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_nrest
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_nrest = hecmat%Iarray(idx_i_nrest)
  end function hecmw_mat_get_nrest
 
  subroutine hecmw_mat_set_nbfgs( hecMAT, nbfgs )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: nbfgs
 
    hecmat%Iarray(idx_i_nbfgs) = nbfgs
  end subroutine hecmw_mat_set_nbfgs
 
  function hecmw_mat_get_nbfgs( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_nbfgs
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_nbfgs = hecmat%Iarray(idx_i_nbfgs)
  end function hecmw_mat_get_nbfgs
 
  subroutine hecmw_mat_set_scaling( hecMAT, scaling )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: scaling
 
    hecmat%Iarray(idx_i_scaling) = scaling
  end subroutine hecmw_mat_set_scaling
 
  function hecmw_mat_get_scaling( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_scaling
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_scaling = hecmat%Iarray(idx_i_scaling)
  end function hecmw_mat_get_scaling
 
  subroutine hecmw_mat_set_penalized( hecMAT, penalized )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: penalized
 
    hecmat%Iarray(idx_i_penalized) = penalized
  end subroutine hecmw_mat_set_penalized
 
  function hecmw_mat_get_penalized( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_penalized
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_penalized = hecmat%Iarray(idx_i_penalized)
  end function hecmw_mat_get_penalized
 
  subroutine hecmw_mat_set_penalized_b( hecMAT, penalized_b )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: penalized_b
 
    hecmat%Iarray(idx_i_penalized_b) = penalized_b
  end subroutine hecmw_mat_set_penalized_b
 
  function hecmw_mat_get_penalized_b( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_penalized_b
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_penalized_b = hecmat%Iarray(idx_i_penalized_b)
  end function hecmw_mat_get_penalized_b
 
  subroutine hecmw_mat_set_mpc_method( hecMAT, mpc_method )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: mpc_method
 
    hecmat%Iarray(idx_i_mpc_method) = mpc_method
  end subroutine hecmw_mat_set_mpc_method
 
  function hecmw_mat_get_mpc_method( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_mpc_method
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_mpc_method = hecmat%Iarray(idx_i_mpc_method)
  end function hecmw_mat_get_mpc_method
 
  function hecmw_mat_get_estcond( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_estcond
    type(hecmwst_matrix) :: hecmat
    hecmw_mat_get_estcond = hecmat%Iarray(idx_i_estcond)
  end function hecmw_mat_get_estcond
 
  subroutine hecmw_mat_set_estcond( hecMAT, estcond )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: estcond
    hecmat%Iarray(idx_i_estcond) = estcond
  end subroutine hecmw_mat_set_estcond
 
  function hecmw_mat_get_recompute_residual( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_recompute_residual
    type(hecmwst_matrix) :: hecmat
    hecmw_mat_get_recompute_residual = hecmat%Iarray(idx_i_recompute_residual)
  end function hecmw_mat_get_recompute_residual
 
  subroutine hecmw_mat_set_recompute_residual( hecMAT, recompute_residual )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: recompute_residual
    hecmat%Iarray(idx_i_recompute_residual) = recompute_residual
  end subroutine hecmw_mat_set_recompute_residual
 
  function hecmw_mat_get_contact_elim( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_contact_elim
    type(hecmwst_matrix) :: hecmat
    hecmw_mat_get_contact_elim = hecmat%Iarray(idx_i_contact_elim)
  end function hecmw_mat_get_contact_elim
 
  subroutine hecmw_mat_set_contact_elim( hecMAT, contact_elim )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: contact_elim
    hecmat%Iarray(idx_i_contact_elim) = contact_elim
  end subroutine hecmw_mat_set_contact_elim
 
  subroutine hecmw_mat_set_iterlog( hecMAT, iterlog )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: iterlog
 
    hecmat%Iarray(idx_i_iterlog) = iterlog
  end subroutine hecmw_mat_set_iterlog
 
  function hecmw_mat_get_iterlog( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_iterlog
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_iterlog = hecmat%Iarray(idx_i_iterlog)
  end function hecmw_mat_get_iterlog
 
  subroutine hecmw_mat_set_timelog( hecMAT, timelog )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: timelog
 
    hecmat%Iarray(idx_i_timelog) = timelog
  end subroutine hecmw_mat_set_timelog
 
  function hecmw_mat_get_timelog( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_timelog
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_timelog = hecmat%Iarray(idx_i_timelog)
  end function hecmw_mat_get_timelog
 
  subroutine hecmw_mat_set_loglevel( hecMAT, loglevel )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: loglevel
 
    hecmat%Iarray(idx_i_loglevel) = loglevel
  end subroutine hecmw_mat_set_loglevel
 
  function hecmw_mat_get_loglevel( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_loglevel
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_loglevel = hecmat%Iarray(idx_i_loglevel)
  end function hecmw_mat_get_loglevel
 
  function hecmw_mat_get_dump( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_dump
    type(hecmwst_matrix) :: hecmat
    hecmw_mat_get_dump = hecmat%Iarray(idx_i_dump)
  end function hecmw_mat_get_dump
 
  subroutine hecmw_mat_set_dump( hecMAT, dump_type )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: dump_type
    hecmat%Iarray(idx_i_dump) = dump_type
  end subroutine hecmw_mat_set_dump
 
  function hecmw_mat_get_dump_exit( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_dump_exit
    type(hecmwst_matrix) :: hecmat
    hecmw_mat_get_dump_exit = hecmat%Iarray(idx_i_dump_exit)
  end function hecmw_mat_get_dump_exit
 
  subroutine hecmw_mat_set_dump_exit( hecMAT, dump_exit )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: dump_exit
    hecmat%Iarray(idx_i_dump_exit) = dump_exit
  end subroutine hecmw_mat_set_dump_exit
 
  function hecmw_mat_get_usejad( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_usejad
    type(hecmwst_matrix) :: hecmat
    hecmw_mat_get_usejad = hecmat%Iarray(idx_i_usejad)
  end function hecmw_mat_get_usejad
 
  subroutine hecmw_mat_set_usejad( hecMAT, usejad )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: usejad
    hecmat%Iarray(idx_i_usejad) = usejad
  end subroutine hecmw_mat_set_usejad
 
  function hecmw_mat_get_matvec_impl( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_matvec_impl
    type(hecmwst_matrix) :: hecmat
    hecmw_mat_get_matvec_impl = hecmat%Iarray(idx_i_matvec_impl)
  end function hecmw_mat_get_matvec_impl
 
  subroutine hecmw_mat_set_matvec_impl( hecMAT, matvec_impl )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: matvec_impl
    hecmat%Iarray(idx_i_matvec_impl) = matvec_impl
  end subroutine hecmw_mat_set_matvec_impl
 
  function hecmw_mat_get_precond_impl( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_precond_impl
    type(hecmwst_matrix) :: hecmat
    hecmw_mat_get_precond_impl = hecmat%Iarray(idx_i_precond_impl)
  end function hecmw_mat_get_precond_impl
 
  subroutine hecmw_mat_set_precond_impl( hecMAT, precond_impl )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: precond_impl
    hecmat%Iarray(idx_i_precond_impl) = precond_impl
  end subroutine hecmw_mat_set_precond_impl
 
  function hecmw_mat_get_ncolor_in( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_ncolor_in
    type(hecmwst_matrix) :: hecmat
    hecmw_mat_get_ncolor_in = hecmat%Iarray(idx_i_ncolor_in)
  end function hecmw_mat_get_ncolor_in
 
  subroutine hecmw_mat_set_ncolor_in( hecMAT, ncolor_in )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: ncolor_in
    hecmat%Iarray(idx_i_ncolor_in) = ncolor_in
  end subroutine hecmw_mat_set_ncolor_in
 
  function hecmw_mat_get_maxrecycle_precond( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_maxrecycle_precond
    type(hecmwst_matrix) :: hecmat
    hecmw_mat_get_maxrecycle_precond = hecmat%Iarray(idx_i_maxrecycle_precond)
  end function hecmw_mat_get_maxrecycle_precond
 
  subroutine hecmw_mat_set_maxrecycle_precond( hecMAT, maxrecycle_precond )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: maxrecycle_precond
    if (maxrecycle_precond > 100) maxrecycle_precond = 100
    hecmat%Iarray(idx_i_maxrecycle_precond) = maxrecycle_precond
  end subroutine hecmw_mat_set_maxrecycle_precond
 
  function hecmw_mat_get_nrecycle_precond( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_nrecycle_precond
    type(hecmwst_matrix) :: hecmat
    hecmw_mat_get_nrecycle_precond = hecmat%Iarray(idx_i_nrecycle_precond)
  end function hecmw_mat_get_nrecycle_precond
 
  subroutine hecmw_mat_reset_nrecycle_precond( hecMAT )
    type(hecmwst_matrix) :: hecmat
    hecmat%Iarray(idx_i_nrecycle_precond) = 0
  end subroutine hecmw_mat_reset_nrecycle_precond
 
  subroutine hecmw_mat_incr_nrecycle_precond( hecMAT )
    type(hecmwst_matrix) :: hecmat
    hecmat%Iarray(idx_i_nrecycle_precond) = hecmat%Iarray(idx_i_nrecycle_precond) + 1
  end subroutine hecmw_mat_incr_nrecycle_precond
 
  function hecmw_mat_get_flag_numfact( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_flag_numfact
    type(hecmwst_matrix) :: hecmat
    hecmw_mat_get_flag_numfact = hecmat%Iarray(idx_i_flag_numfact)
  end function hecmw_mat_get_flag_numfact
 
  subroutine hecmw_mat_set_flag_numfact( hecMAT, flag_numfact )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: flag_numfact
    hecmat%Iarray(idx_i_flag_numfact) = flag_numfact
  end subroutine hecmw_mat_set_flag_numfact
 
  function hecmw_mat_get_flag_symbfact( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_flag_symbfact
    type(hecmwst_matrix) :: hecmat
    hecmw_mat_get_flag_symbfact = hecmat%Iarray(idx_i_flag_symbfact)
  end function hecmw_mat_get_flag_symbfact
 
  subroutine hecmw_mat_set_flag_symbfact( hecMAT, flag_symbfact )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: flag_symbfact
    hecmat%Iarray(idx_i_flag_symbfact) = flag_symbfact
  end subroutine hecmw_mat_set_flag_symbfact
 
  subroutine hecmw_mat_clear_flag_symbfact( hecMAT )
    type(hecmwst_matrix) :: hecmat
    hecmat%Iarray(idx_i_flag_symbfact) = 0
  end subroutine hecmw_mat_clear_flag_symbfact
 
  function hecmw_mat_get_solver_type( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_solver_type
    type(hecmwst_matrix) :: hecmat
    hecmw_mat_get_solver_type = hecmat%Iarray(idx_i_solver_type)
  end function hecmw_mat_get_solver_type
 
  subroutine hecmw_mat_set_solver_type( hecMAT, solver_type )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: solver_type
    hecmat%Iarray(idx_i_solver_type) = solver_type
  end subroutine hecmw_mat_set_solver_type
 
  subroutine hecmw_mat_set_flag_converged( hecMAT, flag_converged )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: flag_converged
    hecmat%Iarray(idx_i_flag_converged) = flag_converged
  end subroutine hecmw_mat_set_flag_converged
 
  function hecmw_mat_get_flag_converged( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_flag_converged
    type(hecmwst_matrix) :: hecmat
    hecmw_mat_get_flag_converged = hecmat%Iarray(idx_i_flag_converged)
  end function hecmw_mat_get_flag_converged
 
  subroutine hecmw_mat_set_flag_diverged( hecMAT, flag_diverged )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: flag_diverged
    hecmat%Iarray(idx_i_flag_diverged) = flag_diverged
  end subroutine hecmw_mat_set_flag_diverged
 
  function hecmw_mat_get_flag_diverged( hecMAT )
    integer(kind=kint) :: hecmw_mat_get_flag_diverged
    type(hecmwst_matrix) :: hecmat
    hecmw_mat_get_flag_diverged = hecmat%Iarray(idx_i_flag_diverged)
  end function hecmw_mat_get_flag_diverged
 
  subroutine hecmw_mat_set_solver_opt( hecMAT, solver_opt )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: solver_opt(:)
    integer(kind=kint) :: nopt
    nopt = idx_i_solver_opt_e - idx_i_solver_opt_s + 1
    hecmat%Iarray(idx_i_solver_opt_s:idx_i_solver_opt_e) = solver_opt(1:nopt)
  end subroutine hecmw_mat_set_solver_opt
 
  subroutine hecmw_mat_get_solver_opt( hecMAT, solver_opt )
    type(hecmwst_matrix) :: hecmat
    integer(kind=kint) :: solver_opt(:)
    integer(kind=kint) :: nopt
    nopt = idx_i_solver_opt_e - idx_i_solver_opt_s + 1
    solver_opt(1:nopt) = hecmat%Iarray(idx_i_solver_opt_s:idx_i_solver_opt_e)
  end subroutine hecmw_mat_get_solver_opt
 
  subroutine hecmw_mat_set_solver_ropt( hecMAT, solver_ropt )
    type(hecmwst_matrix) :: hecmat
    real(kind=kreal) :: solver_ropt(:)
    integer(kind=kint) :: nopt
    nopt = idx_r_solver_opt_e - idx_r_solver_opt_s + 1
    hecmat%Rarray(idx_r_solver_opt_s:idx_r_solver_opt_e) = solver_ropt(1:nopt)
  end subroutine hecmw_mat_set_solver_ropt
 
  subroutine hecmw_mat_get_solver_ropt( hecMAT, solver_ropt )
    type(hecmwst_matrix) :: hecmat
    real(kind=kreal) :: solver_ropt(:)
    integer(kind=kint) :: nopt
    nopt = idx_r_solver_opt_e - idx_r_solver_opt_s + 1
    solver_ropt(1:nopt) = hecmat%Rarray(idx_r_solver_opt_s:idx_r_solver_opt_e)
  end subroutine hecmw_mat_get_solver_ropt
 
  subroutine hecmw_mat_set_resid( hecMAT, resid )
    type(hecmwst_matrix) :: hecmat
    real(kind=kreal) :: resid
 
    hecmat%Rarray(idx_r_resid) = resid
  end subroutine hecmw_mat_set_resid
 
  function hecmw_mat_get_resid( hecMAT )
    real(kind=kreal) :: hecmw_mat_get_resid
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_resid = hecmat%Rarray(idx_r_resid)
  end function hecmw_mat_get_resid
 
  subroutine hecmw_mat_set_sigma_diag( hecMAT, sigma_diag )
    type(hecmwst_matrix) :: hecmat
    real(kind=kreal) :: sigma_diag
 
    if( sigma_diag < 0.d0 ) then
      hecmat%Rarray(idx_r_sigma_diag) = -1.d0
    elseif( sigma_diag < 1.d0 ) then
      hecmat%Rarray(idx_r_sigma_diag) = 1.d0
    elseif( sigma_diag > 2.d0 ) then
      hecmat%Rarray(idx_r_sigma_diag) = 2.d0
    else
      hecmat%Rarray(idx_r_sigma_diag) = sigma_diag
    endif
  end subroutine hecmw_mat_set_sigma_diag
 
  function hecmw_mat_get_sigma_diag( hecMAT )
    real(kind=kreal) :: hecmw_mat_get_sigma_diag
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_sigma_diag = hecmat%Rarray(idx_r_sigma_diag)
  end function hecmw_mat_get_sigma_diag
 
  subroutine hecmw_mat_set_omega( hecMAT, omega )
    type(hecmwst_matrix) :: hecmat
    real(kind=kreal) :: omega
 
    ! the preconditioner is defined only for 0 < omega < 2; anything else, including
    ! values within machine epsilon of 0 or 2, falls back to plain SSOR
    if (omega < 0.d0 .or. omega > 2.d0 .or. abs(omega) < epsilon(1.d0) &
        .or. abs(omega - 2.d0) < epsilon(1.d0)) then
      hecmat%Rarray(idx_r_omega) = 1.d0
    else
      hecmat%Rarray(idx_r_omega) = omega
    endif
  end subroutine hecmw_mat_set_omega
 
  function hecmw_mat_get_omega( hecMAT )
    real(kind=kreal) :: hecmw_mat_get_omega
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_omega = hecmat%Rarray(idx_r_omega)
  end function hecmw_mat_get_omega
 
  subroutine hecmw_mat_set_thresh( hecMAT, thresh )
    type(hecmwst_matrix) :: hecmat
    real(kind=kreal) :: thresh
 
    hecmat%Rarray(idx_r_thresh) = thresh
  end subroutine hecmw_mat_set_thresh
 
  function hecmw_mat_get_thresh( hecMAT )
    real(kind=kreal) :: hecmw_mat_get_thresh
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_thresh = hecmat%Rarray(idx_r_thresh)
  end function hecmw_mat_get_thresh
 
  subroutine hecmw_mat_set_filter( hecMAT, filter )
    type(hecmwst_matrix) :: hecmat
    real(kind=kreal) :: filter
 
    hecmat%Rarray(idx_r_filter) = filter
  end subroutine hecmw_mat_set_filter
 
  function hecmw_mat_get_filter( hecMAT )
    real(kind=kreal) :: hecmw_mat_get_filter
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_filter = hecmat%Rarray(idx_r_filter)
  end function hecmw_mat_get_filter
 
  subroutine hecmw_mat_set_penalty( hecMAT, penalty )
    type(hecmwst_matrix) :: hecmat
    real(kind=kreal) :: penalty
 
    hecmat%Rarray(idx_r_penalty) = penalty
  end subroutine hecmw_mat_set_penalty
 
  function hecmw_mat_get_penalty( hecMAT )
    real(kind=kreal) :: hecmw_mat_get_penalty
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_penalty = hecmat%Rarray(idx_r_penalty)
  end function hecmw_mat_get_penalty
 
  subroutine hecmw_mat_set_penalty_alpha( hecMAT, alpha )
    type(hecmwst_matrix) :: hecmat
    real(kind=kreal) :: alpha
 
    hecmat%Rarray(idx_r_penalty_alpha) = alpha
  end subroutine hecmw_mat_set_penalty_alpha
 
  function hecmw_mat_get_penalty_alpha( hecMAT )
    real(kind=kreal) :: hecmw_mat_get_penalty_alpha
    type(hecmwst_matrix) :: hecmat
 
    hecmw_mat_get_penalty_alpha = hecmat%Rarray(idx_r_penalty_alpha)
  end function hecmw_mat_get_penalty_alpha
 
  function hecmw_mat_diag_max(hecMAT, hecMESH)
    real(kind=kreal) :: hecmw_mat_diag_max
    type (hecmwst_matrix) :: hecmat
    type (hecmwst_local_mesh) :: hecmesh
    integer(kind=kint) :: ndiag, i
 
    hecmw_mat_diag_max = -1.0e20
    ndiag = hecmat%NDOF**2 * hecmat%NP
    do i = 1, ndiag
      if( hecmat%D(i) > hecmw_mat_diag_max ) hecmw_mat_diag_max = hecmat%D(i)
    enddo
    call hecmw_allreduce_r1(hecmesh, hecmw_mat_diag_max, hecmw_max)
  end function hecmw_mat_diag_max
 
  function hecmw_mat_diag(hecMAT) result(diag)
    type(hecmwst_matrix), intent(in), target :: hecmat
    real(kind=kreal), pointer :: diag(:)
    integer(kind=kint) :: i, k, idx, ndof, np
    
    ndof = hecmat%NDOF
    np = hecmat%NP
    allocate(diag(ndof * np))
    
    do i = 1, np
      do k = 1, ndof
        idx = ndof * ndof * (i - 1) + (k-1) * ndof + k
        diag(ndof * (i - 1) + k) = hecmat%D(idx)
      enddo
    enddo
  end function hecmw_mat_diag
 
  subroutine hecmw_mat_recycle_precond_setting( hecMAT )
    type (hecmwst_matrix) :: hecmat
    integer(kind=kint) :: nrecycle, maxrecycle
    if (hecmat%Iarray(idx_i_flag_symbfact) >= 1) then
      hecmat%Iarray(idx_i_flag_numfact)=1
      call hecmw_mat_reset_nrecycle_precond(hecmat)
    elseif (hecmat%Iarray(idx_i_flag_numfact) > 1) then
      call hecmw_mat_reset_nrecycle_precond(hecmat)
      hecmat%Iarray(idx_i_flag_numfact) = 1
    elseif (hecmat%Iarray(idx_i_flag_numfact) == 1) then
      nrecycle = hecmw_mat_get_nrecycle_precond(hecmat)
      maxrecycle = hecmw_mat_get_maxrecycle_precond(hecmat)
      if ( nrecycle < maxrecycle ) then
        hecmat%Iarray(idx_i_flag_numfact) = 0
        call hecmw_mat_incr_nrecycle_precond(hecmat)
      else
        call hecmw_mat_reset_nrecycle_precond(hecmat)
      endif
    endif
  end subroutine hecmw_mat_recycle_precond_setting
 
  subroutine hecmw_mat_substitute( dest, src )
    type (hecmwst_matrix), intent(inout) :: dest
    type (hecmwst_matrix), intent(inout) :: src
    dest%N = src%N
    dest%NP = src%NP
    dest%NPL = src%NPL
    dest%NPU = src%NPU
    dest%NDOF = src%NDOF
    if (associated(src%D)) dest%D => src%D
    if (associated(src%B)) dest%B => src%B
    if (associated(src%X)) dest%X => src%X
    if (associated(src%AL)) dest%AL => src%AL
    if (associated(src%AU)) dest%AU => src%AU
    if (associated(src%indexL)) dest%indexL => src%indexL
    if (associated(src%indexU)) dest%indexU => src%indexU
    if (associated(src%itemL)) dest%itemL => src%itemL
    if (associated(src%itemU)) dest%itemU => src%itemU
    dest%Iarray(:) = src%Iarray(:)
    dest%Rarray(:) = src%Rarray(:)
  end subroutine hecmw_mat_substitute
 
  subroutine hecmw_mat_integrate( hecMAT )
    type (hecmwst_matrix), intent(inout) :: hecmat
#ifdef _OPENACC
    integer(kind=kint) :: i, j, k, nn, pre, pp, js, je
 
    nn = hecmat%NDOF * hecmat%NDOF
    hecmat%NPA = hecmat%NP + hecmat%NPL + hecmat%NPU
    if (associated(hecmat%A))      deallocate(hecmat%A)
    if (associated(hecmat%indexA)) deallocate(hecmat%indexA)
    if (associated(hecmat%itemA))  deallocate(hecmat%itemA)
    allocate (hecmat%A(nn * hecmat%NPA))
    allocate (hecmat%indexA(0:hecmat%NP))
    allocate (hecmat%itemA(hecmat%NPA))
    hecmat%indexA(0) = 0
 
    pre = 0
    pp = 0
    !$acc parallel loop private(i, j, k, pre, pp, jS, jE)
    do i = 1, hecmat%NP
      hecmat%indexA(i) = i + hecmat%indexL(i) + hecmat%indexU(i)
 
      pre = i - 1 + hecmat%indexU(i - 1)
      js= hecmat%indexL(i - 1) + 1
      je= hecmat%indexL(i  )
      do j = js, je
        pp = pre + j
        hecmat%itemA(pp) = hecmat%itemL(j)
        do k = -nn+1, 0
          hecmat%A(nn * pp + k) = hecmat%AL(nn * j + k)
        enddo
      enddo
 
      pp = i + hecmat%indexU(i - 1) + hecmat%indexL(i)
      hecmat%itemA(pp) = i
      do k = -nn+1, 0
        hecmat%A(nn * pp + k) = hecmat%D(nn * i + k)
      enddo
 
      pre = i + hecmat%indexL(i)
      js= hecmat%indexU(i - 1) + 1
      je= hecmat%indexU(i  )
      do j = js, je
        pp = pre + j
        hecmat%itemA(pp) = hecmat%itemU(j)
        do k = -nn+1, 0
          hecmat%A(nn * pp + k) = hecmat%AU(nn * j + k)
        enddo
      enddo
    enddo
    !$acc end parallel
#endif
  end subroutine hecmw_mat_integrate
 
end module hecmw_matrix_misc