module FEMSolverClass
   !Linear soler for FEMDomainClass
   use iso_fortran_env
   use SparseClass
   use FEMDomainClass
   implicit none

   type :: M_Link_Item_
      integer(int32) :: rank_and_rowID_1(1:2) = [0, 0]
      integer(int32) :: rank_and_rowID_2(1:2) = [0, 0]
   end type

   type :: FEMSolver_
      type(FEMDomainp_), allocatable :: femdomains(:)
      integer(int32), allocatable :: DomainIDs(:)

      real(real64), allocatable :: IfaceElemConnectivity(:, :)
      real(real64), allocatable :: IfaceElemDomainID(:, :)

      logical :: initialized = .false.
      logical :: InterfaceExist = .false.
      logical :: debug = .false.

      real(real64), allocatable :: CRS_val(:)
      integer(int32), allocatable :: CRS_Index_Col(:)
      integer(int64), allocatable :: CRS_Index_Row(:)
      real(real64), allocatable :: CRS_RHS(:)

      integer(int32), allocatable :: number_of_element_per_domain(:)
      integer(int32), allocatable :: number_of_point_per_domain(:)

      !> General Eigen Value Problem
      !> [A]{x} = (lambda)[B]{x}

      real(real64), allocatable   :: A_CRS_val(:)
      integer(int32), allocatable :: A_CRS_Index_Col(:)
      integer(int64), allocatable :: A_CRS_Index_Row(:)
      logical                    :: A_empty = .true.

      real(real64), allocatable   :: B_CRS_val(:)
      integer(int32), allocatable :: B_CRS_Index_Col(:)
      integer(int64), allocatable :: B_CRS_Index_Row(:)
      logical                    :: B_empty = .true.

      integer(int32), allocatable :: fix_eig_IDs(:)

      logical, allocatable :: fix_lin_exists(:)
      real(real64), allocatable :: fix_lin_exists_Values(:)

      !integer(int32),allocatable :: fix_lin_IDs(:)
      !real(real64),allocatable :: fix_lin_Values(:)

      real(real64), allocatable :: CRS_x(:)
      integer(int32), allocatable :: CRS_ID_Starts_From(:)

      ! dense matrix
      real(real64), allocatable :: A_dense(:, :)
      integer(int32), allocatable :: Num_nodes_in_Domains(:)

      ! with MPI
      type(MPI_), pointer :: MPI_target => null()
      type(M_Link_Item_), allocatable :: Link_Table(:)

      integer(int32) :: LINK_TABLE_INIT_SIZE = 1000
      integer(int32) :: Link_num = 0

      type(CRS_) :: ILU_MATRIX

      integer(int32) :: itrmax = 100000
      real(real64)   :: er0 = dble(1.0e-10)
      real(real64)   :: relative_er = dble(1.0e-10)

      logical        :: use_LOBPCG = .false.
      integer(int32) :: LOBPCG_MAX_ITR = 100000
      real(real64)   :: LOBPCG_TOL = dble(1.0e-6)
      integer(int32) :: LOBPCG_NUM_MODE = 5

        !!! stack for CRS matrices
      type(CRS_), allocatable :: CRS_matrix_stack(:)
      type(Dictionary_)      :: CRS_matrix_stack_key
      integer(int32) :: LAST_CRS_MATRIX_STACK = 0
      integer(int32) :: MAXSIZE_CRS_MATRIX_STACK = 10

   contains
      !(1) Initialize solver
      procedure, public ::  init => initFEMSolver

      !(2) set Domain info
      procedure, public ::  setDomain => setDomainFEMSolver
      procedure, public ::  setDomains => setDomainFEMSolver

      !(3) setup Ax=b as CRS format
      procedure, pass :: setCRS_CRSobjFEMSolver
      procedure, pass :: setCRSFEMSolver
      procedure, public :: setRHS => setRHSFEMSolver
      generic ::  setCRS => setCRSFEMSolver, setCRS_CRSobjFEMSolver
      procedure, public ::  getCRS => getCRSFEMSolver

      !(4) set Ax=b as CRS format
      procedure, public ::  setValue => setValueFEMSolver
      procedure, public ::  setValues => setValueFEMSolver
      procedure, public :: addMatrixValue => addMatrixValueFEMSolver

      procedure, public ::  setMatrix => setMatrixFEMSolver
      procedure, public ::  setVector => setVectorFEMSolver

      procedure, public ::  keepThisMatrixAs => keepThisMatrixAsFEMSolver

      !(5) fix x=\bar{x} (not implemented yet.)
      ! for Linear Solver
      procedure, public :: fix => fixFEMSolver

      !(5) fix x=\bar{x} (not implemented yet.)
      ! for eigen solver
      procedure, public :: fix_eig => fix_eigFEMSolver

      procedure, pass :: importCRSMatrix_FEMSolver
      generic, public :: import => importCRSMatrix_FEMSolver

      !(6) save matrix
      procedure, public :: saveMatrix => saveMatrixFEMSolver
      procedure, public :: loadMatrix => loadMatrixFEMSolver

      !(7-1) Modal analysis
      procedure, public :: eig => eigFEMSolver

      !(7-2) linear solver
      procedure, pass :: solveFEMSolver
      procedure, pass :: solveFEMSolver_UserDefinedLinearSolver
      procedure, pass :: solveFEMSolver_UserDefinedLinearSolverAsFunc
      generic :: solve => solveFEMSolver, solveFEMSolver_UserDefinedLinearSolver, &
         solveFEMSolver_UserDefinedLinearSolverAsFunc

      !(7-3) condition number
      procedure, public :: conditionNumber => conditionNumberFEMSolver

      !re-zero matrix
      procedure, public :: zeros => zerosFEMSolver
      ! M:diag matrix,  A*M^{-1}
      procedure, public :: matmulDiagMatrix => matmulDiagMatrixFEMSolver
      procedure, public :: diag => diagFEMSolver

      ! others:
      ! Energy-based Overset Mesh (Tomobe et al., under review)
      procedure, public :: setEbOM => setEbOMFEMSolver
      procedure, public :: get_gap_function => get_gap_function_FEMSolver
      procedure, public :: argumented_Lagrangian_RHS => argumented_Lagrangian_RHS_FEMSolver
      procedure, public :: num_EbOFEM_projection_point &
         => num_EbOFEM_projection_point_FEMSolver

      procedure, public :: get_fix_idx => get_fix_idx_FEMSolver
      procedure, public :: get_fix_value => get_fix_value_FEMSolver

      ! FOR MPI
      procedure, public :: MPI_link => MPI_linkFEMSolver
      procedure, public :: MPI_dot_product => MPI_dot_productFEMSolver
      procedure, public :: MPI_matmul => MPI_matmulFEMSolver
      procedure, public :: MPI_BICGSTAB => MPI_BICGSTABFEMSolver

      ! destractor
      procedure, public :: remove => removeFEMSolver

   end type

   interface reverseArray
      module procedure reverseArrayReal64
   end interface reverseArray
   interface reverseVector
      module procedure reverseVectorReal64
   end interface reverseVector
contains

   recursive subroutine initFEMSolver(this, NumDomain, FEMDomains, DomainIDs, DOF, MPI_target, &
                                      NUM_CRS_MATRIX_STACK)
      class(FEMSolver_), intent(inout) :: this
      ! two-way
      integer(int32), optional, intent(in) :: numDomain
      ! optional
      ! you can bypass solver%setDomain and solver%setCRS
      type(FEMDomain_), optional, intent(in) :: FEMDomains(:)
      integer(int32), optional, intent(in) :: DomainIDs(:), DOF, NUM_CRS_MATRIX_STACK

      ! useless
      type(MPI_), optional, target, intent(in) :: MPI_target
      integer(int32) :: i, np
      integeR(int32), allocatable :: default_DomainIDs(:)

      if (present(FEMDomains)) then
         this%number_of_element_per_domain = int(zeros(size(FEMDomains)))
         do i = 1, size(FEMDomains)
            this%number_of_element_per_domain(i) = FEMDomains(i)%ne()
            this%number_of_point_per_domain(i) = FEMDomains(i)%nn()
         end do
      end if

      if (allocated(this%CRS_matrix_stack)) deallocate (this%CRS_matrix_stack)

      np = input(default=10, option=NUM_CRS_MATRIX_STACK)
      allocate (this%CRS_matrix_stack(np))
      call this%CRS_matrix_stack_key%init(np)
      this%LAST_CRS_MATRIX_STACK = 0

      if (present(FEMDomains)) then
         if (present(DomainIDs) .and. present(DOF)) then
            ! bypass mode
            call this%init(NumDomain=size(FEMDomains))
            call this%setDomain(FEMDomains=FEMDomains(:), DomainIDs=DomainIDs)
            call this%setCRS(DOF=DOF)
            return
         elseif (present(DOF)) then
            default_DomainIDs = zeros(size(FEMDomains))
            do i = 1, size(FEMDomains)
               default_DomainIDs(i) = i
            end do

            call this%init(NumDomain=size(FEMDomains))
            call this%setDomain(FEMDomains=FEMDomains(:), DomainIDs=default_DomainIDs)
            call this%setCRS(DOF=DOF)
            return
         else
            print *, "ERROR :: initFEMSolver >> "
            print *, "You are trying to use ByPass-mode,"
            print *, "which requires at least following two arguments"
            print *, "(1) type(FEMDomain_) :: FEMDomains(:)  "
            print *, "(2) Integer(int32)   :: DOF <DEGREE OF FREEDOM> "
            print *, " "
            print *, "and, if you give original domain-ids,"
            print *, "(2) Integer(int32)   :: DomainIDs(:) "
            print *, "is also necessary."
            stop
         end if
      end if

      if (.not. present(numDomain)) then
         print *, "ERROR :: initFEMSolver >> "
         print *, "Please input "
         print *, "integer(int32) :: numDomain "
         print *, "or  "
         print *, "type(FEMDomains_) :: femdomains(:) "
         print *, "integer(int32) :: DOF <DEGREE OF FREEDOM> "
         stop
      end if

      nullify (this%MPI_target)

      if (present(MPI_target)) then
         this%MPI_target => MPI_target
      end if

      this%initialized = .false.
      if (allocated(this%femdomains)) deallocate (this%femdomains)
      if (allocated(this%IfaceElemConnectivity)) then
         deallocate (this%IfaceElemConnectivity)
      end if! IfaceElemConnectivity(:,:)
      if (allocated(this%IfaceElemDomainID)) then
         deallocate (this%IfaceElemDomainID)
      end if! IfaceElemDomainID(:,:)

      this%initialized = .false.
      this%InterfaceExist = .false.

      if (allocated(this%CRS_val)) then
         deallocate (this%CRS_val)
      end if! CRS_val(:)
      if (allocated(this%CRS_Index_Col)) deallocate (this%CRS_Index_Col)!(:)
      if (allocated(this%CRS_Index_Row)) deallocate (this%CRS_Index_Row)!(:)
      if (allocated(this%CRS_RHS)) then
         deallocate (this%CRS_RHS)
      end if! CRS_RHS(:)

      !> General Eigen Value Problem
      !> [A]{x} = (lambda)[B]{x}

      if (allocated(this%A_CRS_val)) deallocate (this%A_CRS_val)!(:)
      if (allocated(this%A_CRS_Index_Col)) deallocate (this%A_CRS_Index_Col)!(:)
      if (allocated(this%A_CRS_Index_Row)) deallocate (this%A_CRS_Index_Row)!(:)
      this%A_empty = .true.

      if (allocated(this%B_CRS_val)) deallocate (this%B_CRS_val)!(:)
      if (allocated(this%B_CRS_Index_Col)) deallocate (this%B_CRS_Index_Col)!(:)
      if (allocated(this%B_CRS_Index_Row)) deallocate (this%B_CRS_Index_Row)!(:)
      this%B_empty = .true.

      if (allocated(this%fix_eig_IDs)) deallocate (this%fix_eig_IDs)!(:)

      if (allocated(this%CRS_x)) then
         deallocate (this%CRS_x)
      end if! CRS_x(:)
      if (allocated(this%CRS_ID_Starts_From)) then
         deallocate (this%CRS_ID_Starts_From)
      end if! CRS_ID_Starts_From(:)

      ! dense matrix
      if (allocated(this%A_dense)) then
         deallocate (this%A_dense)
      end if! A_dense(:,:)
      if (allocated(this%Num_nodes_in_Domains)) deallocate (this%Num_nodes_in_Domains)! Num_nodes_in_Domains(:)

      if (numDomain <= 0) then
         print *, "ERROR :: Number of element should be greater than 1"
         return
      else
         if (allocated(this%femdomains)) then
            deallocate (this%femdomains)
         end if
         allocate (this%femdomains(numDomain))
         allocate (this%DomainIDs(numDomain))
         do i = 1, numDomain
            this%DomainIDs(i) = i
         end do
      end if

!    if(NumInterfaceElement==0)then
!        this%InterfaceExist=.false.
!    else
!        ! check if NumNodePerInterfaceElement exists
!        if(.not. present(NumNodePerInterfaceElement) )then
!            print *, "ERROR :: NumNodePerInterfaceElement should be present."
!            return
!        endif
!        this%InterfaceExist=.true.
!    endif

      this%initialized = .true.

   end subroutine

   recursive subroutine setDomainFEMSolver(this, FEMDomain, FEMDomains, FEMDomainPointers, DomainID, DomainIDs)
      class(FEMSolver_), intent(inout) :: this
      type(FEMDomain_), target, optional, intent(in) :: FEMDomain, FEMDomains(:)
      type(FEMDomainp_), target, optional, intent(in) :: FEMDomainPointers(:)
      integer(int32), optional, intent(in) :: DomainID, DomainIDs(:)
      integer(int32) :: i

      if (present(FEMDomain)) then
         this%number_of_point_per_domain = [FEMDomain%nn()]
      elseif (present(FEMDomains)) then
         this%number_of_point_per_domain = int(zeros(size(FEMDomains)))
         do i = 1, size(FEMDomains)
            this%number_of_point_per_domain(i) = FEMDomains(i)%nn()
         end do
      end if

      if (present(DomainID)) then
         if (DomainID <= 0) then
            print *, "ERROR :: FEMSOlver%setDomain >> DomainID should be >=1"
         end if
      end if

      if (present(DomainIDs)) then
         if (minval(DomainIDs) <= 0) then
            print *, "ERROR :: FEMSOlver%setDomain >> DomainID should be >=1"
         end if
      end if

