hecmw_local_matrix.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_local_matrix
  use hecmw_util
  use hecmw_pair_array
 
  private
  public :: hecmwst_local_matrix
  public :: hecmw_localmat_write
  public :: hecmw_localmat_blocking
  public :: hecmw_localmat_free
  public :: hecmw_localmat_mulvec
  public :: hecmw_trimatmul_ttkt
  public :: hecmw_trimatmul_ttkt_serial
  public :: hecmw_trimatmul_ttkt_mpc
  public :: hecmw_localmat_transpose
  public :: hecmw_localmat_assemble
  public :: hecmw_localmat_add
  public :: hecmw_localmat_init_with_hecmat
  public :: hecmw_localmat_add_hecmat
  public :: hecmw_localmat_multmat
  public :: hecmw_localmat_make_hecmat
 
  type hecmwst_local_matrix
    integer :: nr, nc, nnz, ndof
    integer(kind=kint), pointer :: index(:)
    integer(kind=kint), pointer :: item(:)
    real(kind=kreal), pointer :: a(:)
  end type hecmwst_local_matrix
 
  integer(kind=kint), parameter :: cNCOL_ITEM = 3
  integer(kind=kint), parameter :: cLID = 1
  integer(kind=kint), parameter :: cRANK = 2
  integer(kind=kint), parameter :: cGID = 3
 
  integer(kind=kint), parameter :: DEBUG = 0
  integer(kind=kint), parameter :: DEBUG_MATRIX = 0
  integer(kind=kint), parameter :: TIMER = 0
 
