anchor.diff:      real(kind=kreal) :: Tcomm
   
       Tcomm = 0.d0
  -    call hecmw_matvec_22 (hecMESH, hecMAT, X, R, hecMAT%NP, Tcomm)
  +    call hecmw_matvec_22 (hecMESH, hecMAT, X, R, Tcomm)
       if (present(COMMtime)) COMMtime = COMMtime + Tcomm
  -
  +!$omp parallel default(none),private(i),shared(hecMAT,R,B)
  +!$omp do
       do i = 1, hecMAT%N * 2
         R(i) = B(i) - R(i)
       enddo
  -
  +!$omp end do
  +!$omp end parallel
     end subroutine hecmw_matresid_22
   
     !C
  
positive.code: 
    private
  
    public :: hecmw_matvec_22
    public :: hecmw_matresid_22
    public :: hecmw_rel_resid_L2_22
  
  contains
  
    !C
    !C***
    !C*** hecmw_matvec_22
    !C***
    !C
    subroutine hecmw_matvec_22 (hecMESH, hecMAT, X, Y, n, COMMtime)
      use hecmw_util
      use m_hecmw_comm_f
  
      implicit none
      integer(kind=kint):: n
      real(kind=kreal), dimension(2*n) :: X, Y
      type (hecmwST_matrix)     :: hecMAT
      type (hecmwST_local_mesh) :: hecMESH
      real(kind=kreal), optional :: COMMtime
  
      real(kind=kreal) :: START_TIME, END_TIME
      integer(kind=kint) :: i, j, jS, jE, in
      real(kind=kreal) :: YV1, YV2, X1, X2
  
      START_TIME= HECMW_WTIME()
      call hecmw_update_2_R (hecMESH, X, hecMAT%NP)
      END_TIME= HECMW_WTIME()
      if (present(COMMtime)) COMMtime = COMMtime + END_TIME - START_TIME
  
      do i= 1, hecMAT%N
        X1= X(2*i-1)
        X2= X(2*i  )
        YV1= hecMAT%D(4*i-3)*X1 + hecMAT%D(4*i-2)*X2
        YV2= hecMAT%D(4*i-1)*X1 + hecMAT%D(4*i  )*X2
  
        jS= hecMAT%indexL(i-1) + 1
        jE= hecMAT%indexL(i  )
        do j= jS, jE
          in  = hecMAT%itemL(j)
          X1 = X(2*in-1)
          X2 = X(2*in  )
          YV1= YV1 + hecMAT%AL(4*j-3)*X1 + hecMAT%AL(4*j-2)*X2
          YV2= YV2 + hecMAT%AL(4*j-1)*X1 + hecMAT%AL(4*j  )*X2
        enddo
        jS= hecMAT%indexU(i-1) + 1
        jE= hecMAT%indexU(i  )
        do j= jS, jE
          in  = hecMAT%itemU(j)
          X1 = X(2*in-1)
          X2 = X(2*in  )
          YV1= YV1 + hecMAT%AU(4*j-3)*X1 + hecMAT%AU(4*j-2)*X2
          YV2= YV2 + hecMAT%AU(4*j-1)*X1 + hecMAT%AU(4*j  )*X2
        enddo
        Y(2*i-1)= YV1
        Y(2*i  )= YV2
      enddo
  
    end subroutine hecmw_matvec_22
  
    !C
    !C***
    !C*** hecmw_matresid_22
  