!    if(.not. present(DomainID) .and. .not.present(DomainIDs) )then
!        print *, "ERROR :: DomainID or DomainIDs are to be passed."
!        stop
!    endif
!
!    if(.not. present(FEMDomain) .and. .not.present(FEMDomains) )then
!        print *, "ERROR :: FEMDomain or FEMDomains are to be passed."
!        stop
!    endif

      if (.not. this%initialized) then
         print *, "ERROR :: this%setDomain should be called after this%init()"
      end if

      if (present(FEMDomain) .and. present(DomainID)) then
         if (associated(this%FEMDomains(DomainID)%FEMDomainp)) then
            nullify (this%FEMDomains(DomainID)%FEMDomainp)
         end if
         this%FEMDomains(DomainID)%FEMDomainp => FEMDomain

         return
      end if

      if (present(FEMDomains) .and. present(DomainIDs)) then

         if (maxval(DomainIDs) > size(this%femdomains)) then
            print *, "ERROR :: setFEMDomains >> invalid domain IDs"
            stop
         end if

         do i = 1, size(DomainIDs)
            call this%setDomain( &
               FEMDomain=FEMDomains(i), DomainID=DomainIDs(i))
         end do
         return

      end if

      if (present(FEMDomains) .and. .not. present(DomainIDs)) then
         ! see domain_ptr as id
         do i = 1, size(this%DomainIDs)
            call this%setDomain( &
               FEMDomain=FEMDomains(i), DomainID=this%DomainIDs(i))
         end do
         return
      end if

      if (present(FEMDomainPointers) .and. .not. present(DomainIDs)) then
         ! see domain_ptr as id
         this%DomainIDs = int(zeros(size(FEMDomainPointers)))
         do i = 1, size(FEMDomainPointers)
            this%DomainIDs(i) = i
         end do

         do i = 1, size(FEMDomainPointers)
            call this%setDomain( &
               FEMDomain=FEMDomainPointers(i)%femdomainp, DomainID=i)
         end do
         return
      end if
      print *, "ERROR >> setDomainFEMSolver >> invalid combinations for args"

   end subroutine

   subroutine setCRS_CRSobjFEMSolver(this, CRS)
      class(FEMSolver_), intent(inout) :: this
      type(CRS_), intent(in) :: CRS

      this%CRS_Index_Col = CRS%col_idx
      this%CRS_Index_Row = CRS%row_ptr
      this%CRS_Val = CRS%val

   end subroutine

   subroutine setRHSFEMSolver(this, RHS)
      class(FEMSolver_), intent(inout) :: this
      real(real64), intent(in) :: RHS(:)

      this%CRS_RHS = RHS

   end subroutine

   recursive subroutine setCRSFEMSolver(this, DOF, debug)
      class(FEMSolver_), intent(inout) :: this

      integer(int32), intent(in) :: DOF
      logical, optional, intent(in) :: debug
      integer(int64) :: i, j, k, l, m, n, itr
      integer(int64) :: size_of_global_matrix, offset, node_id, row_id
      integer(int64) :: node_id_p, col_id, kk, ll, drow_offset, buf_1
      integer(int32), allocatable :: Num_nodes_in_Domains(:), num_entry_in_row(:)
      integer(int32), allocatable :: col_local(:), new_col_local(:)
      integer(int32), allocatable :: Overset_CRS_Index_Col(:)
      integer(int64), allocatable :: Overset_CRS_Index_Row(:)
      integer(int32), allocatable :: new_CRS_Index_Col(:)
      integer(int64), allocatable :: new_CRS_Index_Row(:)
      integer(int64) :: row_node_id, col_node_id, row, col, row_domain_id, col_domain_id, &
                        row_DOF, col_DOF, k0
      integer(int32), allocatable :: num_nodes_of_domains(:)
      logical :: debug_mode_on = .true.
      integer(int32) :: DomainID
      type(FEMSolver_) :: single_domain_solver
      type(COO_)       :: COO
      type(CRS_)       :: ZeroMatrix

      if (present(debug)) then
         debug_mode_on = debug
      end if

      if (allocated(this%femdomains) .and. size(this%femdomains) >= 2) then
         ! check interface
         itr = 0

         if (debug_mode_on) then
            print *, "[ok] setCRS started."
         end if

         do domainID = 1, size(this%femdomains)
            if (this%femdomains(domainID)%femdomainp%empty()) then
               cycle
            end if

            call single_domain_solver%init(NumDomain=1)
            call single_domain_solver%setDomain(FEMDomain=this%femdomains(domainID)%femdomainp, &
                                                DomainID=1)
            call single_domain_solver%setCRS(DOF=DOF)

            itr = itr + 1

            if (itr == 1) then
               this%CRS_val = single_domain_solver%CRS_val
               this%CRS_Index_Col = single_domain_solver%CRS_Index_Col
               this%CRS_Index_Row = single_domain_solver%CRS_Index_Row
               this%CRS_RHS = single_domain_solver%CRS_RHS
            else
               n = size(this%CRS_Index_Row) - 1
               single_domain_solver%CRS_Index_Col(:) = single_domain_solver%CRS_Index_Col(:) + n
               this%CRS_Index_Col = this%CRS_Index_Col//single_domain_solver%CRS_Index_Col

               m = this%CRS_Index_Row(size(this%CRS_Index_Row)) - 1
               single_domain_solver%CRS_Index_Row(:) = single_domain_solver%CRS_Index_Row(:) + m

               this%CRS_Index_Row = this%CRS_Index_Row(1:n)//single_domain_solver%CRS_Index_Row

               this%CRS_RHS = this%CRS_RHS//single_domain_solver%CRS_RHS
               this%CRS_val = this%CRS_val//single_domain_solver%CRS_val
            end if
            call single_domain_solver%remove()

         end do

         if (debug_mode_on) then
            print *, "[ok] setCRS >> CRS-allocation done."
         end if
         ! <debug> following values are not unnatural
         !print *, size(this%CRS_Index_Col), maxval(this%CRS_Index_Col),minval(this%CRS_Index_Col)
         !print *, size(this%CRS_Index_Row), maxval(this%CRS_Index_Row),minval(this%CRS_Index_Row)
         !print *, size(this%CRS_RHS), maxval(this%CRS_RHS),minval(this%CRS_RHS)
         !print *, size(this%CRS_Val), maxval(this%CRS_Val),minval(this%CRS_Val)
         !print *, size(this%CRS_Index_Col)/size(this%CRS_Index_Row)

         ! >>>> time consuming part >>>>
         ! >>>> time consuming part >>>>
         ! >>>> time consuming part >>>>
         ! >>>> time consuming part >>>>

         call COO%init(num_row=size(this%CRS_Index_Row) - 1)
         ! next :: interface:
         do domainID = 1, size(this%femdomains)
            if (this%femdomains(domainID)%femdomainp%empty()) then
               cycle
            end if

            if (allocated(this%femdomains(DomainID)%femdomainp%OversetConnect)) then
               ! (a) create interface connectivity

               if (allocated(num_nodes_of_domains)) deallocate (num_nodes_of_domains)
               allocate (num_nodes_of_domains(1:size(this%femdomains)))
               num_nodes_of_domains(1:size(this%femdomains)) = 0
               do i = 1, size(this%femdomains)
                  if (associated(this%femdomains(i)%femdomainp)) then
                     if (.not. this%femdomains(i)%femdomainp%empty()) then
                        num_nodes_of_domains(i) = this%femdomains(i)%femdomainp%nn()
                     end if
                  end if
               end do

               ! Overset_CRS_Index_Col
               ! Overset_CRS_Index_Row
               do i = 1, size(this%femdomains(DomainID)%femdomainp%OversetConnect)
                  if (.not. allocated(this%femdomains(DomainID)%femdomainp%OversetConnect(i)%InterConnect)) then
                     cycle
                  end if
                  do j = 1, size(this%femdomains(DomainID)%femdomainp%OversetConnect(i)%InterConnect)
                     do k = 1, size(this%femdomains(DomainID)%femdomainp%OversetConnect(i)%InterConnect)
                        row_node_id = this%femdomains(DomainID)%femdomainp%OversetConnect(i)%InterConnect(j)
                        col_node_id = this%femdomains(DomainID)%femdomainp%OversetConnect(i)%InterConnect(k)
                        row_domain_id = this%femdomains(DomainID)%femdomainp%OversetConnect(i)%DomainIDs12(j)
                        col_domain_id = this%femdomains(DomainID)%femdomainp%OversetConnect(i)%DomainIDs12(k)
                        do row_DOF = 1, DOF
                           do col_DOF = 1, DOF
                              if (row_domain_id == 1) then
                                 row = (row_node_id - 1)*DOF + row_DOF
                              else
                                 row = sum(num_nodes_of_domains(1:row_domain_id - 1))*DOF + &
                                       (row_node_id - 1)*DOF + row_DOF
                              end if

                              if (col_domain_id == 1) then
                                 col = (col_node_id - 1)*DOF + col_DOF
                              else
                                 col = sum(num_nodes_of_domains(1:col_domain_id - 1))*DOF + &
                                       (col_node_id - 1)*DOF + col_DOF
                              end if

                              call COO%add(row=int(row), col=int(col), val=0.0d0)
                           end do
                        end do
                     end do
                  end do
               end do

               ! (b) and merge them to
               ! this%CRS_Index_Col
               ! this%CRS_Index_Row
               do i = 1, size(this%CRS_Index_Row) - 1
                  do j = this%CRS_Index_row(i), this%CRS_Index_row(i + 1) - 1
                     row = i
                     col = this%CRS_Index_Col(j)
                     call COO%add(row=int(row), col=int(col), val=0.0d0)
                  end do
               end do

               ! COO -> CRSにする作戦，完了．
            end if
         end do
         ! <<<< time consuming part <<<<
         ! <<<< time consuming part <<<<
         ! <<<< time consuming part <<<<
         ! <<<< time consuming part <<<<

         if (present(debug)) then
            if (debug) then
               print *, "[Ready!] Interface Memory "
            end if
         end if

         do i = 1, size(COO%row)
            if (.not. allocated(COO%row(i)%col)) then
               print *, "[CAUTION] >> some of doamins are not overset"
               cycle
            end if
            call heapsort(n=size(COO%row(i)%col), array=COO%row(i)%col)
         end do

         if (present(debug)) then
            if (debug) then
               print *, "[Ready!] Sort ok! "
            end if
         end if

         new_CRS_Index_Row = int(zeros(size(COO%row) + 1))
         new_CRS_Index_Row(1) = 1
         do i = 2, size(COO%row) + 1
            if (.not. allocated(COO%row(i - 1)%col)) then
               new_CRS_Index_Row(i) = new_CRS_Index_Row(i - 1)
               cycle
            end if
            new_CRS_Index_Row(i) = new_CRS_Index_Row(i - 1) + size(COO%row(i - 1)%col)
         end do

         if (present(debug)) then
            if (debug) then
               print *, "[Ready!] Row ok! "
            end if
         end if

         ! this is really heavy.
         k = 0

         ! first, count size:

         new_CRS_Index_Col = COO%getAllCol()

         !do i=1, size(COO%row)
         !    if(allocated(COO%row(i)%col ) )then
         !        k=k+1
         !        if(k==1)then
         !            new_CRS_Index_Col = COO%row(i)%col
         !        else
         !            new_CRS_Index_Col = new_CRS_Index_Col // COO%row(i)%col
         !        endif
         !    endif
         !enddo

         if (present(debug)) then
            if (debug) then
               print *, "[Ready!] Col ok! "
            end if
         end if

         call COO%remove()

         this%CRS_Val = zeros(size(new_CRS_Index_Col))
         this%CRS_Index_Col = new_CRS_Index_Col
         this%CRS_Index_Row = new_CRS_Index_Row

         if (.not. allocated(this%CRS_ID_Starts_From)) then
            allocate (this%CRS_ID_Starts_From(size(this%FEMDomains)))
            this%CRS_ID_Starts_From(1) = 1
            do i = 2, size(this%FEMDomains)
               if (associated(this%FEMDomains(i)%femdomainp)) then
                  this%CRS_ID_Starts_From(i) = this%CRS_ID_Starts_From(i - 1) + this%FEMDomains(i - 1)%femdomainp%nn()*DOF
               else
                  this%CRS_ID_Starts_From(i) = this%CRS_ID_Starts_From(i - 1) + 0
               end if
            end do
         end if

         return
      else

         if (.not. allocated(this%CRS_val)) then
            ! single-domain
            ZeroMatrix = this%femdomains(1)%femdomainp%ZeroMatrix(DOF=DOF)
            ! allocate CRS-formatted Matrix-vector (Ax = b)
            ![A]
            this%CRS_Index_Col = ZeroMatrix%col_idx
            this%CRS_Index_Row = ZeroMatrix%row_ptr
            this%CRS_val = ZeroMatrix%val

            ![b]
            this%CRS_RHS = zeros(size(this%CRS_Index_Row) - 1)
            this%CRS_x = zeros(size(this%CRS_Index_Row) - 1)

            ![x]
            this%CRS_ID_Starts_From = int(ones(1))
            !if(allocated(this%CRS_ID_Starts_From) ) deallocate(this%CRS_ID_Starts_From)