contains
 
  subroutine hecmw_localmat_write(Tmat,iunit)
    implicit none
    type (hecmwst_local_matrix), intent(in) :: tmat
    integer(kind=kint), intent(in) :: iunit
    integer(kind=kint) :: nr, nc, nnz, ndof, ndof2, i, js, je, j, jj
    character(len=64) :: fmt
    nr=tmat%nr
    nc=tmat%nc
    nnz=tmat%nnz
    ndof=tmat%ndof
    ndof2=ndof*ndof
    write(iunit,'(a,4i10)') 'nr, nc, nnz, ndof', nr, nc, nnz, ndof
    write(iunit,'(a)') 'i, j, A'
    write(fmt,'(a,i0,a)') '(',ndof,'f12.3)'
    do i=1,nr
      js=tmat%index(i-1)+1
      je=tmat%index(i)
      do j=js,je
        jj=tmat%item(j)
        if (ndof==1) then
          write(iunit,'(2i10,f12.3)') i, jj, tmat%A(j)
        else
          write(iunit,'(2i10)') i, jj
          write(iunit,fmt) tmat%A((j-1)*ndof2+1:j*ndof2)
        endif
      enddo
    enddo
  end subroutine hecmw_localmat_write
 
  subroutine hecmw_localmat_write_size(Tmat,iunit)
    implicit none
    type (hecmwst_local_matrix), intent(in) :: tmat
    integer(kind=kint), intent(in) :: iunit
    integer(kind=kint) :: nr, nc, nnz, ndof
    nr=tmat%nr
    nc=tmat%nc
    nnz=tmat%nnz
    ndof=tmat%ndof
    write(iunit,'(a,4i10)') 'nr, nc, nnz, ndof', nr, nc, nnz, ndof
  end subroutine hecmw_localmat_write_size
 
  subroutine hecmw_localmat_write_ij(Tmat,iunit)
    implicit none
    type (hecmwst_local_matrix), intent(in) :: tmat
    integer(kind=kint), intent(in) :: iunit
    integer(kind=kint) :: nr, nc, nnz, ndof, i, js, je, j, jj
    nr=tmat%nr
    nc=tmat%nc
    nnz=tmat%nnz
    ndof=tmat%ndof
    write(iunit,'(a,4i10)') 'nr, nc, nnz, ndof', nr, nc, nnz, ndof
    write(iunit,'(a)') 'i, j'
    do i=1,nr
      js=tmat%index(i-1)+1
      je=tmat%index(i)
      do j=js,je
        jj=tmat%item(j)
        write(iunit,'(2i10)') i, jj
      enddo
    enddo
  end subroutine hecmw_localmat_write_ij
 
  subroutine hecmw_localmat_blocking(Tmat, ndof, BTmat)
    implicit none
    type (hecmwst_local_matrix), intent(in) :: tmat
    integer, intent(in) :: ndof
    type (hecmwst_local_matrix), intent(out) :: btmat
    integer, allocatable :: iw(:)
    integer :: ndof2, i, icnt, idof, idx, ls, le, l, j, jb, k, lb0, jdof, ks, ke
    ndof2=ndof*ndof
 
    if (mod(tmat%nr, ndof) /= 0 .or. mod(tmat%nc, ndof) /= 0) then
      write(0,*) tmat%nr, tmat%nc, ndof
      stop 'ERROR: blocking_Tmat failed'
    endif
    btmat%nr=tmat%nr/ndof
    btmat%nc=tmat%nc/ndof
    btmat%ndof=ndof
 
    allocate(iw(btmat%nc))
    allocate(btmat%index(0:btmat%nr))
 
    btmat%index(0)=0
    do i=1,btmat%nr
      icnt=0
      do idof=1,ndof
        idx=(i-1)*ndof+idof
        ls=tmat%index(idx-1)+1
        le=tmat%index(idx)
        lcol: do l=ls,le
          j=tmat%item(l)
          jb=(j-1)/ndof+1
          do k=1,icnt
            if (iw(k)==jb) cycle lcol
          enddo
          icnt=icnt+1
          iw(icnt)=jb
        enddo lcol
      enddo
      btmat%index(i)=btmat%index(i-1)+icnt
    enddo
 
    btmat%nnz=btmat%index(btmat%nr)
    allocate(btmat%item(btmat%nnz))
    allocate(btmat%A(btmat%nnz*ndof2))
    btmat%A=0.d0
 
    do i=1,btmat%nr
      icnt=0
      do idof=1,ndof
        idx=(i-1)*ndof+idof
        ls=tmat%index(idx-1)+1
        le=tmat%index(idx)
        lcol2: do l=ls,le
          j=tmat%item(l)
          jb=(j-1)/ndof+1
          do k=1,icnt
            if (iw(k)==jb) cycle lcol2
          enddo
          icnt=icnt+1
          iw(icnt)=jb
        enddo lcol2
      enddo
      ! if (icnt /= BTmat%index(i)-BTmat%index(i-1)) stop 'ERROR: blocking Tmat'
      ! ! call qsort(iw, 1, icnt)
      lb0=btmat%index(i-1)
      do k=1,icnt
        btmat%item(lb0+k)=iw(k)
      enddo
      do idof=1,ndof
        idx=(i-1)*ndof+idof
        ls=tmat%index(idx-1)+1
        le=tmat%index(idx)
        lcol3: do l=ls,le
          j=tmat%item(l)
          jb=(j-1)/ndof+1
          jdof=mod((j-1), ndof)+1
          ks=btmat%index(i-1)+1
          ke=btmat%index(i)
          do k=ks,ke
            if (btmat%item(k)==jb) then
              btmat%A((k-1)*ndof2+(idof-1)*ndof+jdof)=tmat%A(l)
              cycle lcol3
            endif
          enddo
          stop 'ERROR: something wrong in blocking Tmat'
        enddo lcol3
      enddo
    enddo
  end subroutine hecmw_localmat_blocking
 
  subroutine hecmw_localmat_free(Tmat)
    implicit none
    type (hecmwst_local_matrix), intent(inout) :: tmat
    deallocate(tmat%index)
    if (associated(tmat%item)) deallocate(tmat%item)
    if (associated(tmat%A)) deallocate(tmat%A)
    tmat%nr=0
    tmat%nc=0
    tmat%nnz=0
    tmat%ndof=0
  end subroutine hecmw_localmat_free
 
  subroutine hecmw_trimatmul_ttkt(hecMESH, BTtmat, hecMAT, BTmat, &
      iwS, num_lagrange, hecTKT)
    use hecmw_matrix_misc
    implicit none
    type (hecmwst_local_mesh), intent(inout) :: hecmesh
    type (hecmwst_local_matrix), intent(inout) :: bttmat, btmat
    type (hecmwst_matrix), intent(in) :: hecmat
    integer(kind=kint), intent(in) :: iws(:)
    integer(kind=kint), intent(in) :: num_lagrange
    type (hecmwst_matrix), intent(inout) :: hectkt
    if (hecmesh%n_neighbor_pe == 0) then
      call hecmw_trimatmul_ttkt_serial(hecmesh, bttmat, hecmat, btmat, &
           iws, num_lagrange, hectkt)
    else
      call hecmw_trimatmul_ttkt_parallel(hecmesh, bttmat, hecmat, btmat, &
           iws, num_lagrange, hectkt)
    endif
  end subroutine hecmw_trimatmul_ttkt
 
  subroutine hecmw_trimatmul_ttkt_serial(hecMESH, BTtmat, hecMAT, BTmat, &
      iwS, num_lagrange, hecTKT)
    use hecmw_matrix_misc
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_local_matrix), intent(in) :: bttmat, btmat
    type (hecmwst_matrix), intent(in) :: hecmat
    integer(kind=kint), intent(in) :: iws(:)
    integer(kind=kint), intent(in) :: num_lagrange
    type (hecmwst_matrix), intent(inout) :: hectkt
    type (hecmwst_local_matrix) :: bttkt
    real(kind=kreal) :: num
 
    ! perform three matrices multiplication for elimination
    call trimatmul_ttkt(bttmat, hecmat, btmat, bttkt)
    call debug_write_matrix(bttkt, 'BTtKT(MPC)', debug_matrix)
 
    ! place small numbers where the DOF is eliminated
    !num = hecmw_mat_diag_max(hecMAT, hecMESH) * 1.0d-10
    num = 1.d0
    call place_num_on_diag(bttkt, iws, num_lagrange, num)
    call debug_write_matrix(bttkt, 'BTtKT(MPC) (place 1.0 on slave diag)', debug_matrix)
 
    ! make_new HECMW matrix
    call make_new_hecmat(hecmat, bttkt, hectkt)
    call hecmw_localmat_free(bttkt)
  end subroutine hecmw_trimatmul_ttkt_serial
 
  subroutine hecmw_trimatmul_ttkt_parallel(hecMESH, BTtmat, hecMAT, BTmat, &
      iwS, num_lagrange, hecTKT)
    use hecmw_matrix_misc
    implicit none
    type (hecmwST_local_mesh), intent(inout) :: hecMESH
    type (hecmwST_local_matrix), intent(inout) :: BTtmat, BTmat
    type (hecmwST_matrix), intent(in) :: hecMAT
    integer(kind=kint), intent(in) :: iwS(:)
    integer(kind=kint), intent(in) :: num_lagrange
    type (hecmwST_matrix), intent(inout) :: hecTKT
    type (hecmwST_local_matrix) :: BKmat, BTtKmat, BTtKTmat
    real(kind=kreal) :: num
    real(kind=kreal) :: t0, t1
 
    ! perform three matrices multiplication for elimination
    t0 = hecmw_wtime()
    call hecmw_localmat_init_with_hecmat(bkmat, hecmat)
    call debug_write_matrix(bkmat, 'BKmat (hecMAT)', debug_matrix)
    t1 = hecmw_wtime()
    if (timer >= 1) write(0, '(A,f10.4)') "#### hecmw_trimatmul_TtKT_parallel (1) : ",t1-t0
 
    t0 = hecmw_wtime()
    call hecmw_localmat_multmat(bttmat, bkmat, hecmesh, bttkmat)
    if (debug >= 2) write(0,*) '  DEBUG2: multiply Tt and K done'
    call debug_write_matrix(bttkmat, 'BTtKmat', debug_matrix)
    call hecmw_localmat_free(bkmat)
    t1 = hecmw_wtime()
    if (timer >= 1) write(0, '(A,f10.4)') "#### hecmw_trimatmul_TtKT_parallel (2) : ",t1-t0
 
    t0 = hecmw_wtime()
    call hecmw_localmat_multmat(bttkmat, btmat, hecmesh, bttktmat)
    if (debug >= 2) write(0,*) '  DEBUG2: multiply TtK and T done'
    call debug_write_matrix(bttktmat, 'BTtKTmat', debug_matrix)
    call hecmw_localmat_free(bttkmat)
    t1 = hecmw_wtime()
    if (timer >= 1) write(0, '(A,f10.4)') "#### hecmw_trimatmul_TtKT_parallel (3) : ",t1-t0
 
    t0 = hecmw_wtime()
    ! place small numbers where the DOF is eliminated
    !num = hecmw_mat_diag_max(hecMAT, hecMESH) * 1.0d-10
    num = 1.d0
    call place_num_on_diag(bttktmat, iws, num_lagrange, num)
    if (debug >= 2) then
      write(700+hecmw_comm_get_rank(),*) 'num_lagrange =', num_lagrange
      if (debug >= 3) then
        write(700+hecmw_comm_get_rank(),*) 'iwS(1:num_lagrange)'
        write(700+hecmw_comm_get_rank(),*) iws(1:num_lagrange)
      endif
    endif
    call debug_write_matrix(bttktmat, 'BTtKTmat (place 1.0 on slave diag)', debug_matrix)
    t1 = hecmw_wtime()
    if (timer >= 1) write(0, '(A,f10.4)') "#### hecmw_trimatmul_TtKT_parallel (4) : ",t1-t0
 
    t0 = hecmw_wtime()
    ! make_new HECMW matrix
    call make_new_hecmat(hecmat, bttktmat, hectkt)
    call hecmw_localmat_free(bttktmat)
    t1 = hecmw_wtime()
    if (timer >= 1) write(0, '(A,f10.4)') "#### hecmw_trimatmul_TtKT_parallel (5) : ",t1-t0
  end subroutine hecmw_trimatmul_ttkt_parallel
 
  subroutine trimatmul_ttkt(BTtmat, hecMAT, BTmat, BTtKT)
    implicit none
    type (hecmwST_local_matrix), intent(in) :: BTtmat, BTmat
    type (hecmwST_matrix), intent(in) :: hecMAT
    type (hecmwST_local_matrix), intent(out) :: BTtKT
    integer :: nr, nc, ndof, ndof2, i, icnt, js, je, j, jj, ks, ke, k, kk
    integer :: ls, le, l, ll, m, ms, me, mm
    integer, allocatable :: iw(:)
    real(kind=kreal), pointer :: ttp(:), kp(:), tp(:), ttktp(:)
    ! real(kind=kreal) :: tsym_s, tsym_e, tnum_s, tnum_e
 
    nr=bttmat%nr
    nc=btmat%nc
    ndof=bttmat%ndof
    ndof2=ndof*ndof
 
    bttkt%nr=nr
    bttkt%nc=nc
    bttkt%ndof=ndof
    allocate(bttkt%index(0:nr))
 
    ! tsym_s = hecmw_wtime()
 
    !$omp parallel default(none), &
      !$omp&         private(iw,i,icnt,js,je,j,jj,ks,ke,k,kk,ls,le,l,ll,m), &
      !$omp&         shared(nr,nc,BTtmat,hecMAT,BTmat,BTtKT)
    allocate(iw(nc))
    !$omp do
    do i=1,nr
      icnt=0
      js=bttmat%index(i-1)+1
      je=bttmat%index(i)
      do j=js,je
        jj=bttmat%item(j)
        ! lower
        ks=hecmat%indexL(jj-1)+1
        ke=hecmat%indexL(jj)
        do k=ks,ke
          kk=hecmat%itemL(k)
          ls=btmat%index(kk-1)+1
          le=btmat%index(kk)
          ll1: do l=ls,le
            ll=btmat%item(l)
            do m=1,icnt
              if (iw(m)==ll) cycle ll1
            enddo
            icnt=icnt+1
            iw(icnt)=ll
            !if (i==1) write(0,*) 'l', icnt, jj, kk, ll
          enddo ll1
        enddo
        ! diagonal
        ls=btmat%index(jj-1)+1
        le=btmat%index(jj)
        ll2: do l=ls,le
          ll=btmat%item(l)
          do m=1,icnt
            if (iw(m)==ll) cycle ll2
          enddo
          icnt=icnt+1
          iw(icnt)=ll
          !if (i==1) write(0,*) 'd', icnt, jj, kk, ll
        enddo ll2
        ! upper
        ks=hecmat%indexU(jj-1)+1
        ke=hecmat%indexU(jj)
        do k=ks,ke
          kk=hecmat%itemU(k)
          ls=btmat%index(kk-1)+1
          le=btmat%index(kk)
          ll3: do l=ls,le
            ll=btmat%item(l)
            do m=1,icnt
              if (iw(m)==ll) cycle ll3
            enddo
            icnt=icnt+1
            iw(icnt)=ll
            !if (i==1) write(0,*) 'u', icnt, jj, kk, ll
          enddo ll3
        enddo
      enddo
      if (icnt == 0) icnt=1
      !if (i==1) write(0,*) iw(1:icnt)
      bttkt%index(i)=icnt
    enddo
    !$omp end do
    deallocate(iw)
    !$omp end parallel
 
    ! tsym_e = hecmw_wtime()
    ! write(0,*) 'tsym:',tsym_e-tsym_s
 
    bttkt%index(0)=0
    do i=1,nr
      bttkt%index(i)=bttkt%index(i-1)+bttkt%index(i)
    enddo
    !write(0,*) BTtKT%index(1:n)-BTtKT%index(0:n-1)
 
    bttkt%nnz=bttkt%index(nr)
    allocate(bttkt%item(bttkt%nnz))
    allocate(bttkt%A(bttkt%nnz*ndof2))
    bttkt%item=0
    bttkt%A=0.d0
 
    ! tnum_s = hecmw_wtime()
 
    !$omp parallel default(none), &
      !$omp&         private(i,icnt,js,je,j,jj,ks,ke,k,kk,ls,le,l,ll,m, &
      !$omp&                 ms,me,mm,Ttp,Kp,Tp,TtKTp), &
      !$omp&         shared(nr,nc,BTtmat,hecMAT,BTmat,BTtKT,ndof,ndof2)
    !$omp do
    do i=1,nr
      icnt=0
      ms=bttkt%index(i-1)+1
      !me=BTtKT%index(i)
      js=bttmat%index(i-1)+1
      je=bttmat%index(i)
      do j=js,je
        jj=bttmat%item(j)
        ttp=>bttmat%A((j-1)*ndof2+1:j*ndof2)
        ! lower
        ks=hecmat%indexL(jj-1)+1
        ke=hecmat%indexL(jj)
        do k=ks,ke
          kk=hecmat%itemL(k)
          kp=>hecmat%AL((k-1)*ndof2+1:k*ndof2)
          ls=btmat%index(kk-1)+1
          le=btmat%index(kk)
          do l=ls,le
            ll=btmat%item(l)
            tp=>btmat%A((l-1)*ndof2+1:l*ndof2)
            me=ms-1+icnt
            mm=-1
            do m=ms,me
              if (bttkt%item(m)==ll) mm=m
            enddo
            if (mm<0) then
              icnt=icnt+1
              mm=me+1
              bttkt%item(mm)=ll
              !if (i==1) write(0,*) 'l', mm, jj, kk, ll
            endif
            ttktp=>bttkt%A((mm-1)*ndof2+1:mm*ndof2)
            call blk_trimatmul_add(ndof, ttp, kp, tp, ttktp)
          enddo
        enddo
        ! diagonal
        kp=>hecmat%D((jj-1)*ndof2+1:jj*ndof2)
        ls=btmat%index(jj-1)+1
        le=btmat%index(jj)
        do l=ls,le
          ll=btmat%item(l)
          tp=>btmat%A((l-1)*ndof2+1:l*ndof2)
          me=ms-1+icnt
          mm=-1
          do m=ms,me
            if (bttkt%item(m)==ll) mm=m
          enddo
          if (mm<0) then
            icnt=icnt+1
            mm=me+1
            bttkt%item(mm)=ll
            !if (i==1) write(0,*) 'd', mm, jj, kk, ll
          endif
          ttktp=>bttkt%A((mm-1)*ndof2+1:mm*ndof2)
          call blk_trimatmul_add(ndof, ttp, kp, tp, ttktp)
        enddo
        ! upper
        ks=hecmat%indexU(jj-1)+1
        ke=hecmat%indexU(jj)
        do k=ks,ke
          kk=hecmat%itemU(k)
          kp=>hecmat%AU((k-1)*ndof2+1:k*ndof2)
          ls=btmat%index(kk-1)+1
          le=btmat%index(kk)
          do l=ls,le
            ll=btmat%item(l)
            tp=>btmat%A((l-1)*ndof2+1:l*ndof2)
            me=ms-1+icnt
            mm=-1
            do m=ms,me
              if (bttkt%item(m)==ll) mm=m
            enddo
            if (mm<0) then
              icnt=icnt+1
              mm=me+1
              bttkt%item(mm)=ll
              !if (i==1) write(0,*) 'u', mm, jj, kk, ll
            endif
            ttktp=>bttkt%A((mm-1)*ndof2+1:mm*ndof2)
            call blk_trimatmul_add(ndof, ttp, kp, tp, ttktp)
          enddo
        enddo
      enddo
      if (icnt == 0) then
        icnt=1
        bttkt%item(ms)=i
      endif
      ! error check!
      !write(0,*) BTtKT%item(ms:ms-1+icnt)
      !write(0,*) BTtKT%index(i)-BTtKT%index(i-1), icnt
      if (ms-1+icnt /= bttkt%index(i)) stop 'ERROR: trimatmul'
    enddo
    !$omp end do
    !$omp end parallel
 
    ! tnum_e = hecmw_wtime()
    ! write(0,*) 'tnum:',tnum_e-tnum_s
  end subroutine trimatmul_ttkt
 
  subroutine blk_trimatmul_add(ndof, A, B, C, ABC)
    implicit none
    integer, intent(in) :: ndof
    real(kind=kreal), intent(in) :: a(:), b(:), c(:)
    real(kind=kreal), intent(inout) :: abc(:)
    real(kind=kreal), allocatable :: ab(:)
    integer :: ndof2, i, j, k, i0, j0, ij, ik, jk
 
    ndof2=ndof*ndof
    allocate(ab(ndof2))
    ab=0.d0
 
    do i=1,ndof
      i0=(i-1)*ndof
      do j=1,ndof
        ij=i0+j
        j0=(j-1)*ndof
        do k=1,ndof
          ik=i0+k
          jk=j0+k
          ab(ik)=ab(ik)+a(ij)*b(jk)
        enddo
      enddo
    enddo
 
    do i=1,ndof
      i0=(i-1)*ndof
      do j=1,ndof
        ij=i0+j
        j0=(j-1)*ndof
        do k=1,ndof
          ik=i0+k
          jk=j0+k
          abc(ik)=abc(ik)+ab(ij)*c(jk)
        enddo
      enddo
    enddo
 
    deallocate(ab)
  end subroutine blk_trimatmul_add
 
  subroutine place_num_on_diag(BTtKT, iwS, num_lagrange, num)
    implicit none
    type (hecmwST_local_matrix), intent(inout) :: BTtKT
    integer(kind=kint), intent(in) :: iwS(:)
    integer(kind=kint), intent(in) :: num_lagrange
    real(kind=kreal), intent(in) :: num
    integer(kind=kint) :: ndof, ndof2, ilag, i, idof, js, je, j, jj
    integer(kind=kint) :: nmissing, k, ks, ke
    integer(kind=kint), allocatable :: missing(:), cnt(:)
    integer(kind=kint), pointer :: index(:), item(:)
    real(kind=kreal), pointer :: a(:)
 
    ndof=bttkt%ndof
    ndof2=ndof*ndof
 
    ! check if there are places
    allocate(missing(num_lagrange))
    nmissing = 0
    outer1: do ilag=1,num_lagrange
      i=(iws(ilag)-1)/ndof+1
      idof=mod(iws(ilag)-1, ndof)+1
      js=bttkt%index(i-1)+1
      je=bttkt%index(i)
      do j=js,je
        jj=bttkt%item(j)
        if (jj==i) cycle outer1  ! found place
      enddo
      ! not found
      do k=1,nmissing
        if (missing(k) == i) cycle outer1 ! already marked as missing
      enddo
      nmissing = nmissing + 1
      missing(nmissing) = i
    enddo outer1
 
    ! if not, reallocate
    if (nmissing > 0) then
      allocate(cnt(bttkt%nr))
      allocate(index(0:bttkt%nr))
      do i=1,bttkt%nr
        cnt(i) = bttkt%index(i) - bttkt%index(i-1)
      enddo
      do i=1,nmissing
        cnt(missing(i)) = cnt(missing(i)) + 1
      enddo
      call make_index(bttkt%nr, cnt, index)
      allocate(item(bttkt%nnz + nmissing))
      allocate(a(ndof2 * (bttkt%nnz + nmissing)))
      do i=1,bttkt%nr
        ks=index(i-1)+1
        js=bttkt%index(i-1)+1
        je=bttkt%index(i)
        item(ks:ks+(je-js))=bttkt%item(js:je)
        a(ndof2*(ks-1)+1:ndof2*(ks+(je-js)))=bttkt%A(ndof2*(js-1)+1:ndof2*je)
      enddo
      do i=1,nmissing
        ke=index(missing(i))
        item(ke)=missing(i)
        a(ndof2*(ke-1)+1:ndof2*ke)=0.d0
      enddo
      deallocate(bttkt%index)
      deallocate(bttkt%item)
      deallocate(bttkt%A)
      bttkt%index => index
      bttkt%item => item
      bttkt%A => a
      bttkt%nnz = index(bttkt%nr)
      deallocate(cnt)
    endif
    deallocate(missing)
 
    ! place num
    outer: do ilag=1,num_lagrange
      i=(iws(ilag)-1)/ndof+1
      idof=mod(iws(ilag)-1, ndof)+1
      js=bttkt%index(i-1)+1
      je=bttkt%index(i)
      do j=js,je
        jj=bttkt%item(j)
        if (jj==i) then
          !write(0,*) ilag, i, idof
          bttkt%A((j-1)*ndof2+(idof-1)*ndof+idof)=num
          cycle outer
        endif
      enddo
    enddo outer
  end subroutine place_num_on_diag
 
  subroutine replace_hecmat(hecMAT, BTtKT)
    implicit none
    type (hecmwST_matrix), intent(inout) :: hecMAT
    type (hecmwST_local_matrix), intent(in) :: BTtKT
    integer :: nr, nc, ndof, ndof2, i, nl, nu, js, je, j, jj
    integer :: ksl, ksu, k
 
    nr=bttkt%nr
    nc=bttkt%nc
    ndof=hecmat%NDOF
    ndof2=ndof*ndof
 
    ! free old hecMAT
    if (associated(hecmat%AL)) deallocate(hecmat%AL)
    if (associated(hecmat%AU)) deallocate(hecmat%AU)
    if (associated(hecmat%itemL)) deallocate(hecmat%itemL)
    if (associated(hecmat%itemU)) deallocate(hecmat%itemU)
    hecmat%indexL=0
    hecmat%indexU=0
 
    ! count NPL, NPU
    !$omp parallel default(none),private(i,nl,nu,js,je,j,jj), &
      !$omp&         shared(nr,BTtKT,hecMAT)
    !$omp do
    do i=1,nr
      nl=0
      nu=0
      js=bttkt%index(i-1)+1
      je=bttkt%index(i)
      do j=js,je
        jj=bttkt%item(j)
        if (jj < i) then
          nl=nl+1
        elseif (i < jj) then
          nu=nu+1
        else
          ! diagonal
        endif
      enddo
      hecmat%indexL(i)=nl
      hecmat%indexU(i)=nu
    enddo
    !$omp end do
    !$omp end parallel
 
    hecmat%indexL(0)=0
    hecmat%indexU(0)=0
    do i=1,nc
      hecmat%indexL(i)=hecmat%indexL(i-1)+hecmat%indexL(i)
      hecmat%indexU(i)=hecmat%indexU(i-1)+hecmat%indexU(i)
    enddo
    hecmat%NPL=hecmat%indexL(nc)
    hecmat%NPU=hecmat%indexU(nc)
 
    ! allocate new hecMAT
    allocate(hecmat%itemL(hecmat%NPL), hecmat%itemU(hecmat%NPU))
    allocate(hecmat%AL(hecmat%NPL*ndof2), hecmat%AU(hecmat%NPU*ndof2))
    hecmat%itemL=0
    hecmat%itemU=0
    hecmat%D=0.d0
    hecmat%AL=0.d0
    hecmat%AU=0.d0
 
    ! copy from BTtKT to hecMAT
    !$omp parallel default(none),private(i,nl,nu,js,je,ksl,ksu,j,jj,k), &
      !$omp&  shared(nr,BTtKT,hecMAT,ndof2)
    !$omp do
    do i=1,nr
      nl=0
      nu=0
      js=bttkt%index(i-1)+1
      je=bttkt%index(i)
      ksl=hecmat%indexL(i-1)+1
      ksu=hecmat%indexU(i-1)+1
      do j=js,je
        jj=bttkt%item(j)
        if (jj < i) then
          k=ksl+nl
          hecmat%itemL(k)=jj
          hecmat%AL((k-1)*ndof2+1:k*ndof2)=bttkt%A((j-1)*ndof2+1:j*ndof2)
          nl=nl+1
        elseif (i < jj) then
          k=ksu+nu
          hecmat%itemU(k)=jj
          hecmat%AU((k-1)*ndof2+1:k*ndof2)=bttkt%A((j-1)*ndof2+1:j*ndof2)
          nu=nu+1
        else
          hecmat%D((i-1)*ndof2+1:i*ndof2)=bttkt%A((j-1)*ndof2+1:j*ndof2)
        endif
      enddo
      ! if (ksl+nl /= hecMAT%indexL(i)+1) stop 'ERROR: indexL'
      ! if (ksu+nu /= hecMAT%indexU(i)+1) stop 'ERROR: indexU'
    enddo
    !$omp end do
    !$omp end parallel
 
    ! do i=1,hecMAT%NPL
    !   if (hecMAT%itemL(i) <= 0) stop 'ERROR: negative itemL'
    !   if (hecMAT%itemL(i) > nc) stop 'ERROR: too big itemL'
    ! enddo
    ! do i=1,hecMAT%NPU
    !   if (hecMAT%itemU(i) <= 0) stop 'ERROR: negative itemU'
    !   if (hecMAT%itemU(i) > nc) stop 'ERROR: too big itemU'
    ! enddo
  end subroutine replace_hecmat
 
  subroutine make_new_hecmat(hecMAT, BTtKT, hecTKT)
    implicit none
    type(hecmwst_matrix), intent(in) :: hecMAT
    type(hecmwst_local_matrix), intent(in) :: BTtKT
    type(hecmwst_matrix), intent(inout) :: hecTKT
    integer(kind=kint) :: nr, nc, ndof, ndof2
 
    nr=bttkt%nr
    nc=bttkt%nc
    ndof=bttkt%ndof
    ndof2=ndof*ndof
 
    !write(0,*) 'DEBUG: nr, nc =',nr,nc
 
    ! if (nr /= nc) then
    !   stop 'ERROR: nr /= nc'
    ! endif
    hectkt%N =hecmat%N
    hectkt%NP=nc
    hectkt%NDOF=ndof
 
    if (associated(hectkt%D)) deallocate(hectkt%D)
    allocate(hectkt%D(nc*ndof2))
 
    if (associated(hectkt%indexL)) deallocate(hectkt%indexL)
    if (associated(hectkt%indexU)) deallocate(hectkt%indexU)
    allocate(hectkt%indexL(0:nc))
    allocate(hectkt%indexU(0:nc))
 
    hectkt%Iarray=hecmat%Iarray
    hectkt%Rarray=hecmat%Rarray
 
    call replace_hecmat(hectkt, bttkt)
  end subroutine make_new_hecmat
 
  subroutine hecmw_localmat_mulvec(BTmat, V, TV)
    implicit none
    type (hecmwst_local_matrix), intent(in) :: btmat
    real(kind=kreal), intent(in), target :: v(:)
    real(kind=kreal), intent(out), target :: tv(:)
    real(kind=kreal), pointer :: tvp(:), tp(:), vp(:)
    integer :: nr, ndof, ndof2, i, js, je, j, jj, k, kl0, l
    !!$    real(kind=kreal) :: vnorm
 
    nr=btmat%nr
    ndof=btmat%ndof
    ndof2=ndof*ndof
 
    tv=0.d0
 
    !!$    vnorm=0.d0
    !!$    do i=1,nr*ndof
    !!$      vnorm=vnorm+V(i)**2
    !!$    enddo
    !!$    write(0,*) 'vnorm:', sqrt(vnorm)
 
    !$omp parallel default(none),private(i,TVp,js,je,j,jj,Tp,Vp,k,kl0,l), &
      !$omp&  shared(nr,TV,ndof,BTmat,ndof2,V)
    !$omp do
    do i=1,nr
      tvp=>tv((i-1)*ndof+1:i*ndof)
      js=btmat%index(i-1)+1
      je=btmat%index(i)
      do j=js,je
        jj=btmat%item(j)
        tp=>btmat%A((j-1)*ndof2+1:j*ndof2)
        vp=>v((jj-1)*ndof+1:jj*ndof)
        do k=1,ndof
          kl0=(k-1)*ndof
          do l=1,ndof
            tvp(k)=tvp(k)+tp(kl0+l)*vp(l)
          enddo
        enddo
      enddo
    enddo
    !$omp end do
    !$omp end parallel
  end subroutine hecmw_localmat_mulvec
 
  subroutine hecmw_trimatmul_ttkt_mpc(hecMESH, hecMAT, hecTKT)
    implicit none
    type (hecmwst_local_mesh), intent(inout) :: hecmesh
    type (hecmwst_matrix), intent(in) :: hecmat
    type (hecmwst_matrix), intent(inout) :: hectkt
    type (hecmwst_local_matrix) :: btmat, bttmat
    integer(kind=kint), allocatable :: iws(:)
    integer(kind=kint) :: ndof, n_mpc, i_mpc
    integer(kind=kint) :: i, j, k, kk, ilag
    integer(kind=kint) :: num_lagrange
    real(kind=kreal) :: t0, t1
    t0 = hecmw_wtime()
    ndof=hecmat%NDOF
    n_mpc=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
      n_mpc=n_mpc+1
    enddo outer
    allocate(iws(n_mpc))
    i_mpc=0
    outer2: 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 outer2
      enddo
      i_mpc=i_mpc+1
      k=hecmesh%mpc%mpc_index(i-1)+1
      kk=ndof*(hecmesh%mpc%mpc_item(k)-1)+hecmesh%mpc%mpc_dof(k)
      iws(i_mpc)=kk
    enddo outer2
    if (debug >= 2) then
      write(700+hecmw_comm_get_rank(),*) 'DEBUG: n_mpc, slaves',n_mpc,iws(1:n_mpc)
    endif
    t1 = hecmw_wtime()
    if (timer >= 1) write(0, '(A,f10.4)') "### hecmw_trimatmul_TtKT_mpc (1) : ",t1-t0
    t0 = hecmw_wtime()
    call make_btmat_mpc(hecmesh, ndof, btmat)
    call debug_write_matrix(btmat, 'BTmat(MPC)', debug_matrix)
    t1 = hecmw_wtime()
    if (timer >= 1) write(0, '(A,f10.4)') "### hecmw_trimatmul_TtKT_mpc (2) : ",t1-t0
    t0 = hecmw_wtime()
    call hecmw_localmat_transpose(btmat, bttmat)
    ! if (hecmw_localmat_equal(BTtmat, BTtmat2) == 0) then
    !   write(0,*) 'ERROR: BTtmat2 is incorrect!!!'
    ! else
    !   write(0,*) 'DEBUG: BTtmat2 is correct'
    ! endif
    call debug_write_matrix(bttmat, 'BTtmat(MPC)', debug_matrix)
 
    if (debug >= 3) then
      write(700+hecmw_comm_get_rank(),*) 'hecMESH%node_ID before trimatmul_TtKT'
      do i=hecmesh%nn_internal+1, hecmesh%n_node
        write(700+hecmw_comm_get_rank(),*) i,hecmesh%node_ID(2*i-1),hecmesh%node_ID(2*i),hecmesh%global_node_ID(i)
      enddo
    endif
    t1 = hecmw_wtime()
    if (timer >= 1) write(0, '(A,f10.4)') "### hecmw_trimatmul_TtKT_mpc (3) : ",t1-t0
    t0 = hecmw_wtime()
    call hecmw_trimatmul_ttkt(hecmesh, bttmat, hecmat, btmat, iws, n_mpc, hectkt)
    t1 = hecmw_wtime()
    if (timer >= 1) write(0, '(A,f10.4)') "### hecmw_trimatmul_TtKT_mpc (4) : ",t1-t0
    t0 = hecmw_wtime()
    if (debug >= 3) then
      write(700+hecmw_comm_get_rank(),*) 'hecMESH%node_ID after trimatmul_TtKT'
      do i=hecmesh%nn_internal+1, hecmesh%n_node
        write(700+hecmw_comm_get_rank(),*) i,hecmesh%node_ID(2*i-1),hecmesh%node_ID(2*i),hecmesh%global_node_ID(i)
      enddo
    endif
 
    if (associated(hectkt%B)) deallocate(hectkt%B)
    if (associated(hectkt%X)) deallocate(hectkt%X)
    num_lagrange = size(hecmat%B) - hecmat%NP*ndof
    allocate(hectkt%B(ndof*hectkt%NP + num_lagrange))
    allocate(hectkt%X(ndof*hectkt%NP + num_lagrange))
    hectkt%B(:) = 0.d0
    hectkt%X(:) = 0.d0
    do i=1, ndof*hecmat%NP
      hectkt%B(i) = hecmat%B(i)
      hectkt%X(i) = hecmat%X(i)
    enddo
    do i=1, num_lagrange
      hectkt%B(ndof*hectkt%NP+i) = hecmat%B(ndof*hecmat%NP+i)
      hectkt%X(ndof*hectkt%NP+i) = hecmat%X(ndof*hecmat%NP+i)
    enddo
    do ilag=1,n_mpc
      hectkt%X(iws(ilag)) = 0.d0
    enddo
 
    call hecmw_localmat_free(btmat)
    call hecmw_localmat_free(bttmat)
    ! call hecmw_localmat_free(BTtmat2)
    deallocate(iws)
    t1 = hecmw_wtime()
    if (timer >= 1) write(0, '(A,f10.4)') "### hecmw_trimatmul_TtKT_mpc (5) : ",t1-t0
  end subroutine hecmw_trimatmul_ttkt_mpc
 
  subroutine make_btmat_mpc(hecMESH, ndof, BTmat)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    integer(kind=kint), intent(in) :: ndof
    type (hecmwst_local_matrix), intent(out) :: btmat
    type (hecmwst_local_matrix) :: tmat
    integer(kind=kint) :: n_mpc
    integer(kind=kint) :: i,j,k,js,jj,kk
    n_mpc=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
      n_mpc=n_mpc+1
    enddo outer
    tmat%nr=hecmesh%n_node*ndof
    tmat%nc=tmat%nr
    tmat%ndof=1
    allocate(tmat%index(0:tmat%nr))
    ! count nonzero in each row
    tmat%index(1:tmat%nr)=1
    outer2: 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 outer2
      enddo
      k=hecmesh%mpc%mpc_index(i-1)+1
      kk=ndof*(hecmesh%mpc%mpc_item(k)-1)+hecmesh%mpc%mpc_dof(k)
      tmat%index(kk)=hecmesh%mpc%mpc_index(i)-hecmesh%mpc%mpc_index(i-1)-1
    enddo outer2
    ! index
    tmat%index(0)=0
    do i=1,tmat%nr
      tmat%index(i)=tmat%index(i-1)+tmat%index(i)
    enddo
    tmat%nnz=tmat%index(tmat%nr)
    allocate(tmat%item(tmat%nnz), tmat%A(tmat%nnz))
    ! diag
    do i=1,tmat%nr
      js=tmat%index(i-1)+1
      tmat%item(js)=i
      tmat%A(js)=1.d0
    enddo
    ! others
    outer3: 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 outer3
      enddo
      k=hecmesh%mpc%mpc_index(i-1)+1
      kk=ndof*(hecmesh%mpc%mpc_item(k)-1)+hecmesh%mpc%mpc_dof(k)
      js=tmat%index(kk-1)+1
      do j= hecmesh%mpc%mpc_index(i-1) + 2, hecmesh%mpc%mpc_index(i)
        jj = ndof * (hecmesh%mpc%mpc_item(j) - 1) + hecmesh%mpc%mpc_dof(j)
        tmat%item(js)=jj
        tmat%A(js)=-hecmesh%mpc%mpc_val(j)
        js=js+1
      enddo
    enddo outer3
    !call debug_write_matrix(Tmat, 'Tmat(MPC)', DEBUG_MATRIX)
    call hecmw_localmat_blocking(tmat, ndof, btmat)
    call hecmw_localmat_free(tmat)
  end subroutine make_btmat_mpc
 
 
  subroutine hecmw_localmat_transpose(Tmat, Ttmat)
    implicit none
    type (hecmwst_local_matrix), intent(in) :: tmat
    type (hecmwst_local_matrix), intent(out) :: ttmat
    integer(kind=kint), allocatable :: iw(:)
    integer(kind=kint) :: i, j, jj, ndof, ndof2, k, idof, jdof
    allocate(iw(tmat%nc))
    iw = 0
    do i = 1, tmat%nr
      do j = tmat%index(i-1)+1, tmat%index(i)
        jj = tmat%item(j)
        iw(jj) = iw(jj) + 1
      enddo
    enddo
    ttmat%nr = tmat%nc
    ttmat%nc = tmat%nr
    ttmat%nnz = tmat%nnz
    ttmat%ndof = tmat%ndof
    ndof = tmat%ndof
    ndof2 = ndof * ndof
    allocate(ttmat%index(0:ttmat%nr))
    allocate(ttmat%item(ttmat%nnz))
    allocate(ttmat%A(ttmat%nnz*ndof2))
    ttmat%index(0) = 0
    do i = 1, ttmat%nr
      ttmat%index(i) = ttmat%index(i-1) + iw(i)
      iw(i) = ttmat%index(i-1) + 1
    enddo
    do i = 1, tmat%nr
      do j = tmat%index(i-1)+1, tmat%index(i)
        jj = tmat%item(j)
        k = iw(jj)
        ttmat%item( k ) = i
        do idof = 1, ndof
          do jdof = 1, ndof
            ttmat%A((k-1)*ndof2+(idof-1)*ndof+jdof) = &
              tmat%A((j-1)*ndof2+(jdof-1)*ndof+idof)
          enddo
        enddo
        iw(jj) = k + 1
      enddo
    enddo
  end subroutine hecmw_localmat_transpose
 
  function hecmw_localmat_equal(Tmat1, Tmat2)
    implicit none
    type (hecmwst_local_matrix), intent(in) :: tmat1, tmat2
    integer(kind=kint) :: hecmw_localmat_equal
    integer(kind=kint) :: i, j, k0, k, ndof, ndof2
    hecmw_localmat_equal = 0
    if (tmat1%nr /= tmat2%nr) return
    if (tmat1%nc /= tmat2%nc) return
    if (tmat1%nnz /= tmat2%nnz) return
    if (tmat1%ndof /= tmat2%ndof) return
    ndof = tmat1%ndof
    ndof2 = ndof * ndof
    do i = 1, tmat1%nr
      if (tmat1%index(i) /= tmat2%index(i)) return
      do j = tmat1%index(i-1)+1, tmat1%index(i)
        if (tmat1%item(j) /= tmat2%item(j)) return
        k0 = (j-1)*ndof2
        do k = 1, ndof2
          if (tmat1%A(k0+k) /= tmat2%A(k0+k)) return
        enddo
      enddo
    enddo
    hecmw_localmat_equal = 1
  end function hecmw_localmat_equal
 