negative.code:     hecmw_rel_resid_L2_22 = sqrt(rnorm2 / bnorm2)
  
    end function hecmw_rel_resid_L2_22
    !C
    !C***
    !C*** hecmw_Tvec_22
    !C***
    !C
    subroutine hecmw_Tvec_22 (hecMESH, X, Y, COMMtime)
      use hecmw_util
      use m_hecmw_comm_f
      implicit none
      type (hecmwST_local_mesh), intent(in) :: hecMESH
      real(kind=kreal), intent(in) :: X(:)
      real(kind=kreal), intent(out) :: Y(:)
      real(kind=kreal), intent(inout) :: COMMtime
  
      real(kind=kreal) :: START_TIME, END_TIME
      integer(kind=kint) :: i, j, jj, k, kk
  
      START_TIME= HECMW_WTIME()
      call hecmw_update_3_R (hecMESH, X, hecMESH%n_node)
      END_TIME= HECMW_WTIME()
      COMMtime = COMMtime + END_TIME - START_TIME
  
  !$omp parallel default(none),private(i,k,kk,j,jj),shared(hecMESH,X,Y)
  !$omp do
      do i= 1, hecMESH%nn_internal * hecMESH%n_dof
        Y(i)= X(i)
      enddo
  !$omp end do
  
  !$omp do
      do i= 1, hecMESH%mpc%n_mpc
        k = hecMESH%mpc%mpc_index(i-1) + 1
        kk = 2 * (hecMESH%mpc%mpc_item(k) - 1) + hecMESH%mpc%mpc_dof(k)
        Y(kk) = 0.d0
        do j= hecMESH%mpc%mpc_index(i-1) + 2, hecMESH%mpc%mpc_index(i)
          jj = 2 * (hecMESH%mpc%mpc_item(j) - 1) + hecMESH%mpc%mpc_dof(j)
          Y(kk) = Y(kk) - hecMESH%mpc%mpc_val(j) * X(jj)
        enddo
      enddo
  !$omp end do
  !$omp end parallel
  
    end subroutine hecmw_Tvec_22
  
    !C
    !C***
    !C*** hecmw_Ttvec_22
    !C***
    !C
    subroutine hecmw_Ttvec_22 (hecMESH, X, Y, COMMtime)
      use hecmw_util
      use m_hecmw_comm_f
      implicit none
      type (hecmwST_local_mesh), intent(in) :: hecMESH
      real(kind=kreal), intent(in) :: X(:)
      real(kind=kreal), intent(out) :: Y(:)
      real(kind=kreal), intent(inout) :: COMMtime
  
      real(kind=kreal) :: START_TIME, END_TIME
      integer(kind=kint) :: i, j, jj, k, kk
  
      START_TIME= HECMW_WTIME()
      call hecmw_update_2_R (hecMESH, X, hecMESH%n_node)
      END_TIME= HECMW_WTIME()
      COMMtime = COMMtime + END_TIME - START_TIME
  
  !$omp parallel default(none),private(i,k,kk,j,jj),shared(hecMESH,X,Y)
  !$omp do
      do i= 1, hecMESH%nn_internal * hecMESH%n_dof
        Y(i)= X(i)
      enddo
  !$omp end do
  
  !$omp do
      do i= 1, hecMESH%mpc%n_mpc
        k = hecMESH%mpc%mpc_index(i-1) + 1
        kk = 2 * (hecMESH%mpc%mpc_item(k) - 1) + hecMESH%mpc%mpc_dof(k)
        Y(kk) = 0.d0
        do j= hecMESH%mpc%mpc_index(i-1) + 2, hecMESH%mpc%mpc_index(i)
          jj = 2 * (hecMESH%mpc%mpc_item(j) - 1) + hecMESH%mpc%mpc_dof(j)
          Y(jj) = Y(jj) - hecMESH%mpc%mpc_val(j) * X(kk)
        enddo
      enddo
  !$omp end do
  !$omp end parallel
  
    end subroutine hecmw_Ttvec_22
  
    !C
    !C***
    !C*** hecmw_TtmatTvec_22
    !C***
    !C
    subroutine hecmw_TtmatTvec_22 (hecMESH, hecMAT, X, Y, W, COMMtime)
      use hecmw_util
      implicit none
      type (hecmwST_local_mesh), intent(in) :: hecMESH
      type (hecmwST_matrix), intent(in)     :: hecMAT
      real(kind=kreal), intent(in) :: X(:)
      real(kind=kreal), intent(out) :: Y(:), W(:)
      real(kind=kreal), intent(inout) :: COMMtime
  
      call hecmw_Tvec_22(hecMESH, X, Y, COMMtime)
      call hecmw_matvec_22_inner(hecMESH, hecMAT, Y, W, COMMtime)
      call hecmw_Ttvec_22(hecMESH, W, Y, COMMtime)
  
    end subroutine hecmw_TtmatTvec_22
  
  
    !C
    !C***
    !C*** hecmw_trans_b_22
    !C***
    !C
    subroutine hecmw_trans_b_22(hecMESH, hecMAT, B, BT, COMMtime)
      use hecmw_util
      implicit none
      type (hecmwST_local_mesh), intent(in) :: hecMESH
      type (hecmwST_matrix), intent(in)     :: hecMAT
      real(kind=kreal), intent(in) :: B(:)
      real(kind=kreal), intent(out), target :: BT(:)
      real(kind=kreal), intent(inout) :: COMMtime
  
      real(kind=kreal), allocatable :: W(:)
      real(kind=kreal), pointer :: XG(:)
      integer(kind=kint) :: i, k, kk
  
      allocate(W(hecMESH%n_node * 3))
  
      !C===
      !C +---------------------------+
      !C | {bt}= [T']({b} - [A]{xg}) |
      !C +---------------------------+
      !C===
      XG => BT
      XG = 0.d0
  
  !$omp parallel default(none),private(i,k,kk),shared(hecMESH,XG)
  !$omp do
      !C-- Generate {xg} from mpc_const
      do i = 1, hecMESH%mpc%n_mpc
        k = hecMESH%mpc%mpc_index(i-1) + 1
        kk = 2 * hecMESH%mpc%mpc_item(k) + hecMESH%mpc%mpc_dof(k) - 2
        XG(kk) = hecMESH%mpc%mpc_const(i)
      enddo
  !$omp end do
  !$omp end parallel
  
      !C-- {w} = {b} - [A]{xg}
      call hecmw_matresid_22 (hecMESH, hecMAT, XG, B, W, COMMtime)
  
      !C-- {bt} = [T'] {w}
      call hecmw_Ttvec_22(hecMESH, W, BT, COMMtime)
  
      deallocate(W)
    end subroutine hecmw_trans_b_22
  
    !C
    !C***
    !C*** hecmw_tback_x_22
    !C***
    !C
    subroutine hecmw_tback_x_22(hecMESH, X, COMMtime)
      use hecmw_util
      implicit none
      type (hecmwST_local_mesh), intent(in) :: hecMESH
      real(kind=kreal), intent(inout) :: X(:)
      real(kind=kreal) :: COMMtime
  
      real(kind=kreal), allocatable :: W(:)
      integer(kind=kint) :: i, k, kk
  
      allocate(W(hecMESH%n_node * 2))
  
      !C-- {tx} = [T]{x}
      call hecmw_Tvec_22(hecMESH, X, W, COMMtime)
  
      !C-- {x} = {tx} + {xg}
  
  !$omp parallel default(none),private(i,k,kk),shared(hecMESH,X,W)
  !$omp do
      do i= 1, hecMESH%nn_internal * 2
        X(i)= W(i)
      enddo
  !$omp end do
  
  !$omp do
      do i = 1, hecMESH%mpc%n_mpc
        k = hecMESH%mpc%mpc_index(i-1) + 1
        kk = 2 * hecMESH%mpc%mpc_item(k) + hecMESH%mpc%mpc_dof(k) - 2
        X(kk) = X(kk) + hecMESH%mpc%mpc_const(i)
      enddo
  !$omp end do
  !$omp end parallel
  
      deallocate(W)
    end subroutine hecmw_tback_x_22
  
    !C
    !C***
    !C*** hecmw_mat_diag_sr_22
    !C***
    !C
    subroutine hecmw_mat_diag_sr_22(hecMESH, hecMAT, COMMtime)
      use hecmw_util
      use m_hecmw_comm_f
      implicit none
      type (hecmwST_local_mesh), intent(in) :: hecMESH
      type (hecmwST_matrix), intent(inout), target :: hecMAT
      real(kind=kreal), intent(inout), optional :: COMMtime
      real(kind=kreal), allocatable :: W(:,:)
      real(kind=kreal), pointer :: D(:)
      integer(kind=kint) :: ip
      real(kind=kreal) :: START_TIME, END_TIME, Tcomm
      allocate(W(2*hecMAT%NP,2))
      D => hecMAT%D
      do ip= 1, hecMAT%N
        W(2*ip-1,1)= D(4*ip-3); W(2*ip-1,2)= D(4*ip-2);
        W(2*ip-0,1)= D(4*ip-1); W(2*ip-0,2)= D(4*ip-0);
      enddo
      START_TIME= HECMW_WTIME()
      call hecmw_update_2_R (hecMESH, W(:,1), hecMAT%NP)
      call hecmw_update_2_R (hecMESH, W(:,2), hecMAT%NP)
      END_TIME= HECMW_WTIME()
      if (present(COMMtime)) COMMtime = COMMtime + END_TIME - START_TIME
      do ip= hecMAT%N+1, hecMAT%NP
        D(4*ip-3)= W(2*ip-1,1); D(4*ip-2)= W(2*ip-1,2);
        D(4*ip-1)= W(2*ip-0,1); D(4*ip-0)= W(2*ip-0,2);
      enddo
      deallocate(W)
    end subroutine hecmw_mat_diag_sr_22
  end module hecmw_solver_las_22
  