!            !通し番号を振る
!            !First, For Domains
!            !DomainID -> NodeID -> DOF(x-y-z, etc.)
!            ! count number of global unknowns
!            size_of_global_matrix = 0
!            Num_nodes_in_Domains = int(zeros(size(this%FEMDomains)) )
!            do i=1,size(this%FEMDomains)
!                if(associated(this%FEMDomains(i)%femdomainp ) )then
!                    size_of_global_matrix = size_of_global_matrix &
!                        + this%FEMDomains(i)%femdomainp%nn()*DOF
!                    Num_nodes_in_Domains(i)=this%FEMDomains(i)%femdomainp%nn()
!                endif
!            enddo
!
!            allocate(this%CRS_ID_Starts_From(size(this%FEMDomains) ))
!            this%CRS_ID_Starts_From(1) = 1
!            do i=2,size(this%FEMDomains)
!                if(associated(this%FEMDomains(i)%femdomainp ) )then
!                    this%CRS_ID_Starts_From(i) = this%CRS_ID_Starts_From(i-1) + this%FEMDomains(i-1)%femdomainp%nn()*DOF
!                else
!                    this%CRS_ID_Starts_From(i) = this%CRS_ID_Starts_From(i-1) + 0
!                endif
!            enddo
!            !
!            this%CRS_Index_Row = int(zeros(size_of_global_matrix+1))
!            this%CRS_RHS       = int(zeros(size_of_global_matrix))
!            this%CRS_x         = int(zeros(size_of_global_matrix))
!
!            num_entry_in_row   = int(zeros(size_of_global_matrix))
!
!
!            ! First, create CRS-Index-Row
!            !print *, "! First, create CRS-Index-Row"
!            !DomainID -> NodeID -> DOF(x-y-z, etc.)
!            !CRSだが，重複を許し大目に見積もる
!            ! 本当のCRS_INdex_Rowではない．あくまで，各Rowに最大いくつのcolumnが非ゼロとなりうるか．
!
!            do i=1,size(Num_nodes_in_Domains)
!
!                if(i==1)then
!                    if(associated(this%FEMDomains(i)%femdomainp ))then
!
!                        offset=0
!                        do j=1,this%FEMDomains(i)%femdomainp%ne()
!                            do k=1,this%FEMDomains(i)%femdomainp%nne()
!                                do l=1,DOF
!                                    node_id = this%FEMDomains(i)%femdomainp%mesh%elemnod(j,k)
!                                    this%CRS_Index_Row(DOF*Offset + DOF*(node_id-1)+l ) = &
!                                    this%CRS_Index_Row(DOF*Offset + DOF*(node_id-1)+l ) &
!                                        + (this%FEMDomains(i)%femdomainp%nne())*DOF
!                                    !最大でも，この数までの未知数としか係数行列を持たない
!                                enddo
!                            enddo
!                        enddo
!                    endif
!
!                else
!                    if(associated(this%FEMDomains(i)%femdomainp ))then
!                        offset=sum(Num_nodes_in_Domains(1:i-1))
!                        do j=1,this%FEMDomains(i)%femdomainp%ne()
!                            do k=1,this%FEMDomains(i)%femdomainp%nne()
!                                do l=1,DOF
!                                    node_id = this%FEMDomains(i)%femdomainp%mesh%elemnod(j,k)
!                                    this%CRS_Index_Row(DOF*Offset + DOF*(node_id-1)+l ) = &
!                                    this%CRS_Index_Row(DOF*Offset + DOF*(node_id-1)+l ) &
!                                        + (this%FEMDomains(i)%femdomainp%nne())*DOF
!                                    !最大でも，この数までの未知数としか係数行列を持たない
!                                enddo
!                            enddo
!                        enddo
!                    endif
!
!                endif
!            enddo
!
!
!            !this%CRS_Index_Rowに，あと，Interfaceのconnectivityのぶんを足す（あとで）
!
!
!            ! CRS-Index-colを作成
!            !print *, "! CRS-Index-colをAllocate"
!            this%CRS_Index_Col = int(zeros(sum(this%CRS_Index_Row(:)) ))
!            this%CRS_Val = zeros(  sum(this%CRS_Index_Row(:)) )
!
!            ! 本当のCRS_INdex_Rowにする．
!            buf_1 = this%CRS_Index_Row(1)
!            do i=2,size(this%CRS_Index_Row)
!                this%CRS_Index_Row(i) =this%CRS_Index_Row(i-1)+ this%CRS_Index_Row(i)
!            enddo
!            do i=size(this%CRS_Index_Row),2,-1
!                this%CRS_Index_Row(i) =this%CRS_Index_Row(i-1)+1
!            enddo
!            this%CRS_Index_Row(1)=1
!
!            !print *, "! CRS-Index-colを作成"
!            num_entry_in_row(:) = 0
!            do i=1,size(Num_nodes_in_Domains)
!
!                if(i==1)then
!                    if(associated(this%FEMDomains(i)%femdomainp ))then
!                        offset=0
!                        do j=1,this%FEMDomains(i)%femdomainp%ne()
!                            do k=1,this%FEMDomains(i)%femdomainp%nne()
!                                do l=1,DOF
!                                    node_id = this%FEMDomains(i)%femdomainp%mesh%elemnod(j,k)
!                                    row_id = DOF*Offset+DOF*(node_id-1)+l
!                                    do kk=1,this%FEMDomains(i)%femdomainp%nne()
!                                        do ll=1,DOF
!                                            node_id_p=this%FEMDomains(i)%femdomainp%mesh%elemnod(j,kk)
!                                            col_id =DOF*Offset+DOF*(node_id_p-1)+ll
!                                            num_entry_in_row(row_id) = num_entry_in_row(row_id) + 1
!                                            drow_offset = this%CRS_Index_Row(row_id)-1
!
!                                            this%CRS_Index_Col( drow_offset  &
!                                            + num_entry_in_row(row_id)) = col_id
!                                        enddo
!                                    enddo
!                                enddo
!                            enddo
!                        enddo
!                    endif
!
!                else
!                    if(associated(this%FEMDomains(i)%femdomainp ))then
!                        offset=sum(Num_nodes_in_Domains(1:i-1))
!                        do j=1,this%FEMDomains(i)%femdomainp%ne()
!                            do k=1,this%FEMDomains(i)%femdomainp%nne()
!                                do l=1,DOF
!                                    node_id = this%FEMDomains(i)%femdomainp%mesh%elemnod(j,k)
!                                    row_id = DOF*Offset+DOF*(node_id-1)+l
!                                    do kk=1,this%FEMDomains(i)%femdomainp%nne()
!                                        do ll=1,DOF
!                                            node_id_p=this%FEMDomains(i)%femdomainp%mesh%elemnod(j,kk)
!                                            col_id =DOF*Offset+DOF*(node_id_p-1)+ll
!                                            num_entry_in_row(row_id) = num_entry_in_row(row_id) + 1
!                                            drow_offset = this%CRS_Index_Row(row_id)-1
!
!                                            this%CRS_Index_Col( drow_offset  &
!                                            + num_entry_in_row(row_id)) = col_id
!                                        enddo
!                                    enddo
!                                enddo
!                            enddo
!                        enddo
!                    endif
!
!                endif
!            enddo
!
!            ! crs_index_colに被っているものがあるので，それを省く
!            !print *,"! crs_index_colに被っているものがあるので，それを省く"
!
!            num_entry_in_row(:) = 0
!
!            !$OMP parallel do default(shared) private(col_local,new_col_local)
!            do i=1,size(this%CRS_Index_Row)-1
!
!                col_local = int(zeros(this%CRS_Index_Row(i+1)-this%CRS_Index_Row(i)))
!                col_local(1:this%CRS_Index_Row(i+1)-this%CRS_Index_Row(i) ) = &
!                this%CRS_Index_Col(this%CRS_Index_Row(i):this%CRS_Index_Row(i+1)-1 )
!
!                new_col_local = RemoveIF(col_local,equal_to=0)
!                new_col_local = RemoveOverlap(new_col_local)
!
!
!                col_local(:)  = 0
!                col_local(1:size(new_col_local) ) = new_col_local(:)
!                this%CRS_Index_Col(this%CRS_Index_Row(i):this%CRS_Index_Row(i+1)-1 ) = &
!                col_local(1:this%CRS_Index_Row(i+1)-this%CRS_Index_Row(i) )
!                num_entry_in_row(i) = size(new_col_local)
!            enddo
!            !$OMP end parallel do
!
!            !print *, "Final"
!
!            this%CRS_Index_Col = RemoveIF(this%CRS_Index_Col,equal_to=0)
!
!
!            buf_1 = num_entry_in_row(size(num_entry_in_row)  )
!            do i=2,size(num_entry_in_row)
!                num_entry_in_row(i) =num_entry_in_row(i-1)+ num_entry_in_row(i)
!            enddo
!            do i=size(num_entry_in_row),2,-1
!                num_entry_in_row(i) =num_entry_in_row(i-1)+1
!            enddo
!            num_entry_in_row(1)=1
!
!            this%CRS_Index_Row(1:size(num_entry_in_row )) = num_entry_in_row(:)
!            this%CRS_Index_Row(size(num_entry_in_row )+1) = num_entry_in_row(size(num_entry_in_row) ) + buf_1
!
!            !do i=1,size(num_entry_in_row)
!            !    this%CRS_Index_Row(i) = sum(num_entry_in_row(1:i)) -num_entry_in_row(i) +  1
!            !enddo
!
!
!            ! then, this%CRS_Index_Row and this%CRS_Index_Col are filled.
!            this%CRS_Val = zeros(size(this%CRS_Index_Col))
!
         end if

      end if

   end subroutine
! #############################################################################

   subroutine setMatrixFEMSolver(this, DomainID, ElementID, DOF, Matrix, as_Dense)
      class(FEMSolver_), intent(inout) :: this
      integer(int32), intent(in) :: DomainID, ElementID, DOF
      real(real64), intent(in) :: Matrix(:, :)
      logical, optional, intent(in) :: as_Dense

      call this%setValue(DomainID=DomainID, ElementID=ElementID, DOF=DOF, Matrix=Matrix, as_dense=as_dense)

   end subroutine

! #############################################################################

   subroutine setVectorFEMSolver(this, DomainID, ElementID, DOF, Vector)
      class(FEMSolver_), intent(inout) :: this
      integer(int32), intent(in) :: DomainID, ElementID, DOF
      real(real64), intent(in) :: Vector(:)

      call this%setValue(DomainID=DomainID, ElementID=ElementID, DOF=DOF, Vector=Vector)

   end subroutine

! #############################################################################
   subroutine addMatrixValueFEMSolver(this, row_id, col_id, SingleValue, as_Dense)
      class(FEMSolver_), intent(inout) :: this
      integer(int32), intent(in) :: row_id, col_id
      real(real64), intent(in) :: SingleValue
      logical, optional, intent(in) :: as_Dense
      logical :: successfully_done = .false.
      integer(int64) :: i, j

      if (present(As_Dense)) then
         if (As_Dense) then
            this%A_dense(row_id, col_id) = SingleValue
            return
         end if
      end if

      if (.not. allocated(this%CRS_val)) then
         print *, "ERROR ::addMatrixValueFEMSolver >>  .not. allocated(this%CRS_val)"
         stop
      end if

      do i = this%CRS_Index_row(row_id), this%CRS_Index_row(row_id + 1) - 1
         if (this%CRS_Index_col(i) == col_id) then
            ! add
            this%CRS_Val(i) = this%CRS_Val(i) + SingleValue
            successfully_done = .true.
         end if
      end do
      if (.not. successfully_done) then
         print *, "ERROR ::addMatrixValueFEMSolver >>  the address is not allocated in CRS"
         print *, "row", row_id
         print *, "col", col_id
         stop
      end if

   end subroutine
! #############################################################################

! #############################################################################
   subroutine setValueFEMSolver(this, DomainID, ElementID, DOF, Matrix, Vector, as_Dense)
      class(FEMSolver_), intent(inout) :: this
      integer(int32), intent(in) :: DomainID, ElementID, DOF
      real(real64), optional, intent(in) :: Matrix(:, :), Vector(:)
      logical, optional, intent(in) :: as_Dense

      integer(int32) :: row_id, col_id, CRS_id, row_node_id, col_node_id
      integer(int32) :: i, j, ii, jj, eRow_id, eCol_id, k, id, nne, l_row_id, l_col_id
      integer(int32) :: g_row_id, g_col_id, g_node_id_row, g_node_id_col, offset
      integer(int32), allocatable ::  local_col_ids(:)

      if (present(as_Dense)) then
         if (as_Dense) then
            ! store as dense matrix

            if (.not. allocated(this%A_dense)) then
               this%Num_nodes_in_Domains = int(zeros(size(this%femdomains)))
               do i = 1, sizE(this%FEMDomains)
                  if (associated(this%femdomains(i)%femdomainp)) then
                     this%Num_nodes_in_Domains(i) = this%femdomains(i)%femdomainp%nn()
                  end if
               end do
               this%A_dense = zeros(sum(this%Num_nodes_in_Domains)*DOF, sum(this%Num_nodes_in_Domains)*DOF)
            end if

            do i = 1, this%femdomains(DomainID)%femdomainp%nne() ! row
               do j = 1, DOF
                  do ii = 1, this%femdomains(DomainID)%femdomainp%nne() ! col
                     do jj = 1, DOF
                        nne = this%femdomains(DomainID)%femdomainp%nne()

                        l_row_id = DOF*(i - 1) + j
                        l_col_id = DOF*(ii - 1) + jj

                        g_node_id_row = this%femdomains(DomainID)%femdomainp%mesh%elemnod(ElementID, i)
                        g_node_id_col = this%femdomains(DomainID)%femdomainp%mesh%elemnod(ElementID, ii)

                        if (DomainID == 1) then
                           offset = 0
                        else
                           offset = sum(this%Num_nodes_in_Domains(1:DomainID - 1))
                        end if

                        g_row_id = Offset*DOF + (g_node_id_row - 1)*DOF + j
                        g_col_id = Offset*DOF + (g_node_id_col - 1)*DOF + jj
                        this%A_dense(g_row_id, g_col_id) = this%A_dense(g_row_id, g_col_id) + Matrix(l_row_id, l_col_id)
                     end do
                  end do
               end do
            end do
            return
         end if
      end if

      if (.not. allocated(this%CRS_val)) then
         call this%setCRS(DOF=DOF)
      end if

      ! bugなし?
      if (present(Matrix)) then
         do i = 1, this%femdomains(DomainID)%femdomainp%nne()
            do j = 1, DOF
               row_node_id = this%femdomains(DomainID)%femdomainp%mesh%elemnod(ElementID, i)
               row_id = this%CRS_ID_Starts_From(DomainID) - 1 + DOF*(row_node_id - 1) + j

               do ii = 1, this%femdomains(DomainID)%femdomainp%nne()
                  do jj = 1, DOF
                     col_node_id = this%femdomains(DomainID)%femdomainp%mesh%elemnod(ElementID, ii)
                     col_id = this%CRS_ID_Starts_From(DomainID) - 1 + DOF*(col_node_id - 1) + jj
                     local_col_ids = zeros(this%CRS_Index_Row(row_id + 1) - this%CRS_Index_Row(row_id))
                     local_col_ids(:) = this%CRS_Index_Col(this%CRS_Index_Row(row_id):this%CRS_Index_Row(row_id + 1) - 1)

                     id = -1
                     do k = 1, size(local_col_ids)
                        if (local_col_ids(k) == col_id) then
                           id = k
                           exit
                        end if
                     end do
                     if (id == -1) then
                        print *, "[ERROR] ::setValueFEMSolver (fill-in request) No memory is allocated in CRS format for "
                        print *, "Domain  : ", DomainID
                        print *, "Element : ", ElementID
                        print *, "Nodes    : ", i, ii
                        print *, "Dims    : ", j, jj
                        stop
                     end if
                     CRS_id = this%CRS_Index_Row(row_id) - 1 + id
                     eRow_id = DOF*(i - 1) + j
                     eCol_id = DOF*(ii - 1) + jj
                        !!$OMP atomic
                     this%CRS_val(CRS_id) = this%CRS_val(CRS_id) + Matrix(eRow_id, eCol_id)
                  end do
               end do
            end do
         end do
      end if

      if (present(Vector)) then
         do i = 1, this%femdomains(DomainID)%femdomainp%nne()
            do j = 1, DOF
               row_node_id = this%femdomains(DomainID)%femdomainp%mesh%elemnod(ElementID, i)
               row_id = this%CRS_ID_Starts_From(DomainID) - 1 + DOF*(row_node_id - 1) + j
               eRow_id = DOF*(i - 1) + j
                !!$OMP atomic
               this%CRS_RHS(Row_id) = this%CRS_RHS(Row_id) + Vector(eRow_id)
            end do
         end do
      end if

   end subroutine
