hecmw_mpc_prepost.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_mpc_prepost
  use hecmw_util
  use m_hecmw_comm_f
  use hecmw_matrix_misc
  use hecmw_matrix_ass
  use hecmw_local_matrix
  use hecmw_solver_las
  implicit none
 
  private
  public :: hecmw_mpc_mat_init
  public :: hecmw_mpc_mat_init_explicit
  public :: hecmw_mpc_mat_finalize
  public :: hecmw_mpc_mat_finalize_explicit
  public :: hecmw_mpc_mat_ass
  public :: hecmw_mpc_trans_rhs
  public :: hecmw_mpc_tback_sol
  public :: hecmw_mpc_trans_mass
  public :: hecmw_mpc_tback_eigvec
  public :: hecmw_mpc_mark_slave
 
  integer, parameter :: DEBUG = 0
  logical, parameter :: DEBUG_VECTOR = .false.
 
contains
 
  !C
  !C***
  !C*** hecmw_mpc_mat_init
  !C***
  !C
  subroutine hecmw_mpc_mat_init(hecMESH, hecMAT, hecMESHmpc, hecMATmpc, conMAT, conMATmpc)
    implicit none
    type (hecmwst_local_mesh), intent(inout), target :: hecmesh
    type (hecmwst_matrix), intent(inout), target :: hecmat
    type (hecmwst_local_mesh), pointer :: hecmeshmpc
    type (hecmwst_matrix), pointer :: hecmatmpc
    type (hecmwst_matrix), intent(in), target, optional :: conmat
    type (hecmwst_matrix), pointer, optional :: conmatmpc
    integer(kind=kint) :: totalmpc, mpc_method, solver_type
 
    totalmpc = hecmesh%mpc%n_mpc
    call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
 
    if (totalmpc == 0) then
      hecmeshmpc => hecmesh
      hecmatmpc => hecmat
      if (present(conmat).and.present(conmatmpc)) conmatmpc => conmat
      return
    endif
 
    call hecmw_mpc_scale(hecmesh)
 
    mpc_method = hecmw_mat_get_mpc_method(hecmat)
    if (mpc_method == 2 .and. hecmesh%my_rank == 0) then
      write(*,*) 'WARNING: MPCMETHOD=2 (MPCCG) has been removed; falling back to the default'
    endif
    if (mpc_method /= 1 .and. mpc_method /= 3) then
      solver_type = hecmw_mat_get_solver_type(hecmat)
      if (solver_type > 1) then  ! DIRECT SOLVER
        mpc_method = 1           !   default: penalty
      else                       ! ITERATIVE SOLVER
        mpc_method = 3           !   default: elimination
      endif
      call hecmw_mat_set_mpc_method(hecmat, mpc_method)
    endif
 
    select case (mpc_method)
    case (1)  ! penalty
      hecmeshmpc => hecmesh
      hecmatmpc => hecmat
      if (present(conmat).and.present(conmatmpc)) conmatmpc => conmat
    case (3)  ! elimination
      allocate(hecmeshmpc)
      call hecmw_mpc_mesh_copy(hecmesh, hecmeshmpc)
      allocate(hecmatmpc)
      call hecmw_mat_init(hecmatmpc)
      if (present(conmat).and.present(conmatmpc)) then
        allocate(conmatmpc)
        call hecmw_mat_init(conmatmpc)
      endif
    end select
 
  end subroutine hecmw_mpc_mat_init
 
  !C
  !C***
  !C*** hecmw_mpc_mat_init_explicit
  !C***
  !C
  subroutine hecmw_mpc_mat_init_explicit(hecMESH, hecMAT, hecMATmpc)
    implicit none
    type (hecmwst_local_mesh), intent(inout), target :: hecmesh
    type (hecmwst_matrix), intent(inout), target :: hecmat
    type (hecmwst_matrix), pointer :: hecmatmpc
    integer(kind=kint) :: totalmpc, mpc_method
 
    totalmpc = hecmesh%mpc%n_mpc
    call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
 
    if (totalmpc == 0) then
      hecmatmpc => hecmat
      return
    endif
 
    call hecmw_mpc_scale(hecmesh)
 
    ! Force MPC_METHOD=3
    mpc_method = 3
    call hecmw_mat_set_mpc_method(hecmat, mpc_method)
 
    allocate(hecmatmpc)
    call hecmw_mat_init(hecmatmpc)
 
    hecmatmpc%N = hecmat%N
    hecmatmpc%NP = hecmat%NP
    hecmatmpc%NDOF = hecmat%NDOF
    allocate(hecmatmpc%B(size(hecmat%B)))
    allocate(hecmatmpc%X(size(hecmat%X)))
  end subroutine hecmw_mpc_mat_init_explicit
 
  !C
  !C***
  !C*** hecmw_mpc_mat_finalize
  !C***
  !C
  subroutine hecmw_mpc_mat_finalize(hecMESH, hecMAT, hecMESHmpc, hecMATmpc, conMATmpc)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(in) :: hecmat
    type (hecmwst_local_mesh), pointer :: hecmeshmpc
    type (hecmwst_matrix), pointer :: hecmatmpc
    type (hecmwst_matrix), pointer, optional :: conmatmpc
    integer(kind=kint) :: totalmpc, mpc_method
 
    totalmpc = hecmesh%mpc%n_mpc
    call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
 
    if (totalmpc == 0) then
      nullify(hecmeshmpc)
      nullify(hecmatmpc)
      if (present(conmatmpc)) nullify(conmatmpc)
      return
    endif
 
    mpc_method = hecmw_mat_get_mpc_method(hecmat)
 
    select case (mpc_method)
    case (1)  ! penalty
      nullify(hecmeshmpc)
      nullify(hecmatmpc)
      if (present(conmatmpc)) nullify(conmatmpc)
    case (3) ! elimination
      call hecmw_mpc_mesh_free(hecmeshmpc)
      deallocate(hecmeshmpc)
      nullify(hecmeshmpc)
      call hecmw_mat_finalize(hecmatmpc)
      deallocate(hecmatmpc)
      nullify(hecmatmpc)
      if (present(conmatmpc)) then
        call hecmw_mat_finalize(conmatmpc)
        deallocate(conmatmpc)
        nullify(conmatmpc)
      endif
    end select
 
  end subroutine hecmw_mpc_mat_finalize
 
  !C
  !C***
  !C*** hecmw_mpc_mat_finalize_explicit
  !C***
  !C
  subroutine hecmw_mpc_mat_finalize_explicit(hecMESH, hecMAT, hecMATmpc)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(in) :: hecmat
    type (hecmwst_matrix), pointer :: hecmatmpc
    integer(kind=kint) :: totalmpc, mpc_method
 
    totalmpc = hecmesh%mpc%n_mpc
    call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
 
    if (totalmpc == 0) then
      nullify(hecmatmpc)
      return
    endif
 
    mpc_method = hecmw_mat_get_mpc_method(hecmat)
 
    select case (mpc_method)
    case (1)  ! penalty
      nullify(hecmatmpc)
    case (3) ! elimination
      call hecmw_mat_finalize(hecmatmpc)
      deallocate(hecmatmpc)
      nullify(hecmatmpc)
    end select
 
  end subroutine hecmw_mpc_mat_finalize_explicit
 
  !C
  !C***
  !C*** hecmw_mpc_mat_ass
  !C***
  !C
  subroutine hecmw_mpc_mat_ass(hecMESH, hecMAT, hecMESHmpc, hecMATmpc, conMAT, conMATmpc, hecLagMAT)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(inout) :: hecmat
    type (hecmwst_local_mesh), pointer :: hecmeshmpc
    type (hecmwst_matrix), pointer :: hecmatmpc
    type (hecmwst_matrix), intent(in), optional :: conmat
    type (hecmwst_matrix), pointer, optional :: conmatmpc
    type (hecmwst_matrix_lagrange), intent(inout), optional :: heclagmat
    integer(kind=kint) :: totalmpc, mpc_method
 
    totalmpc = hecmesh%mpc%n_mpc
    call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
 
    if (totalmpc == 0) return
 
    mpc_method = hecmw_mat_get_mpc_method(hecmat)
 
    select case (mpc_method)
    case (1)  ! penalty
      !if (hecMESH%my_rank.eq.0) write(0,*) "MPC Method: Penalty"
      call hecmw_mat_ass_equation ( hecmesh, hecmat )
    case (3)  ! elimination
      !if (hecMESH%my_rank.eq.0) write(0,*) "MPC Method: Elimination"
      call hecmw_trimatmul_ttkt_mpc(hecmeshmpc, hecmat, hecmatmpc)
      if (present(conmat).and.present(conmatmpc).and.present(heclagmat)) then
        call hecmw_trimatmul_ttkt_mpc(hecmeshmpc, conmat, conmatmpc)
        call resize_heclagmat(conmat%NP, conmatmpc%NP, conmat%NDOF, heclagmat)
      endif
    end select
 
  end subroutine hecmw_mpc_mat_ass
 
 
  subroutine resize_heclagmat(NP_orig, NP_new, ndof, hecLagMAT)
    integer(kind=kint), intent(in) :: np_orig, np_new, ndof
    type (hecmwst_matrix_lagrange), intent(inout) :: heclagmat
    integer(kind=kint), pointer :: itemp(:)
 
    if (heclagmat%num_lagrange == 0) return
 
    allocate(itemp(0:np_new))
    itemp(0:np_orig) = heclagmat%indexU_lagrange(0:np_orig)
    itemp(np_orig+1:np_new) = heclagmat%indexU_lagrange(np_orig)
 
    deallocate(heclagmat%indexU_lagrange)
    heclagmat%indexU_lagrange => itemp
 
  end subroutine resize_heclagmat
 
  !C
  !C***
  !C*** hecmw_mpc_trans_rhs
  !C***
  !C
  subroutine hecmw_mpc_trans_rhs(hecMESH, hecMAT, hecMATmpc)
    implicit none
    type (hecmwst_local_mesh), intent(inout) :: hecmesh
    type (hecmwst_matrix), intent(inout) :: hecmat
    type (hecmwst_matrix), pointer :: hecmatmpc
    real(kind=kreal) :: time_dumm
    integer(kind=kint) :: totalmpc, mpc_method
 
    totalmpc = hecmesh%mpc%n_mpc
    call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
 
    if (totalmpc == 0) return
 
    mpc_method = hecmw_mat_get_mpc_method(hecmat)
 
    select case (mpc_method)
    case (1)  ! penalty
      call hecmw_mat_ass_equation_rhs ( hecmesh, hecmatmpc )
    case (3) ! elimination
      time_dumm = 0.d0
      call hecmw_trans_b(hecmesh, hecmat, hecmat%B, hecmatmpc%B, time_dumm)
      hecmatmpc%Iarray=hecmat%Iarray
      hecmatmpc%Rarray=hecmat%Rarray
      ! the elimination T^t A T preserves symmetry; direct solvers read this flag,
      ! which hecmw_mpc_mat_init defaulted to .true. on the newly created hecMATmpc
      hecmatmpc%symmetric=hecmat%symmetric
    end select
 
  end subroutine hecmw_mpc_trans_rhs
 
  !C
  !C***
  !C*** hecmw_mpc_tback_sol
  !C***
  !C
  subroutine hecmw_mpc_tback_sol(hecMESH, hecMAT, hecMATmpc)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(inout) :: hecmat
    type (hecmwst_matrix), pointer :: hecmatmpc
    real(kind=kreal) :: time_dumm
    integer(kind=kint) :: totalmpc, mpc_method, i
    integer(kind=kint) :: npndof, npndof_mpc, num_lagrange
 
    totalmpc = hecmesh%mpc%n_mpc
    call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
 
    if (totalmpc == 0) return
 
    mpc_method = hecmw_mat_get_mpc_method(hecmat)
 
    select case (mpc_method)
    case (1)  ! penalty
      ! do nothing
    case (3)  ! elimination
      npndof = hecmat%NP * hecmat%NDOF
      do i = 1, npndof
        hecmat%X(i) = hecmatmpc%X(i)
      enddo
      time_dumm = 0.d0
      call hecmw_tback_x(hecmesh, hecmat%NDOF, hecmat%X, time_dumm)
      num_lagrange = size(hecmat%X) - npndof
      npndof_mpc = hecmatmpc%NP * hecmatmpc%NDOF
      do i = 1, num_lagrange
        hecmat%X(npndof+i) = hecmatmpc%X(npndof_mpc+i)
      enddo
      hecmat%Iarray=hecmatmpc%Iarray
      hecmat%Rarray=hecmatmpc%Rarray
    end select
  end subroutine hecmw_mpc_tback_sol
 
  !C
  !C***
  !C*** hecmw_mpc_trans_mass
  !C***
  !C
  subroutine hecmw_mpc_trans_mass(hecMESH, hecMAT, hecMATmpc, mass)
    implicit none
    type (hecmwst_local_mesh), intent(inout) :: hecmesh
    type (hecmwst_matrix), intent(inout) :: hecmat
    type (hecmwst_matrix), pointer :: hecmatmpc
    real(kind=kreal), pointer :: mass(:)
 
    real(kind=kreal), allocatable :: w(:), mtmp(:)
    real(kind=kreal) :: time_dumm
    integer(kind=kint) :: totalmpc, mpc_method, i
    integer(kind=kint) :: npndof, npndof_mpc
 
    totalmpc = hecmesh%mpc%n_mpc
    call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
 
    if (totalmpc == 0) return
 
    mpc_method = hecmw_mat_get_mpc_method(hecmat)
 
    select case (mpc_method)
    case (1)  ! penalty
      ! do nothing
    case (3) ! elimination
      npndof = hecmat%NP * hecmat%NDOF
      npndof_mpc = hecmatmpc%NP * hecmatmpc%NDOF
      allocate(w(npndof))
      allocate(mtmp(npndof))
      !C-- lumped mass of the reduced system = row sums of [T'][M][T]
      !C-- {w} = [T]{1}
      mtmp(:) = 1.d0
      time_dumm = 0.d0
      call hecmw_tvec(hecmesh, hecmat%NDOF, mtmp, w, time_dumm)
      !C-- {w} = [M]{w}
      do i = 1, npndof
        w(i) = mass(i) * w(i)
      enddo
      !C-- {Mt} = [T']{w}
      call hecmw_ttvec(hecmesh, hecmat%NDOF, w, mtmp, time_dumm)
      !C-- external nodes added by the reduction get zero mass, as the RHS does
      deallocate(mass)
      allocate(mass(npndof_mpc))
      mass(:) = 0.d0
      do i = 1, npndof
        mass(i) = mtmp(i)
      enddo
      deallocate(w)
      deallocate(mtmp)
    end select
 
  end subroutine hecmw_mpc_trans_mass
 
  !C
  !C***
  !C*** hecmw_mpc_tback_eigvec
  !C***
  !C
  subroutine hecmw_mpc_tback_eigvec(hecMESH, hecMAT, neig, eigvec)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(inout) :: hecmat
    integer(kind=kint), intent(in) :: neig
    real(kind=kreal), intent(inout) :: eigvec(:,:)
 
    real(kind=kreal) :: time_dumm
    integer(kind=kint) :: totalmpc, mpc_method, i
 
    totalmpc = hecmesh%mpc%n_mpc
    call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
 
    if (totalmpc == 0) return
 
    mpc_method = hecmw_mat_get_mpc_method(hecmat)
 
    select case (mpc_method)
    case (1)  ! penalty
      ! do nothing
    case (3)  ! elimination
      time_dumm = 0.d0
      do i = 1, neig
        call hecmw_tback_x(hecmesh, hecmat%NDOF, eigvec(:,i), time_dumm)
        !!! need normalization???
      enddo
    end select
  end subroutine hecmw_mpc_tback_eigvec
 
  !C
  !C***
  !C*** hecmw_mpc_mark_slave
  !C***
  !C
  subroutine hecmw_mpc_mark_slave(hecMESH, hecMAT, mark)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(inout) :: hecmat
    integer(kind=kint), intent(out) :: mark(:)
 
    integer(kind=kint) :: ndof, i, j, k, kk
 
    ndof = hecmat%NDOF
    mark(:) = 0
    outer: do i = 1, hecmesh%mpc%n_mpc
      do j = hecmesh%mpc%mpc_index(i-1)+1, hecmesh%mpc%mpc_index(i)
        if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer
      enddo
      k = hecmesh%mpc%mpc_index(i-1)+1
      kk = ndof * (hecmesh%mpc%mpc_item(k) - 1) + hecmesh%mpc%mpc_dof(k)
      mark(kk) = 1
    enddo outer
  end subroutine hecmw_mpc_mark_slave
 
  !C
  !C***
  !C*** hecmw_mpc_scale
  !C***
  !C
  subroutine hecmw_mpc_scale(hecMESH)
    implicit none
    type (hecmwst_local_mesh), intent(inout) :: hecmesh
    integer(kind=kint) :: i, j, k
    real(kind=kreal) :: wval
 
    !$omp parallel default(none),private(i,j,k,WVAL),shared(hecMESH)
    !$omp do
    do i = 1, hecmesh%mpc%n_mpc
      k = hecmesh%mpc%mpc_index(i-1)+1
      wval = 1.d0 / hecmesh%mpc%mpc_val(k)
      hecmesh%mpc%mpc_val(k) = 1.d0
      do j = hecmesh%mpc%mpc_index(i-1)+2, hecmesh%mpc%mpc_index(i)
        hecmesh%mpc%mpc_val(j) = hecmesh%mpc%mpc_val(j) * wval
      enddo
      hecmesh%mpc%mpc_const(i) = hecmesh%mpc%mpc_const(i) * wval
    enddo
    !$omp end do
    !$omp end parallel
 
  end subroutine hecmw_mpc_scale
 
 
  !C
  !C***
  !C*** hecmw_trans_b
  !C***
  !C
  subroutine hecmw_trans_b(hecMESH, hecMAT, B, BT, COMMtime)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(in) :: hecmat
    real(kind=kreal), intent(in) :: b(:)
    real(kind=kreal), intent(out), target :: bt(:)
    real(kind=kreal), intent(inout) :: commtime
 
    real(kind=kreal), allocatable :: w(:)
    real(kind=kreal), pointer :: xg(:)
    integer(kind=kint) :: ndof, i, j, k, kk
 
    ndof = hecmat%NDOF
 
    call debug_write_vector(b, 'original RHS', 'B', ndof, hecmat%N, hecmat%NP, .true.)
 
    allocate(w(hecmesh%n_node * ndof))
 
    !C===
    !C +---------------------------+
    !C | {bt}= [T']({b} - [A]{xg}) |
    !C +---------------------------+
    !C===
    xg => bt
    do i = 1, hecmat%N * ndof
      xg(i) = 0.d0
    enddo
 
    !C-- Generate {xg} from mpc_const
    !$omp parallel default(none),private(i,k,kk),shared(hecMESH,XG),firstprivate(ndof)
    !$omp do
    outer: do i = 1, hecmesh%mpc%n_mpc
      do j = hecmesh%mpc%mpc_index(i-1)+1, hecmesh%mpc%mpc_index(i)
        if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer
      enddo
      k = hecmesh%mpc%mpc_index(i-1) + 1
      kk = ndof * (hecmesh%mpc%mpc_item(k) - 1) + hecmesh%mpc%mpc_dof(k)
      xg(kk) = hecmesh%mpc%mpc_const(i)
    enddo outer
    !$omp end do
    !$omp end parallel
 
    !C-- {w} = {b} - [A]{xg}
    call hecmw_matresid(hecmesh, hecmat, xg, b, w, commtime)
 
    !C-- {bt} = [T'] {w}
    call hecmw_ttvec(hecmesh, ndof, w, bt, commtime)
 
    deallocate(w)
 
    call debug_write_vector(bt, 'transformed RHS', 'BT', ndof, hecmat%N, hecmat%NP, .true.)
  end subroutine hecmw_trans_b
 
 
  !C
  !C***
  !C*** hecmw_tback_x
  !C***
  !C
  subroutine hecmw_tback_x(hecMESH, ndof, X, COMMtime)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    integer(kind=kint), intent(in) :: ndof
    real(kind=kreal), intent(inout) :: x(:)
    real(kind=kreal), intent(inout) :: commtime
 
    real(kind=kreal), allocatable :: w(:)
    integer(kind=kint) :: i, j, k, kk
 
    call debug_write_vector(x, 'solution for transformed eqn', 'X', ndof, hecmesh%nn_internal, &
        hecmesh%n_node, .true.)
 
    allocate(w(hecmesh%n_node * ndof))
 
    !C-- {tx} = [T]{x}
    call hecmw_tvec(hecmesh, ndof, x, w, commtime)
 
    !C-- {x} = {tx} + {xg}
    !$omp parallel default(none),private(i,k,kk),shared(hecMESH,X,W),firstprivate(ndof)
    !$omp do
    do i= 1, hecmesh%nn_internal * ndof
      x(i)= w(i)
    enddo
    !$omp end do
 
    !$omp do
    outer: do i = 1, hecmesh%mpc%n_mpc
      do j = hecmesh%mpc%mpc_index(i-1)+1, hecmesh%mpc%mpc_index(i)
        if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer
      enddo
      k = hecmesh%mpc%mpc_index(i-1) + 1
      kk = ndof * (hecmesh%mpc%mpc_item(k) - 1) + hecmesh%mpc%mpc_dof(k)
      x(kk) = x(kk) + hecmesh%mpc%mpc_const(i)
    enddo outer
    !$omp end do
    !$omp end parallel
 
    deallocate(w)
 
    call hecmw_update_r(hecmesh, x, hecmesh%n_node, ndof)
    call debug_write_vector(x, 'recovered solution', 'X', ndof, hecmesh%nn_internal, &
        hecmesh%n_node, .true.)
  end subroutine hecmw_tback_x
 
  subroutine hecmw_mpc_mesh_copy(src, dst)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: src
    type (hecmwst_local_mesh), intent(out) :: dst
    dst%zero          = src%zero
    dst%MPI_COMM      = src%MPI_COMM
    dst%PETOT         = src%PETOT
    dst%PEsmpTOT      = src%PEsmpTOT
    dst%my_rank       = src%my_rank
    dst%n_subdomain   = src%n_subdomain
    dst%n_node        = src%n_node
    dst%nn_internal   = src%nn_internal
    dst%n_elem        = src%n_elem
    dst%ne_internal   = src%ne_internal
    dst%n_elem_type   = src%n_elem_type
    dst%n_dof         = src%n_dof
    dst%n_neighbor_pe = src%n_neighbor_pe
    if (src%n_neighbor_pe > 0) then
      allocate(dst%neighbor_pe(dst%n_neighbor_pe))
      dst%neighbor_pe(:) = src%neighbor_pe(:)
      allocate(dst%import_index(0:dst%n_neighbor_pe))
      dst%import_index(:)= src%import_index(:)
      allocate(dst%export_index(0:dst%n_neighbor_pe))
      dst%export_index(:)= src%export_index(:)
      allocate(dst%import_item(dst%import_index(dst%n_neighbor_pe)))
      dst%import_item(:) = src%import_item(:)
      allocate(dst%export_item(dst%export_index(dst%n_neighbor_pe)))
      dst%export_item(:) = src%export_item(:)
    endif
    allocate(dst%global_node_ID(dst%n_node))
    dst%global_node_ID(1:dst%n_node) = src%global_node_ID(1:dst%n_node)
    allocate(dst%node_ID(2*dst%n_node))
    dst%node_ID(1:2*dst%n_node) = src%node_ID(1:2*dst%n_node)
    allocate(dst%elem_type_item(dst%n_elem_type))
    dst%elem_type_item(:) = src%elem_type_item(:)
    !
    dst%mpc%n_mpc =  src%mpc%n_mpc
    dst%mpc%mpc_index => src%mpc%mpc_index
    dst%mpc%mpc_item => src%mpc%mpc_item
    dst%mpc%mpc_dof => src%mpc%mpc_dof
    dst%mpc%mpc_val => src%mpc%mpc_val
    dst%mpc%mpc_const => src%mpc%mpc_const
    !
    dst%node_group%n_grp = src%node_group%n_grp
    dst%node_group%n_bc = src%node_group%n_bc
    dst%node_group%grp_name => src%node_group%grp_name
    dst%node_group%grp_index => src%node_group%grp_index
    dst%node_group%grp_item => src%node_group%grp_item
    dst%node_group%bc_grp_ID => src%node_group%bc_grp_ID
    dst%node_group%bc_grp_type => src%node_group%bc_grp_type
    dst%node_group%bc_grp_index => src%node_group%bc_grp_index
    dst%node_group%bc_grp_dof => src%node_group%bc_grp_dof
    dst%node_group%bc_grp_val => src%node_group%bc_grp_val
    !
    dst%node      => src%node
    dst%elem_type_index => src%elem_type_index
    dst%elem_node_index => src%elem_node_index
    dst%elem_node_item  => src%elem_node_item
  end subroutine hecmw_mpc_mesh_copy
 
  subroutine hecmw_mpc_mesh_free(hecMESH)
    implicit none
    type (hecmwst_local_mesh), intent(inout) :: hecmesh
    if (hecmesh%n_neighbor_pe > 0) then
      deallocate(hecmesh%neighbor_pe)
      deallocate(hecmesh%import_index)
      deallocate(hecmesh%export_index)
      deallocate(hecmesh%import_item)
      deallocate(hecmesh%export_item)
    endif
    deallocate(hecmesh%global_node_ID)
    deallocate(hecmesh%node_ID)
    deallocate(hecmesh%elem_type_item)
  end subroutine hecmw_mpc_mesh_free
 
  subroutine debug_write_vector(Vec, label, name, ndof, N, &
      NP, write_ext, slaves)
    real(kind=kreal),   intent(in) :: vec(:) 
    character(len=*),   intent(in) :: label
    character(len=*),   intent(in) :: name
    integer(kind=kint), intent(in) :: ndof
    integer(kind=kint), intent(in) :: n
    integer(kind=kint), intent(in), optional :: np
    logical,            intent(in), optional :: write_ext
    integer(kind=kint), intent(in), optional :: slaves(:)
    !
    integer(kind=kint) :: iunit
    character(len=128) :: fmt
 
    if (.not. debug_vector) return
 
    write(fmt,'(a,i0,a)') '(',ndof,'f12.3)'
 
    iunit = 1000 + hecmw_comm_get_rank()
    write(iunit,*) trim(label),'------------------------------------------------------------'
    write(iunit,*) 'size of ',trim(name),size(vec)
    write(iunit,*) trim(name),': 1-',n*ndof
    write(iunit,fmt) vec(1:n*ndof)
    if (present(write_ext) .and. present(np)) then
      if (write_ext) then
        write(iunit,*) trim(name),'(external): ',n*ndof+1,'-',np*ndof
        write(iunit,fmt) vec(n*ndof+1:np*ndof)
      endif
    endif
    if (present(slaves)) then
      if (size(slaves) > 0) then
        write(iunit,*) trim(name),'(slave):',slaves(:)
        write(iunit,fmt) vec(slaves(:))
      endif
    endif
  end subroutine debug_write_vector
end module hecmw_mpc_prepost