!!!
!!! Subroutines for parallel contact analysis with iterative linear solver
!!!
 
  subroutine hecmw_localmat_assemble(BTmat, hecMESH, hecMESHnew)
    implicit none
    type (hecmwst_local_matrix), intent(inout) :: btmat
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_local_mesh), intent(inout) :: hecmeshnew
    integer(kind=kint) :: nn_int, np, ndof, ndof2, nr_ext, nnz_ext
    integer(kind=kint), allocatable :: exp_rows_index(:), exp_cols_index(:)
    integer(kind=kint), allocatable :: exp_rows_item(:,:), exp_cols_item(:,:)
    type (hecmwst_local_matrix), allocatable :: bt_ext(:)
    type (hecmwst_local_matrix) :: bt_int
    type (hecmwst_local_matrix) :: btnew
    ! some checks
    if (debug >= 1) write(0,*) 'DEBUG: nr,nc,nnz,ndof',btmat%nr,btmat%nc,btmat%nnz,btmat%ndof
    if (btmat%nr /= hecmesh%n_node) stop 'ERROR: invalid size in hecmw_localmat_assemble'
    !
    nn_int = hecmesh%nn_internal
    np = hecmesh%n_node
    ndof = btmat%ndof
    ndof2 = ndof*ndof
    !
    nr_ext = np - nn_int
    nnz_ext = btmat%index(np) - btmat%index(nn_int)
    !
    call prepare_bt_ext(btmat, hecmesh, exp_rows_index, exp_rows_item, bt_ext)
    if (debug >= 1) write(0,*) 'DEBUG: prepare_BT_ext done'
    !
    call prepare_column_info(hecmesh, bt_ext, exp_cols_index, exp_cols_item)
    if (debug >= 1) write(0,*) 'DEBUG: prepare_column info done'
    !
    call send_bt_ext_and_recv_bt_int(hecmesh, exp_rows_index, exp_rows_item, bt_ext, &
         exp_cols_index, exp_cols_item, bt_int, hecmeshnew)
    if (debug >= 1) write(0,*) 'DEBUG: send BT_ext and recv BT_int done'
    !
    !write(0,*) 'BTmat%ndof,BT_int%ndof',BTmat%ndof,BT_int%ndof
    call hecmw_localmat_add(btmat, bt_int, btnew)
    if (debug >= 1) write(0,*) 'DEBUG: localmat_add done'
    !
    call hecmw_localmat_free(btmat)
    call hecmw_localmat_free(bt_int)
    !
    btmat%nr = btnew%nr
    btmat%nc = btnew%nc
    btmat%nnz = btnew%nnz
    btmat%ndof = btnew%ndof
    btmat%index => btnew%index
    btmat%item => btnew%item
    btmat%A => btnew%A
    !
    ! hecMESH%n_node = hecMESHnew%n_node
    ! hecMESH%n_neighbor_pe = hecMESHnew%n_neighbor_pe
    ! deallocate(hecMESH%neighbor_pe)
    ! deallocate(hecMESH%import_index)
    ! deallocate(hecMESH%export_index)
    ! deallocate(hecMESH%import_item)
    ! deallocate(hecMESH%export_item)
    ! deallocate(hecMESH%node_ID)
    ! deallocate(hecMESH%global_node_ID)
    ! hecMESH%neighbor_pe => hecMESHnew%neighbor_pe
    ! hecMESH%import_index => hecMESHnew%import_index
    ! hecMESH%export_index => hecMESHnew%export_index
    ! hecMESH%import_item => hecMESHnew%import_item
    ! hecMESH%export_item => hecMESHnew%export_item
    ! hecMESH%node_ID => hecMESHnew%node_ID
    ! hecMESH%global_node_ID => hecMESHnew%global_node_ID
    !
    if (debug >= 1) write(0,*) 'DEBUG: update BTmat and hecMESH done'
  end subroutine hecmw_localmat_assemble
 
  subroutine prepare_bt_ext(BTmat, hecMESH, exp_rows_index, exp_rows_item, BT_ext)
    implicit none
    type (hecmwst_local_matrix), intent(in) :: btmat
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    integer(kind=kint), allocatable, intent(out) :: exp_rows_index(:)
    integer(kind=kint), allocatable, intent(out) :: exp_rows_item(:,:)
    type (hecmwst_local_matrix), allocatable, intent(out) :: bt_ext(:)
    integer(kind=kint), allocatable :: incl_nz(:), exp_cols_per_row(:), exp_rows_per_rank(:)
    integer(kind=kint) :: nn_int
    logical, parameter :: flg_check_nonzero_numerically = .true.
    nn_int = hecmesh%nn_internal
    !
    if (flg_check_nonzero_numerically) then
      ! efficient for assembling conMAT which has same non-zero profile as hecMAT
      ! but only dofs related to cntact are actually non-zero
      call check_external_nz_blocks(btmat, nn_int, incl_nz)
    else
      ! probably good enough to assemble T or T^t
      call incl_all_external_nz_blocks(btmat, nn_int, incl_nz)
    endif
    !
    call count_ext_rows_with_nz(btmat, nn_int, incl_nz, exp_cols_per_row)
    !
    call count_exp_rows_per_rank(hecmesh, exp_cols_per_row, exp_rows_per_rank)
    !
    allocate(exp_rows_index(0:hecmesh%n_neighbor_pe))
    call make_index(hecmesh%n_neighbor_pe, exp_rows_per_rank, exp_rows_index)
    !write(0,*) 'exp_rows_index',exp_rows_index(:)
    !
    deallocate(exp_rows_per_rank)
    !
    call make_exp_rows_item(hecmesh, exp_cols_per_row, exp_rows_index, exp_rows_item)
    !
    deallocate(exp_cols_per_row)
    !
    allocate(bt_ext(hecmesh%n_neighbor_pe))
    call extract_bt_ext(hecmesh, btmat, incl_nz, exp_rows_index, exp_rows_item, bt_ext)
    !
    deallocate(incl_nz)
  end subroutine prepare_bt_ext
 
  subroutine check_external_nz_blocks(BTmat, nn_internal, incl_nz)
    implicit none
    type (hecmwst_local_matrix), intent(in) :: btmat
    integer(kind=kint), intent(in) :: nn_internal
    integer(kind=kint), allocatable, intent(out) :: incl_nz(:)
    integer(kind=kint) :: ndof2, i0, nnz_ext, i, k, nnz_blk
    if (nn_internal > btmat%nr) stop 'ERROR: invalid nn_internal'
    ndof2 = btmat%ndof ** 2
    i0 = btmat%index(nn_internal)
    nnz_ext = btmat%index(btmat%nr) - i0
    allocate(incl_nz(nnz_ext))
    nnz_blk = 0
    do i = 1, nnz_ext
      incl_nz(i) = 0
      do k = 1, ndof2
        if (btmat%A(ndof2*(i0+i-1)+k) /= 0.0d0) then
          incl_nz(i) = 1
          nnz_blk = nnz_blk + 1
          exit
        endif
      enddo
    enddo
    if (debug >= 1) write(0,*) 'DEBUG: nnz_blk',nnz_blk
  end subroutine check_external_nz_blocks
 
  subroutine incl_all_external_nz_blocks(BTmat, nn_internal, incl_nz)
    implicit none
    type (hecmwst_local_matrix), intent(in) :: btmat
    integer(kind=kint), intent(in) :: nn_internal
    integer(kind=kint), allocatable, intent(out) :: incl_nz(:)
    integer(kind=kint) :: i0, nnz_ext
    if (nn_internal > btmat%nr) stop 'ERROR: invalid nn_internal'
    i0 = btmat%index(nn_internal)
    nnz_ext = btmat%index(btmat%nr) - i0
    allocate(incl_nz(nnz_ext))
    incl_nz(1:nnz_ext) = 1
  end subroutine incl_all_external_nz_blocks
 
  subroutine count_ext_rows_with_nz(BTmat, nn_internal, incl_nz, exp_cols_per_row)
    implicit none
    type (hecmwst_local_matrix), intent(in) :: btmat
    integer(kind=kint), intent(in) :: nn_internal
    integer(kind=kint), intent(in) :: incl_nz(:)
    integer(kind=kint), allocatable, intent(out) :: exp_cols_per_row(:)
    integer(kind=kint) :: nr_ext, nnz_int, i, irow, js, je, j, jcol
    nr_ext = btmat%nr - nn_internal
    nnz_int = btmat%index(nn_internal)
    allocate(exp_cols_per_row(nr_ext))
    exp_cols_per_row(:) = 0
    do i = 1, nr_ext
      irow = nn_internal+i
      js = btmat%index(irow-1)+1
      je = btmat%index(irow)
      do j = js, je
        jcol = btmat%item(j)
        if (incl_nz(j-nnz_int) == 1) exp_cols_per_row(i) = exp_cols_per_row(i) + 1
      enddo
    enddo
    !write(0,*) 'exp_cols_per_row',exp_cols_per_row(:)
  end subroutine count_ext_rows_with_nz
 
  subroutine count_exp_rows_per_rank(hecMESH, exp_cols_per_row, exp_rows_per_rank)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    integer(kind=kint), intent(in) :: exp_cols_per_row(:)
    integer(kind=kint), allocatable, intent(out) :: exp_rows_per_rank(:)
    integer(kind=kint) :: nn_int, np, nr_ext, i, irow, exp_rank, idom
    allocate(exp_rows_per_rank(hecmesh%n_neighbor_pe))
    exp_rows_per_rank(1:hecmesh%n_neighbor_pe) = 0
    nn_int = hecmesh%nn_internal
    np = hecmesh%n_node
    nr_ext = np - nn_int
    do i = 1, nr_ext
      if (exp_cols_per_row(i) > 0) then
        irow = nn_int + i
        exp_rank = hecmesh%node_ID(2*irow)
        call rank_to_idom(hecmesh, exp_rank, idom)
        exp_rows_per_rank(idom) = exp_rows_per_rank(idom) + 1
      endif
    enddo
    !write(0,*) 'exp_rows_per_rank',exp_rows_per_rank(:)
  end subroutine count_exp_rows_per_rank
 
  subroutine rank_to_idom(hecMESH, rank, idom)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    integer(kind=kint), intent(in) :: rank
    integer(kind=kint), intent(out) :: idom
    integer(kind=kint) :: i
    do i = 1, hecmesh%n_neighbor_pe
      if (hecmesh%neighbor_pe(i) == rank) then
        idom = i
        return
      endif
    enddo
    stop 'ERROR: exp_rank not found in neighbor_pe'
  end subroutine rank_to_idom
 
  subroutine make_index(len, cnt, index)
    implicit none
    integer(kind=kint), intent(in) :: len
    integer(kind=kint), intent(in) :: cnt(len)
    integer(kind=kint), intent(out) :: index(0:)
    integer(kind=kint) :: i
    ! write(0,*) 'make_index: len',len
    index(0) = 0
    do i = 1,len
      index(i) = index(i-1) + cnt(i)
    enddo
  end subroutine make_index
 
  subroutine make_exp_rows_item(hecMESH, exp_cols_per_row, exp_rows_index, exp_rows_item)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    integer(kind=kint), intent(in) :: exp_cols_per_row(:)
    integer(kind=kint), allocatable, intent(in) :: exp_rows_index(:)
    integer(kind=kint), allocatable, intent(out) :: exp_rows_item(:,:)
    integer(kind=kint), allocatable :: cnt(:)
    integer(kind=kint) :: nn_int, np, nr_ext, i, irow, exp_rank, idom, idx
    allocate(exp_rows_item(2,exp_rows_index(hecmesh%n_neighbor_pe)))
    allocate(cnt(hecmesh%n_neighbor_pe))
    cnt(:) = 0
    nn_int = hecmesh%nn_internal
    np = hecmesh%n_node
    nr_ext = np - nn_int
    do i = 1, nr_ext
      if (exp_cols_per_row(i) > 0) then
        irow = nn_int + i
        exp_rank = hecmesh%node_ID(2*irow)
        call rank_to_idom(hecmesh, exp_rank, idom)
        cnt(idom) = cnt(idom) + 1
        idx = exp_rows_index(idom-1) + cnt(idom)
        exp_rows_item(1,idx) = irow
        exp_rows_item(2,idx) = exp_cols_per_row(i)
      endif
    enddo
    !write(0,*) 'cnt',cnt(:)
    do idom = 1, hecmesh%n_neighbor_pe
      if (cnt(idom) /= exp_rows_index(idom)-exp_rows_index(idom-1)) stop 'ERROR: make exp_rows_item'
    enddo
    !write(0,*) 'exp_rows_item(1,:)',exp_rows_item(1,:)
    !write(0,*) 'exp_rows_item(2,:)',exp_rows_item(2,:)
  end subroutine make_exp_rows_item
 
  subroutine extract_bt_ext(hecMESH, BTmat, incl_nz, exp_rows_index, exp_rows_item, BT_ext)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_local_matrix), intent(in) :: btmat
    integer(kind=kint), intent(in) :: incl_nz(:)
    integer(kind=kint), allocatable, intent(in) :: exp_rows_index(:)
    integer(kind=kint), intent(in) :: exp_rows_item(:,:)
    type (hecmwst_local_matrix), allocatable, intent(out) :: bt_ext(:)
    integer(kind=kint) :: ndof, ndof2, nn_int, nnz_int, idom, j, idx, ncol, cnt, jrow, ks, ke, k, kcol
    allocate(bt_ext(hecmesh%n_neighbor_pe))
    ndof = btmat%ndof
    ndof2 = ndof * ndof
    nn_int = hecmesh%nn_internal
    nnz_int = btmat%index(nn_int)
    do idom = 1, hecmesh%n_neighbor_pe
      bt_ext(idom)%nr = exp_rows_index(idom) - exp_rows_index(idom-1)
      bt_ext(idom)%nc = btmat%nc
      bt_ext(idom)%nnz = 0
      bt_ext(idom)%ndof = ndof
      allocate(bt_ext(idom)%index(0:bt_ext(idom)%nr))
      bt_ext(idom)%index(0) = 0
      do j = 1, bt_ext(idom)%nr
        idx = exp_rows_index(idom-1) + j
        ncol = exp_rows_item(2,idx)
        bt_ext(idom)%index(j) = bt_ext(idom)%index(j-1) + ncol
      enddo
      bt_ext(idom)%nnz = bt_ext(idom)%index(bt_ext(idom)%nr)
      if (debug >= 1) write(0,*) 'DEBUG: idom,nr,nc,nnz,ndof', &
           idom,bt_ext(idom)%nr,bt_ext(idom)%nc,bt_ext(idom)%nnz,bt_ext(idom)%ndof
      allocate(bt_ext(idom)%item(bt_ext(idom)%nnz))
      allocate(bt_ext(idom)%A(bt_ext(idom)%nnz * ndof2))
      cnt = 0
      do j = 1, bt_ext(idom)%nr
        idx = exp_rows_index(idom-1) + j
        jrow = exp_rows_item(1,idx)
        if (jrow < 1 .or. btmat%nr < jrow) stop 'ERROR: extract BT_ext: jrow'
        ks = btmat%index(jrow-1)+1
        ke = btmat%index(jrow)
        do k = ks, ke
          kcol = btmat%item(k)
          if (incl_nz(k-nnz_int) == 0) cycle
          cnt = cnt + 1
          bt_ext(idom)%item(cnt) = kcol
          bt_ext(idom)%A(ndof2*(cnt-1)+1:ndof2*cnt) = btmat%A(ndof2*(k-1)+1:ndof2*k)
        enddo
        if (cnt /= bt_ext(idom)%index(j)) stop 'ERROR: extract BT_ext'
      enddo
      ! write(label,'(a,i0,a)') 'BT_ext(',idom,')'
      ! call debug_write_matrix(Bt_ext(idom), label, DEBUG_MATRIX)
    enddo
  end subroutine extract_bt_ext
 
  subroutine prepare_column_info(hecMESH, BT_ext, exp_cols_index, exp_cols_item)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_local_matrix), intent(in) :: bt_ext(:)
    integer(kind=kint), allocatable, intent(out) :: exp_cols_index(:)
    integer(kind=kint), allocatable, intent(out) :: exp_cols_item(:,:)
    !
    call make_exp_cols_index(hecmesh%n_neighbor_pe, bt_ext, exp_cols_index)
    if (debug >= 2) write(0,*) '  DEBUG2: make exp_cols_index done'
    if (debug >= 3) write(0,*) '    DEBUG3: exp_cols_index', exp_cols_index(0:hecmesh%n_neighbor_pe)
    !
    ! (col ID, rank, global ID)
    !
    call make_exp_cols_item(hecmesh, bt_ext, exp_cols_index, exp_cols_item)
    if (debug >= 2) write(0,*) '  DEBUG2: make exp_cols_item done'
    ! if (DEBUG >= 3) write(0,*) '    DEBUG3: exp_cols_item', exp_cols_item(1:cNCOL_ITEM,1:exp_cols_index(hecMESH%n_neighbor_pe))
  end subroutine prepare_column_info
 
  subroutine make_exp_cols_index(nnb, BT_ext, exp_cols_index)
    implicit none
    integer(kind=kint), intent(in) :: nnb
    type (hecmwst_local_matrix), intent(in) :: bt_ext(:)
    integer(kind=kint), allocatable, intent(out) :: exp_cols_index(:)
    integer(kind=kint) :: idom
    allocate(exp_cols_index(0:nnb))
    exp_cols_index(0) = 0
    do idom = 1, nnb
      exp_cols_index(idom) = exp_cols_index(idom-1) + bt_ext(idom)%nnz
    enddo
  end subroutine make_exp_cols_index
 
  subroutine make_exp_cols_item(hecMESH, BT_ext, exp_cols_index, exp_cols_item)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_local_matrix), intent(in) :: bt_ext(:)
    integer(kind=kint), allocatable, intent(in) :: exp_cols_index(:)
    integer(kind=kint), allocatable, intent(out) :: exp_cols_item(:,:)
    integer(kind=kint) :: cnt, idom, j, jcol
    allocate(exp_cols_item(cncol_item,exp_cols_index(hecmesh%n_neighbor_pe)))
    cnt = 0
    do idom = 1, hecmesh%n_neighbor_pe
      do j = 1, bt_ext(idom)%nnz
        cnt = cnt + 1
        jcol = bt_ext(idom)%item(j)
        ! if (DEBUG >= 3) write(0,*) '    DEBUG3: idom,j,cnt,jcol,nn_internal,n_node',&
        !      idom,j,cnt,jcol,hecMESH%nn_internal,hecMESH%n_node
        ! if (DEBUG >= 3) write(0,*) '    DEBUG3: size of exp_cols_item',size(exp_cols_item)
        ! if (DEBUG >= 3) write(0,*) '    DEBUG3: size of node_ID',size(hecMESH%node_ID)
        ! if (DEBUG >= 3) write(0,*) '    DEBUG3: size of global_node_ID',size(hecMESH%global_node_ID)
        exp_cols_item(clid,cnt) = hecmesh%node_ID(2*jcol-1)
        exp_cols_item(crank,cnt) = hecmesh%node_ID(2*jcol)
        if (cncol_item >= 3) exp_cols_item(cgid,cnt) = hecmesh%global_node_ID(jcol)
        ! if (DEBUG >= 3) write(0,*) '    DEBUG3: lid,rank(,gid)',exp_cols_item(1:cNCOL_ITEM,cnt)
      enddo
      if (cnt /= exp_cols_index(idom)) stop 'ERROR: make exp_cols_item'
    enddo
  end subroutine make_exp_cols_item
 
  subroutine send_bt_ext_and_recv_bt_int(hecMESH, exp_rows_index, exp_rows_item, BT_ext, &
       exp_cols_index, exp_cols_item, BT_int, hecMESHnew)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    integer(kind=kint), allocatable, intent(inout) :: exp_rows_index(:), exp_cols_index(:)
    integer(kind=kint), allocatable, intent(inout) :: exp_rows_item(:,:), exp_cols_item(:,:)
    type (hecmwst_local_matrix), allocatable, intent(inout) :: bt_ext(:)
    type (hecmwst_local_matrix), intent(out) :: bt_int
    type (hecmwst_local_mesh), intent(inout) :: hecmeshnew
    integer(kind=kint), allocatable :: imp_rows_index(:), imp_cols_index(:)
    integer(kind=kint), allocatable :: imp_rows_item(:,:), imp_cols_item(:,:)
    real(kind=kreal), allocatable :: imp_vals_item(:)
    integer(kind=kint), allocatable :: map(:), add_nodes(:,:)
    integer(kind=kint) :: ndof, ndof2, idom, n_add_node, i0
    if (hecmesh%n_neighbor_pe == 0) return
    ndof = bt_ext(1)%ndof
    ndof2 = ndof*ndof
    !
    call convert_rowid_to_remote_localid(hecmesh, exp_rows_index(hecmesh%n_neighbor_pe), exp_rows_item)
    if (debug >= 2) write(0,*) '  DEBUG2: convert rowID to remote localID done'
    !
    call send_recv_bt_ext_nr_nnz(hecmesh, bt_ext, imp_rows_index, imp_cols_index)
    if (debug >= 2) write(0,*) '  DEBUG2: send recv BT_ext nr and nnz done'
    !
    call send_recv_bt_ext_contents(hecmesh, bt_ext, &
         exp_rows_index, exp_cols_index, exp_rows_item, exp_cols_item, &
         imp_rows_index, imp_cols_index, &
         imp_rows_item, imp_cols_item, imp_vals_item)
    if (debug >= 2) write(0,*) '  DEBUG2: send recv BT_ext contents done'
    !
    do idom = 1, hecmesh%n_neighbor_pe
      call hecmw_localmat_free(bt_ext(idom))
    enddo
    deallocate(bt_ext)
    !
    call allocate_bt_int(hecmesh, ndof, imp_rows_index, imp_rows_item, bt_int)
    if (debug >= 2) write(0,*) '  DEBUG2: allocate BT_int done'
    !
    ! call copy_mesh(hecMESH, hecMESHnew)
    ! if (DEBUG >= 2) write(0,*) '  DEBUG2: copy mesh done'
    !
    call map_imported_cols(hecmeshnew, imp_cols_index(hecmesh%n_neighbor_pe), &
         imp_cols_item, n_add_node, add_nodes, map, i0)
    if (debug >= 2) write(0,*) '  DEBUG2: map imported cols done'
    !
    call update_comm_table(hecmeshnew, n_add_node, add_nodes, i0)
    if (debug >= 2) write(0,*) '  DEBUG2: update comm_table done'
    !
    bt_int%nc = hecmeshnew%n_node
    !
    call copy_vals_to_bt_int(hecmesh%n_neighbor_pe, imp_rows_index, imp_cols_index, &
         imp_rows_item, map, ndof2, imp_vals_item, bt_int)
    if (debug >= 2) write(0,*) '  DEBUG2: copy vals to BT_int done'
    !
    deallocate(imp_rows_index)
    deallocate(imp_cols_index)
    deallocate(imp_rows_item)
    deallocate(imp_cols_item)
    deallocate(imp_vals_item)
    deallocate(map)
    !
    call sort_and_uniq_rows(bt_int)
    if (debug >= 2) write(0,*) '  DEBUG2: sort and uniq rows of BT_int done'
  end subroutine send_bt_ext_and_recv_bt_int
 
  subroutine convert_rowid_to_remote_localid(hecMESH, len, exp_rows_item)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    integer(kind=kint), intent(in) :: len
    integer(kind=kint), intent(out) :: exp_rows_item(:,:)
    integer(kind=kint) :: i
    do i = 1, len
      exp_rows_item(1,i) = hecmesh%node_ID(2 * exp_rows_item(1,i) - 1)
    enddo
  end subroutine convert_rowid_to_remote_localid
 
  subroutine send_recv_bt_ext_nr_nnz(hecMESH, BT_ext, imp_rows_index, imp_cols_index)
    use m_hecmw_comm_f
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_local_matrix), intent(in) :: bt_ext(:)
    integer(kind=kint), allocatable, intent(out) :: imp_rows_index(:), imp_cols_index(:)
    integer(kind=kint) :: nnb, idom, irank, tag, recvbuf(2)
    integer(kind=kint), allocatable :: sendbuf(:,:)
    integer(kind=kint), allocatable :: requests(:)
    integer(kind=kint), allocatable :: statuses(:,:)
    nnb = hecmesh%n_neighbor_pe
    allocate(imp_rows_index(0:nnb))
    allocate(imp_cols_index(0:nnb))
    allocate(requests(nnb))
    allocate(statuses(hecmw_status_size, nnb))
    allocate(sendbuf(2,nnb))
    do idom = 1, nnb
      irank = hecmesh%neighbor_pe(idom)
      ! nr = exp_rows_per_rank(idom)
      sendbuf(1,idom) = bt_ext(idom)%nr
      sendbuf(2,idom) = bt_ext(idom)%nnz
      tag=2001
      call hecmw_isend_int(sendbuf(1,idom), 2, irank, tag, hecmesh%MPI_COMM, &
           requests(idom))
    enddo
    imp_rows_index(0) = 0
    imp_cols_index(0) = 0
    do idom = 1, nnb
      irank = hecmesh%neighbor_pe(idom)
      tag = 2001
      call hecmw_recv_int(recvbuf, 2, irank, tag, &
           hecmesh%MPI_COMM, statuses(:,1))
      imp_rows_index(idom) = imp_rows_index(idom-1) + recvbuf(1)
      imp_cols_index(idom) = imp_cols_index(idom-1) + recvbuf(2)
    enddo
    call hecmw_waitall(nnb, requests, statuses)
    deallocate(requests)
    deallocate(statuses)
    deallocate(sendbuf)
  end subroutine send_recv_bt_ext_nr_nnz
 
  subroutine send_recv_bt_ext_contents(hecMESH, BT_ext, &
       exp_rows_index, exp_cols_index, exp_rows_item, exp_cols_item, &
       imp_rows_index, imp_cols_index, &
       imp_rows_item, imp_cols_item, imp_vals_item)
    use m_hecmw_comm_f
    implicit none
    type (hecmwST_local_mesh), intent(in) :: hecMESH
    type (hecmwST_local_matrix), intent(in) :: BT_ext(:)
    integer(kind=kint), allocatable, intent(inout) :: exp_rows_index(:), exp_cols_index(:)
    integer(kind=kint), allocatable, intent(inout) :: exp_rows_item(:,:), exp_cols_item(:,:)
    integer(kind=kint), allocatable, intent(in) :: imp_rows_index(:), imp_cols_index(:)
    integer(kind=kint), allocatable, intent(out) :: imp_rows_item(:,:), imp_cols_item(:,:)
    real(kind=kreal), allocatable, intent(out) :: imp_vals_item(:)
    integer(kind=kint) :: nnb, ndof2, n_send, idom, irank, tag, nr, nnz
    integer(kind=kint), allocatable :: requests(:)
    integer(kind=kint), allocatable :: statuses(:,:)
    nnb = hecmesh%n_neighbor_pe
    if (nnb < 1) return
    ndof2 = bt_ext(1)%ndof ** 2
    allocate(imp_rows_item(2,imp_rows_index(nnb)))
    allocate(imp_cols_item(cncol_item,imp_cols_index(nnb)))
    allocate(imp_vals_item(ndof2*imp_cols_index(nnb)))
    allocate(requests(3*nnb))
    allocate(statuses(hecmw_status_size, 3*nnb))
    n_send = 0
    do idom = 1, nnb
      irank = hecmesh%neighbor_pe(idom)
      if (bt_ext(idom)%nr > 0) then
        n_send = n_send + 1
        tag = 2002
        call hecmw_isend_int(exp_rows_item(1,exp_rows_index(idom-1)+1), &
             2*bt_ext(idom)%nr, irank, tag, hecmesh%MPI_COMM, &
             requests(n_send))
        n_send = n_send + 1
        tag = 2003
        call hecmw_isend_int(exp_cols_item(1,exp_cols_index(idom-1)+1), &
             cncol_item*bt_ext(idom)%nnz, irank, tag, hecmesh%MPI_COMM, &
             requests(n_send))
        n_send = n_send + 1
        tag = 2004
        call hecmw_isend_r(bt_ext(idom)%A, ndof2*bt_ext(idom)%nnz, irank, &
             tag, hecmesh%MPI_COMM, requests(n_send))
      endif
    enddo
    do idom = 1, nnb
      irank = hecmesh%neighbor_pe(idom)
      nr = imp_rows_index(idom) - imp_rows_index(idom-1)
      nnz = imp_cols_index(idom) - imp_cols_index(idom-1)
      if (nr > 0) then
        tag = 2002
        call hecmw_recv_int(imp_rows_item(1,imp_rows_index(idom-1)+1), &
             2*nr, irank, tag, hecmesh%MPI_COMM, statuses(:,1))
        tag = 2003
        call hecmw_recv_int(imp_cols_item(1,imp_cols_index(idom-1)+1), &
             cncol_item*nnz, irank, tag, hecmesh%MPI_COMM, statuses(:,1))
        tag = 2004
        call hecmw_recv_r(imp_vals_item(ndof2*imp_cols_index(idom-1)+1), &
             ndof2*nnz, irank, tag, hecmesh%MPI_COMM, statuses(:,1))
      endif
    enddo
    call hecmw_waitall(n_send, requests, statuses)
    deallocate(exp_rows_index)
    deallocate(exp_rows_item)
    deallocate(exp_cols_index)
    deallocate(exp_cols_item)
  end subroutine send_recv_bt_ext_contents
 
  subroutine allocate_bt_int(hecMESH, ndof, imp_rows_index, imp_rows_item, BT_int)
    implicit none
    type (hecmwST_local_mesh), intent(in) :: hecMESH
    integer(kind=kint), intent(in) :: ndof
    integer(kind=kint), allocatable, intent(in) :: imp_rows_index(:), imp_rows_item(:,:)
    type (hecmwST_local_matrix), intent(out) :: BT_int
    integer(kind=kint), allocatable :: cnt(:)
    integer(kind=kint) :: idom, is, ie, i, irow, ncol, ndof2
    ndof2 = ndof*ndof
    bt_int%nr = hecmesh%nn_internal
    bt_int%nc = hecmesh%n_node
    bt_int%nnz = 0
    bt_int%ndof = ndof
    allocate(cnt(bt_int%nr))
    cnt(:) = 0
    do idom = 1, hecmesh%n_neighbor_pe
      is = imp_rows_index(idom-1)+1
      ie = imp_rows_index(idom)
      do i = is, ie
        irow = imp_rows_item(1,i)
        ncol = imp_rows_item(2,i)
        if (irow < 1 .or. bt_int%nr < irow) stop 'ERROR: allocate BT_int'
        cnt(irow) = cnt(irow) + ncol   !!! might include duplicate cols
      enddo
    enddo
    !
    allocate(bt_int%index(0:bt_int%nr))
    call make_index(bt_int%nr, cnt, bt_int%index)
    !
    bt_int%nnz = bt_int%index(bt_int%nr)
    allocate(bt_int%item(bt_int%nnz))
    allocate(bt_int%A(bt_int%nnz * ndof2))
    bt_int%A(:) = 0.d0
  end subroutine allocate_bt_int
 
  subroutine copy_mesh(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_dof = src%n_dof
    dst%n_neighbor_pe = src%n_neighbor_pe
    allocate(dst%neighbor_pe(dst%n_neighbor_pe))
    dst%neighbor_pe(:) = src%neighbor_pe(:)
    allocate(dst%import_index(0:dst%n_neighbor_pe))
    allocate(dst%export_index(0:dst%n_neighbor_pe))
    dst%import_index(:)= src%import_index(:)
    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(:)
    allocate(dst%node_ID(2*dst%n_node))
    dst%node_ID(1:2*dst%n_node) = src%node_ID(1:2*src%n_node)
    allocate(dst%global_node_ID(dst%n_node))
    dst%global_node_ID(1:dst%n_node) = src%global_node_ID(1:src%n_node)
    dst%mpc%n_mpc = 0
    dst%node => src%node
  end subroutine copy_mesh
 
  subroutine map_imported_cols(hecMESHnew, ncols, cols, n_add_node, add_nodes, map, i0)
    implicit none
    type (hecmwST_local_mesh), intent(inout) :: hecMESHnew
    integer(kind=kint), intent(in) :: ncols
    integer(kind=kint), intent(in) :: cols(cNCOL_ITEM,ncols)
    integer(kind=kint), allocatable, intent(out) :: map(:)
    integer(kind=kint), intent(out) :: n_add_node
    integer(kind=kint), allocatable, intent(out) :: add_nodes(:,:)
    integer(kind=kint), intent(out) :: i0
    allocate(map(ncols))
    !
    call map_present_nodes(hecmeshnew, ncols, cols, map, n_add_node)
    !
    ! add nodes == unmapped nodes
    !
    call extract_add_nodes(ncols, cols, map, n_add_node, add_nodes)
    !
    call append_nodes(hecmeshnew, n_add_node, add_nodes, i0)
    !
    call map_additional_nodes(ncols, cols, n_add_node, add_nodes, i0, map)
  end subroutine map_imported_cols
 
  subroutine map_present_nodes(hecMESH, ncols, cols, map, n_add_node)
    implicit none
    type (hecmwST_local_mesh), intent(in) :: hecMESH
    integer(kind=kint), intent(in) :: ncols
    integer(kind=kint), intent(in) :: cols(cNCOL_ITEM,ncols)
    integer(kind=kint), intent(out) :: map(ncols)
    integer(kind=kint), intent(out) :: n_add_node
    integer(kind=kint) :: i, j, lid, rank, llid, n_ext_node, idx
    integer(kind=kint), allocatable :: ext_node(:)
    type (hecmwST_pair_array) :: parray
    !
    call hecmw_pair_array_init(parray, hecmesh%n_node - hecmesh%nn_internal)
    do i = hecmesh%nn_internal + 1, hecmesh%n_node
      call hecmw_pair_array_append(parray, i, hecmesh%node_ID(2*i-1), hecmesh%node_ID(2*i))
    enddo
    call hecmw_pair_array_sort(parray)
    !
    n_add_node = 0
    n_ext_node = 0
    allocate(ext_node(ncols))
    !$omp parallel default(none), &
      !$omp& private(i,lid,rank,llid,idx,j), &
      !$omp& shared(ncols,hecMESH,cols,map,n_ext_node,ext_node,parray), &
      !$omp& reduction(+:n_add_node)
    !$omp do
    do i = 1, ncols
      lid = cols(clid,i)
      rank = cols(crank,i)
      !   check rank
      if (rank == hecmesh%my_rank) then  !   internal: set mapping
        map(i) = lid
      else                               !   external
        !$omp atomic capture
        n_ext_node = n_ext_node + 1
        idx = n_ext_node
        !$omp end atomic
        ext_node(idx) = i
      endif
    enddo
    !$omp end do
    !$omp do
    do j = 1, n_ext_node
      i = ext_node(j)
      lid = cols(clid,i)
      rank = cols(crank,i)
      !     search node_ID in external nodes
      llid = hecmw_pair_array_find_id(parray, lid, rank)
      if (llid > 0) then  !     found: set mapping
        map(i) = llid
      else                !     not found
        map(i) = -1
        n_add_node = n_add_node + 1
      endif
    enddo
    !$omp end do
    !$omp end parallel
    deallocate(ext_node)
    !
    call hecmw_pair_array_finalize(parray)
  end subroutine map_present_nodes
 
  subroutine extract_add_nodes(ncols, cols, map, n_add_node, add_nodes)
    implicit none
    integer(kind=kint), intent(in) :: ncols
    integer(kind=kint), intent(in) :: cols(cNCOL_ITEM,ncols), map(ncols)
    integer(kind=kint), intent(inout) :: n_add_node
    integer(kind=kint), allocatable, intent(out) :: add_nodes(:,:)
    integer(kind=kint) :: cnt, i
    allocate(add_nodes(cncol_item,n_add_node))
    cnt = 0
    do i = 1, ncols
      if (map(i) == -1) then
        cnt = cnt + 1
        add_nodes(1:cncol_item,cnt) = cols(1:cncol_item,i)
      endif
    enddo
    if (cnt /= n_add_node) stop 'ERROR: extract add_nodes'
    call sort_and_uniq_add_nodes(n_add_node, add_nodes)
  end subroutine extract_add_nodes
 
  subroutine sort_and_uniq_add_nodes(n_add_node, add_nodes)
    implicit none
    integer(kind=kint), intent(inout) :: n_add_node
    integer(kind=kint), intent(inout) :: add_nodes(cNCOL_ITEM,n_add_node)
    integer(kind=kint) :: ndup
    call sort_add_nodes(add_nodes, 1, n_add_node)
    call uniq_add_nodes(add_nodes, n_add_node, ndup)
    n_add_node = n_add_node - ndup
  end subroutine sort_and_uniq_add_nodes
 
  recursive subroutine sort_add_nodes(add_nodes, id1, id2)
    implicit none
    integer(kind=kint), intent(inout) :: add_nodes(:,:)
    integer(kind=kint), intent(in) :: id1, id2
    integer(kind=kint) :: center, left, right
    integer(kind=kint) :: pivot(cNCOL_ITEM), tmp(cNCOL_ITEM)
    if (id1 >= id2) return
    center = (id1 + id2) / 2
    pivot(1:cncol_item) = add_nodes(1:cncol_item,center)
    left = id1
    right = id2
    do
      do while ((add_nodes(crank,left) < pivot(crank)) .or. &
           (add_nodes(crank,left) == pivot(crank) .and. add_nodes(clid,left) < pivot(clid)))
        left = left + 1
      enddo
      do while ((pivot(crank) < add_nodes(crank,right)) .or. &
           (pivot(crank) == add_nodes(crank,right) .and. pivot(clid) < add_nodes(clid,right)))
        right = right - 1
      enddo
      if (left >= right) exit
      tmp(1:cncol_item) = add_nodes(1:cncol_item,left)
      add_nodes(1:cncol_item,left) = add_nodes(1:cncol_item,right)
      add_nodes(1:cncol_item,right) = tmp(1:cncol_item)
      left = left + 1
      right = right - 1
    enddo
    if (id1 < left-1) call sort_add_nodes(add_nodes, id1, left-1)
    if (right+1 < id2) call sort_add_nodes(add_nodes, right+1, id2)
    return
  end subroutine sort_add_nodes
 
  subroutine uniq_add_nodes(add_nodes, len, ndup)
    implicit none
    integer(kind=kint), intent(inout) :: add_nodes(:,:)
    integer(kind=kint), intent(in) :: len
    integer(kind=kint), intent(out) :: ndup
    integer(kind=kint) :: i
    ndup = 0
    do i = 2,len
      if (add_nodes(clid,i) == add_nodes(clid,i-1-ndup) .and. &
           add_nodes(crank,i) == add_nodes(crank,i-1-ndup)) then
        ndup = ndup + 1
      else if (ndup > 0) then
        add_nodes(1:cncol_item,i-ndup) = add_nodes(1:cncol_item,i)
      endif
    enddo
  end subroutine uniq_add_nodes
 
  subroutine search_add_nodes(n_add_node, add_nodes, rank, lid, idx)
    implicit none
    integer(kind=kint), intent(in) :: n_add_node
    integer(kind=kint), intent(in) :: add_nodes(cNCOL_ITEM,n_add_node)
    integer(kind=kint), intent(in) :: rank
    integer(kind=kint), intent(in) :: lid
    integer(kind=kint), intent(out) :: idx
    integer(kind=kint) :: left, right, center
    left = 1
    right = n_add_node
    do while (left <= right)
      center = (left + right) / 2
      if ((rank == add_nodes(crank,center)) .and. (lid == add_nodes(clid,center))) then
        idx = center
        return
      else if ((rank < add_nodes(crank,center)) .or. &
           (rank == add_nodes(crank,center) .and. lid < add_nodes(clid,center))) then
        right = center - 1
      else if ((add_nodes(crank,center) < rank) .or. &
           (add_nodes(crank,center) == rank .and. add_nodes(clid,center) < lid)) then
        left = center + 1
      endif
    end do
    idx = -1
  end subroutine search_add_nodes
 
  subroutine append_nodes(hecMESHnew, n_add_node, add_nodes, i0)
    implicit none
    type (hecmwST_local_mesh), intent(inout) :: hecMESHnew
    integer(kind=kint), intent(in) :: n_add_node
    integer(kind=kint), intent(in) :: add_nodes(:,:)
    integer(kind=kint), intent(out) :: i0
    integer(kind=kint) :: n_node, i, ii
    integer(kind=kint), pointer :: node_ID(:), global_node_ID(:)
    i0 = hecmeshnew%n_node
    n_node = hecmeshnew%n_node + n_add_node
    allocate(node_id(2*n_node))
    allocate(global_node_id(n_node))
    do i = 1, hecmeshnew%n_node
      node_id(2*i-1) = hecmeshnew%node_ID(2*i-1)
      node_id(2*i  ) = hecmeshnew%node_ID(2*i  )
      global_node_id(i) = hecmeshnew%global_node_ID(i)
    enddo
    do i = 1, n_add_node
      ii = hecmeshnew%n_node + i
      node_id(2*ii-1) = add_nodes(clid,i)
      node_id(2*ii  ) = add_nodes(crank,i)
      if (cncol_item >= 3) then
        global_node_id(ii) = add_nodes(cgid,i)
      else
        global_node_id(ii) = -1
      endif
    enddo
    deallocate(hecmeshnew%node_ID)
    deallocate(hecmeshnew%global_node_ID)
    hecmeshnew%n_node = n_node
    hecmeshnew%node_ID => node_id
    hecmeshnew%global_node_ID => global_node_id
  end subroutine append_nodes
 
  subroutine map_additional_nodes(ncols, cols, n_add_node, add_nodes, i0, map)
    implicit none
    integer(kind=kint), intent(in) :: ncols
    integer(kind=kint), intent(in) :: cols(cNCOL_ITEM,ncols)
    integer(kind=kint), intent(in) :: n_add_node
    integer(kind=kint), intent(in) :: add_nodes(cNCOL_ITEM,n_add_node)
    integer(kind=kint), intent(in) :: i0
    integer(kind=kint), intent(inout) :: map(ncols)
    integer(kind=kint) :: i, j
    do i = 1, ncols
      if (map(i) > 0) cycle
      call search_add_nodes(n_add_node, add_nodes, cols(crank,i), cols(clid,i), j)
      if (j == -1) stop 'ERROR: map_additional_nodes'
      map(i) = i0 + j
    enddo
  end subroutine map_additional_nodes
 
  subroutine update_comm_table(hecMESHnew, n_add_node, add_nodes, i0)
    use m_hecmw_comm_f
    implicit none
    type (hecmwST_local_mesh), intent(inout) :: hecMESHnew
    integer(kind=kint), intent(in) :: n_add_node
    integer(kind=kint), allocatable, intent(inout) :: add_nodes(:,:)
    integer(kind=kint), intent(in) :: i0
    integer(kind=kint), allocatable :: n_add_imp(:), add_imp_index(:)
    integer(kind=kint), allocatable :: add_imp_item_remote(:), add_imp_item_local(:)
    integer(kind=kint), allocatable :: n_add_exp(:), add_exp_index(:), add_exp_item(:)
    integer(kind=kint), allocatable :: n_new_imp(:), n_new_exp(:)
    integer(kind=kint) :: npe, nnb, comm, new_nnb
    integer(kind=kint), pointer :: nbpe(:), new_nbpe(:)
    integer(kind=kint), pointer :: import_index(:), export_index(:), import_item(:), export_item(:)
    integer(kind=kint), pointer :: new_import_index(:), new_export_index(:)
    integer(kind=kint), pointer :: new_import_item(:), new_export_item(:)
    npe = hecmeshnew%PETOT
    nnb = hecmeshnew%n_neighbor_pe
    comm = hecmeshnew%MPI_COMM
    nbpe => hecmeshnew%neighbor_pe
    import_index => hecmeshnew%import_index
    export_index => hecmeshnew%export_index
    import_item => hecmeshnew%import_item
    export_item => hecmeshnew%export_item
    !
    call count_add_imp_per_rank(n_add_node, add_nodes, npe, n_add_imp)
    if (debug >= 3) write(0,*) '    DEBUG3: count add_imp per rank done'
    !
    allocate(add_imp_index(0:npe))
    call make_index(npe, n_add_imp, add_imp_index)
    if (debug >= 3) write(0,*) '    DEBUG3: make add_imp_index done'
    !
    call make_add_imp_item(n_add_node, add_nodes, npe, i0, add_imp_index, &
         add_imp_item_remote, add_imp_item_local)
    if (debug >= 3) write(0,*) '    DEBUG3: make add_imp_item done'
    !
    deallocate(add_nodes)
    !
    ! all_to_all n_add_imp -> n_add_exp
    !
    allocate(n_add_exp(npe))
    call hecmw_alltoall_int(n_add_imp, 1, n_add_exp, 1, comm)
    if (debug >= 3) write(0,*) '    DEBUG3: alltoall n_add_imp to n_add_exp done'
    !
    allocate(add_exp_index(0:npe))
    call make_index(npe, n_add_exp, add_exp_index)
    if (debug >= 3) write(0,*) '    DEBUG3: make add_exp_index done'
    !
    call send_recv_add_imp_exp_item(npe, add_imp_index, add_imp_item_remote, &
         add_exp_index, add_exp_item, comm)
    if (debug >= 3) write(0,*) '    DEBUG3: send recv add_imp/exp_item done'
    !
    ! count new import
    !
    call count_new_comm_nodes(npe, nnb, nbpe, import_index, n_add_imp, n_new_imp)
    if (debug >= 3) write(0,*) '    DEBUG3: count new comm_nodes (import) done'
    !
    ! count new export
    !
    call count_new_comm_nodes(npe, nnb, nbpe, export_index, n_add_exp, n_new_exp)
    if (debug >= 3) write(0,*) '    DEBUG3: count new comm_nodes (export) done'
    !
    call update_neighbor_pe(npe, n_new_imp, n_new_exp, new_nnb, new_nbpe)
    if (debug >= 3) write(0,*) '    DEBUG3: update neighbor_pe done'
    !
    ! merge import table: import
    !
    call merge_comm_table(npe, nnb, nbpe, import_index, import_item, &
         new_nnb, new_nbpe, add_imp_index, add_imp_item_local, n_add_imp, n_new_imp, &
         new_import_index, new_import_item)
    if (debug >= 3) write(0,*) '    DEBUG3: merge comm_table (import) done'
    !
    deallocate(n_add_imp)
    deallocate(add_imp_index)
    deallocate(add_imp_item_remote, add_imp_item_local)
    deallocate(n_new_imp)
    !
    ! merge export table: export
    !
    call merge_comm_table(npe, nnb, nbpe, export_index, export_item, &
         new_nnb, new_nbpe, add_exp_index, add_exp_item, n_add_exp, n_new_exp, &
         new_export_index, new_export_item)
    if (debug >= 3) write(0,*) '    DEBUG3: merge comm_table (export) done'
    !
    deallocate(n_add_exp)
    deallocate(add_exp_index)
    deallocate(add_exp_item)
    deallocate(n_new_exp)
    !
    deallocate(nbpe)
    deallocate(import_index,import_item)
    deallocate(export_index,export_item)
    hecmeshnew%n_neighbor_pe = new_nnb
    hecmeshnew%neighbor_pe => new_nbpe
    hecmeshnew%import_index => new_import_index
    hecmeshnew%export_index => new_export_index
    hecmeshnew%import_item => new_import_item
    hecmeshnew%export_item => new_export_item
  end subroutine update_comm_table
 
  subroutine count_add_imp_per_rank(n_add_node, add_nodes, npe, n_add_imp)
    implicit none
    integer(kind=kint), intent(in) :: n_add_node
    integer(kind=kint), intent(in) :: add_nodes(cNCOL_ITEM,n_add_node)
    integer(kind=kint), intent(in) :: npe
    integer(kind=kint), allocatable, intent(out) :: n_add_imp(:)
    integer(kind=kint) :: i, rank
    allocate(n_add_imp(npe))
    n_add_imp(:) = 0
    do i = 1, n_add_node
      rank = add_nodes(crank,i)
      n_add_imp(rank+1) = n_add_imp(rank+1) + 1
    enddo
  end subroutine count_add_imp_per_rank
 
  subroutine make_add_imp_item(n_add_node, add_nodes, npe, i0, add_imp_index, &
       add_imp_item_remote, add_imp_item_local)
    implicit none
    integer(kind=kint), intent(in) :: n_add_node
    integer(kind=kint), intent(in) :: add_nodes(cNCOL_ITEM,n_add_node)
    integer(kind=kint), intent(in) :: npe, i0
    integer(kind=kint), allocatable, intent(in) :: add_imp_index(:)
    integer(kind=kint), allocatable, intent(out) :: add_imp_item_remote(:), add_imp_item_local(:)
    integer(kind=kint), allocatable :: cnt(:)
    integer(kind=kint) :: i, lid, rank, ipe
    allocate(add_imp_item_remote(add_imp_index(npe)))
    allocate(add_imp_item_local(add_imp_index(npe)))
    allocate(cnt(npe))
    cnt(:) = 0
    do i = 1, n_add_node
      lid = add_nodes(clid,i)
      rank = add_nodes(crank,i)
      ipe = rank + 1
      cnt(ipe) = cnt(ipe) + 1
      add_imp_item_remote(add_imp_index(ipe-1) + cnt(ipe)) = lid
      add_imp_item_local(add_imp_index(ipe-1) + cnt(ipe)) = i0 + i
    enddo
    deallocate(cnt)
  end subroutine make_add_imp_item
 
  subroutine send_recv_add_imp_exp_item(npe, add_imp_index, add_imp_item_remote, &
       add_exp_index, add_exp_item, mpi_comm)
    use m_hecmw_comm_f
    implicit none
    integer(kind=kint), intent(in) :: npe
    integer(kind=kint), allocatable, intent(in) :: add_imp_index(:), add_imp_item_remote(:)
    integer(kind=kint), allocatable, intent(in) :: add_exp_index(:)
    integer(kind=kint), allocatable, intent(out) :: add_exp_item(:)
    integer(kind=kint), intent(in) :: mpi_comm
    integer(kind=kint) :: n_send, i, irank, is, ie, len, tag
    integer(kind=kint), allocatable :: requests(:)
    integer(kind=kint), allocatable :: statuses(:,:)
    allocate(add_exp_item(add_exp_index(npe)))
    allocate(requests(npe))
    allocate(statuses(hecmw_status_size, npe))
    n_send = 0
    do i = 1, npe
      irank = i-1
      is = add_imp_index(i-1)+1
      ie = add_imp_index(i)
      len = ie - is + 1
      if (len == 0) cycle
      tag = 4001
      n_send = n_send + 1
      call hecmw_isend_int(add_imp_item_remote(is:ie), len, irank, tag, &
           mpi_comm, requests(n_send))
    enddo
    !
    do i = 1, npe
      irank = i-1
      is = add_exp_index(i-1)+1
      ie = add_exp_index(i)
      len = ie - is + 1
      if (len == 0) cycle
      tag = 4001
      call hecmw_recv_int(add_exp_item(is:ie), len, irank, tag, &
           mpi_comm, statuses(:,1))
    enddo
    call hecmw_waitall(n_send, requests, statuses)
  end subroutine send_recv_add_imp_exp_item
 
  subroutine count_new_comm_nodes(npe, org_nnb, org_nbpe, org_index, n_add, n_new)
    implicit none
    integer(kind=kint), intent(in) :: npe, org_nnb
    !integer(kind=kint), intent(in) :: org_nbpe(org_nnb), org_index(0:org_nnb), n_add(npe)
    integer(kind=kint), pointer, intent(in) :: org_nbpe(:), org_index(:)
    integer(kind=kint), intent(in) ::  n_add(:)
    integer(kind=kint), allocatable, intent(out) :: n_new(:)
    integer(kind=kint) :: i, irank, n_org
    allocate(n_new(npe))
    n_new(:) = n_add(:)
    do i = 1, org_nnb
      irank = org_nbpe(i)
      n_org = org_index(i) - org_index(i-1)
      n_new(irank+1) = n_new(irank+1) + n_org
    enddo
  end subroutine count_new_comm_nodes
 
  subroutine update_neighbor_pe(npe, n_new_imp, n_new_exp, &
       new_nnb, new_nbpe)
    implicit none
    integer(kind=kint), intent(in) :: npe
    integer(kind=kint), intent(in) :: n_new_imp(npe), n_new_exp(npe)
    integer(kind=kint), intent(out) :: new_nnb
    integer(kind=kint), pointer, intent(out) :: new_nbpe(:)
    integer(kind=kint) :: i
    new_nnb = 0
    do i = 1, npe
      if (n_new_imp(i) > 0 .or. n_new_exp(i) > 0) new_nnb = new_nnb+1
    enddo
    allocate(new_nbpe(new_nnb))
    new_nnb = 0
    do i = 1, npe
      if (n_new_imp(i) > 0 .or. n_new_exp(i) > 0) then
        new_nnb = new_nnb+1
        new_nbpe(new_nnb) = i-1
      endif
    enddo
  end subroutine update_neighbor_pe
 
  subroutine merge_comm_table(npe, org_nnb, org_nbpe, org_index, org_item, &
       new_nnb, new_nbpe, add_index, add_item, n_add, n_new, new_index, new_item)
    implicit none
    integer(kind=kint), intent(in) :: npe, org_nnb
    !integer(kind=kint), intent(in) :: org_nbpe(org_nnb), org_index(0:org_nnb), org_item(:)
    integer(kind=kint), pointer, intent(in) :: org_nbpe(:), org_index(:), org_item(:)
    integer(kind=kint), intent(in) :: new_nnb
    !integer(kind=kint), intent(in) :: new_nbpe(new_nnb), add_index(0:npe), add_item(:)
    integer(kind=kint), pointer, intent(in) :: new_nbpe(:)
    integer(kind=kint), allocatable, intent(in) :: add_index(:), add_item(:)
    integer(kind=kint), intent(in) :: n_add(npe), n_new(npe)
    integer(kind=kint), pointer, intent(out) :: new_index(:), new_item(:)
    integer(kind=kint), allocatable :: cnt(:)
    integer(kind=kint) :: i, irank, j, jrank, i0, j0, len
    ! if (associated(new_index)) deallocate(new_index)
    ! if (associated(new_item)) deallocate(new_item)
    allocate(new_index(0:new_nnb))
    new_index(0) = 0
    do i = 1, new_nnb
      irank = new_nbpe(i)
      new_index(i) = new_index(i-1) + n_new(irank+1)
    enddo
    allocate(new_item(new_index(new_nnb)))
    allocate(cnt(npe))
    cnt(:) = 0
    j = 1
    jrank = new_nbpe(j)
    do i = 1, org_nnb
      if (org_index(i) - org_index(i-1) == 0) cycle
      irank = org_nbpe(i)
      do while (jrank < irank)
        j = j + 1
        if (j > new_nnb) exit
        jrank = new_nbpe(j)
      enddo
      if (jrank /= irank) stop 'ERROR: merging comm table: org into new'
      i0 = org_index(i-1)
      len = org_index(i) - i0
      j0 = new_index(j-1)
      new_item(j0+1:j0+len) = org_item(i0+1:i0+len)
      cnt(jrank+1) = len
    enddo
    j = 1
    jrank = new_nbpe(j)
    do i = 1, npe
      if (n_add(i) == 0) cycle
      irank = i-1
      do while (jrank < irank)
        j = j + 1
        jrank = new_nbpe(j)
      enddo
      if (jrank /= irank) stop 'ERROR: merging comm table: add into new'
      i0 = add_index(i-1)
      len = add_index(i) - i0
      j0 = new_index(j-1) + cnt(jrank+1)
      new_item(j0+1:j0+len) = add_item(i0+1:i0+len)
      cnt(jrank+1) = cnt(jrank+1) + len
      if (cnt(jrank+1) /= new_index(j)-new_index(j-1)) stop 'ERROR: merging comm table'
    enddo
    deallocate(cnt)
  end subroutine merge_comm_table
 
  subroutine copy_vals_to_bt_int(nnb, imp_rows_index, imp_cols_index, &
       imp_rows_item, map, ndof2, imp_vals_item, BT_int)
    implicit none
    integer(kind=kint), intent(in) :: nnb
    integer(kind=kint), allocatable, intent(in) :: imp_rows_index(:), imp_cols_index(:)
    integer(kind=kint), intent(in) :: imp_rows_item(:,:), map(:)
    integer(kind=kint), intent(in) :: ndof2
    real(kind=kreal), intent(in) :: imp_vals_item(:)
    type (hecmwST_local_matrix), intent(inout) :: BT_int
    integer(kind=kint), allocatable :: cnt(:)
    integer(kind=kint) :: idom, is, ie, ic0, i, irow, ncol, j0, j
    allocate(cnt(bt_int%nr))
    cnt(:) = 0
    do idom = 1, nnb
      is = imp_rows_index(idom-1)+1
      ie = imp_rows_index(idom)
      ic0 = imp_cols_index(idom-1)
      do i = is, ie
        irow = imp_rows_item(1,i)
        ncol = imp_rows_item(2,i)
        if (irow < 1 .or. bt_int%nr < irow) stop 'ERROR: copy vals to BT_int: irow'
        j0 = bt_int%index(irow-1) + cnt(irow)
        do j = 1, ncol
          bt_int%item(j0+j) = map(ic0+j)
          bt_int%A(ndof2*(j0+j-1)+1:ndof2*(j0+j)) = imp_vals_item(ndof2*(ic0+j-1)+1:ndof2*(ic0+j))
        enddo
        cnt(irow) = cnt(irow) + ncol
        ic0 = ic0 + ncol
      enddo
      if (ic0 /= imp_cols_index(idom)) stop 'ERROR: copy vals to BT_int: ic0'
    enddo
    deallocate(cnt)
  end subroutine copy_vals_to_bt_int
 
  subroutine sort_and_uniq_rows(BTmat)
    use hecmw_array_util
    implicit none
    type (hecmwST_local_matrix), intent(inout) :: BTmat
    integer(kind=kint) :: nr, ndof, ndof2
    integer(kind=kint) :: irow, is, ie, is_new, ie_new, i, i_new
    integer(kind=kint) :: ndup, ndup_tot
    integer(kind=kint) :: js, je, js_new, je_new
    integer(kind=kint) :: new_nnz
    integer(kind=kint), allocatable :: cnt(:)
    integer(kind=kint), pointer :: sort_item(:), new_index(:), new_item(:)
    real(kind=kreal), pointer :: new_a(:)
    logical :: sorted
    real(kind=kreal) :: t0, t1
    t0 = hecmw_wtime()
    nr = btmat%nr
    ! check if already sorted
    sorted = .true.
    outer: do irow = 1, nr
      is = btmat%index(irow-1)+1
      ie = btmat%index(irow)
      do i = is, ie-1
        if (btmat%item(i) >= btmat%item(i+1)) then
          sorted = .false.
          exit outer
        endif
      enddo
    end do outer
    t1 = hecmw_wtime()
    if (timer >= 4) write(0, '(A,f10.4,L2)') "####### sort_and_uniq_rows (1) : ",t1-t0,sorted
    t0 = hecmw_wtime()
    if (sorted) return
    ! perform sort
    ndof = btmat%ndof
    ndof2 = ndof*ndof
    ! duplicate item array (sort_item)
    allocate(sort_item(btmat%nnz))
    do i = 1, btmat%nnz
      sort_item(i) = btmat%item(i)
    enddo
    ! sort and uniq item for each row
    allocate(cnt(nr))
    ndup_tot = 0
    !$omp parallel do default(none), &
      !$omp& schedule(dynamic,1), &
      !$omp& private(irow,is,ie,ndup), &
      !$omp& shared(nr,BTmat,sort_item,cnt), &
      !$omp& reduction(+:ndup_tot)
    do irow = 1, nr
      is = btmat%index(irow-1)+1
      ie = btmat%index(irow)
      call hecmw_qsort_int_array(sort_item, is, ie)
      call hecmw_uniq_int_array(sort_item, is, ie, ndup)
      cnt(irow) = (ie-is+1) - ndup
      ndup_tot = ndup_tot + ndup
    enddo
    !$omp end parallel do
    t1 = hecmw_wtime()
    if (timer >= 4) write(0, '(A,f10.4,I5)') "####### sort_and_uniq_rows (2) : ",t1-t0,ndup_tot
    t0 = hecmw_wtime()
    ! make new index and item array (new_index, new_item)
    if (ndup_tot == 0) then
      new_index => btmat%index
      new_nnz = btmat%nnz
      new_item => sort_item
    else
      allocate(new_index(0:nr))
      call make_index(nr, cnt, new_index)
      new_nnz = new_index(nr)
      allocate(new_item(new_nnz))
      do irow = 1, nr
        is = btmat%index(irow-1)+1
        ie = is+cnt(irow)-1
        is_new = new_index(irow-1)+1
        ie_new = is_new+cnt(irow)-1
        new_item(is_new:ie_new) = sort_item(is:ie)
      enddo
      deallocate(sort_item)
    endif
    deallocate(cnt)
    t1 = hecmw_wtime()
    if (timer >= 4) write(0, '(A,f10.4)') "####### sort_and_uniq_rows (3) : ",t1-t0
    t0 = hecmw_wtime()
    ! allocate and clear value array (new_A)
    allocate(new_a(ndof2*new_nnz))
    new_a(:) = 0.d0
    ! copy/add value from old A to new A
    !$omp parallel do default(none), &
      !$omp& schedule(dynamic,1), &
      !$omp& private(irow,is,ie,is_new,ie_new,i,i_new,js,je,js_new,je_new), &
      !$omp& shared(nr,BTmat,new_index,new_item,ndof2,new_A)
    do irow = 1, nr
      is = btmat%index(irow-1)+1
      ie = btmat%index(irow)
      is_new = new_index(irow-1)+1
      ie_new = new_index(irow)
      ! for each item in row
      do i = is, ie
        ! find place in new item
        call hecmw_bsearch_int_array(new_item, is_new, ie_new, btmat%item(i), i_new)
        if (i_new == -1) stop 'ERROR: sort_and_uniq_rows'
        js = ndof2*(i-1)+1
        je = ndof2*i
        js_new = ndof2*(i_new-1)+1
        je_new = ndof2*i_new
        new_a(js_new:je_new) = new_a(js_new:je_new) + btmat%A(js:je)
      enddo
    enddo
    !$omp end parallel do
    t1 = hecmw_wtime()
    if (timer >= 4) write(0, '(A,f10.4)') "####### sort_and_uniq_rows (4) : ",t1-t0
    t0 = hecmw_wtime()
    ! deallocate/update nnz, index, item, A
    if (ndup_tot == 0) then
      deallocate(btmat%item)
      btmat%item => new_item
      deallocate(btmat%A)
      btmat%A => new_a
    else
      btmat%nnz = new_nnz
      deallocate(btmat%index)
      btmat%index => new_index
      deallocate(btmat%item)
      btmat%item => new_item
      deallocate(btmat%A)
      btmat%A => new_a
    endif
  end subroutine sort_and_uniq_rows
 
  subroutine hecmw_localmat_add(Amat, Bmat, Cmat)
    implicit none
    type (hecmwst_local_matrix), intent(in) :: amat
    type (hecmwst_local_matrix), intent(in) :: bmat
    type (hecmwst_local_matrix), intent(out) :: cmat
    integer(kind=kint) :: ndof, ndof2, nr, nc, i, icnt, js, je, j, jcol, idx, i0, k
    integer(kind=kint), allocatable :: iw(:)
    if (amat%ndof /= bmat%ndof) stop 'ERROR: hecmw_localmat_add: non-matching ndof'
    ndof = amat%ndof
    ndof2 = ndof*ndof
    nr = min(amat%nr, bmat%nr)
    nc = max(amat%nc, bmat%nc)
    cmat%ndof = ndof
    cmat%nr = nr
    cmat%nc = nc
    cmat%nnz = 0
    allocate(cmat%index(0:nr))
    cmat%index(0) = 0
    allocate(iw(nc))
    do i = 1, nr
      icnt = 0
      ! Amat
      js = amat%index(i-1)+1
      je = amat%index(i)
      do j = js, je
        jcol = amat%item(j)
        icnt = icnt + 1
        iw(icnt) = jcol
      enddo
      ! Bmat
      js = bmat%index(i-1)+1
      je = bmat%index(i)
      lj1: do j = js, je
        jcol = bmat%item(j)
        do k = 1, icnt
          if (iw(k) == jcol) cycle lj1
        enddo
        icnt = icnt + 1
        iw(icnt) = jcol
      enddo lj1
      cmat%index(i) = cmat%index(i-1) + icnt
    enddo
    cmat%nnz = cmat%index(nr)
    allocate(cmat%item(cmat%nnz))
    allocate(cmat%A(ndof2*cmat%nnz))
    do i = 1, nr
      i0 = cmat%index(i-1)
      icnt = 0
      ! Amat
      js = amat%index(i-1)+1
      je = amat%index(i)
      do j = js, je
        jcol = amat%item(j)
        icnt = icnt + 1
        idx = i0 + icnt
        cmat%item(idx) = jcol
        cmat%A(ndof2*(idx-1)+1:ndof2*idx) = amat%A(ndof2*(j-1)+1:ndof2*j)
      enddo
      ! Bmat
      js = bmat%index(i-1)+1
      je = bmat%index(i)
      lj2: do j = js, je
        jcol = bmat%item(j)
        do k = 1, icnt
          idx = i0 + k
          if (cmat%item(idx) == jcol) then
            cmat%A(ndof2*(idx-1)+1:ndof2*idx) = &
                 cmat%A(ndof2*(idx-1)+1:ndof2*idx) + bmat%A(ndof2*(j-1)+1:ndof2*j)
            cycle lj2
          endif
        enddo
        icnt = icnt + 1
        idx = i0 + icnt
        cmat%item(idx) = jcol
        cmat%A(ndof2*(idx-1)+1:ndof2*idx) = bmat%A(ndof2*(j-1)+1:ndof2*j)
      enddo lj2
      if (i0 + icnt /= cmat%index(i)) stop 'ERROR: merge localmat'
    enddo
    call sort_and_uniq_rows(cmat)
  end subroutine hecmw_localmat_add
 
  ! subroutine hecmw_localmat_add(Amat, Bmat, Cmat)
  !   implicit none
  !   type (hecmwST_local_matrix), intent(in) :: Amat
  !   type (hecmwST_local_matrix), intent(in) :: Bmat
  !   type (hecmwST_local_matrix), intent(out) :: Cmat
  !   integer(kind=kint) :: ndof, ndof2, nr, nc, i, js, je, j, jcol, nnz_row, idx, ks, ke, k, kcol
  !   if (Amat%ndof /= Bmat%ndof) stop 'ERROR: hecmw_localmat_add: non-matching ndof'
  !   ndof = Amat%ndof
  !   ndof2 = ndof*ndof
  !   nr = min(Amat%nr, Bmat%nr)
  !   nc = max(Amat%nc, Bmat%nc)
  !   Cmat%ndof = ndof
  !   Cmat%nr = nr
  !   Cmat%nc = nc
  !   Cmat%nnz = Amat%index(nr) + Bmat%index(nr)
  !   allocate(Cmat%index(0:nr))
  !   allocate(Cmat%item(Cmat%nnz))
  !   allocate(Cmat%A(ndof2 * Cmat%nnz))
  !   Cmat%index(0) = 0
  !   idx = 0
  !   do i = 1, nr
  !     ! Amat
  !     js = Amat%index(i-1)+1
  !     je = Amat%index(i)
  !     do j = js, je
  !       idx = idx + 1
  !       Cmat%item(idx) = Amat%item(j)
  !       Cmat%A(ndof2*(idx-1)+1:ndof2*idx) = Amat%A(ndof2*(j-1)+1:ndof2*j)
  !     enddo
  !     ! Bmat
  !     js = Bmat%index(i-1)+1
  !     je = Bmat%index(i)
  !     do j = js, je
  !       idx = idx + 1
  !       Cmat%item(idx) = Bmat%item(j)
  !       Cmat%A(ndof2*(idx-1)+1:ndof2*idx) = Bmat%A(ndof2*(j-1)+1:ndof2*j)
  !     enddo
  !     Cmat%index(i) = idx
  !   enddo
  !   if (Cmat%index(nr) /= Cmat%nnz) stop 'ERROR: merge localmat'
  !   call sort_and_uniq_rows(Cmat)
  ! end subroutine hecmw_localmat_add
 
  subroutine hecmw_localmat_init_with_hecmat(BKmat, hecMAT, num_lagrange)
    implicit none
    type (hecmwst_local_matrix), intent(inout) :: bkmat
    type (hecmwst_matrix), intent(in) :: hecmat
    integer(kind=kint), optional, intent(in) :: num_lagrange
    integer(kind=kint) :: ndof, ndof2, i, idx, idx2, js, je, j, k
    integer(kind=kint), allocatable :: incl_nz(:), cnt(:)
    logical :: check_nonzero
    check_nonzero = .false.
    !check_nonzero = .true.  !!! always checking nonzero seems to be faster
    !
    ndof = hecmat%NDOF
    ndof2 = ndof*ndof
    ! nr, nc, nnz
    bkmat%nr = hecmat%NP
    bkmat%nc = hecmat%NP
    bkmat%ndof = ndof
    !
    if (present(num_lagrange)) then       !!! TEMPORARY (DUE TO WRONG conMAT WHEN num_lagrange==0) !!!
      check_nonzero = .true.
    endif
    !
    if (check_nonzero) then
      allocate(incl_nz(hecmat%NPL + hecmat%NPU + hecmat%NP))
      allocate(cnt(bkmat%nr))
      incl_nz(:) = 0
      !$omp parallel default(none), &
        !$omp& private(i,idx,js,je,j,k), &
        !$omp& shared(BKmat,hecMAT,cnt,ndof2,incl_nz)
      !$omp do
      do i = 1, bkmat%nr
        idx = hecmat%indexL(i-1) + (i-1) + hecmat%indexU(i-1)
        cnt(i) = 0
        ! lower
        js = hecmat%indexL(i-1)+1
        je = hecmat%indexL(i)
        do j = js, je
          idx = idx + 1
          do k = 1, ndof2
            if (hecmat%AL(ndof2*(j-1)+k) /= 0.0d0) then
              incl_nz(idx) = 1
              cnt(i) = cnt(i) + 1
              exit
            endif
          enddo
        enddo
        ! diag
        idx = idx + 1
        do k = 1, ndof2
          if (hecmat%D(ndof2*(i-1)+k) /= 0.0d0) then
            incl_nz(idx) = 1
            cnt(i) = cnt(i) + 1
            exit
          endif
        enddo
        ! upper
        js = hecmat%indexU(i-1)+1
        je = hecmat%indexU(i)
        do j = js, je
          idx = idx + 1
          do k = 1, ndof2
            if (hecmat%AU(ndof2*(j-1)+k) /= 0.0d0) then
              incl_nz(idx) = 1
              cnt(i) = cnt(i) + 1
              exit
            endif
          enddo
        enddo
        if (idx /= hecmat%indexL(i) + i + hecmat%indexU(i)) stop 'ERROR: hecmw_localmat_init_with_hecmat: count'
      enddo
      !$omp end do
      !$omp end parallel
      ! index
      allocate(bkmat%index(0:bkmat%nr))
      call make_index(bkmat%nr, cnt, bkmat%index)
      deallocate(cnt)
      bkmat%nnz = bkmat%index(bkmat%nr)
      ! item, A
      allocate(bkmat%item(bkmat%nnz))
      allocate(bkmat%A(ndof2 * bkmat%nnz))
      !$omp parallel default(none), &
        !$omp& private(i,idx,idx2,js,je,j), &
        !$omp& shared(BKmat,hecMAT,ndof2,incl_nz)
      !$omp do
      do i = 1, bkmat%nr
        idx = hecmat%indexL(i-1) + (i-1) + hecmat%indexU(i-1)
        idx2 = bkmat%index(i-1)
        ! lower
        js = hecmat%indexL(i-1)+1
        je = hecmat%indexL(i)
        do j = js, je
          idx = idx + 1
          if (incl_nz(idx) == 1) then
            idx2 = idx2 + 1
            bkmat%item(idx2) = hecmat%itemL(j)
            bkmat%A(ndof2*(idx2-1)+1:ndof2*idx2) = hecmat%AL(ndof2*(j-1)+1:ndof2*j)
          endif
        enddo
        ! diag
        idx = idx + 1
        if (incl_nz(idx) == 1) then
          idx2 = idx2 + 1
          bkmat%item(idx2) = i
          bkmat%A(ndof2*(idx2-1)+1:ndof2*idx2) = hecmat%D(ndof2*(i-1)+1:ndof2*i)
        endif
        ! upper
        js = hecmat%indexU(i-1)+1
        je = hecmat%indexU(i)
        do j = js, je
          idx = idx + 1
          if (incl_nz(idx) == 1) then
            idx2 = idx2 + 1
            bkmat%item(idx2) = hecmat%itemU(j)
            bkmat%A(ndof2*(idx2-1)+1:ndof2*idx2) = hecmat%AU(ndof2*(j-1)+1:ndof2*j)
          endif
        enddo
        if (idx /= hecmat%indexL(i) + i + hecmat%indexU(i)) stop 'ERROR: hecmw_localmat_init_with_hecmat: copy'
        if (idx2 /= bkmat%index(i)) stop 'ERROR: hecmw_localmat_init_with_hecmat: index'
      enddo
      !$omp end do
      !$omp end parallel
      deallocate(incl_nz)
    else
      bkmat%nnz = hecmat%NPL + hecmat%NP + hecmat%NPU
      allocate(bkmat%index(0:bkmat%nr))
      allocate(bkmat%item(bkmat%nnz))
      allocate(bkmat%A(ndof2 * bkmat%nnz))
      bkmat%index(0) = 0
      !$omp parallel do default(none), &
        !$omp& private(i,idx,js,je,j), &
        !$omp& shared(BKmat,hecMAT,ndof2)
      do i = 1, bkmat%nr
        idx = hecmat%indexL(i-1) + (i-1) + hecmat%indexU(i-1)
        ! lower
        js = hecmat%indexL(i-1)+1
        je = hecmat%indexL(i)
        do j = js, je
          idx = idx + 1
          bkmat%item(idx) = hecmat%itemL(j)
          bkmat%A(ndof2*(idx-1)+1:ndof2*idx) = hecmat%AL(ndof2*(j-1)+1:ndof2*j)
        enddo
        ! diag
        idx = idx + 1
        bkmat%item(idx) = i
        bkmat%A(ndof2*(idx-1)+1:ndof2*idx) = hecmat%D(ndof2*(i-1)+1:ndof2*i)
        ! upper
        js = hecmat%indexU(i-1)+1
        je = hecmat%indexU(i)
        do j = js, je
          idx = idx + 1
          bkmat%item(idx) = hecmat%itemU(j)
          bkmat%A(ndof2*(idx-1)+1:ndof2*idx) = hecmat%AU(ndof2*(j-1)+1:ndof2*j)
        enddo
        bkmat%index(i) = idx
        if (idx /= hecmat%indexL(i) + i + hecmat%indexU(i)) stop 'ERROR: hecmw_localmat_init_with_hecmat: copy'
      enddo
      !$omp end parallel do
    endif
  end subroutine hecmw_localmat_init_with_hecmat
 
  subroutine hecmw_localmat_add_hecmat(BKmat, hecMAT)
    implicit none
    type (hecmwst_local_matrix), intent(inout) :: bkmat
    type (hecmwst_matrix), intent(in) :: hecmat
    type (hecmwst_local_matrix) :: w1mat, w2mat
    !! Should Be Simple If Non-Zero Profile Is Kept !!
    call hecmw_localmat_init_with_hecmat(w1mat, hecmat)
    call debug_write_matrix(w1mat, 'BKmat (hecMAT)', debug_matrix)
    call hecmw_localmat_add(bkmat, w1mat, w2mat)
    call hecmw_localmat_free(bkmat)
    call hecmw_localmat_free(w1mat)
    bkmat%nr = w2mat%nr
    bkmat%nc = w2mat%nc
    bkmat%nnz = w2mat%nnz
    bkmat%ndof = w2mat%ndof
    bkmat%index => w2mat%index
    bkmat%item => w2mat%item
    bkmat%A => w2mat%A
  end subroutine hecmw_localmat_add_hecmat
 
  subroutine hecmw_localmat_multmat(BKmat, BTmat, hecMESH, BKTmat)
    implicit none
    type (hecmwst_local_matrix), intent(in) :: bkmat
    type (hecmwst_local_matrix), intent(inout) :: btmat
    type (hecmwst_local_mesh), intent(inout) :: hecmesh
    type (hecmwst_local_matrix), intent(out) :: bktmat
    type (hecmwst_matrix_comm) :: heccomm
    type (hecmwst_local_mesh) :: hecmeshnew
    type (hecmwst_local_matrix), allocatable :: bt_exp(:)
    type (hecmwst_local_matrix) :: bt_imp, bt_all
    integer(kind=kint), allocatable :: exp_cols_index(:)
    integer(kind=kint), allocatable :: exp_cols_item(:,:)
    real(kind=kreal) :: t0, t1
    t0 = hecmw_wtime()
    !
    if (hecmesh%PETOT > 1) then
      call make_comm_table(bkmat, hecmesh, heccomm)
      if (debug >= 1) write(0,*) 'DEBUG: hecmw_localmat_multmat: make_comm_table done'
      t1 = hecmw_wtime()
      if (timer >= 2) write(0,'(A,f10.4)') '##### hecmw_localmat_multmat (1) : ',t1-t0
      t0 = hecmw_wtime()
      !
      if (btmat%nr > hecmesh%nn_internal) then
        ! consider only internal part of BTmat
        if (debug >= 1) write(0,'(A)') 'DEBUG: hecmw_localmat_multmat: ignore external part of BTmat'
        btmat%nr = hecmesh%nn_internal
        btmat%nnz = btmat%index(btmat%nr)
      endif
      !
      call extract_bt_exp(btmat, heccomm, bt_exp)
      if (debug >= 1) write(0,*) 'DEBUG: hecmw_localmat_multmat: extract_BT_exp done'
      t1 = hecmw_wtime()
      if (timer >= 2) write(0,'(A,f10.4)') '##### hecmw_localmat_multmat (2) : ',t1-t0
      t0 = hecmw_wtime()
      !
      call prepare_column_info(hecmesh, bt_exp, exp_cols_index, exp_cols_item)
      if (debug >= 1) write(0,*) 'DEBUG: hecmw_localmat_multmat: prepare column info done'
      t1 = hecmw_wtime()
      if (timer >= 2) write(0,'(A,f10.4)') '##### hecmw_localmat_multmat (3) : ',t1-t0
      t0 = hecmw_wtime()
      !
      call send_bt_exp_and_recv_bt_imp(hecmesh, heccomm, bt_exp, exp_cols_index, exp_cols_item, bt_imp, hecmeshnew)
      if (debug >= 1) write(0,*) 'DEBUG: hecmw_localmat_multmat: send BT_exp and recv BT_imp done'
      t1 = hecmw_wtime()
      if (timer >= 2) write(0,'(A,f10.4)') '##### hecmw_localmat_multmat (4) : ',t1-t0
      t0 = hecmw_wtime()
      call free_comm_table(heccomm)
      !
      call concat_btmat_and_bt_imp(btmat, bt_imp, bt_all)
      if (debug >= 1) write(0,*) 'DEBUG: hecmw_localmat_multmat: concat BTmat and BT_imp into BT_all done'
      t1 = hecmw_wtime()
      if (timer >= 2) write(0,'(A,f10.4)') '##### hecmw_localmat_multmat (5) : ',t1-t0
      t0 = hecmw_wtime()
      call hecmw_localmat_free(bt_imp)
      !
      call multiply_mat_mat(bkmat, bt_all, bktmat)
      if (debug >= 1) write(0,*) 'DEBUG: hecmw_localmat_multmat: multiply BKmat and BT_all into BKTmat done'
      t1 = hecmw_wtime()
      if (timer >= 2) write(0,'(A,f10.4)') '##### hecmw_localmat_multmat (6) : ',t1-t0
      t0 = hecmw_wtime()
      call hecmw_localmat_free(bt_all)
      !
      if (hecmesh%n_neighbor_pe > 0) then
        hecmesh%n_node = hecmeshnew%n_node
        hecmesh%n_neighbor_pe = hecmeshnew%n_neighbor_pe
        deallocate(hecmesh%neighbor_pe)
        deallocate(hecmesh%import_index)
        deallocate(hecmesh%export_index)
        deallocate(hecmesh%import_item)
        deallocate(hecmesh%export_item)
        deallocate(hecmesh%node_ID)
        deallocate(hecmesh%global_node_ID)
        hecmesh%neighbor_pe => hecmeshnew%neighbor_pe
        hecmesh%import_index => hecmeshnew%import_index
        hecmesh%export_index => hecmeshnew%export_index
        hecmesh%import_item => hecmeshnew%import_item
        hecmesh%export_item => hecmeshnew%export_item
        hecmesh%node_ID => hecmeshnew%node_ID
        hecmesh%global_node_ID => hecmeshnew%global_node_ID
        if (debug >= 1) write(0,*) 'DEBUG: hecmw_localmat_multmat: update hecMESH done'
        t1 = hecmw_wtime()
        if (timer >= 2) write(0,'(A,f10.4)') '##### hecmw_localmat_multmat (7) : ',t1-t0
      endif
    else
      call multiply_mat_mat(bkmat, btmat, bktmat)
      if (debug >= 1) write(0,*) 'DEBUG: hecmw_localmat_multmat: multiply BKmat and BTmat into BKTmat done'
      t1 = hecmw_wtime()
      if (timer >= 2) write(0,'(A,f10.4)') '##### hecmw_localmat_multmat : ',t1-t0
    endif
  end subroutine hecmw_localmat_multmat
 
  subroutine make_comm_table(BKmat, hecMESH, hecCOMM)
    use m_hecmw_comm_f
    implicit none
    type (hecmwst_local_matrix), intent(in) :: bkmat
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix_comm), intent(out) :: heccomm
    integer(kind=kint) :: nn_int, nn_ext, nnb, i, icol, irank, idom, idx, n_send, tag, js, je, len
    integer(kind=kint), allocatable :: is_nz_col(:), imp_cnt(:), exp_cnt(:), import_item_remote(:)
    integer(kind=kint), allocatable :: requests(:), statuses(:,:)
    heccomm%zero = hecmesh%zero
    heccomm%HECMW_COMM = hecmesh%MPI_COMM
    heccomm%PETOT = hecmesh%PETOT
    heccomm%PEsmpTOT = hecmesh%PEsmpTOT
    heccomm%my_rank = hecmesh%my_rank
    heccomm%errnof = hecmesh%errnof
    heccomm%n_subdomain = hecmesh%n_subdomain
    heccomm%n_neighbor_pe = hecmesh%n_neighbor_pe
    allocate(heccomm%neighbor_pe(heccomm%n_neighbor_pe))
    heccomm%neighbor_pe(:) = hecmesh%neighbor_pe(:)
    !
    nn_int = hecmesh%nn_internal
    nn_ext = hecmesh%n_node - hecmesh%nn_internal
    nnb = heccomm%n_neighbor_pe
    !
    ! check_external_nz_cols (by profile (not number))
    allocate(is_nz_col(nn_ext))
    is_nz_col(:) = 0
    do i = 1, bkmat%index(nn_int)
      icol = bkmat%item(i)
      if (icol > nn_int) is_nz_col(icol - nn_int) = 1
    enddo
    !
    ! count_nz_cols_per_rank
    allocate(imp_cnt(nnb))
    imp_cnt(:) = 0
    do i = 1, nn_ext
      if (is_nz_col(i) == 1) then
        irank = hecmesh%node_ID(2*(nn_int+i))
        call rank_to_idom(hecmesh, irank, idom)
        imp_cnt(idom) = imp_cnt(idom) + 1
      endif
    enddo
    if (debug >= 3) write(0,*) '    DEBUG3: imp_cnt',imp_cnt(:)
    !
    ! make_index
    allocate(heccomm%import_index(0:nnb))
    call make_index(nnb, imp_cnt, heccomm%import_index)
    if (debug >= 3) write(0,*) '    DEBUG3: import_index',heccomm%import_index(:)
    !
    ! fill item
    allocate(heccomm%import_item(heccomm%import_index(nnb)))
    imp_cnt(:) = 0
    do i = 1, nn_ext
      if (is_nz_col(i) == 1) then
        irank = hecmesh%node_ID(2*(nn_int+i))
        call rank_to_idom(hecmesh, irank, idom)
        imp_cnt(idom) = imp_cnt(idom) + 1
        idx = heccomm%import_index(idom-1)+imp_cnt(idom)
        heccomm%import_item(idx) = nn_int+i
      endif
    enddo
    if (debug >= 3) write(0,*) '    DEBUG3: import_item',heccomm%import_item(:)
    !
    allocate(import_item_remote(heccomm%import_index(nnb)))
    do i = 1, heccomm%import_index(nnb)
      import_item_remote(i) = hecmesh%node_ID(2*heccomm%import_item(i)-1)
    enddo
    if (debug >= 3) write(0,*) '    DEBUG3: import_item_remote',import_item_remote(:)
    !
    allocate(requests(2*nnb))
    allocate(statuses(hecmw_status_size, 2*nnb))
    !
    ! send/recv
    n_send = 0
    do idom = 1, nnb
      irank = heccomm%neighbor_pe(idom)
      n_send = n_send + 1
      tag = 6001
      call hecmw_isend_int(imp_cnt(idom), 1, irank, tag, heccomm%HECMW_COMM, requests(n_send))
      if (imp_cnt(idom) > 0) then
        js = heccomm%import_index(idom-1)+1
        je = heccomm%import_index(idom)
        len = je-js+1
        n_send = n_send + 1
        tag = 6002
        call hecmw_isend_int(import_item_remote(js:je), len, irank, tag, &
             heccomm%HECMW_COMM, requests(n_send))
      endif
    enddo
    !
    ! index
    allocate(exp_cnt(nnb))
    do idom = 1, nnb
      irank = heccomm%neighbor_pe(idom)
      tag = 6001
      call hecmw_recv_int(exp_cnt(idom), 1, irank, tag, heccomm%HECMW_COMM, statuses(:,1))
    enddo
    allocate(heccomm%export_index(0:nnb))
    call make_index(nnb, exp_cnt, heccomm%export_index)
    if (debug >= 3) write(0,*) '    DEBUG3: export_index',heccomm%export_index(:)
    !
    ! item
    allocate(heccomm%export_item(heccomm%export_index(nnb)))
    do idom = 1, nnb
      if (exp_cnt(idom) <= 0) cycle
      irank = heccomm%neighbor_pe(idom)
      js = heccomm%export_index(idom-1)+1
      je = heccomm%export_index(idom)
      len = je-js+1
      tag = 6002
      call hecmw_recv_int(heccomm%export_item(js:je), len, irank, tag, &
           heccomm%HECMW_COMM, statuses(:,1))
    enddo
    if (debug >= 3) write(0,*) '    DEBUG3: export_item',heccomm%export_item(:)
    call hecmw_waitall(n_send, requests, statuses)
    !
    deallocate(imp_cnt)
    deallocate(exp_cnt)
    deallocate(import_item_remote)
  end subroutine make_comm_table
 
  subroutine free_comm_table(hecCOMM)
    implicit none
    type (hecmwST_matrix_comm), intent(inout) :: hecCOMM
    deallocate(heccomm%neighbor_pe)
    deallocate(heccomm%import_index)
    deallocate(heccomm%import_item)
    deallocate(heccomm%export_index)
    deallocate(heccomm%export_item)
  end subroutine free_comm_table
 
  subroutine extract_bt_exp(BTmat, hecCOMM, BT_exp)
    implicit none
    type (hecmwST_local_matrix), intent(in) :: BTmat
    type (hecmwST_matrix_comm), intent(in) :: hecCOMM
    type (hecmwST_local_matrix), allocatable, intent(out) :: BT_exp(:)
    integer(kind=kint) :: ndof, ndof2, idom, idx_0, idx_n, j, jrow, nnz_row, idx, ks, ke, k
    if (heccomm%n_neighbor_pe == 0) return
    allocate(bt_exp(heccomm%n_neighbor_pe))
    ndof = btmat%ndof
    ndof2 = ndof * ndof
    do idom = 1, heccomm%n_neighbor_pe
      idx_0 = heccomm%export_index(idom-1)
      idx_n = heccomm%export_index(idom)
      bt_exp(idom)%nr = idx_n - idx_0
      bt_exp(idom)%nc = btmat%nc
      bt_exp(idom)%nnz = 0
      bt_exp(idom)%ndof = ndof
      allocate(bt_exp(idom)%index(0:bt_exp(idom)%nr))
      bt_exp(idom)%index(0) = 0
      do j = 1, bt_exp(idom)%nr
        jrow = heccomm%export_item(idx_0 + j)
        nnz_row = btmat%index(jrow) - btmat%index(jrow-1)
        bt_exp(idom)%index(j) = bt_exp(idom)%index(j-1) + nnz_row
      enddo
      bt_exp(idom)%nnz = bt_exp(idom)%index(bt_exp(idom)%nr)
      allocate(bt_exp(idom)%item(bt_exp(idom)%nnz))
      allocate(bt_exp(idom)%A(ndof2 * bt_exp(idom)%nnz))
      idx = 0
      do j = 1, bt_exp(idom)%nr
        jrow = heccomm%export_item(idx_0 + j)
        ks = btmat%index(jrow-1) + 1
        ke = btmat%index(jrow)
        do k = ks, ke
          idx = idx + 1
          bt_exp(idom)%item(idx) = btmat%item(k)
          bt_exp(idom)%A(ndof2*(idx-1)+1:ndof2*idx) = btmat%A(ndof2*(k-1)+1:ndof2*k)
        enddo
        if (idx /= bt_exp(idom)%index(j)) stop 'ERROR: extract BT_exp'
      enddo
    enddo
  end subroutine extract_bt_exp
 
  subroutine send_bt_exp_and_recv_bt_imp(hecMESH, hecCOMM, BT_exp, exp_cols_index, exp_cols_item, BT_imp, hecMESHnew)
    use m_hecmw_comm_f
    implicit none
    type (hecmwST_local_mesh), intent(in) :: hecMESH
    type (hecmwST_matrix_comm), intent(in) :: hecCOMM
    type (hecmwST_local_matrix), allocatable, intent(inout) :: BT_exp(:)
    integer(kind=kint), allocatable, intent(inout) :: exp_cols_index(:)
    integer(kind=kint), allocatable, intent(inout) :: exp_cols_item(:,:)
    type (hecmwST_local_matrix), intent(out) :: BT_imp
    type (hecmwST_local_mesh), intent(inout) :: hecMESHnew
    integer(kind=kint), allocatable :: nnz_imp(:), cnt(:), index_imp(:)
    integer(kind=kint), allocatable :: imp_cols_index(:)
    integer(kind=kint), allocatable :: imp_cols_item(:,:)
    real(kind=kreal), allocatable :: imp_vals_item(:)
    integer(kind=kint) :: nnb, ndof, ndof2, idom, irank, nr, n_send, tag, idx_0, idx_n, j, jj, nnz
    integer(kind=kint), allocatable :: requests(:)
    integer(kind=kint), allocatable :: statuses(:,:)
    integer(kind=kint), allocatable :: map(:), add_nodes(:,:)
    integer(kind=kint) :: n_add_node, i0
    nnb = heccomm%n_neighbor_pe
    if (nnb == 0) then
      bt_imp%nr = 0
      bt_imp%nc = 0
      bt_imp%nnz = 0
      bt_imp%ndof = 0
      allocate(bt_imp%index(0:0))
      bt_imp%index(0) = 0
      return
    endif
    ndof = bt_exp(1)%ndof
    ndof2 = ndof*ndof
    allocate(requests(nnb*3))
    allocate(statuses(hecmw_status_size, nnb*3))
    n_send = 0
    do idom = 1, nnb
      irank = heccomm%neighbor_pe(idom)
      nr = bt_exp(idom)%nr
      if (nr == 0) cycle
      n_send = n_send + 1
      tag = 3001
      call hecmw_isend_int(bt_exp(idom)%index(0:bt_exp(idom)%nr), bt_exp(idom)%nr + 1, &
           irank, tag, heccomm%HECMW_COMM, requests(n_send))
      if (bt_exp(idom)%nnz == 0) cycle
      n_send = n_send + 1
      tag = 3002
      call hecmw_isend_int(exp_cols_item(1,exp_cols_index(idom-1)+1), &
           cncol_item * bt_exp(idom)%nnz, irank, tag, heccomm%HECMW_COMM, requests(n_send))
      n_send = n_send + 1
      tag = 3003
      call hecmw_isend_r(bt_exp(idom)%A, ndof2 * bt_exp(idom)%nnz, &
           irank, tag, heccomm%HECMW_COMM, requests(n_send))
    enddo
    !
    ! BT_imp%nr = hecCOMM%import_index(nnb)
    bt_imp%nr = hecmesh%n_node - hecmesh%nn_internal
    bt_imp%nc = 0  !!! TEMPORARY
    bt_imp%nnz = 0
    bt_imp%ndof = ndof
    !
    allocate(nnz_imp(nnb))
    allocate(cnt(bt_imp%nr))
    !
    cnt(:) = 0
    do idom = 1, nnb
      irank = heccomm%neighbor_pe(idom)
      idx_0 = heccomm%import_index(idom-1)
      idx_n = heccomm%import_index(idom)
      nr = idx_n - idx_0
      if (nr == 0) then
        nnz_imp(idom) = 0
        cycle
      endif
      allocate(index_imp(0:nr))
      tag = 3001
      call hecmw_recv_int(index_imp(0:nr), nr+1, irank, tag, &
           heccomm%HECMW_COMM, statuses(:,1))
      nnz_imp(idom) = index_imp(nr)
      do j = 1, nr
        jj = heccomm%import_item(idx_0 + j) - hecmesh%nn_internal
        if (jj < 1 .or. bt_imp%nr < jj) stop 'ERROR: jj out of range'
        if (cnt(jj) /= 0) stop import rows?'
        cnt(jj) = index_imp(j) - index_imp(j-1)
      enddo
      deallocate(index_imp)
    enddo
    !
    allocate(imp_cols_index(0:nnb))
    call make_index(nnb, nnz_imp, imp_cols_index)
    deallocate(nnz_imp)
    !
    allocate(BT_imp%index(0:BT_imp%nr))
    call make_index(BT_imp%nr, cnt, BT_imp%index)
    deallocate(cnt)
    !
    BT_imp%nnz = BT_imp%index(BT_imp%nr)
    if (BT_imp%nnz /= imp_cols_index(nnb)) &
         stop 'error: total num of nonzero of bt_imp'
    !
    allocate(imp_cols_item(cNCOL_ITEM, BT_imp%nnz))
    allocate(imp_vals_item(ndof2 * BT_imp%nnz))
    !
    do idom = 1, nnb
      irank = hecCOMM%neighbor_pe(idom)
      idx_0 = imp_cols_index(idom-1)
      idx_n = imp_cols_index(idom)
      nnz = idx_n - idx_0
      if (nnz == 0) cycle
      tag = 3002
      call HECMW_RECV_INT(imp_cols_item(1, idx_0 + 1), cNCOL_ITEM * nnz, &
           irank, tag, hecCOMM%HECMW_COMM, statuses(:,1))
      tag = 3003
      call HECMW_RECV_R(imp_vals_item(ndof2*idx_0 + 1), ndof2 * nnz, &
           irank, tag, hecCOMM%HECMW_COMM, statuses(:,1))
    enddo
    call HECMW_Waitall(n_send, requests, statuses)
    if (DEBUG >= 2) write(0,*) '  debug2: send bt_imp and recv into temporary data done'
    !
    deallocate(requests)
    deallocate(statuses)
    !
    do idom = 1, nnb
      call hecmw_localmat_free(BT_exp(idom))
    enddo
    deallocate(BT_exp)
    deallocate(exp_cols_index)
    deallocate(exp_cols_item)
    !
    call copy_mesh(hecMESH, hecMESHnew)
    !
    call map_imported_cols(hecMESHnew, imp_cols_index(nnb), imp_cols_item, n_add_node, add_nodes, map, i0)
    if (DEBUG >= 2) write(0,*) '  debug2: map imported cols done'
    !
    call update_comm_table(hecMESHnew, n_add_node, add_nodes, i0)
    if (DEBUG >= 2) write(0,*) '  debug2: update comm_table done'
    !
    BT_imp%nc = hecMESHnew%n_node
    !
    allocate(BT_imp%item(BT_imp%nnz))
    allocate(BT_imp%A(ndof2 * BT_imp%nnz))
    call copy_vals_to_BT_imp(hecCOMM, hecMESH%nn_internal, imp_cols_index, map, imp_vals_item, BT_imp)
    if (DEBUG >= 2) write(0,*) '  debug2: copy vals to bt_imp done'
    !
    deallocate(imp_cols_index)
    deallocate(imp_cols_item)
    deallocate(imp_vals_item)
    deallocate(map)
  end subroutine send_BT_exp_and_recv_BT_imp
 
  subroutine copy_vals_to_BT_imp(hecCOMM, nn_internal, imp_cols_index, map, imp_vals_item, BT_imp)
    implicit none
    type (hecmwST_matrix_comm), intent(in) :: hecCOMM
    integer(kind=kint), intent(in) :: nn_internal
    integer(kind=kint), allocatable, intent(in) :: imp_cols_index(:)
    integer(kind=kint), intent(in) :: map(:)
    real(kind=kreal), intent(in) :: imp_vals_item(:)
    type (hecmwST_local_matrix), intent(inout) :: BT_imp
    integer(kind=kint) :: nnb, ndof2, idx, idom, idx_0, idx_n, nr, j, jrow, ks, ke, k
    nnb = hecCOMM%n_neighbor_pe
    ndof2 = BT_imp%ndof ** 2
    idx = 0
    do idom = 1, nnb
      idx_0 = hecCOMM%import_index(idom-1)
      idx_n = hecCOMM%import_index(idom)
      nr = idx_n - idx_0
      if (nr == 0) cycle
      do j = 1, nr
        jrow = hecCOMM%import_item(idx_0 + j) - nn_internal
        ks = BT_imp%index(jrow-1)+1
        ke = BT_imp%index(jrow)
        do k = ks, ke
          idx = idx + 1
          BT_imp%item(k) = map(idx)
          BT_imp%A(ndof2*(k-1)+1:ndof2*k) = imp_vals_item(ndof2*(idx-1)+1:ndof2*idx)
        enddo
      enddo
      if (idx /= imp_cols_index(idom)) stop 'error: copy vals to bt_imp'
    enddo
  end subroutine copy_vals_to_BT_imp
 
  subroutine concat_BTmat_and_BT_imp(BTmat, BT_imp, BT_all)
    implicit none
    type (hecmwST_local_matrix), intent(in) :: BTmat
    type (hecmwST_local_matrix), intent(in) :: BT_imp
    type (hecmwST_local_matrix), intent(out) :: BT_all
    integer(kind=kint) :: ndof, ndof2, i, ii
    ndof = BTmat%ndof