! ###################################################################
   subroutine fixFEMSolver(this, DomainID, IDs, FixValue, FixValues)
      class(FEMSolver_), intent(inout) :: this

      real(real64), optional, intent(in) :: FixValue
      real(real64), optional, intent(in) :: FixValues(:)
      integer(int32), optional, intent(in) :: DomainID
      integer(int32), intent(in) :: IDs(:)

      integer(int32), allocatable :: buf(:)
      integer(int32), allocatable :: buf_real(:), ned(:)
      integer(int32) :: i, n, offset

      offset = 0
      if (present(DomainID)) then
         ned = this%number_of_point_per_domain
         if (DomainID >= 2) then
            do i = 1, DomainID - 1
               offset = offset + ned(i)
            end do
         end if
      end if

      ! fix unknowns for linear solvers
      ! fix IDs(:)-th amplitudes as zero (Fixed boundary)
      ! only create list
      if (.not. allocated(this%fix_lin_exists)) then
         n = size(this%CRS_RHS)
         allocate (this%fix_lin_exists(n))

         this%fix_lin_exists_values = zeros(n)
         this%fix_lin_exists(:) = .false.

         do i = 1, size(IDs)
            if (IDs(i) < 1) cycle

            if (IDs(i) > size(this%fix_lin_exists)) cycle
            if (IDs(i) <= 0) cycle

            this%fix_lin_exists(IDs(i)) = .true.
            if (present(FixValue)) then
               this%fix_lin_exists_values(IDs(i) + offset) = FixValue
            elseif (present(FixValues)) then
               this%fix_lin_exists_values(IDs(i) + offset) = FixValues(i)
            else
               this%fix_lin_exists_values(IDs(i) + offset) = 0.0d0
            end if

         end do
      else

         do i = 1, size(IDs)
            if (IDs(i) < 1) cycle
            this%fix_lin_exists(IDs(i)) = .true.
            if (present(FixValue)) then
               this%fix_lin_exists_values(IDs(i) + offset) = FixValue
            elseif (present(FixValues)) then
               this%fix_lin_exists_values(IDs(i) + offset) = FixValues(i)
            else
               this%fix_lin_exists_values(IDs(i) + offset) = 0.0d0
            end if
         end do
      end if

   end subroutine

!
   function diagFEMSolver(this, cell_centered) result(diag_vector)
      class(FEMSolver_), intent(in) :: this
      logical, optional, intent(in)  :: cell_centered
      real(real64), allocatable :: diag_vector(:)

      integer(int32) :: row, col, id

      if (present(cell_centered)) then
         if (cell_centered) then

            ! diagonal components of CRS matrix
            if (allocated(this%CRS_val)) then
               diag_vector = zeros(size(this%CRS_Index_row) - 1)
               do row = 1, size(this%CRS_Index_row) - 1
                  do id = this%CRS_Index_row(row), this%CRS_Index_row(row + 1) - 1
                     diag_vector(row) = diag_vector(row) + this%CRS_val(id)
                  end do
               end do
            end if
         end if
      end if

      ! diagonal components of CRS matrix
      if (allocated(this%CRS_val)) then
         diag_vector = zeros(size(this%CRS_Index_row) - 1)
         do row = 1, size(this%CRS_Index_row) - 1
            do id = this%CRS_Index_row(row), this%CRS_Index_row(row + 1) - 1
               col = this%CRS_Index_col(id)
               if (col == row) then
                  diag_vector(row) = this%CRS_val(id)
               end if
            end do
         end do
      end if

   end function

! ####################################################

!
   subroutine matmulDiagMatrixFEMSolver(this, diagMat)
      class(FEMSolver_), intent(inout) :: this
      real(real64), intent(in) :: diagMat(:)
      integer(int32) :: n

      integer(int32) :: row, col, id

      n = size(diagMat)
!> diag is diagonal component of n x n matrix

      ! diagonal components of CRS matrix
      if (allocated(this%CRS_val)) then
         do row = 1, size(this%CRS_Index_row) - 1
            do id = this%CRS_Index_row(row), this%CRS_Index_row(row + 1) - 1
               col = this%CRS_Index_col(id)
               if (col == row) then
                  this%CRS_val(id) = this%CRS_val(id)*diagMat(row)
               end if
            end do
         end do
      end if

   end subroutine
! ###################################################################
   subroutine loadMatrixFEMSolver(this, from)
      class(FEMSolver_), intent(inout) :: this
      character(1), intent(in) :: from

      ! overwrite this%CRS_*
      select case (from)
      case ("A")
         this%CRS_val = this%A_CRS_val
         this%CRS_Index_Col = this%A_CRS_Index_Col
         this%CRS_Index_Row = this%A_CRS_Index_Row
      case ("B")
         this%CRS_val = this%B_CRS_val
         this%CRS_Index_Col = this%B_CRS_Index_Col
         this%CRS_Index_Row = this%B_CRS_Index_Row
      case default
         print *, "ERROR :: loadMatrixFEMSolver >> from should be A or B."
         stop
      end select

   end subroutine

! ###################################################################
   subroutine saveMatrixFEMSolver(this, name, CRS_as_dense, if_dense_exists, zero_or_nonzero)
      class(FEMSolver_), intent(in) :: this
      character(*), intent(in) :: name
      logical, optional, intent(in) :: CRS_as_dense, if_dense_exists, zero_or_nonzero
      integer(int32) :: i, j, k, n
      real(real64), allocatable :: row_vector(:)
      type(IO_)::f

      if (present(if_dense_exists)) then
         if (if_dense_exists) then
            if (.not. allocated(this%A_dense)) return
            call f%open(name + "_dense.csv", "w")
            n = size(this%A_dense, 1)

            do i = 1, n
               do j = 1, n - 1
                  write (f%fh, '(A)', advance='no') str(this%A_dense(i, j)) + ","
               end do
               write (f%fh, '(A)', advance='yes') str(this%A_dense(i, n))
            end do
            call f%close()
         end if
         return
      end if

      if (.not. allocated(this%CRS_val)) then
         print *, "[ERROR] saveMatrixFEMSolver >> (.not.allocated(this%CRS_val) )"
         return
      end if

      if (present(CRS_as_dense)) then
         if (CRS_as_dense) then
            call f%open(name + "_dense.csv", "w")
            n = size(this%CRS_Index_Row) - 1
            do i = 1, n
               row_vector = zeros(n)
               do j = this%CRS_Index_Row(i), this%CRS_Index_Row(i + 1) - 1
                  row_vector(this%CRS_Index_Col(j)) = this%CRS_val(j)
               end do
               do j = 1, size(row_vector) - 1
                  if (present(zero_or_nonzero)) then
                     if (zero_or_nonzero) then
                        if (row_vector(j) == 0.0d0) then
                           write (f%fh, '(A)', advance='no') "0"
                        else
                           write (f%fh, '(A)', advance='no') "*"
                        end if
                        cycle
                     end if
                  end if
                  write (f%fh, '(A)', advance='no') str(row_vector(j)) + ","
               end do
               if (present(zero_or_nonzero)) then
                  if (zero_or_nonzero) then
                     if (row_vector(n) == 0.0d0) then
                        write (f%fh, '(A)', advance='yes') "0"
                     else
                        write (f%fh, '(A)', advance='yes') "*"
                     end if
                     cycle
                  end if
               end if
               write (f%fh, '(A)', advance='yes') str(row_vector(n))
            end do
            call f%close()
         end if
      else

         call f%open(name + "_data.txt", "w")
         call f%write(this%CRS_val)
         call f%close()

         call f%open(name + "_indices.txt", "w")
         call f%write(this%CRS_Index_Col)
         call f%close()

         call f%open(name + "_indptr.txt", "w")
         call f%write(int(this%CRS_Index_Row))
         call f%close()

      end if

   end subroutine

! ###################################################################
   subroutine zerosFEMSolver(this)
      class(FEMSolver_), intent(inout)::this

      this%CRS_val(:) = 0.0d0
      if (allocated(this%CRS_RHS)) then
         this%CRS_RHS(:) = 0.0d0
      end if

      if (allocated(this%A_dense)) then
         this%A_dense(:, :) = 0.0d0
      end if

   end subroutine

! ###################################################################

!function eigFEMSolver(this,num_eigen,tol,eigen_value,as_dense) result(eig_vec)
!    class(FEMSolver_),intent(in)::this
!    real(real64),allocatable :: eig_vec(:,:),dense_mat(:,:)
!    real(real64),optional,allocatable,intent(inout) :: eigen_value(:)
!    real(real64),intent(in) :: tol
!    integer(int32),optional,intent(in) :: num_eigen
!    integer(int32) :: ndim
!    logical,optional,intent(in) :: as_dense
!    !> default =>> get eigen vectors of this%CRS
!    !> eigens(:,0) are eigen values
!    !> eigens(:,n) are n-th eigen vectors
!    if(present(as_Dense))then
!        if(as_Dense)then
!            call to_Dense(this%CRS_val,this%CRS_index_col,this%CRS_index_row,&
!                dense_mat)
!            dense_mat = 0.50d0*(dense_mat + transpose(dense_mat) )
!            call eigenValueAndVector(A=dense_mat,&
!                lambda=eigen_value,x=eig_vec,tol=tol)
!            return
!        endif
!    endif
!
!    ndim = size(this%CRS_Index_Row) - 1
!
!    eigen_value = zeros(num_eigen)
!    eig_vec =  LOBPCG_sparse(&
!        A_val=this%CRS_val,&
!        A_col=this%CRS_index_col,&
!        A_rowptr=this%CRS_index_row,&
!        lambda_min=eigen_value,&
!        tolerance=tol)
!
!
!end function

! ###################################################################
   subroutine LanczosMethod(this, eigen_value, Eigen_vectors, max_itr)
      clasS(FEMSolver_), intent(inout) :: this
      real(real64), allocatable :: eigen_value(:), Eigen_vectors(:, :)
      real(real64), allocatable :: w(:)
      real(real64)::alpha, beta
      integer(int32), intent(in) :: max_itr
      integer(int32) :: num_dim
      integer(int32) :: i, j

      num_dim = size(this%CRS_index_row) - 1
      !http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/cgi-bin/wiki/index.php?%CF%A2%CE%A91%BC%A1%CA%FD%C4%F8%BC%B0%A1%A7Lanczos%CB%A1

      eigen_value = zeros(num_dim)
      eigen_vectors = eyesMatrix(num_dim, num_dim)
      print *, "Lanczos method is not implemented."
      do i = 1, max_itr

      end do

   end subroutine
! ###################################################################

   recursive subroutine eigFEMSolver(this, num_eigen, eigen_value, eigen_vectors)
      ! solve Ku = \lambda M x by LAPACK
      clasS(FEMSolver_), intent(inout) :: this
      integer(int32), optional, intent(in)::num_eigen
      !logical,optional,intent(in) :: Lanczos

      !>>>>>>>>>>>>>> INPUT
      integer(int32) :: ITYPE = 1   ! A*x = (lambda)*B*x
      character(1) :: JOBZ = 'V' ! Compute eigenvalues and eigenvectors.
      character(1) :: UPLO = 'U' ! Upper triangles of A and B are stored;
      !<<<<<<<<<<<<<< INPUT

      integer(int32) :: N ! order of matrix
      real(real64), allocatable :: AP(:)
      real(real64), allocatable :: BP(:)
      real(real64), allocatable :: W(:)
      real(real64), allocatable :: Z(:, :), M(:)
      real(real64), allocatable :: WORK(:), ID(:)
      real(real64), allocatable, intent(inout) :: eigen_value(:)
      real(real64), allocatable, intent(inout) :: eigen_vectors(:, :)
      integer(int32), allocatable :: IWORK(:), IDS(:)
      integer(int32) :: LDZ
      integer(int32) :: LWORK
      integer(int32) :: LIWORK
      integer(int32) :: INFO
      integer(int32) :: from, to, k, j, i
      integer(int32), allocatable :: new_id_from_old_id(:)
      real(real64), allocatable :: dense_mat(:, :)
      !logical :: use_lanczos
      type(IO_) :: f
      type(CRS_) :: crs, A, B

      !use_lanczos = .false.
      !if(present(Lanczos) )then
      !    use_lanczos = Lanczos
      !endif

      A = this%getCRS("A")
      B = this%getCRS("B")

      if (allocated(this%fix_eig_IDs)) then
         ! amplitudes are zero@ this%fix_eig_IDs
         ! remove from problem [A][U] = w[B][U]
         ! sort before it

         if (size(this%fix_eig_IDs) >= 1) then
            ! first, for [A]
            call heapsort(n=size(this%fix_eig_IDs), array=this%fix_eig_IDs)
            call reduce_crs_matrix(CRS_val=A%val, CRS_col=A%col_idx, &
                                   CRS_rowptr=A%row_ptr, remove_IDs=this%fix_eig_IDs)
            call reduce_crs_matrix(CRS_val=B%val, CRS_col=B%col_idx, &
                                   CRS_rowptr=B%row_ptr, remove_IDs=this%fix_eig_IDs)
         end if
      end if

      !>>>>>>>>>>>>>> INPUT
      N = size(A%row_ptr) - 1
      LDZ = input(default=N, option=num_eigen)
      LWORK = 1 + 6*N + 2*N**2
      LIWORK = 3 + 5*N
      !<<<<<<<<<<<<<< INPUT

      if (this%use_LOBPCG) then
         print *, ">> Solver :: LOBPCG"
         call LOBPCG( &
            A=A, &
            B=B, &
            X=Z, lambda=W, &
            m=input(default=this%LOBPCG_NUM_MODE, option=num_eigen), &
            MAX_ITR=this%LOBPCG_MAX_ITR, &
            TOL=this%LOBPCG_TOL, &
            debug=this%debug)

      else
         !>>>>>>>>>>>>>>  INPUT/OUTPUT
         AP = zeros(N*(N + 1)/2)
         BP = zeros(N*(N + 1)/2)
         ! Upper triangle matrix
         AP = UpperTriangularMatrix(CRS_val=A%val, CRS_col=A%col_idx, &
                                    CRS_rowptr=A%row_ptr)
         BP = UpperTriangularMatrix(CRS_val=B%val, CRS_col=B%col_idx, &
                                    CRS_rowptr=B%row_ptr)
         !<<<<<<<<<<<<<< INPUT/OUTPUT

         !>>>>>>>>>>>>>> INPUT
         N = size(A%row_ptr) - 1
         LDZ = input(default=N, option=num_eigen)
         LWORK = 1 + 6*N + 2*N**2
         LIWORK = 3 + 5*N
         !<<<<<<<<<<<<<< INPUT

         !>>>>>>>>>>>>>>  OUTPUT
         W = zeros(N)
         Z = zeros(LDZ, N)
         WORK = zeros(LWORK)
         IWORK = zeros(LIWORK)
         INFO = 0
         !<<<<<<<<<<<<<< OUTPUT

         print *, ">> Solver :: LAPACK/DSPGVD"

         call DSPGVD(ITYPE, JOBZ, UPLO, N, AP, BP, W, Z, LDZ, WORK, &
                     LWORK, IWORK, LIWORK, INFO)
      end if
      eigen_value = w
      if (allocated(this%fix_eig_IDs)) then
         ! U(this%fix_eig_IDs(i),: ) = 0.0d0
         if (size(this%fix_eig_IDs) >= 1) then
            new_id_from_old_id = zeros(N)

            !new_id_from_old_id(Dirichlet境界を除いた固有ベクトルzのj番成分が，もとの何番目に対応するか)

            k = 0
            do j = 1, this%fix_eig_IDs(1) - 1
               k = k + 1
               new_id_from_old_id(k) = k
            end do

            do i = 2, size(this%fix_eig_IDs)
               from = this%fix_eig_IDs(i - 1) + 1
               to = this%fix_eig_IDs(i) - 1
               do j = from, to
                  k = k + 1
                  if (k > size(new_id_from_old_id)) cycle
                  new_id_from_old_id(k) = j
               end do
            end do

            do j = this%fix_eig_IDs(size(this%fix_eig_IDs)) + 1, N + size(this%fix_eig_IDs)
               k = k + 1
               if (k > size(new_id_from_old_id)) cycle
               new_id_from_old_id(k) = j
            end do

            eigen_vectors = zeros(size(Z, 1) + size(this%fix_eig_IDs), size(Z, 2))
            do i = 1, size(Z, 2)
               do j = 1, size(new_id_from_old_id, 1)
                  eigen_vectors(new_id_from_old_id(j), i) = Z(j, i)
               end do
            end do

         else
            eigen_vectors = Z
         end if
      else
         eigen_vectors = Z
      end if

   end subroutine
! ###################################################################

   subroutine keepThisMatrixAsFEMSolver(this, As)
      class(FEMSolver_), intent(inout) :: this
      character(*), intent(in) :: As ! any length
      integer(int32) :: i
      !character(1),intent(in) :: As ! [A] or [B]

!    if(As == "A")then
!        this%A_CRS_Index_%CRS_Index_Col
!        this%A_CRS_Index_Row = this%CRS_Index_Row
!        this%A_CRS_val       = this%CRS_val
!        this%A_empty         = .false.
!        return
!    endif
!
!    if(As == "B")then
!        this%B_CRS_Index_Col = this%CRS_Index_Col
!        this%B_CRS_Index_Row = this%CRS_Index_Row
!        this%B_CRS_val       = this%CRS_val
!        this%B_empty         = .false.
!        return
!    endif

      ! for other cases:
      if (.not. allocated(this%CRS_matrix_stack)) then
         allocate (this%CRS_matrix_stack(this%MAXSIZE_CRS_MATRIX_STACK))
         allocate (this%CRS_matrix_stack_key%pages(this%MAXSIZE_CRS_MATRIX_STACK))
         this%LAST_CRS_MATRIX_STACK = 0
      end if
      if (this%LAST_CRS_MATRIX_STACK == 0) then
         this%CRS_matrix_stack(1)%row_ptr = this%CRS_Index_Row
         this%CRS_matrix_stack(1)%col_idx = this%CRS_Index_Col
         this%CRS_matrix_stack(1)%val = this%CRS_Val
         this%CRS_matrix_stack_key%pages(1)%charvalue = As
         this%LAST_CRS_MATRIX_STACK = 1
      else
         do i = 1, this%LAST_CRS_MATRIX_STACK
            if (this%CRS_matrix_stack_key%pages(i)%charvalue == As) then
               this%CRS_matrix_stack(i)%row_ptr = this%CRS_Index_Row
               this%CRS_matrix_stack(i)%col_idx = this%CRS_Index_Col
               this%CRS_matrix_stack(i)%val = this%CRS_Val
               return
            else
               cycle
            end if
         end do
         this%LAST_CRS_MATRIX_STACK = this%LAST_CRS_MATRIX_STACK + 1

         if (size(this%CRS_matrix_stack) < this%LAST_CRS_MATRIX_STACK) then
            print *, "[ERROR] keepThisMatrixAsFEMSolver size(this%CRS_matrix_stack)<this%LAST_CRS_MATRIX_STACK)"
            stop
         else
            i = this%LAST_CRS_MATRIX_STACK
            this%CRS_matrix_stack(i)%row_ptr = this%CRS_Index_Row
            this%CRS_matrix_stack(i)%col_idx = this%CRS_Index_Col
            this%CRS_matrix_stack(i)%val = this%CRS_Val
            this%CRS_matrix_stack_key%pages(i)%charvalue = As
         end if

      end if

   end subroutine

! ################################################################
   subroutine fix_eigFEMSolver(this, IDs)
      class(FEMSolver_), intent(inout) :: this
      integer(int32), intent(in) :: IDs(:)
      integer(int32), allocatable :: buf(:)

      ! fix IDs(:)-th amplitudes as zero (Fixed boundary)
      ! only create list
      if (.not. allocated(this%fix_eig_IDs)) then
         this%fix_eig_IDs = IDs
      elseif (size(this%fix_eig_IDs) == 0) then
         this%fix_eig_IDs = IDs
      else

         buf = this%fix_eig_IDs
         this%fix_eig_IDs = zeros(size(buf) + size(IDs))
         this%fix_eig_IDs(1:size(buf)) = buf(:)
         this%fix_eig_IDs(size(buf) + 1:) = IDs(:)
      end if

      this%fix_eig_IDs = unique(this%fix_eig_IDs)

   end subroutine
! #####################################################

! #####################################################
   function solveFEMSolver_UserDefinedLinearSolver(this, LinearSolver, x0) result(x)
      class(FEMSolver_), intent(inout) :: this
      real(real64), optional, intent(in) :: x0(:)
      real(real64), allocatable :: x(:)

      ! CRS formatted Linear solver
      interface
         subroutine LinearSolver(row_ptr, col_idx, val, rhs, x)
            use iso_fortran_env
            implicit none
            real(real64), intent(in) :: val(:), rhs(:)
            real(real64), intent(inout) :: x(:)
            integer(int32), intent(in) :: col_idx(:)
            integer(int64), intent(in) :: row_ptr(:)

         end subroutine
      end interface

      if (present(x0)) then
         x = x0
      else
         x = zeros(size(this%CRS_index_row) - 1)
      end if

      ! 外部ソルバの利用(ただし，subroutineのみ)
      call LinearSolver(this%CRS_index_row, this%CRS_Index_Col, this%CRS_val, &
                        this%CRS_RHS, x)

      this%CRS_x = x

   end function

! #####################################################

! #####################################################
   function solveFEMSolver_UserDefinedLinearSolverAsFunc(this, LinearSolver, x0) result(x)
      class(FEMSolver_), intent(inout) :: this
      real(real64), intent(in) :: x0(:)
      real(real64), allocatable :: x(:)

      ! CRS formatted Linear solver
      interface
         function LinearSolver(row_ptr, col_idx, val, rhs, x0) result(x)
            use iso_fortran_env
            implicit none
            real(real64), intent(in) :: val(:), rhs(:)
            real(real64), intent(in) :: x0(:)
            real(real64), allocatable :: x(:)
            integer(int32), intent(in) :: col_idx(:)
            integer(int64), intent(in) :: row_ptr(:)

         end function
      end interface

      ! 外部ソルバの利用(ただし，subroutineのみ)
      x = LinearSolver(this%CRS_index_row, this%CRS_Index_Col, this%CRS_val, &
                       this%CRS_RHS, x0)

      this%CRS_x = x

   end function

! #####################################################

! #####################################################
   function solveFEMSolver(this, algorithm, preconditioning, x0) result(x)
      class(FEMSolver_), intent(inout) :: this
      real(real64), allocatable :: x(:), dense_mat(:, :), fix_value(:)
      real(real64), optional, intent(in) :: x0(:)
      integer(int64) :: i, j, ElementID, col, row_ptr, col_row_fix
      logical, allocatable :: need_fix(:)
      character(*), optional, intent(in) :: algorithm, preconditioning

      type(IO_) :: f

      if (allocated(this%fix_lin_exists)) then
         need_fix = this%fix_lin_exists
         fix_value = this%fix_lin_exists_values

         do i = 1, size(this%CRS_Index_Row) - 1 !すべての行
            do col = this%CRS_index_row(i), this%CRS_index_row(i + 1) - 1

               if (need_fix(this%CRS_index_col(col))) then
                  this%CRS_RHS(i) = this%CRS_RHS(i) &
                                    - this%CRS_val(col)*fix_value(this%CRS_index_col(col)) ! 移項
               end if
            end do
         end do
      end if

      if (allocated(this%fix_lin_exists)) then

         ! 右辺ベクトルに強制値を導入
         ! for each boundary conditioned-node
         do i = 1, size(this%CRS_RHS)
            if (this%fix_lin_exists(i)) then
               this%CRS_RHS(i) = this%fix_lin_exists_values(i)
            end if
         end do

         do i = 1, size(this%CRS_Index_row) - 1
            do j = this%CRS_Index_row(i), this%CRS_Index_row(i + 1) - 1
               if (this%fix_lin_exists(i)) then
                  this%CRS_val(j) = 0.0d0
               end if

               if (this%fix_lin_exists(this%CRS_Index_Col(j))) then
                  this%CRS_val(j) = 0.0d0
               end if
            end do
         end do

         do i = 1, size(this%CRS_Index_row) - 1
            do j = this%CRS_Index_row(i), this%CRS_Index_row(i + 1) - 1
               if (this%fix_lin_exists(i) .and. &
                   this%CRS_Index_Col(j) == i) then
                  this%CRS_val(j) = 1.0d0
               else
                  cycle
               end if

            end do
         end do

      end if

      if (this%debug) then
         print *, "[ok] b.c. loaded"
      end if

      if (present(x0)) then
         x = x0
      else
         x = zeros(size(this%CRS_RHS))
      end if

      if (present(preconditioning)) then
         if (preconditioning == "PointJacobi") then

            if (this%debug) then
               print *, "PBiCGSTAB (PointJacobi)"
            end if
            call JacobiPreconditionerCRS(val=this%CRS_val, row_ptr=this%CRS_index_row, &
                                         col_idx=this%CRS_Index_col, rhs=this%CRS_RHS)

         elseif (preconditioning == "incompleteLU" .or. preconditioning == "ILU") then

            if (this%debug) then
               print *, "PBiCGSTAB (ILU(0))"
            end if
            call bicgstab_CRS_ILU(this%CRS_val, this%CRS_index_row, this%CRS_index_col, &
                                  x, this%CRS_RHS, this%itrmax, this%er0, this%relative_er, this%debug, &
                                  this%ILU_MATRIX)
            
         else
            print *, "[Warning!] :: FEMSolver :: invalid preconditioning"
            stop
         end if
      else

         if (associated(this%mpi_target)) then
            call this%MPI_BICGSTAB(x)
         else
            if (present(algorithm)) then
               if (algorithm == "GaussJordan") then
                  call gauss_jordan_crs(this%CRS_val, this%CRS_index_row, this%CRS_index_col, &
                                        x, this%CRS_RHS, size(this%CRS_RHS))
                  
               else
                  call bicgstab_CRS_2(this%CRS_val, this%CRS_index_row, this%CRS_index_col, &
                                      x, this%CRS_RHS, this%itrmax, this%er0, this%relative_er, this%debug)
                  
               end if
            else
               call bicgstab_CRS_2(this%CRS_val, this%CRS_index_row, this%CRS_index_col, &
                                   x, this%CRS_RHS, this%itrmax, this%er0, this%relative_er, this%debug)
               
            end if
         end if
   
      end if
      
      if (allocated(this%fix_lin_exists)) then

         ! 右辺ベクトルに強制値を導入
         ! for each boundary conditioned-node
         do i = 1, size(this%CRS_RHS)
            if (this%fix_lin_exists(i)) then
               x(i) = this%fix_lin_exists_values(i)
            end if
         end do
      endif

   end function