.and.    if (BT_imp%nr > 0  BT_imp%ndof /= ndof) stop 'error: concat btmat and bt_imp: ndof'
    ndof2 = ndof*ndof
    BT_all%nr = BTmat%nr + BT_imp%nr
    BT_all%nc = max(BTmat%nc, BT_imp%nc)
    BT_all%nnz = BTmat%nnz + BT_imp%nnz
    BT_all%ndof = ndof
    allocate(BT_all%index(0:BT_all%nr))
    allocate(BT_all%item(BT_all%nnz))
    allocate(BT_all%A(ndof2 * BT_all%nnz))
    BT_all%index(0) = 0
    do i = 1, BTmat%nr
      BT_all%index(i) = BTmat%index(i)
    enddo
    do i = 1, BT_imp%nr
      BT_all%index(BTmat%nr+i) = BT_all%index(BTmat%nr+i-1) + &
           BT_imp%index(i) - BT_imp%index(i-1)
    enddo
    do i = 1, BTmat%nnz
      BT_all%item(i) = BTmat%item(i)
      BT_all%A(ndof2*(i-1)+1:ndof2*i) = BTmat%A(ndof2*(i-1)+1:ndof2*i)
    enddo
    do i = 1, BT_imp%nnz
      ii = BTmat%nnz + i
      BT_all%item(ii) = BT_imp%item(i)
      BT_all%A(ndof2*(ii-1)+1:ndof2*ii) = BT_imp%A(ndof2*(i-1)+1:ndof2*i)
    enddo
  end subroutine concat_BTmat_and_BT_imp
 
  subroutine multiply_mat_mat(Amat, Bmat, Cmat)
    implicit none
    type (hecmwST_local_matrix), intent(in) :: Amat
    type (hecmwST_local_matrix), intent(in) :: Bmat
    type (hecmwST_local_matrix), intent(out) :: Cmat
    integer(kind=kint) :: ndof, ndof2, nr, nc, nnz, i, icnt
    integer(kind=kint) :: js, je, j, jj, ks, ke, k, kk, l, ll, l0
    integer(kind=kint), allocatable :: iw(:)
    real(kind=kreal), pointer :: Ap(:), Bp(:), Cp(:)
    real(kind=kreal) :: t0, t1
    t0 = hecmw_wtime()
    if (Amat%ndof /= Bmat%ndof) stop 'error: multiply_mat_mat: unmatching ndof'
    ndof = Amat%ndof
    ndof2 = ndof*ndof
    nr = Amat%nr
    nc = Bmat%nc
    if (Amat%nc /= Bmat%nr) then
      write(0,*) 'amat: nr, nc = ', Amat%nr, Amat%nc
      write(0,*) 'bmat: nr, nc = ', Bmat%nr, Bmat%nc
      stop 'error: multiply_mat_mat: unmatching size'
    endif
    Cmat%ndof = ndof
    Cmat%nr = nr
    Cmat%nc = nc
    allocate(Cmat%index(0:nr))
    Cmat%index(0) = 0
    !$omp parallel default(none), &
      !$omp& private(iw,i,icnt,js,je,j,jj,ks,ke,k,kk,l), &
      !$omp& shared(nr,nc,Amat,Bmat,Cmat)
    allocate(iw(nc))
    !$omp do
    do i = 1, nr
      icnt = 0
      js = Amat%index(i-1)+1
      je = Amat%index(i)
      do j = js, je
        jj = Amat%item(j)
        ks = Bmat%index(jj-1)+1
        ke = Bmat%index(jj)
        kl1: do k = ks, ke
          kk = Bmat%item(k)
          do l = 1, icnt
            if (iw(l) == kk) cycle kl1
          enddo
          icnt = icnt + 1
          iw(icnt) = kk
        enddo kl1
      enddo
      Cmat%index(i) = icnt
    enddo
    !$omp end do
    deallocate(iw)
    !$omp end parallel
    do i = 1, nr
      Cmat%index(i) = Cmat%index(i-1) + Cmat%index(i)
    enddo
    nnz = Cmat%index(nr)
    Cmat%nnz = nnz
    !write(0,*) 'nnz',nnz
    t1 = hecmw_wtime()
    if (TIMER >= 3) write(0, '(a,f10.4)') "###### multiply_mat_mat (1) : ",t1-t0
    t0 = hecmw_wtime()
    allocate(Cmat%item(nnz))
    allocate(Cmat%A(ndof2 * nnz))
    Cmat%A(:) = 0.0d0
    !$omp parallel default(none), &
      !$omp& private(i,icnt,l0,js,je,j,jj,Ap,ks,ke,k,kk,Bp,ll,l,Cp), &
      !$omp& shared(nr,Cmat,Amat,Bmat,ndof2,ndof)
    !$omp do
    do i = 1, nr
      icnt = 0
      l0 = Cmat%index(i-1)
      ! item
      js = Amat%index(i-1)+1
      je = Amat%index(i)
      do j = js, je
        jj = Amat%item(j)
        Ap => Amat%A(ndof2*(j-1)+1:ndof2*j)
        ks = Bmat%index(jj-1)+1
        ke = Bmat%index(jj)
        do k = ks, ke
          kk = Bmat%item(k)
          Bp => Bmat%A(ndof2*(k-1)+1:ndof2*k)
          ll = -1
          do l = 1, icnt
            if (Cmat%item(l0+l) == kk) then
              ll = l0 + l
              exit
            endif
          enddo
          if (ll < 0) then
            icnt = icnt + 1
            ll = l0 + icnt
            Cmat%item(ll) = kk
          endif
          Cp => Cmat%A(ndof2*(ll-1)+1:ndof2*ll)
          call blk_matmul_add(ndof, Ap, Bp, Cp)
        enddo
      enddo
      !write(0,*) 'l0,icnt,index(i)',Cmat%index(i-1),icnt,Cmat%index(i)
      if (l0+icnt /= Cmat%index(i)) stop 'error: multiply_mat_mat: unknown error'
    enddo
    !$omp end do
    !$omp end parallel
    t1 = hecmw_wtime()
    if (TIMER >= 3) write(0, '(a,f10.4)') "###### multiply_mat_mat (2) : ",t1-t0
    t0 = hecmw_wtime()
    call sort_and_uniq_rows(Cmat)
    t1 = hecmw_wtime()
    if (TIMER >= 3) write(0, '(a,f10.4)') "###### multiply_mat_mat (3) : ",t1-t0
  end subroutine multiply_mat_mat
 
  subroutine blk_matmul_add(ndof, A, B, AB)
    implicit none
    integer, intent(in) :: ndof
    real(kind=kreal), intent(in) :: A(:), B(:)
    real(kind=kreal), intent(inout) :: AB(:)
    integer :: ndof2, i, j, k, i0, j0, ij, ik, jk
    ndof2=ndof*ndof
    do i=1,ndof
      i0=(i-1)*ndof
      do j=1,ndof
        ij=i0+j
        j0=(j-1)*ndof
        do k=1,ndof
          ik=i0+k
          jk=j0+k
          !$omp atomic
          AB(ik)=AB(ik)+A(ij)*B(jk)
        enddo
      enddo
    enddo
  end subroutine blk_matmul_add
 
  subroutine hecmw_localmat_make_hecmat(hecMAT, BTtKTmat, hecTKT)
    implicit none
    type (hecmwST_matrix), intent(in) :: hecMAT
    type (hecmwST_local_matrix), intent(in) :: BTtKTmat
    type (hecmwST_matrix), intent(inout) :: hecTKT
    call make_new_hecmat(hecMAT, BTtKTmat, hecTKT)
  end subroutine hecmw_localmat_make_hecmat
 
  !> \brief Debug write matrix
  !>
  subroutine debug_write_matrix(Mat, label, level)
    type(hecmwST_local_matrix), intent(in) :: Mat   !< matrix
    character(len=*),           intent(in) :: label !< label for matrix
    integer(kind=kint),         intent(in) :: level !< debug level
    !
    integer(kind=kint) :: iunit
 
    if (level <= 0) return
 
    iunit = 700 + hecmw_comm_get_rank()
    write(iunit,'(a,a)') trim(label),'============================================================'
    if (level == 1) then
      call hecmw_localmat_write_size(Mat, iunit)
    else if (level == 2) then
      call hecmw_localmat_write_ij(Mat, iunit)
    else
      call hecmw_localmat_write(Mat, iunit)
    endif
  end subroutine debug_write_matrix
 
end module hecmw_local_matrix