! #####################################################
   subroutine removeFEMSolver(this, only_matrices)
      class(FEMSolver_), intent(inout) :: this
      logical, optional, intent(in) :: only_matrices

      if (allocated(this%CRS_val)) deallocate (this%CRS_val)!(:)
      if (allocated(this%CRS_Index_Col)) deallocate (this%CRS_Index_Col)!(:)
      if (allocated(this%CRS_Index_Row)) deallocate (this%CRS_Index_Row)!(:)
      if (allocated(this%CRS_RHS)) deallocate (this%CRS_RHS)!(:)

      !> General Eigen Value Problem
      !> [A]{x} = (lambda)[B]{x}

      if (allocated(this%A_CRS_val)) deallocate (this%A_CRS_val)!(:)
      if (allocated(this%A_CRS_Index_Col)) deallocate (this%A_CRS_Index_Col)!(:)
      if (allocated(this%A_CRS_Index_Row)) deallocate (this%A_CRS_Index_Row)!(:)
      this%A_empty = .true.

      if (allocated(this%B_CRS_val)) deallocate (this%B_CRS_val)!(:)
      if (allocated(this%B_CRS_Index_Col)) deallocate (this%B_CRS_Index_Col)!(:)
      if (allocated(this%B_CRS_Index_Row)) deallocate (this%B_CRS_Index_Row)!(:)
      this%B_empty = .true.

      if (present(only_matrices)) then
         if (only_matrices) then
            return
         end if
      end if
      if (allocated(this%femdomains)) deallocate (this%femdomains)!(:)
      if (allocated(this%DomainIDs)) deallocate (this%DomainIDs)!(:)

      if (allocated(this%IfaceElemConnectivity)) deallocate (this%IfaceElemConnectivity)!(:,:)
      if (allocated(this%IfaceElemDomainID)) deallocate (this%IfaceElemDomainID)!(:,:)

      this%initialized = .false.
      this%InterfaceExist = .false.
      this%debug = .false.

      if (allocated(this%fix_eig_IDs)) deallocate (this%fix_eig_IDs)!(:)

      if (allocated(this%fix_lin_exists)) deallocate (this%fix_lin_exists)!(:)
      if (allocated(this%fix_lin_exists_Values)) deallocate (this%fix_lin_exists_Values)!(:)

      ! if(allocated( this%fix_lin_IDs)) deallocate(this%fix_lin_IDs)!(:)
      ! if(allocated( this%fix_lin_Values)) deallocate(this%fix_lin_Values)!(:)

      if (allocated(this%CRS_x)) deallocate (this%CRS_x)!(:)
      if (allocated(this%CRS_ID_Starts_From)) deallocate (this%CRS_ID_Starts_From)!(:)

      ! dense matrix
      if (allocated(this%A_dense)) deallocate (this%A_dense)!(:,:)
      if (allocated(this%Num_nodes_in_Domains)) deallocate (this%Num_nodes_in_Domains)!(:)

      ! with MPI
      if (associated(this%MPI_target)) nullify (this%MPI_target)
      if (allocated(this%Link_Table)) deallocate (this%Link_Table)

      this%LINK_TABLE_INIT_SIZE = 1000
      this%Link_num = 0
      this%itrmax = 100000
      this%er0 = dble(1.0e-10)
      this%relative_er = dble(1.0e-10)

   end subroutine
! #####################################################
   subroutine setEbOMFEMSolver(this, penalty, DOF)
      class(FEMSolver_), intent(inout) :: this

      real(real64), intent(in) :: penalty
      integer(int32), intent(in) :: DOF
      integer(int32) :: DomainID, ElementID, GaussPointID, myDomainID, pairDomainID
      integer(int32) :: i, j, k, m, n, row_id, col_id, row_Domain_id, col_domain_id
      real(real64) :: singleValue
      real(real64), allocatable :: A_ij(:, :), position(:)
      type(ShapeFunction_) :: sf

      ! it was
      ! do DomainID=1,2
      do DomainID = 1, size(this%FEMDomains)
         do i = 1, this%FEMDomains(DomainID)%femdomainp &
            %numOversetElements()

            if (this%FEMDomains(DomainID)%femdomainp &
                %OversetConnect(i)%active) then

               ElementID = this%FEMDomains(DomainID)%femdomainp &
                           %OversetConnect(i)%ElementID
               GaussPointID = this%FEMDomains(DomainID)%femdomainp &
                              %OversetConnect(i)%GaussPointID

               position = this%FEMDomains(DomainID)%femdomainp &
                          %OversetConnect(i)%position

               myDomainID = this%FEMDomains(DomainID)%femdomainp &
                            %OversetConnect(i)%DomainIDs12(1)
               n = size(this%FEMDomains(DomainID)%femdomainp &
                        %OversetConnect(i)%DomainIDs12)
               pairDomainID = this%FEMDomains(DomainID)%femdomainp &
                              %OversetConnect(i)%DomainIDs12(n)

               if (GaussPointID > 0) then
                  ! GPP

                  sf = this%FEMDomains(myDomainID)%femdomainp%mesh%getShapeFunction(ElementID, GaussPointID)

                  sf%ElementID = ElementID

                  A_ij = penalty*this%FEMDomains(pairDomainID)%femdomainp%connectMatrix( &
                         position=position, &
                         DOF=DOF, &
                         shapefunction=sf)
               else
                  ! P2P
                  A_ij = penalty*this%FEMDomains(pairDomainID)%femdomainp &
                         %connectMatrix(position=position, DOF=DOF)
               end if

               ! assemble them
               do j = 1, size(this%FEMDomains(DomainID)%femdomainp &
                              %OversetConnect(i)%interConnect)
                  do k = 1, size(this%FEMDomains(DomainID)%femdomainp &
                                 %OversetConnect(i)%interConnect)
                     do m = 1, DOF
                        do n = 1, DOF

                           row_Domain_id = this%FEMDomains(DomainID)%femdomainp &
                                           %OversetConnect(i)%DomainIDs12(j)
                           col_Domain_id = this%FEMDomains(DomainID)%femdomainp &
                                           %OversetConnect(i)%DomainIDs12(k)

                           row_id = &
                              +this%FEMDomains(DomainID)%femdomainp &
                              %OversetConnect(i)%interConnect(j)
                           col_id = &
                              +this%FEMDomains(DomainID)%femdomainp &
                              %OversetConnect(i)%interConnect(k)

                           row_id = this%CRS_ID_Starts_From(row_Domain_id) - 1 + (row_id - 1)*DOF + m
                           col_id = this%CRS_ID_Starts_From(col_Domain_id) - 1 + (col_id - 1)*DOF + n

                           singleValue = A_ij((j - 1)*DOF + m, (k - 1)*DOF + n)

                           call this%addMatrixValue( &
                              row_id=row_id, &
                              col_id=col_id, &
                              singleValue=singleValue &
                              )
                        end do
                     end do
                  end do
               end do
               !call obj%solver%assemble(&
               !    connectivity=InterConnect,&
               !    DOF=domain2%nd() ,&
               !    eMatrix=A_ij,&
               !    DomainIDs=this%FEMDomains(DomainID)%OversetConnect(i)%DomainIDs12)

               !call solver%setMatrix(DomainID=DomainID,ElementID=ElementID,DOF=3,&
               !    Matrix=this%FEMDomains(DomainID)%StiffnessMatrix(ElementID=ElementID,E=300000.0d0, v=0.330d0),&
               !    )
            end if
         end do
      end do

   end subroutine
! #####################################################
   function num_EbOFEM_projection_point_FEMSolver(this) result(ret)
      class(FEMSolver_), intent(in) :: this
      integer(int32) :: ret, DomainID, i

      ! count number of projection points
      ret = 0
      do DomainID = 1, size(this%FEMDomains)
         do i = 1, this%FEMDomains(DomainID)%femdomainp &
            %numOversetElements()
            if (this%FEMDomains(DomainID)%femdomainp%OversetConnect(i)%active) then
               ret = ret + 1
            end if
         end do
      end do

   end function

! #####################################################
   function get_gap_function_FEMSolver(this, DOF, X) result(gap_function) ! projection point-wise vector
      class(FEMSolver_), intent(in) :: this
      integer(int32), intent(in) :: DOF
      real(real64), intent(in) :: X(:)

      real(real64), allocatable :: gap_function(:)

      ! compute g(x) = \Sigma [N]{x}
      integer(int32) :: DomainID, ElementID, GaussPointID, myDomainID, pairDomainID
      integer(int32) :: i, j, k, m, n, row_id, col_id, row_Domain_id, col_domain_id, itr, node_id
      real(real64) :: singleValue
      real(real64), allocatable :: N_I(:), position(:), Cvec(:)
      type(ShapeFunction_) :: sf

      ! projection point-wise vector
      gap_function = zeros(this%num_EbOFEM_projection_point())

      ! do DomainID=1,2
      itr = 0
      do DomainID = 1, size(this%FEMDomains)
         do i = 1, this%FEMDomains(DomainID)%femdomainp &
            %numOversetElements()

            if (this%FEMDomains(DomainID)%femdomainp &
                %OversetConnect(i)%active) then

               ElementID = this%FEMDomains(DomainID)%femdomainp &
                           %OversetConnect(i)%ElementID
               GaussPointID = this%FEMDomains(DomainID)%femdomainp &
                              %OversetConnect(i)%GaussPointID

               position = this%FEMDomains(DomainID)%femdomainp &
                          %OversetConnect(i)%position

               myDomainID = this%FEMDomains(DomainID)%femdomainp &
                            %OversetConnect(i)%DomainIDs12(1)
               n = size(this%FEMDomains(DomainID)%femdomainp &
                        %OversetConnect(i)%DomainIDs12)
               pairDomainID = this%FEMDomains(DomainID)%femdomainp &
                              %OversetConnect(i)%DomainIDs12(n)

               itr = itr + 1
               if (GaussPointID > 0) then
                  ! GPP

                  sf = this%FEMDomains(myDomainID)%femdomainp%mesh%getShapeFunction(ElementID, GaussPointID)

                  sf%ElementID = ElementID

                  N_I = this%FEMDomains(pairDomainID)%femdomainp%connectVector( &
                        position=position, &
                        DOF=DOF, &
                        shapefunction=sf)
                  Cvec = zeros( &
                         size(this%FEMDomains(DomainID)%femdomainp%OversetConnect(i)%InterConnect)*DOF &
                         )

                  do k = 1, size(this%FEMDomains(DomainID)%femdomainp &
                                 %OversetConnect(i)%InterConnect)
                     node_id = this%FEMDomains(DomainID)%femdomainp%OversetConnect(i)%InterConnect(k)
                     node_id = node_id + this%CRS_ID_Starts_From( &
                               this%FEMDomains(DomainID)%femdomainp%OversetConnect(i)%DomainIDs12(k) &
                               ) - 1
                     do n = 1, DOF
                        Cvec(DOF*(k - 1) + n) = X(DOF*(node_id - 1) + n)
                     end do
                  end do

                  gap_function(itr) = dot_product(N_I, Cvec)

               else
                  ! P2P
                  N_I = this%FEMDomains(pairDomainID)%femdomainp &
                        %connectVector(position=position, DOF=DOF)
                  Cvec = zeros( &
                         size(this%FEMDomains(DomainID)%femdomainp%OversetConnect(i)%InterConnect)*DOF &
                         )
                  do k = 1, size(this%FEMDomains(DomainID)%femdomainp &
                                 %OversetConnect(i)%InterConnect)
                     node_id = this%FEMDomains(DomainID)%femdomainp%OversetConnect(i)%InterConnect(k)
                     node_id = node_id + this%CRS_ID_Starts_From( &
                               this%FEMDomains(DomainID)%femdomainp%OversetConnect(i)%DomainIDs12(k) &
                               ) - 1
                     do n = 1, DOF
                        Cvec(DOF*(k - 1) + n) = X(DOF*(node_id - 1) + n)
                     end do
                  end do

                  gap_function(itr) = dot_product(N_I, Cvec)

               end if

            end if
         end do
      end do

   end function
! #####################################################
   function argumented_Lagrangian_RHS_FEMSolver(this, DOF, lambda) result(aL_RHS)
      class(FEMSolver_), intent(in) :: this
      real(real64), intent(in) :: lambda(:)
      real(real64), allocatable :: aL_RHS(:)

      integer(int32), intent(in) :: DOF

      real(real64), allocatable :: gap_function(:)

      ! compute g(x) = \Sigma [N]{x}
      integer(int32) :: DomainID, ElementID, GaussPointID, myDomainID, pairDomainID
      integer(int32) :: i, j, k, m, n, row_id, col_id, row_Domain_id, col_domain_id, itr, node_id
      real(real64) :: singleValue
      real(real64), allocatable :: N_I(:), position(:), Cvec(:)
      type(ShapeFunction_) :: sf

      aL_RHS = this%CRS_RHS
      aL_RHS(:) = 0.0d0

      ! do DomainID=1,2
      itr = 0
      do DomainID = 1, size(this%FEMDomains)
         do i = 1, this%FEMDomains(DomainID)%femdomainp &
            %numOversetElements()

            if (this%FEMDomains(DomainID)%femdomainp &
                %OversetConnect(i)%active) then

               ElementID = this%FEMDomains(DomainID)%femdomainp &
                           %OversetConnect(i)%ElementID
               GaussPointID = this%FEMDomains(DomainID)%femdomainp &
                              %OversetConnect(i)%GaussPointID

               position = this%FEMDomains(DomainID)%femdomainp &
                          %OversetConnect(i)%position

               myDomainID = this%FEMDomains(DomainID)%femdomainp &
                            %OversetConnect(i)%DomainIDs12(1)
               n = size(this%FEMDomains(DomainID)%femdomainp &
                        %OversetConnect(i)%DomainIDs12)
               pairDomainID = this%FEMDomains(DomainID)%femdomainp &
                              %OversetConnect(i)%DomainIDs12(n)

               itr = itr + 1
               if (GaussPointID > 0) then
                  ! GPP

                  sf = this%FEMDomains(myDomainID)%femdomainp%mesh%getShapeFunction(ElementID, GaussPointID)

                  sf%ElementID = ElementID

                  N_I = this%FEMDomains(pairDomainID)%femdomainp%connectVector( &
                        position=position, &
                        DOF=DOF, &
                        shapefunction=sf)

                  do k = 1, size(this%FEMDomains(DomainID)%femdomainp &
                                 %OversetConnect(i)%InterConnect)
                     node_id = this%FEMDomains(DomainID)%femdomainp%OversetConnect(i)%InterConnect(k)
                     node_id = node_id + this%CRS_ID_Starts_From( &
                               this%FEMDomains(DomainID)%femdomainp%OversetConnect(i)%DomainIDs12(k) &
                               ) - 1
                     do n = 1, DOF
                        aL_RHS(DOF*(node_id - 1) + n) = aL_RHS(DOF*(node_id - 1) + n) + lambda(i)*N_I(DOF*(k - 1) + n)
                     end do
                  end do

               else
                  ! P2P
                  N_I = this%FEMDomains(pairDomainID)%femdomainp &
                        %connectVector(position=position, DOF=DOF)

                  do k = 1, size(this%FEMDomains(DomainID)%femdomainp &
                                 %OversetConnect(i)%InterConnect)
                     node_id = this%FEMDomains(DomainID)%femdomainp%OversetConnect(i)%InterConnect(k)
                     node_id = node_id + this%CRS_ID_Starts_From( &
                               this%FEMDomains(DomainID)%femdomainp%OversetConnect(i)%DomainIDs12(k) &
                               ) - 1
                     do n = 1, DOF
                        aL_RHS(DOF*(node_id - 1) + n) = aL_RHS(DOF*(node_id - 1) + n) + lambda(i)*N_I(DOF*(k - 1) + n)
                     end do
                  end do

               end if

            end if
         end do
      end do

   end function
! #####################################################
   subroutine bicgstab_CRS_2(a, ptr_i, index_j, x, b, itrmax, er, relative_er, debug)
      integer(int32), intent(inout) :: index_j(:), itrmax
      integer(int64), intent(inout) :: ptr_i(:)
      real(real64), intent(inout) :: a(:), b(:), er
      real(real64), optional, intent(in) :: relative_er
      real(real64), intent(inout) :: x(:)
      logical, optional, intent(in) :: debug
      logical :: speak = .false.
      integer(int32) itr, i, j, n
      real(real64) alp, bet, c1, c2, c3, ev, vv, rr, er0, init_rr, re_er0
      real(real64), allocatable:: r(:), r0(:), p(:), y(:), e(:), v(:), pa(:), ax(:)

      er0 = er
      if (present(debug)) then
         speak = debug
      end if

      n = size(b)
      if (speak) print *, "BiCGSTAB STARTED >> DOF:", n
      allocate (r(n), r0(n), p(n), y(n), e(n), v(n))

      r(:) = b(:)
      if (speak) print *, "BiCGSTAB >> [1] initialize"

      call sub_crs_matvec(CRS_value=a, CRS_col=index_j, &
                          CRS_row_ptr=ptr_i, old_vector=x, new_vector=ax)
      r = b - ax

      if (speak) print *, "BiCGSTAB >> [2] dp1"

      c1 = dot_product(r, r)
      !call omp_dot_product(r,r,c1)

      init_rr = c1
      if (speak) print *, "BiCGSTAB >> [2] init_rr", c1
      !if(speak) print *, "BiCGSTAB >>      |r|^2 = ",init_rr

      !if (c1 < er0) return

      p(:) = r(:)
      r0(:) = r(:)

      do itr = 1, itrmax
         if (speak) print *, "BiCGSTAB >> ["//str(itr)//"] initialize"
         c1 = dot_product(r0, r)
         !call omp_dot_product(r0,r,c1)

         y(:) = 0.0d0
         call sub_crs_matvec(CRS_value=a, CRS_col=index_j, &
                             CRS_row_ptr=ptr_i, old_vector=p, new_vector=y)

         c2 = dot_product(r0, y)
         !call omp_dot_product(r0,y,c2)

         alp = c1/c2
         e(:) = r(:) - alp*y(:)
         v(:) = 0.0d0
         call sub_crs_matvec(CRS_value=a, CRS_col=index_j, &
                             CRS_row_ptr=ptr_i, old_vector=e, new_vector=v)

         if (speak) print *, "BiCGSTAB >> ["//str(itr)//"] half"

         ev = dot_product(e, v)
         vv = dot_product(v, v)
         !call omp_dot_product(e,v,ev)
         !call omp_dot_product(v,v,vv)

         if (vv == 0.0d0) stop "Bicgstab devide by zero"
         c3 = ev/vv
         if (speak) print *, "BiCGSTAB >> c3 = ev/vv", c3
         x(:) = x(:) + alp*p(:) + c3*e(:)
         r(:) = e(:) - c3*v(:)

         rr = dot_product(r, r)
         !call omp_dot_product(r,r,rr)

         if (itr == 1) then
            re_er0 = rr
         end if

         if (speak) then
            print *, sqrt(rr)
         end if

         if (present(relative_er)) then
            if (sqrt(rr/re_er0) < relative_er) then
               exit
            end if
         end if
         !    write(*,*) 'itr, er =', itr,rr
         if (sqrt(rr) < er0) exit

         c1 = dot_product(r0, r)
         !call omp_dot_product(r0,r,c1)

         bet = c1/(c2*c3)
         if ((c2*c3) == 0.0d0) stop "Bicgstab devide by zero"
         p(:) = r(:) + bet*(p(:) - c3*y(:))
      end do
   end subroutine
!===============================================================

! #####################################################
   subroutine MPI_BiCGSTABFEMSolver(this, x)
      class(FEMSolver_), intent(inout) :: this
      !integer(int32), intent(inout) :: ptr_i(:),index_j(:), itrmax
      !real(real64), intent(inout) :: a(:), b(:), er
      !real(real64),optional,intent(in) :: relative_er
      ! MPI
      !type(M_Link_Item_),intent(in) :: link_table(:)
      !type(MPI_) :: mpi_target

      real(real64), allocatable, intent(inout) :: x(:)

      logical :: speak = .false.
      integer(int32) itr, i, j, n
      real(real64) alp, bet, c1, c2, c3, ev, vv, rr, er0, init_rr, re_er0
      real(real64), allocatable:: r(:), r0(:), p(:), y(:), e(:), v(:), pa(:), ax(:)

      er0 = this%er0
      speak = this%debug
      if (.not. allocated(this%CRS_RHS)) then
         print *, "[ERROR] >> bicgstab_CRS_MPIFEMSolver :: detected .not.allocated(this%CRS_RHS)"
         stop
      end if
      n = size(this%CRS_RHS)
      if (speak) print *, "BiCGSTAB STARTED >> DOF:", n

      if (.not. allocated(x)) then
         x = zeros(n)
      end if

      allocate (r(n), r0(n), p(n), y(n), e(n), v(n))

      r(:) = this%CRS_RHS(:)

      if (speak) print *, "BiCGSTAB >> [1] initialize"

      ! matrix-vector multiplication
      !ax = crs_matvec(CRS_value=a,CRS_col=index_j,&
      !    CRS_row_ptr=ptr_i,old_vector=x)
      ax = this%MPI_matmul(b=x)
      r = this%CRS_RHS - ax

      if (speak) print *, "BiCGSTAB >> [2] dp1"

      ! dot_product
      c1 = this%mpi_dot_product(r, r)

      init_rr = c1
      !if(speak) print *, "BiCGSTAB >>      |r|^2 = ",init_rr

      if (c1 < er0) return

      p(:) = r(:)
      r0(:) = r(:)

      do itr = 1, this%itrmax
         if (speak) print *, "BiCGSTAB >> ["//str(itr)//"] initialize"
         c1 = this%mpi_dot_product(r0, r)

         !y = crs_matvec(CRS_value=a,CRS_col=index_j,&
         !CRS_row_ptr=ptr_i,old_vector=p)
         y = this%MPI_matmul(b=p)

         c2 = this%mpi_dot_product(r0, y)

         alp = c1/c2

         e(:) = r(:) - alp*y(:)

         !v = crs_matvec(CRS_value=a,CRS_col=index_j,&
         !CRS_row_ptr=ptr_i,old_vector=e)
         v = this%MPI_matmul(b=e)

         if (speak) print *, "BiCGSTAB >> ["//str(itr)//"] half"

         ev = this%mpi_dot_product(e, v)
         vv = this%mpi_dot_product(v, v)

         if (vv == 0.0d0) stop "Bicgstab devide by zero"

         c3 = ev/vv

         x(:) = x(:) + alp*p(:) + c3*e(:)

         r(:) = e(:) - c3*v(:)

         rr = this%mpi_dot_product(r, r)

         if (itr == 1) then
            re_er0 = rr
         end if

         if (speak) then
            print *, sqrt(rr)
         end if
         if (sqrt(rr/re_er0) < this%relative_er) then
            exit
         end if
         !    write(*,*) 'itr, er =', itr,rr
         if (sqrt(rr) < er0) exit
         c1 = this%mpi_dot_product(r0, r)

         bet = c1/(c2*c3)

         if ((c2*c3) == 0.0d0) stop "Bicgstab devide by zero"

         p(:) = r(:) + bet*(p(:) - c3*y(:))
      end do
   end subroutine
!===============================================================

!===================================================================================
   subroutine gauss_jordan_crs(a0_val, a0_row_ptr, a0_col_idx, x, b, n)
      integer(int32), intent(in) :: n
      real(real64), intent(in) :: a0_val(n), b(n)
      integer(int32), intent(in) ::  a0_col_idx(:)
      integer(int64), intent(in) :: a0_row_ptr(:)
      real(real64), intent(inout) :: x(n)
      real(real64) :: a_ik
      integer(int64) i, j, k, m, nn, mm

      integer(int32) :: i_1, i_2, i_3, i_4

      real(real64) ar, am, t, a_val(n), w(n)

      nn = size(b)

      ! Gauss-Jordan method with CRS format
      !a(:,:)= a0(:,:)
      a_val = a0_val

      !x(:) = b(:)
      x(:) = b(:)

      do k = 1, n
         m = k

         !am = abs(a(k,k))
         am = abs(getCRSval(a_val, a0_row_ptr, a0_col_idx, row=int(k), col=int(k)))

         !do i = k+1, n
         !    if (abs(a(i,k)) > am) then
         !        am = abs(a(i,k))
         !        m = i
         !    endif
         !enddo

         !do i = k+1, n
         !    a_ik = abs(getCRSval(a_val, a0_row_ptr, a0_col_idx, row=i, col=k))
         !    if (abs(a(i,k)) > am) then
         !        am = abs( a_ik )
         !        m = i
         !    endif
         !enddo

         if (am == 0.0d0) stop ' A is singular '

         print *, "Under construction :: gauss_jordan_crs"
         stop

!        if ( k /= m) then
!            w(k:n) = a(k, k:n)
!            a(k,k:n) = a(m, k:n) ! [*] a_val and pointers are to be reallocated.
!            a(m, k:n) =w(k:n)    ! [*] a_val and pointers are to be reallocated.
!            t = x(k)
!            x(k) = x(m)
!            x(m) = t
!        endif
!
!
!       ! gauss_jordan
!       if (a(k, k) == 0.0d0)  stop  'devide by zero3gauss_jordan_pv'
!          ar = 1.0d0 / a(k,k)
!       a(k,k) = 1.0d0
!       a(k,k+1:n) = ar * a(k, k+1:n) ! [*] a_val and pointers are to be reallocated.
!       x(k) = ar * x(k)
!       do i= 1, n
!         if (i /= k) then
!           a(i, k+1:n) = a(i, K+1:n) - a(i,k) * a(k, k+1:n) ! [*] a_val and pointers are to be reallocated.
!           x(i) = x(i) - a(i,k) * x(k)
!           a(i,k) = 0.0d0 ! [*] a_val and pointers are to be reallocated.
!         endif
!       enddo
      end do

   end subroutine
!===========================================================================

   function getCRSval(val, row_ptr, col_idx, row, col) result(ret)
      real(real64), intent(in) :: val(:)
      integer(int32), intent(in) :: col_idx(:), col, row
      integer(int64), intent(in) :: row_ptr(:)
      integer(int64) :: j, k
      real(real64) :: ret

      ret = 0.0d0

      do j = row_ptr(row), row_ptr(row + 1) - 1
         if (col_idx(j) == col) then
            ret = val(j)
         end if
      end do

   end function
!===========================================================================

!===========================================================================
   subroutine JacobiPreconditionerCRS(val, row_ptr, col_idx, rhs)
      real(real64), intent(inout) :: val(:), rhs(:)
      integer(int32), intent(inout) :: col_idx(:)
      integer(int64), intent(inout) :: row_ptr(:)
      real(real64) :: A_k_k, A_k_k_inv
      integer(int64) :: i, j

    !!$OMP parallel
    !!$OMP do
      do i = 1, size(row_ptr) - 1
         A_k_k = getCRSval(val=val, row_ptr=row_ptr, col_idx=col_idx, row=int(i), col=int(i))
         if (A_k_k == 0.0d0) cycle
         A_k_k_inv = 1.0d0/A_k_k
         rhs(i) = rhs(i)*A_k_k_inv
         do j = row_ptr(i), row_ptr(i + 1) - 1
            !!$OMP atomic
            val(j) = val(j)*A_k_k_inv
         end do
      end do
    !!$OMP end do
    !!$OMP end parallel
   end subroutine
!===========================================================================

!===========================================================================
!subroutine incompleteLUCRS(val,row_ptr,col_idx,rhs)
!    real(real64),intent(inout) :: val(:),rhs(:)
!    integer(int32),intent(inout) :: row_ptr(:),col_idx(:)
!    real(real64),allocatable :: diag_vec(:)
!    real(real64) :: A_k_k, A_k_k_inv,A_i_k, A_k_j
!    integer(int32) :: i,j,n,m,k,col,col_id,col_id_2,col_2
!    type(CRS_) :: crs
!    print *, "incompleteLUCRS >> bug exists"
!    stop
!
!    ! ILU(0)
!    call crs%init(val=val,row_ptr=row_ptr,col_idx=col_idx)
!    call crs%ILU(0,rhs=rhs)

!    do row=2,n
!        do col_id = row_ptr(row),row_ptr(row+1)-1
!            col = col_idx(col_id)
!            if(1<= col .and. col <= row - 1)then
!                val(col_id) = val(col_id) / diag_vec(row)
!                A_i_k = val(col_id)
!            endif
!
!            do col_id_2 = row_ptr(row),row_ptr(row+1)-1
!                col_2 = col_idx(col_id_2)
!                if(col + 1<= col_2 .and. col_2 <= row - 1)then
!                    val(col_id_2) = val(col_id_2) - A_i_k*A_k_j
!                endif
!            enddo
!        enddo
!    enddo

!    !$OMP parallel
!    !$OMP do
!    do k=1,size(row_ptr) - 1
!        A_k_k = getCRSval(val=val, row_ptr=row_ptr, col_idx=col_idx, row=k, col=k)
!        A_k_k_inv = 1.0d0/A_k_k
!
!        do i=1,size(row_ptr) - 1
!            do n=row_ptr(i), row_ptr(i+1)-1
!                j = col_idx(n)
!                ! a_ij := val(j)
!                if(j > k )then
!                    A_i_k=0.0d0
!                    do m=row_ptr(i), row_ptr(i+1)-1
!                        col = col_idx(m)
!                        if(col==k )then
!                            A_i_k = val(col)
!                        else
!                            cycle
!                        endif
!                    enddo
!                    if(A_i_k==0.0d0)then
!                        cycle
!                    endif
!                    !A_i_k = getCRSval(val=val, row_ptr=row_ptr, col_idx=col_idx, row=i, col=k)
!                    A_k_j = 0.0d0
!                    do m=row_ptr(k), row_ptr(k+1)-1
!                        col = col_idx(m)
!                        if(col==j )then
!                            A_k_j = val(col)
!                        else
!                            cycle
!                        endif
!                    enddo
!                    if(A_k_j==0.0d0)then
!                        cycle
!                    endif
!
!                    !A_k_j = getCRSval(val=val, row_ptr=row_ptr, col_idx=col_idx, row=k, col=j)
!                    !!$OMP atomic
!                    val(j) = val(j) - A_i_k*A_k_k_inv*A_k_j
!
!                endif
!            enddo
!        enddo
!    enddo
!    !$OMP end do
!    !$OMP end parallel
!end subroutine
!===========================================================================
!
!subroutine incompleteCholesky(val,row_ptr,col_idx,rhs)
!    real(real64),intent(inout) :: val(:),rhs(:)
!    integer(int32),intent(inout) :: row_ptr(:),col_idx(:)
!    real(real64) :: A_k_k, A_k_k_inv,A_i_k, A_k_j
!    integer(int32) :: i,j,n,m,k,col
!
!    n = size(a,1);
!
!        do k=1,n
!                !a(k,k) = sqrt(a(k,k))
!        a_k_k  = getCRSval(val=val, row_ptr=row_ptr, col_idx=col_idx, row=k, col=k)
!
!                do i=(k+1),n
!                    if (a(i,k)/=0)then
!                        a(i,k) = a(i,k)/a(k,k)
!                    endif
!                enddo
!
!                do j=(k+1),n
!                    do i=j,n
!                        if (a(i,j)/=0)then
!                            a(i,j) = a(i,j)-a(i,k)*a(j,k)
!                        endif
!                    enddo
!                enddo
!        enddo
!
!    do i=1,n
!        do j=i+1,n
!            a(i,j) = 0
!        enddo
!    enddo
!
!end subroutine
!

   subroutine reverseVectorReal64(A)
      real(real64), intent(inout) :: A(:)
      real(real64) :: buf
      integer(int32) :: i, n

      n = size(A)
      do i = 1, n/2
         buf = A(i)
         A(i) = A(n - i + 1)
         A(n - i + 1) = buf
      end do

   end subroutine

   subroutine reverseArrayReal64(A)
      real(real64), intent(inout) :: A(:, :)
      real(real64), allocatable :: buf(:)
      integer(int32) :: i, n

      ! see columns as vectors, and reverse order
      n = size(A, 1)
      buf = zeros(n)
      do i = 1, n
         call reverseVectorReal64(A(i, :))
      end do

   end subroutine

   function sortByIDreal64ColisVector(vectors, ID) result(new_vectors)
      real(real64), intent(in) :: vectors(:, :)
      integer(int32), intent(in) :: ID(:)
      real(real64), allocatable :: new_vectors(:, :)
      integer(int32) :: i
      new_vectors = zeros(size(vectors, 1), size(vectors, 2))
      do i = 1, size(ID)
         new_vectors(:, i) = vectors(:, ID(i))
      end do

   end function

   pure function getCRSFEMSolver(this, name) result(ret)
      class(FEMSolver_), intent(in) :: this
      character(*), optional, intent(in) :: name
      type(CRS_) :: ret
      integer(int32) :: i

      if (present(name)) then
!        select case(name)
!            case("A","a","K","StiffnessMatrix")
!                ret%col_idx = this%A_CRS_Index_Col
!                ret%row_ptr = this%A_CRS_Index_Row
!                ret%val     = this%A_CRS_val
!                return
!
!            case("B","b","M","MassMatrix")
!                ret%col_idx = this%B_CRS_Index_Col
!                ret%row_ptr = this%B_CRS_Index_Row
!                ret%val     = this%B_CRS_val
!                return
!
!            case default
         do i = 1, this%LAST_CRS_MATRIX_STACK
            if (this%CRS_matrix_stack_key%pages(i)%charvalue == name) then
               ret = this%CRS_matrix_stack(i)
               return
            else
               cycle
            end if
         end do
!                stop
!
!        end select
      else
         ret%col_idx = this%CRS_Index_Col
         ret%row_ptr = this%CRS_Index_Row
         ret%val = this%CRS_val
      end if
   end function

   subroutine MPI_linkFEMSolver(this, rank_and_rowID_1, rank_and_rowID_2)
      class(FEMSolver_), intent(inout) :: this
      integer(int32), intent(in) :: rank_and_rowID_1(1:2), rank_and_rowID_2(1:2)
      type(M_Link_Item_), allocatable :: buf_table(:)
      integer(int32) :: i

      ! if not allocated, allocate
      if (.not. allocated(this%Link_Table)) then
         allocate (this%Link_Table(this%LINK_TABLE_INIT_SIZE))
         this%Link_num = 0
      end if

      ! if not related, return
      if (rank_and_rowID_1(1) /= this%MPI_target%myrank &
          .and. rank_and_rowID_2(1) /= this%MPI_target%myrank) then
         return
      end if

      ! related. but table is full
      if (this%Link_num >= size(this%Link_Table)) then
         buf_table = this%Link_Table
         deallocate (this%Link_Table)
         allocate (this%Link_Table(2*size(buf_table)))
         do i = 1, size(buf_table)
            this%Link_Table(i) = buf_table(i)
         end do
      end if

      ! related and buffer has enough space
      this%Link_num = this%Link_num + 1

      if (rank_and_rowID_1(1) == this%mpi_target%myrank) then
         this%Link_Table(this%Link_num)%rank_and_rowID_1 = rank_and_rowID_1
         this%Link_Table(this%Link_num)%rank_and_rowID_2 = rank_and_rowID_2
      elseif (rank_and_rowID_2(1) == this%mpi_target%myrank) then
         this%Link_Table(this%Link_num)%rank_and_rowID_1 = rank_and_rowID_2
         this%Link_Table(this%Link_num)%rank_and_rowID_2 = rank_and_rowID_1
      else
         ! do nothing
         return
      end if

   end subroutine

   function MPI_dot_productFEMSolver(this, a, b) result(dp)
      class(FEMSolver_), intent(in) :: this
      real(real64), intent(in) :: a(:), b(:)
      integer(int32), allocatable :: num_link(:)
      integer(int32) :: i, j, n, m
      real(real64) :: dp, my_dp, sendobj(1), recvobj(1)

      if (.not. associated(this%MPI_target)) then
         print *, "ERROR ::MPI_dot_productFEMSolver >> please associate this%mpi_target "
         print *, "this%mpi_target => your_mpi_object"
         stop
      end if

      n = size(a)
      num_link = int(eyes(n))
      if (allocated(this%Link_Table)) then
         do i = 1, this%Link_num

            m = this%Link_Table(i)%rank_and_rowID_1(2)
            if (this%Link_Table(i)%rank_and_rowID_2(1) > this%MPI_target%petot - 1) cycle
            if (m > n) cycle
            num_link(m) = num_link(m) + 1
         end do
      end if

      my_dp = 0.0d0
      !$OMP parallel
      !$OMP do reduction(+:my_dp)
      do i = 1, n
         my_dp = my_dp + a(i)*b(i)/dble(num_link(i))
      end do
      !$OMP end do
      !$OMP end parallel

      sendobj(1) = my_dp
      recvobj(1) = 0.0d0
      call this%MPI_target%AllReduce(sendobj=sendobj, recvobj=recvobj, count=1, sum=.true.)
      dp = recvobj(1)

   end function

   function MPI_matmulFEMSolver(this, A, b) result(my_c)
      class(FEMSolver_), intent(in) :: this
      type(CRS_), optional, intent(in) :: A
      real(real64), intent(in) :: b(:)
      real(real64), allocatable :: sendbuf(:), recvbuf(:), my_c(:)
      integer(int32), allocatable :: send_recv_rank(:)
      integer(int32) :: i, j, n, cor_rank, my_row_id

      if (.not. associated(this%MPI_target)) then
         print *, "ERROR ::MPI_matmulFEMSolver >> please associate this%mpi_target "
         print *, "this%mpi_target => your_mpi_object"
         stop
      end if

!    ! create sendbuf and recvbuf
!    sendbuf  = zeros(this%Link_num)
!    send_recv_rank = int(zeros(this%Link_num))
!    do i=1,this%Link_num
!        my_row_id  = this%Link_Table(i)%rank_and_rowID_1(2)
!        cor_rank   = this%Link_Table(i)%rank_and_rowID_2(1)
!        sendbuf(i)  = b_copy(my_row_id)
!        send_recv_rank(i) = cor_rank
!    enddo
!    recvbuf  = zeros(this%Link_num)
!
!    ! MPI COMMUNICATION
!    ! ISEND :: >>> NON-BLOCKING
!    call this%mpi_target%isend_irecv(sendobj=sendbuf,recvobj=recvbuf,send_recv_rank=send_recv_rank)
!
!
!    do i=1,this%Link_num
!        my_row_id  = this%Link_Table(i)%rank_and_rowID_1(2)
!        cor_rank   = this%Link_Table(i)%rank_and_rowID_2(1)
!        b_copy(my_row_id) = b_copy(my_row_id) + recvbuf(i)
!    enddo
      if (present(A)) then
         my_c = A%matmul(b)
      else
         my_c = zeros(size(b))
         call sub_crs_matvec(CRS_value=this%CRS_val, CRS_col=this%CRS_Index_Col, &
                             CRS_row_ptr=this%CRS_Index_Row, old_vector=b, new_vector=my_c)
      end if

      ! create sendbuf and recvbuf

      sendbuf = zeros(this%Link_num)
      send_recv_rank = int(zeros(this%Link_num))
      recvbuf = zeros(this%Link_num)

      do i = 1, this%Link_num
         my_row_id = this%Link_Table(i)%rank_and_rowID_1(2)
         cor_rank = this%Link_Table(i)%rank_and_rowID_2(1)
         sendbuf(i) = my_c(my_row_id)
         send_recv_rank(i) = cor_rank
      end do
      recvbuf = zeros(this%Link_num)

      call this%mpi_target%isend_irecv(sendobj=sendbuf, recvobj=recvbuf, send_recv_rank=send_recv_rank)

      do i = 1, this%Link_num
         my_row_id = this%Link_Table(i)%rank_and_rowID_1(2)
         cor_rank = this%Link_Table(i)%rank_and_rowID_2(1)
         my_c(my_row_id) = my_c(my_row_id) + recvbuf(i)
      end do

   end function

! ################################################
   function conditionNumberFEMSolver(this) result(RCOND)
      class(FEMSolver_), intent(in) :: this
      integer(int32) :: N, LDA, INFO
      real(real64), allocatable :: A(:, :)
      real(real64) :: ANORM, RCOND
      real(real64), allocatable :: WORK(:)
      integer(int32), allocatable :: IWORK(:)
      character(1) :: NORM
      type(CRS_) :: crs

      N = size(this%CRS_Index_Row) - 1
      LDA = N
      CRS = this%getCRS()
      A = CRS%to_Dense()
      NORM = "1"
      ANORM = sqrt(dot_product(this%CRS_val, this%CRS_val))

      allocate (WORK(4*N), IWORK(N))
      call dgecon(NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)

      if (INFO < 0) then
         print *, "ERROR >> DEGCON failed."
      end if

   end function
! ################################################

! ################################################
   function get_fix_idx_FEMSolver(this) result(fix_idx)
      class(FEMSolver_), intent(in) :: this
      integer(int32), allocatable :: fix_idx(:)

      integer(int32) :: i, num_fix_idx

      if (.not. allocated(this%fix_lin_exists)) then
         allocate (fix_idx(0))
         return
      end if

      num_fix_idx = 0
      do i = 1, size(this%fix_lin_exists)
         if (this%fix_lin_exists(i)) then
            num_fix_idx = num_fix_idx + 1
         end if
      end do
      fix_idx = int(zeros(num_fix_idx))
      num_fix_idx = 0
      do i = 1, size(this%fix_lin_exists)
         if (this%fix_lin_exists(i)) then
            num_fix_idx = num_fix_idx + 1
            fix_idx(num_fix_idx) = i
         end if
      end do

   end function
! ################################################

! ################################################
   function get_fix_value_FEMSolver(this) result(fix_value)
      class(FEMSolver_), intent(in) :: this
      integer(int32), allocatable :: fix_value(:)

      integer(int32) :: i, num_fix_value

      if (.not. allocated(this%fix_lin_exists)) then
         allocate (fix_value(0))
         return
      end if

      num_fix_value = 0
      do i = 1, size(this%fix_lin_exists)
         if (this%fix_lin_exists(i)) then
            num_fix_value = num_fix_value + 1
         end if
      end do
      fix_value = zeros(num_fix_value)

      num_fix_value = 0
      do i = 1, size(this%fix_lin_exists)
         if (this%fix_lin_exists(i)) then
            num_fix_value = num_fix_value + 1
            fix_value(num_fix_value) = this%fix_lin_exists_Values(i)
         end if
      end do

   end function
! ################################################

   subroutine importCRSMatrix_FEMSolver(this, CRS)
      class(FEMSolver_), intent(inout) :: this
      type(CRS_), intent(in) :: CRS

      this%CRS_val = CRS%val
      this%CRS_Index_Col = CRS%col_idx
      this%CRS_Index_Row = CRS%row_ptr

   end subroutine

end module
