module SparseClass
   use iso_c_binding
   use ArrayClass
   use MathClass
   use RandomClass
   use RangeClass
   implicit none

!> Please edit "src/SparseClass/SparseClass.c"
!> for
!>  'c_dot_product'
!>  'crs_spmv_real32'
!>  'crs_spmv_real64'

   interface
      Subroutine c_dot_product(a, b, n, ret) bind(C, Name='c_dot_product')
         import
         integer(C_size_t), value :: n
         real(c_double), intent(in) :: a(n), b(n)
         real(c_double), intent(out) :: ret(1)
      End Subroutine
   End Interface

   interface
      Subroutine crs_spmv_real32(val, row_ptr, col_idx, old_vector, new_vector, n, col_size) &
         Bind(C, Name="crs_spmv_real32")
         import
         Integer(C_size_t), Value :: n, col_size
         real(C_float), Intent(In) :: val(n)
         integer(C_int64_t), Intent(In) :: row_ptr(n + 1)
         integer(C_int), Intent(In) :: col_idx(col_size)
         real(C_float), Intent(In)  :: old_vector(n)
         real(C_float), Intent(InOut) :: new_vector(n)

      End Subroutine
   End interface

   interface
      Subroutine crs_spmv_real64(val, row_ptr, col_idx, old_vector, new_vector, n, col_size) &
         Bind(C, Name="crs_spmv_real64")
         use iso_c_binding
         import
         Integer(C_size_t), Value :: n, col_size
         real(C_double), Intent(In) :: val(n)
         integer(C_int64_t), Intent(In) :: row_ptr(n + 1)
         integer(C_int), Intent(In) :: col_idx(col_size)
         real(C_double), Intent(In)  :: old_vector(n)
         real(C_double), Intent(InOut) :: new_vector(n)

      End Subroutine
   End interface

!  interface
!  Subroutine crs_spmv_real64_opencl(val,row_ptr,col_idx,old_vector,new_vector,n,col_size) &
!        Bind(C,Name="crs_spmv_real64_opencl")
!        use iso_c_binding
!    import
!    Integer(C_size_t),Value :: n,col_size
!    real(C_double),Intent(In) :: val(n)
!    integer(C_int64_t),Intent(In) :: row_ptr(n+1)
!    integer(C_int),Intent(In) :: col_idx(col_size)
!    real(C_double),Intent(In)  :: old_vector(n)
!    real(C_double),Intent(InOut) :: new_vector(n)
!
!  End Subroutine
!End interface

   interface
      Subroutine c_sparse_matvec(row_ptr, col_idx, val, x, n, n_col, ret) bind(C, Name='c_sparse_matvec')
         import
         integer(C_size_t), value :: n, n_col
         integer(c_int), intent(in) :: row_ptr(n), col_idx(n_col)
         real(c_double), intent(in) :: val(n - 1), x(n - 1)
         real(c_double), intent(out) :: ret(n - 1)
      End Subroutine
   End Interface

   interface sinc
      module procedure sinc_complex64, sinc_real64, sinc_real64_vec
   end interface

   type :: COO_Row_
      real(real64), allocatable :: val(:)
      complex(real64), allocatable :: val_complex64(:)
      integer(int32), allocatable :: col(:)
   end type

   type :: COO_
      type(COO_Row_), allocatable :: row(:)
   contains
      procedure, public :: init => initCOO
      procedure, public :: update => updateCOO
      procedure, public :: set => updateCOO

      procedure, pass :: addCOO
      procedure, pass :: add_complexCOO
      generic :: add => addCOO, add_complexCOO

      procedure, public :: getDenseMatrix => getDenseMatrixCOO
      procedure, public :: to_dense => getDenseMatrixCOO
      procedure, public :: remove => removeCOO
      procedure, public :: getAllCol => getAllColCOO
      procedure, public :: DOF => DOFCOO
      
      procedure, public :: to_CRS => to_CRSCOO
      procedure, public :: get => getCOO
      procedure, public :: ne => neCOO
      procedure, public :: maxval => maxvalCOO
      procedure, public :: random => randomCOO ! create random sparse matrix

      ! typical matrix
      procedure, public :: eyes => eyesCOO
      procedure, public :: poisson => poissonCOO

      !procedure,public ::getAllCol_as_row_obj => getAllCol_as_row_objCOO
   end type

   type :: CCS_
      integer(int32), allocatable :: col_ptr(:)
      integer(int32), allocatable :: row_idx(:)
      real(real64), allocatable :: val(:)
   contains
      procedure, public :: get_column => get_column_CCS
   end type

   type :: CRS_
      integer(int32), allocatable :: col_idx(:)

      ! Destructive
      integer(int64), allocatable :: row_ptr(:)
      real(real64), allocatable :: val(:)

      real(real32), allocatable :: val_real32(:)
      complex(real64), allocatable :: val_complex64(:)
      integer(int32) :: dtype = real64

   contains
      procedure, public :: init => initCRS
      procedure, public :: eyes => eyesCRS
      procedure, public :: poisson => poissonCRS
      procedure, public :: Lanczos => LanczosCRS

      procedure, pass :: matmulCRS
      !procedure,pass :: matmul_opencl_CRS

      procedure, pass :: matmul_c_CRS
      procedure, pass :: matmul_c_real32_CRS

      procedure, pass :: matmul_real32_CRS
      procedure, pass :: matmul_complex_CRS
      generic :: matmul => matmulCRS, matmul_complex_CRS, matmul_real32_CRS, matmul_c_CRS, &
         matmul_c_real32_CRS!,matmul_opencl_CRS

      procedure, public :: to_real32 => to_real32_CRS

      procedure, public :: SpMV => spmv_as_subroutine_CRS
      procedure, public :: eig => eigCRS
      procedure, public :: to_dense => to_denseCRS
      procedure, public :: DOF => DOFCRS

      procedure, public :: remove => removeCRS
      procedure, public :: size => sizeCRS
      procedure, public :: shape => shapeCRS
      procedure, public :: update => updateCRS
      procedure, public :: add => addCRS

      ! 対象行列とみなし，埋まっていない要素をすべて埋める．
      procedure, public :: fillSymmetric => fillSymmetric_CRS

      procedure, public :: get => getCRS
      procedure, public :: is_nonzero => is_nonzeroCRS
      procedure, public :: diag => diagCRS
      procedure, public :: lumped => lumpedCRS

      procedure, public :: divide_by => divide_by_CRS
      ! subroutine version
      procedure, public :: divide_by_vector => divide_by_vector_CRS
      procedure, public :: mult_by => mult_by_CRS

      procedure, public :: to_CCS => to_CCSCRS

      procedure, public :: load => loadCRS

      procedure, public :: ILU => ILUCRS
      procedure, public :: ILU_matvec => ILU_matvecCRS
      procedure, public :: BICGSTAB => BICGSTAB_CRSSparse

      procedure, pass :: tensor_exponential_crs
      procedure, pass :: tensor_sqrt_crs
      procedure, pass :: tensor_exp_sqrt_crs
      procedure, pass :: tensor_log_crs
      procedure, public :: tensor_log_crs_modified => tensor_log_crs_modified_Sparse

      procedure, pass :: fixCRS
      procedure, pass :: fix_zeroCRS

      procedure, pass :: tensor_exponential_complex64_crs
      procedure, pass :: tensor_exp_sqrt_complex64_crs

      procedure, pass :: tensor_cos_sqrt_complex64_crs
      procedure, pass :: tensor_cos_sqrt_real64_crs
      procedure, pass :: tensor_cos_sqrt_LPF_real64_crs

      procedure, pass :: tensor_sinc_sqrt_complex64_crs
      procedure, pass :: tensor_sinc_sqrt_real64_crs
      procedure, pass :: tensor_t_sinc_sqrt_LPF_real64_crs

      procedure, pass :: tensor_sqrt_complex64_crs
      procedure, pass :: tensor_log_complex64_crs
      procedure, pass :: fix_complex64_CRS
      procedure, pass :: tensor_d1_wave_kernel_complex64_crs
      procedure, pass :: tensor_wave_kernel_complex_64_crs
      procedure, pass :: tensor_wave_kernel_real_64_crs
      procedure, pass :: tensor_wave_kernel_LPF_real_64_crs
      procedure, pass :: tensor_wave_kernel_RHS_real_64_crs
      procedure, pass :: tensor_wave_kernel_RHS_LPF_real_64_crs
      procedure, pass :: tensor_wave_kernel_RHS_complex_64_crs

      procedure, pass :: tensor_cos_sqrt_cos_sqrt_complex64_crs
      procedure, pass :: tensor_cos_sqrt_sinc_sqrt_complex64_crs
      procedure, pass :: tensor_sinc_sqrt_cos_sqrt_complex64_crs
      procedure, pass :: tensor_sinc_sqrt_sinc_sqrt_complex64_crs
      
      generic, public :: tensor_exponential => tensor_exponential_complex64_crs, tensor_exponential_crs
      generic, public :: tensor_exp => tensor_exponential_complex64_crs, tensor_exponential_crs
      generic, public :: tensor_exp_sqrt => tensor_exp_sqrt_complex64_crs, tensor_exp_sqrt_crs

      generic, public :: tensor_cos_sqrt => tensor_cos_sqrt_complex64_crs, &
         tensor_cos_sqrt_real64_crs
      generic, public :: tensor_sinc_sqrt => tensor_sinc_sqrt_complex64_crs, &
         tensor_sinc_sqrt_real64_crs

      generic, public :: tensor_cos_sqrt_LPF => tensor_cos_sqrt_LPF_real64_crs
      generic, public :: tensor_t_sinc_sqrt_LPF => tensor_t_sinc_sqrt_LPF_real64_crs

      generic, public :: tensor_sqrt => tensor_sqrt_complex64_crs, tensor_sqrt_crs
      generic, public :: tensor_log => tensor_log_complex64_crs, tensor_log_crs

      !>>> not verified.
      generic, public :: tensor_cos_sqrt_cos_sqrt => tensor_cos_sqrt_cos_sqrt_complex64_crs
      generic, public :: tensor_cos_sqrt_sinc_sqrt => tensor_cos_sqrt_sinc_sqrt_complex64_crs
      generic, public :: tensor_sinc_sqrt_cos_sqrt => tensor_sinc_sqrt_cos_sqrt_complex64_crs
      generic, public :: tensor_sinc_sqrt_sinc_sqrt => tensor_sinc_sqrt_sinc_sqrt_complex64_crs
      !<<< not verified.

      generic, public :: tensor_wave_kernel => tensor_wave_kernel_complex_64_crs, &
         tensor_wave_kernel_real_64_crs
      generic, public :: tensor_wave_kernel_LPF => tensor_wave_kernel_LPF_real_64_crs

      generic, public :: tensor_wave_kernel_RHS => tensor_wave_kernel_RHS_real_64_crs, &
         tensor_wave_kernel_RHS_complex_64_crs, tensor_wave_kernel_RHS_LPF_real_64_crs

      generic, public :: tensor_d1_wave_kernel => tensor_d1_wave_kernel_complex64_crs
      generic, public :: fix => fix_complex64_CRS, fixCRS, fix_zeroCRS
         
      ! assemble large matrix by combining small matrices
      ! >> If you want to do this, use type(BCRS_), which is for block CRS format.
      !procedure,public :: assemble => assemble_CRS <- this is/will not be implemented.
   end type

   type :: BCRS_
      type(CRS_),allocatable :: CRS(:,:)
   contains
      procedure,public :: set       => setBCRS
      procedure,public :: add       => addBCRS
      procedure,public :: showShape => showShapeBCRS
      procedure,public :: shape     => shapeBCRS
      procedure,public :: row_range => row_range_BCRS
      procedure,public :: col_range => col_range_BCRS
      procedure,public :: matmul    => matmulBCRS
      procedure,public :: to_dense  => to_dense_BCRS
      procedure,public :: exp       => expBCRS
   end type 

   public :: operator(+)
   public :: operator(-)
   public :: operator(*)
   !public :: operator(/)

   interface operator(+)
      module procedure addCRS_and_CRS
   end interface

   interface operator(-)
      module procedure diffCRS_and_CRS
   end interface

   interface LAPACK_EIG
      module procedure LAPACK_EIG_DENSE, LAPACK_EIG_SPARSE
   end interface

   interface operator(*)
      module procedure multReal64_and_CRS, multCRS_and_Real64, multReal64_and_BCRS, multBCRS_and_Real64
   end interface

   interface matmul
      module procedure :: matmul_CRS_CRS, matmul_CRS_vec
   end interface

   interface LOBPCG
      module procedure LOBPCG_CRS, LOBPCG_SINGLE_CRS
   end interface LOBPCG

   interface to_diag
      module procedure to_diag_vector_to_CRS
   end interface

   interface allocated
      module procedure allocated_CRS
   end interface

   interface speye
      module procedure speyeCOO
   end interface

   interface imag_I
      module procedure imag_ICOO
   end interface

!    interface operator(/)
!        module procedure divide_by_diagonal_vector_CRS, divide_by_scalar_CRS
!    end interface

   interface to_COO
      module procedure to_COO_from_DenseMatrix, to_COO_from_ArrayObject
   end interface

   interface to_CRS
      module procedure to_CRS_from_DenseMatrix, to_CRS_from_ArrayObject
   end interface

   type :: Sparse_
      type(COO_) :: COO
      type(CRS_) :: CRS
      type(CCS_) :: CCS
   end type
contains

   subroutine initCOO(this, num_row)
      class(COO_), intent(inout) :: this
      integer(int32), intent(in) :: num_row
      allocate (this%row(num_row))
   end subroutine

   function to_CRSCOO(this, remove_coo) result(CRS_version)
      class(COO_), intent(inout) :: this
      logical, optional, intent(in) :: remove_coo
      type(CRS_) :: CRS_version
      integer(int32) :: i, idx
      integer(int64) :: n, j

      CRS_version%col_idx = this%getAllCol()

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

         if (present(remove_coo)) then
            if (remove_coo) then
               deallocate (this%row(i - 1)%col)
            end if
         end if

      end do

      !CRS_version%val
      n = size(CRS_version%col_idx)
      allocate (CRS_version%val(n))
      CRS_version%val(:) = 0.0d0

      do i = 1, size(this%row)
         if (.not. allocated(this%row(i)%val)) then
            cycle
         end if
         idx = 0
         do j = CRS_version%row_ptr(i), CRS_version%row_ptr(i + 1) - 1
            idx = idx + 1
            CRS_version%val(j) = this%row(i)%val(idx)
         end do

         if (present(remove_coo)) then
            if (remove_coo) then
               deallocate (this%row(i)%val)
            end if
         end if

      end do

      call this%remove()

   end function

   subroutine updateCOO(this, row, col, val)
      class(COO_), intent(inout) :: this
      integer(int32), intent(in)   :: row
      integer(int32), intent(in)   :: col
      real(real64), intent(in)   :: val
      integer(int32) :: i, col_id

      if (row > this%DOF() .or. row < 1) return
      if (col > this%DOF() .or. col < 1) return

      if (.not. allocated(this%row)) then
         print *, "ERROR :: initCOO"
         print *, "Please call [COO]%init() before this operation."
         return
      end if

      if (.not. allocated(this%row(row)%val)) then
         this%row(row)%val = [val]
         this%row(row)%col = [col]
      else

         ! check duplication
         if (minval(abs(this%row(row)%col(:) - col)) == 0) then
            ! duplication
            do i = 1, size(this%row(row)%col)
               if (this%row(row)%col(i) == col) then
                  col_id = i
                  exit
               end if
            end do
            this%row(row)%val(col_id) = val
         else
            this%row(row)%val = this%row(row)%val//[val]
            this%row(row)%col = this%row(row)%col//[col]
         end if
      end if

   end subroutine

   subroutine addCOO(this, row, col, val)
      class(COO_), intent(inout) :: this
      type(COO_) :: buf
      integer(int32), intent(in)   :: row
      integer(int32), intent(in)   :: col
      real(real64), intent(in)   :: val
      integer(int32) :: i, col_id,n

      !if (row > this%DOF() .or. row < 1) return
      !if (col > this%DOF() .or. col < 1) return

      if ( row < 1) return
      if (col < 1) return

      if(row > this%DOF())then
         ! add row
         buf = this
         n = this%DOF()
         call this%remove()
         call this%init(row)
         this%row(1:n) = buf%row(1:n)
      endif

      if (.not. allocated(this%row)) then
         print *, "ERROR :: initCOO"
         print *, "Please call [COO]%init() before this operation."
         return
      end if

      if (.not. allocated(this%row(row)%val)) then
         this%row(row)%val = [val]
         this%row(row)%col = [col]
      else

         ! check duplication
         if (minval(abs(this%row(row)%col(:) - col)) == 0) then
            ! duplication
            do i = 1, size(this%row(row)%col)
               if (this%row(row)%col(i) == col) then
                  col_id = i
                  exit
               end if
            end do
            this%row(row)%val(col_id) = this%row(row)%val(col_id) + val
         else
            this%row(row)%val = this%row(row)%val//[val]
            this%row(row)%col = this%row(row)%col//[col]
         end if
      end if

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

! ###############################################
   subroutine add_complexCOO(this, row, col, val)
      class(COO_), intent(inout) :: this
      integer(int32), intent(in)   :: row
      integer(int32), intent(in)   :: col
      complex(real64), intent(in)   :: val
      integer(int32) :: i, col_id

      if (row > this%DOF() .or. row < 1) return
      if (col > this%DOF() .or. col < 1) return

      if (.not. allocated(this%row)) then
         print *, "ERROR :: initCOO"
         print *, "Please call [COO]%init() before this operation."
         return
      end if

      if (.not. allocated(this%row(row)%val_complex64)) then
         this%row(row)%val_complex64 = [val]
         this%row(row)%col = [col]
      else

         ! check duplication
         if (minval(abs(this%row(row)%col(:) - col)) == 0) then
            ! duplication
            do i = 1, size(this%row(row)%col)
               if (this%row(row)%col(i) == col) then
                  col_id = i
                  exit
               end if
            end do
            this%row(row)%val_complex64(col_id) = this%row(row)%val_complex64(col_id) + val
         else
            this%row(row)%val_complex64 = this%row(row)%val_complex64//[val]
            this%row(row)%col = this%row(row)%col//[col]
         end if
      end if

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

! ###############################################
   function getDenseMatrixCOO(this) result(dense_matrix)
      class(COO_), intent(in) :: this
      real(real64), allocatable :: dense_matrix(:, :)
      integer(int32) :: i, j

      dense_matrix = zeros(size(this%row), size(this%row))

      do i = 1, size(this%row)
         print *, i
         if (allocated(this%row(i)%col)) then
            do j = 1, size(this%row(i)%col)
               dense_matrix(i, this%row(i)%col(j)) = this%row(i)%val(j)
            end do
         end if
      end do

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

! ###############################################
   subroutine removeCOO(this)
      class(COO_), intent(inout) :: this

      if (allocated(this%row)) deallocate (this%row)

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

   function getAllColCOO(this) result(cols)
      class(COO_), intent(in) :: this
      type(COO_ROW_) :: row
      integer(int32), allocatable :: cols(:)
      integer(int32) :: num_total_cols, i

      num_total_cols = 0
      do i = 1, size(this%row)
         if (allocated(this%row(i)%col)) then
            num_total_cols = num_total_cols + size(this%row(i)%col)
         else
            cycle
         end if
      end do
      allocate (cols(num_total_cols))
      cols(:) = 0
      num_total_cols = 0
      do i = 1, size(this%row)
         if (allocated(this%row(i)%col)) then
            cols(num_total_cols + 1:num_total_cols + size(this%row(i)%col)) &
               = this%row(i)%col(1:size(this%row(i)%col))
            num_total_cols = num_total_cols + size(this%row(i)%col)
         else
            cycle
         end if
      end do

   end function

!recursive function getAllCol_as_row_objCOO(this,row_old) result(row_new)
!    class(COO_),intent(in) :: this
!    type(COO_Row_),optional,intent(in) :: row_old(:)
!    type(COO_Row_) :: row_new
!    type(COO_Row_),allocatable :: row_new_buf(:)
!
!    integer(int32) :: num_row,half_row,num_elem
!
!    print *, "hello"
!    if(.not.present(row_old) )then
!        print *, "hello"
!        row_new = this%getAllCol_as_row_obj(this%row)
!
!    else
!        if(size(row_old) ==1 )then
!            if(allocated(row_old(1)%col) )then
!                row_new%col = row_old(1)%col
!            else
!                allocate(row_new%col(0) )
!            endif
!        elseif(size(row_old) ==2 )then
!            if(allocated(row_old(1)%col) .and. allocated(row_old(1)%col) )then
!                row_new%col = row_old(1)%col
!                row_new%col = row_new%col // row_old(2)%col
!            elseif(allocated(row_old(1)%col))then
!                row_new%col = row_old(1)%col
!            elseif(allocated(row_old(2)%col))then
!                row_new%col = row_old(2)%Col
!            else
!                allocate(row_new%col(0) )
!            endif
!        else
!            num_row = size(row_old)
!            half_row = num_row/2
!            print *, half_row, num_row
!            allocate(row_new_buf(2) )
!            row_new_buf(1) = this%getAllCol_as_row_obj(this%row(1:half_row) )
!            row_new_buf(2) = this%getAllCol_as_row_obj(this%row(half_row:) )
!            row_new%col = row_new_buf(1)%col // row_new_buf(2)%col
!        endif
!    endif
!
!end function
!
   subroutine LanczosCRS(this, DiagonalVector, subDiagonalVector, V)
      class(CRS_), intent(in) :: this

      real(real64), allocatable, intent(inout) :: V(:, :)!,T(:,:)
      real(real64), allocatable, intent(inout) :: DiagonalVector(:)
      real(real64), allocatable, intent(inout) :: subDiagonalVector(:)
      real(real64), allocatable :: v1(:), w1(:), vj(:), wj(:)
      real(real64) :: alp1, beta_j, alp_j
      type(Random_) :: random
      integer(int32) :: i, N, j, m

      ! https://en.wikipedia.org/wiki/Lanczos_algorithm

      !N = size(A,1)
      N = size(this%row_ptr) - 1
      V = zeros(N, N)
      !T = zeros(N,N)
      DiagonalVector = zeros(N)
      subDiagonalVector = zeros(N - 1)

      v1 = 2.0d0*random%randn(N) - 1.0d0
      v1 = v1/norm(v1)
      w1 = this%matmul(v1)
      alp1 = dot_product(w1, v1)
      w1 = w1 - alp1*v1
      beta_j = norm(w1)

      !T(1,1) = alp1
      V(:, 1) = v1
      DiagonalVector(1) = alp1

      do i = 2, N
         beta_j = norm(w1)
         if (beta_j /= 0.0d0) then
            vj = w1/beta_j
         else
            vj = 2.0d0*random%randn(N) - 1.0d0
            stop
         end if
         wj = this%matmul(vj)
         alp_j = dot_product(wj, vj)
         wj = wj - alp_j*vj - beta_j*v1
         v1 = vj
         w1 = wj

         !T(i,i) = alp_j
         !T(i,i-1) = beta_j
         !T(i-1,i) = beta_j
         V(:, i) = vj
         DiagonalVector(i) = alp_j
         subDiagonalVector(i - 1) = beta_j
      end do

   end subroutine

   function matmulCRS(CRS, old_vector, run_mode, cache_size) result(new_vector)
      class(CRS_), intent(in) :: CRS
      real(real64), intent(in)  :: Old_vector(:)
      real(real64), allocatable :: new_vector(:)
      integer(int32), optional, intent(in) :: run_mode
      integer(int32) :: mode_id
      integer(int32), optional, intent(in) :: cache_size

      mode_id = input(default=0, option=run_mode)
      if (mode_id == 1) then
         ! elemental mode
         new_vector = crs_matvec_generic_elemental_SparseClass( &
                      CRS_value=CRS%val, &
                      CRS_col=CRS%col_idx, &
                      CRS_row_ptr=CRS%row_ptr, &
                      old_vector=old_vector)
      else
         new_vector = crs_matvec_generic_SparseClass( &
                      CRS_value=CRS%val, &
                      CRS_col=CRS%col_idx, &
                      CRS_row_ptr=CRS%row_ptr, &
                      old_vector=old_vector, &
                      cache_size=cache_size)
      end if
   end function

   function matmul_c_CRS(CRS, old_vector, c_routine) result(new_vector)

      class(CRS_), intent(in) :: CRS
      real(real64), intent(in)  :: Old_vector(:)
      real(real64), allocatable :: new_vector(:)
      character(*), intent(in) :: c_routine
      integer(int64) :: n, col_size

      n = crs%size()
      col_size = size(crs%col_idx)
      allocate (new_vector(crs%size()))
      new_vector(:) = 0.0d0
      call crs_spmv_real64(val=crs%val, row_ptr=crs%row_ptr, &
                           col_idx=crs%col_idx, &
                           old_vector=old_vector, new_vector=new_vector, n=n, col_size=col_size)

   end function

!  function matmul_opencl_CRS(CRS,old_vector,opencl) result(new_vector)
!
!    class(CRS_),intent(in) :: CRS
!    real(real64),intent(in)  :: Old_vector(:)
!    real(real64),allocatable :: new_vector(:)
!    logical,intent(in) :: opencl
!    integer(int64) :: n,col_size
!
!    n = crs%size()
!    col_size=size(crs%col_idx)
!    allocate(new_vector(crs%size() ) )
!    new_vector(:) = 0.0d0
!    call crs_spmv_real64_opencl(val=crs%val,row_ptr=crs%row_ptr,&
!        col_idx=crs%col_idx,&
!        old_vector=old_vector,new_vector=new_vector,n=n,col_size=col_size)
!
!  end function

   function matmul_c_real32_CRS(CRS, old_vector, c_routine) result(new_vector)

      class(CRS_), intent(in) :: CRS
      real(real32), intent(in)  :: Old_vector(:)
      real(real32), allocatable :: new_vector(:)
      character(*), intent(in) :: c_routine
      integer(int64) :: n, col_size

      n = crs%size()
      col_size = size(crs%col_idx)
      allocate (new_vector(crs%size()))
      new_vector(:) = 0.0d0
      call crs_spmv_real32(val=real(crs%val), row_ptr=crs%row_ptr, &
                           col_idx=crs%col_idx, &
                           old_vector=old_vector, new_vector=new_vector, n=n, col_size=col_size)

   end function

   function matmul_real32_CRS(CRS, old_vector) result(new_vector)
      class(CRS_), intent(in) :: CRS
      real(real32), intent(in)  :: Old_vector(:)
      real(real32), allocatable :: new_vector(:)

      if (CRS%dtype /= real32) then
         print *, "ERROR :: matmul_real32_CRS >> vector type is real32, but &
 &            CRS%dtype /= real32. Plase use this%to_real32() before calling this."
         return
      end if

      new_vector = crs_matvec_generic_real32_SparseClass( &
                   CRS_value=CRS%val_real32, &
                   CRS_col=CRS%col_idx, &
                   CRS_row_ptr=CRS%row_ptr, &
                   old_vector=old_vector)

   end function

   function matmul_complex_CRS(CRS, old_vector) result(new_vector)
      class(CRS_), intent(in) :: CRS
      complex(real64), intent(in)  :: Old_vector(:)
      complex(real64), allocatable :: new_vector(:)
      integer(int32) :: n, row
      integer(int64) :: CRS_id

      n = size(CRS%row_ptr) - 1
      !if (size(old_vector) /= n) then
      !   print *, "ERROR crs_matvec :: inconsistent size for old_vector"
      !   return
      !end if

      new_vector = zeros(n)

      ! accerelation
      !$OMP parallel default(shared) private(CRS_id,row)
      !$OMP do reduction(+:new_vector)
      do row = 1, n
         do CRS_id = CRS%row_ptr(row), CRS%row_ptr(row + 1) - 1
            new_vector(row) = new_vector(row) + CRS%val(CRS_id)*old_vector(CRS%col_idx(CRS_id))
         end do
      end do
      !$OMP end do
      !$OMP end parallel

!    new_vector = crs_matvec_generic_complex_SparseClass(&
!      CRS_value= dcmplx(CRS%val),&
!      CRS_col=CRS%col_idx,&
!      CRS_row_ptr=CRS%row_ptr,&
!      old_vector=old_vector)

   end function

   function crs_matvec_generic_SparseClass(CRS_value, CRS_col, CRS_row_ptr, old_vector, &
                                           cache_size) result(new_vector)
      real(real64), intent(in)  :: CRS_value(:), Old_vector(:)
      integeR(int32), intent(in):: CRS_col(:)
      integeR(int64), intent(in):: CRS_row_ptr(:)
      integer(int32), optional, intent(in) :: cache_size

      real(real64), allocatable :: new_vector(:)
      integer(int32) :: i, j, n, gid, lid, row, col
      integer(int64) :: CRS_id

      integer(int64) :: block_idx, num_block, mem_size, num_row_in_block
      integer(int64), allocatable :: row_idx_from(:), row_idx_to(:)

      !> x_i = A_ij b_j

      n = size(CRS_row_ptr) - 1
      !if (size(old_vector) /= n) then
      !   print *, "ERROR crs_matvec :: inconsistent size for old_vector"
      !   return
      !end if

      new_vector = zeros(n)

      ! accerelation
!    do concurrent (row=1:n)
!        new_vector(row) = new_vector(row) + dot_product( &
!            CRS_value(CRS_row_ptr(row):CRS_row_ptr(row+1)-1),  &
!            old_vector(CRS_col(CRS_row_ptr(row):CRS_row_ptr(row+1)-1) )&
!            )
!    enddo

      !v2.0

      !v2.0
!    !if(.not.present(cache_size) )then
      !$OMP parallel default(shared)
      !$OMP do reduction(+:new_vector)
      do row = 1, n
         new_vector(row) = new_vector(row) + dot_product( &
                           CRS_value(CRS_row_ptr(row):CRS_row_ptr(row + 1) - 1), &
                           old_vector(CRS_col(CRS_row_ptr(row):CRS_row_ptr(row + 1) - 1)) &
                           )
      end do
      !$OMP end do
      !$OMP end parallel
!    else
!        ! This does not work
!
!        ! (1) 行列およびベクトルのメモリ消費量を算出する．
!        mem_size = sizeof(CRS_value) + sizeof(CRS_col)+sizeof(CRS_row_ptr)&
!            +sizeof(new_vector)+sizeof(old_vector)
!        ! (2) メモリ消費量/キャッシュサイズ=ブロック数を算出する．
!
!        ! (3) 1~nを，1~k, k+1~2*k, 2k+1~3*k,...とブロックに分割し，
!        ! 　　j番目の開始idxをrow_idx_from(j), 終了idxをrow_idx_to(j)へ格納
!        if(mem_size <= cache_size)then
!            num_block = 1
!            row_idx_from = [1]
!            row_idx_to   = [n]
!
!        else
!            num_block = mem_size/cache_size
!            if(mod(n,num_block)==0 )then
!                allocate(row_idx_from(num_block) )
!                allocate(row_idx_to(num_block) )
!                num_row_in_block = n/num_block
!                do i=1,num_block
!                    row_idx_from(i) = (i-1)*num_row_in_block+1
!                    row_idx_to(i)   = i*num_row_in_block
!                enddo
!            else
!                allocate(row_idx_from(num_block+1) )
!                allocate(row_idx_to(num_block+1) )
!                num_row_in_block = n/num_block
!                do i=1,num_block
!                    row_idx_from(i) = (i-1)*num_row_in_block+1
!                    row_idx_to(i)   = i*num_row_in_block
!                enddo
!                row_idx_from(num_block+1) = num_block*num_row_in_block+1
!                row_idx_to(num_block+1)   = n
!            endif
!        endif
!
!        do block_idx = 1,num_block
!            !$OMP parallel default(shared)
!            !$OMP do reduction(+:new_vector)
!            do row = row_idx_from(block_idx), row_idx_to(block_idx)
!                new_vector(row) = new_vector(row) + dot_product( &
!                    CRS_value(CRS_row_ptr(row):CRS_row_ptr(row+1)-1),  &
!                    old_vector(CRS_col(CRS_row_ptr(row):CRS_row_ptr(row+1)-1) )&
!                    )
!            enddo
!            !$OMP end do
!            !$OMP end parallel
!        enddo
!    endif

      !<v1.0>
!    !$OMP parallel default(shared) private(CRS_id,col)
!    !$OMP do reduction(+:new_vector)
!    do row = 1, n
!        do concurrent(CRS_id=CRS_row_ptr(row):CRS_row_ptr(row+1)-1)
!            new_vector(row) = new_vector(row) + CRS_value(CRS_id)*old_vector(CRS_col(CRS_id))
!        enddo
!    enddo
!    !$OMP end do
!    !$OMP end parallel
!
   end function
! ###################################################################

   function crs_matvec_generic_elemental_SparseClass(CRS_value, CRS_col, CRS_row_ptr, old_vector) result(new_vector)
      real(real64), intent(in)  :: CRS_value(:), Old_vector(:)
      integeR(int32), intent(in):: CRS_col(:)
      integeR(int64), intent(in):: CRS_row_ptr(:)

      real(real64), allocatable :: new_vector(:)
      integer(int32) :: i, j, n, gid, lid, row, col
      integer(int64) :: CRS_id

      !> x_i = A_ij b_j

      n = size(CRS_row_ptr) - 1
      !if (size(old_vector) /= n) then
      !   print *, "ERROR crs_matvec :: inconsistent size for old_vector"
      !   return
      !end if

      new_vector = zeros(n)

      ! accerelation
!    do concurrent (row=1:n)
!        new_vector(row) = new_vector(row) + dot_product( &
!            CRS_value(CRS_row_ptr(row):CRS_row_ptr(row+1)-1),  &
!            old_vector(CRS_col(CRS_row_ptr(row):CRS_row_ptr(row+1)-1) )&
!            )
!    enddo

      !v2.0

      !v2.0

      !$OMP parallel default(shared)
      !$OMP do reduction(+:new_vector)
      do row = 1, n
         new_vector(row) = new_vector(row) + dot_product( &
                           CRS_value(CRS_row_ptr(row):CRS_row_ptr(row + 1) - 1), &
                           old_vector(CRS_col(CRS_row_ptr(row):CRS_row_ptr(row + 1) - 1)) &
                           )
      end do
      !$OMP end do
      !$OMP end parallel

      !<v1.0>
!    !$OMP parallel default(shared) private(CRS_id,col)
!    !$OMP do reduction(+:new_vector)
!    do row = 1, n
!        do concurrent(CRS_id=CRS_row_ptr(row):CRS_row_ptr(row+1)-1)
!            new_vector(row) = new_vector(row) + CRS_value(CRS_id)*old_vector(CRS_col(CRS_id))
!        enddo
!    enddo
!    !$OMP end do
!    !$OMP end parallel
!
   end function
! ###################################################################

   subroutine spmv_as_subroutine_CRS(this, old_vector, new_vector)
      class(CRS_), intent(in) :: this
      real(real64), intent(in)  :: old_vector(:)
      real(real64), intent(out) :: new_vector(:)
      integer(int32) :: row
      new_vector(:) = 0.0d0

      !$OMP parallel default(shared)
      !$OMP do reduction(+:new_vector)
      do row = 1, size(new_vector)
         new_vector(row) = new_vector(row) + dot_product( &
                           this%val(this%row_ptr(row):this%row_ptr(row + 1) - 1), &
                           old_vector(this%col_idx(this%row_ptr(row):this%row_ptr(row + 1) - 1)) &
                           )
      end do
      !$OMP end do
      !$OMP end parallel

   end subroutine

   function crs_matvec_generic_real32_SparseClass(CRS_value, CRS_col, CRS_row_ptr, old_vector) result(new_vector)
      real(real32), intent(in)  :: CRS_value(:), Old_vector(:)
      integeR(int32), intent(in):: CRS_col(:)
      integeR(int64), intent(in):: CRS_row_ptr(:)

      real(real32), allocatable :: new_vector(:)
      integer(int32) :: i, j, n, gid, lid, row, col
      integer(int64) :: CRS_id

      !> x_i = A_ij b_j

      n = size(CRS_row_ptr) - 1
      !if (size(old_vector) /= n) then
      !   print *, "ERROR crs_matvec :: inconsistent size for old_vector"
      !   return
      !end if

      new_vector = zeros(n)

      ! accerelation
!    do concurrent (row=1:n)
!        new_vector(row) = new_vector(row) + dot_product( &
!            CRS_value(CRS_row_ptr(row):CRS_row_ptr(row+1)-1),  &
!            old_vector(CRS_col(CRS_row_ptr(row):CRS_row_ptr(row+1)-1) )&
!            )
!    enddo

      !v2.0
      !$OMP parallel default(shared)
      !$OMP do reduction(+:new_vector)
      do row = 1, n
         new_vector(row) = new_vector(row) + dot_product( &
                           CRS_value(CRS_row_ptr(row):CRS_row_ptr(row + 1) - 1), &
                           old_vector(CRS_col(CRS_row_ptr(row):CRS_row_ptr(row + 1) - 1)) &
                           )
      end do
      !$OMP end do
      !$OMP end parallel

      !<v1.0>
!    !$OMP parallel default(shared) private(CRS_id,col)
!    !$OMP do reduction(+:new_vector)
!    do row = 1, n
!        do concurrent(CRS_id=CRS_row_ptr(row):CRS_row_ptr(row+1)-1)
!            new_vector(row) = new_vector(row) + CRS_value(CRS_id)*old_vector(CRS_col(CRS_id))
!        enddo
!    enddo
!    !$OMP end do
!    !$OMP end parallel
!
   end function
! ###################################################################

!function crs_opencl_matvec(CRS_row_ptr,CRS_col,CRS_value,old_vector) result(new_vector)
!    integer(int32), intent(in)  :: CRS_row_ptr(:),CRS_col(:)
!    real(real64),intent(in) :: CRS_value(:), old_vector(:)
!
!
!end function
! ###################################################################

   function crs_matvec_generic_complex_SparseClass(CRS_value, CRS_col, CRS_row_ptr, old_vector) result(new_vector)
      complex(real64), intent(in)  :: CRS_value(:), Old_vector(:)
      integeR(int32), intent(in):: CRS_col(:)
      integeR(int64), intent(in):: CRS_row_ptr(:)

      complex(real64), allocatable :: new_vector(:)
      integer(int32) :: i, j, n, gid, lid, row, col
      integer(int64) :: CRS_id
      !> x_i = A_ij b_j

      n = size(CRS_row_ptr) - 1
      !if (size(old_vector) /= n) then
      !   print *, "ERROR crs_matvec :: inconsistent size for old_vector"
      !   return
      !end if

      new_vector = zeros(n)

! accerelation

      !$OMP parallel default(shared) private(CRS_id,row,col)
      !$OMP do reduction(+:new_vector)
      do row = 1, n
         do CRS_id = CRS_row_ptr(row), CRS_row_ptr(row + 1) - 1
            new_vector(row) = new_vector(row) + CRS_value(CRS_id)*old_vector(CRS_col(CRS_id))
         end do
      end do
      !$OMP end do
      !$OMP end parallel

!    !$OMP parallel do default(shared) private(CRS_id,col)
!    do row = 1 , n
!        do CRS_id = CRS_row_ptr(row), CRS_row_ptr(row+1)-1
!            col = CRS_col(CRS_id)
!            !$OMP atomic
!            new_vector(row) = new_vector(row) + CRS_value(CRS_id)*old_vector(col)
!        enddo
!    enddo
!    !$OMP end parallel do

   end function
! ###################################################################
!function crs_opencl_matvec(CRS_row_ptr,CRS_col,CRS_value,old_vector) result(new_vector)
!    integer(int32), intent(in)  :: CRS_row_ptr(:),CRS_col(:)
!    real(real64),intent(in) :: CRS_value(:), old_vector(:)
!
!
!end function
! ###################################################################

   subroutine eigCRS(this, Eigen_vectors, eigen_values)
      class(CRS_), intent(in) :: this
      real(real64), allocatable :: A(:, :), V(:, :)
      real(real64), allocatable :: D(:), E(:)
      real(real64), allocatable :: WORK(:)
      integer(int32), allocatable :: IWORK(:), IFAIL(:)
      real(real64), allocatable, intent(inout) :: Eigen_vectors(:, :), eigen_values(:)
      integer(int32) :: num_eigen, INFO
      real(real64)  :: VL, VU, ABSTOL
      integer(int32):: IL, IU
      character(1) :: JOBZ = "V"
      character(1) :: RANGE = "A"

      call this%Lanczos(DiagonalVector=D, subDiagonalVector=E, V=V)

      !call print("--")
      !call print(V)

      num_eigen = size(D)
      eigen_values = zeros(size(D))
      eigen_vectors = zeros(size(D), num_eigen)
      WORK = zeros(5*size(D))
      IWORK = int(zeros(5*size(D)))
      IFAIL = int(zeros(size(D)))
      ABSTOL = dble(1.0e-14)

      ! LAPACK
      call DSTEVX(JOBZ, RANGE, size(D), &
                  D, E, VL, VU, IL, IU, ABSTOL, &
                  num_eigen, eigen_values, eigen_vectors, size(D), &
                  WORK, IWORK, IFAIL, INFO)

      !eigen_Vectors = matmul(V,matmul(eigen_values, transpose(V) ) )
      eigen_Vectors = matmul(V, eigen_vectors)
      do i_i = 1, 3
         eigen_vectors(:, i_i) = eigen_vectors(:, i_i)/norm(eigen_vectors(:, i_i))
         eigen_vectors(:, i_i) = eigen_vectors(:, i_i)/eigen_vectors(3, i_i)
      end do

   end subroutine

   subroutine initCRS(this, val, col_idx, row_ptr)
      class(CRS_), intent(inout) :: this
      integer(int32), intent(in) :: col_idx(:)
      integer(int64), intent(in) ::  row_ptr(:)
      real(real64), intent(in)   :: val(:)

      this%val = val
      this%col_idx = col_idx
      this%row_ptr = row_ptr

   end subroutine

   function to_denseCRS(this) result(dense_mat)
      class(CRS_), intent(in) :: this
      real(real64), allocatable  :: dense_mat(:, :)
      integer(int32) :: i, n, m, row, col
      integeR(int64) :: j

      n = size(this%row_ptr) - 1
      m = maxval(this%col_idx)
      allocate (dense_mat(n, m))
      dense_mat(:, :) = 0.0d0
      do row = 1, n
         do j = this%row_ptr(row), this%row_ptr(row + 1) - 1
            col = this%col_idx(j)
            dense_mat(row, col) = this%val(j)
         end do
      end do

   end function

   pure function addCRS_and_CRS(CRS1, CRS2) result(CRS_ret)
      type(CRS_), intent(in) :: CRS1, CRS2
      type(CRS_) :: CRS_ret
      integer(int32) :: i, j, row
      integer(int64) :: col_2, col_1

      ! sum : CRS_ret = CRS1 + CRS2
      CRS_ret = CRS1

      ! ignore fill-in
    !!$OMP parallel do private(col_2,col_1)
      do row = 1, size(CRS2%row_ptr) - 1
         ! for each row
         do col_2 = CRS2%row_ptr(row), CRS2%row_ptr(row + 1) - 1
            ! search same col
            do col_1 = CRS1%row_ptr(row), CRS1%row_ptr(row + 1) - 1
               if (CRS1%col_idx(col_1) == CRS2%col_idx(col_2)) then
                    !!$OMP atomic
                  CRS_ret%val(col_1) = CRS_ret%val(col_1) + CRS2%val(col_2)
                  exit
               end if
            end do
         end do
      end do
    !!$OMP end parallel do

   end function

   pure function diffCRS_and_CRS(CRS1, CRS2) result(CRS_ret)
      type(CRS_), intent(in) :: CRS1, CRS2
      type(CRS_) :: CRS_ret
      integer(int32) :: i, j, row
      integer(int64) :: col_2, col_1

      ! sum : CRS_ret = CRS1 + CRS2
      CRS_ret = CRS1

      ! ignore fill-in
      do row = 1, size(CRS2%row_ptr) - 1
         ! for each row
         do col_2 = CRS2%row_ptr(row), CRS2%row_ptr(row + 1) - 1
            ! search same col
            do col_1 = CRS1%row_ptr(row), CRS1%row_ptr(row + 1) - 1
               if (CRS1%col_idx(col_1) == CRS2%col_idx(col_2)) then
                  CRS_ret%val(col_1) = CRS_ret%val(col_1) - CRS2%val(col_2)
                  exit
               end if
            end do
         end do
      end do

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

! #############################################
   function multReal64_and_CRS(scalar64, CRS1) result(CRS_ret)
      real(real64), intent(in) :: scalar64
      type(CRS_), intent(in) :: CRS1
      type(CRS_) :: CRS_ret
      integer(int32) :: i

      CRS_ret = CRS1
      !$OMP parallel do
      do i = 1, size(CRS_ret%val)
         CRS_ret%val(i) = scalar64*CRS_ret%val(i)
      end do
      !$OMP end parallel do

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


! #############################################
   function multReal64_and_BCRS(scalar64, BCRS1) result(BCRS_ret)
      real(real64), intent(in) :: scalar64
      type(BCRS_), intent(in) :: BCRS1
      type(BCRS_) :: BCRS_ret
      integer(int32) :: i,j

      BCRS_ret = BCRS1
      do i = 1, size(BCRS_ret%CRS,1)
         do j = 1, size(BCRS_ret%CRS,1)
            if(.not.allocated(BCRS_ret%CRS(i,j)%val)) cycle
            BCRS_ret%CRS(i,j) = scalar64*BCRS_ret%CRS(i,j)
         enddo
      end do
      
   end function
! #############################################


   ! #############################################
   function multBCRS_and_Real64(BCRS1,scalar64) result(BCRS_ret)
      real(real64), intent(in) :: scalar64
      type(BCRS_), intent(in) :: BCRS1
      type(BCRS_) :: BCRS_ret
      integer(int32) :: i,j

      BCRS_ret = BCRS1
      do i = 1, size(BCRS_ret%CRS,1)
         do j = 1, size(BCRS_ret%CRS,1)
            if(.not.allocated(BCRS_ret%CRS(i,j)%val)) cycle
            BCRS_ret%CRS(i,j) = scalar64*BCRS_ret%CRS(i,j)
         enddo
      end do
      
   end function
! #############################################

! #############################################
   function multCRS_and_Real64(CRS1, scalar64) result(CRS_ret)
      type(CRS_), intent(in) :: CRS1
      real(real64), intent(in) :: scalar64
      type(CRS_) :: CRS_ret
      integer(int32)::i

      CRS_ret = CRS1

      !$OMP parallel do
      do i = 1, size(CRS_ret%val)
         CRS_ret%val = scalar64*CRS_ret%val
      end do
      !$OMP end parallel do
   end function

   pure function DOFCOO(this) result(Degree_of_freedom)
      class(COO_), intent(in) :: this
      integer(int32) :: Degree_of_freedom

      if (allocated(this%row)) then
         Degree_of_freedom = size(this%row)
      else
         Degree_of_freedom = 0
      end if
   end function

   pure function DOFCRS(this) result(Degree_of_freedom)
      class(CRS_), intent(in) :: this
      integer(int32) :: Degree_of_freedom

      if (allocated(this%row_ptr)) then
         Degree_of_freedom = size(this%row_ptr) - 1
      else
         Degree_of_freedom = 0
      end if
   end function

! #######################################################
   pure function getCOO(this, row, col) result(ret)
      class(COO_), intent(in) :: this
      integer(int32), intent(in) :: row, col
      real(real64) :: ret
      integer(int32) :: i

      ret = 0.0d0
      if (.not. allocated(this%row(row)%col)) then
         return
      end if

      do i = 1, size(this%row(row)%col, 1)
         if (this%row(row)%col(i) == col) then
            ret = this%row(row)%val(i)
         end if
      end do

   end function

   pure function neCOO(this) result(ret)
      class(COO_), intent(in) :: this
      real(real64) :: ret
      integer(int32) :: i

      ! get number of entity
      ret = 0.0d0
      do i = 1, size(this%row)
         if (allocated(this%row(i)%col)) then
            ret = ret + size(this%row(i)%col)
         end if
      end do

   end function

   pure function maxvalCOO(this) result(ret)
      class(COO_), intent(in) :: this
      real(real64) :: ret
      real(real64), allocatable :: val(:)
      integer(int32) :: i, itr

      ! get number of entity
      ret = 0.0d0
      itr = 0
      do i = 1, size(this%row)
         if (allocated(this%row(i)%col)) then
            if (itr == 0) then
               val = this%row(i)%val(:)
               ret = maxval(val)
               itr = itr + 1
            else
               val = this%row(i)%val(:)
               if (ret < maxval(val)) then
                  ret = maxval(val)
               else
                  cycle
               end if
            end if
         end if
      end do

   end function

! ##################################################
   function sizeCRS(this) result(n)
      class(CRS_), intent(in) :: this
      integer(int32) :: n

      if (allocated(this%row_ptr)) then
         n = size(this%row_ptr) - 1
      else
         n = 0
      end if

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


! ##################################################
   function shapeCRS(this) result(ret)
      class(CRS_), intent(in) :: this
      integer(int32) :: ret(1:2)

      if(.not.allocated(this))then
         ret(:) = 0
         return
      endif

      if (allocated(this%row_ptr)) then
         ret(1) = size(this%row_ptr) - 1
      else
         ret(1) = 0
      end if

      ret(2) = maxval(this%col_idx)

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

! ##################################################
   subroutine updateCRS(this, row, col, val)
      class(CRS_), intent(inout) :: this
      integer(int32), intent(in) :: row, col
      real(real64), intent(in) :: val
      integer(int64) :: i, j

      ! update but ignore fill-in
      if (col > maxval(this%col_idx(this%row_ptr(row):this%row_ptr(row + 1) - 1))) return
      if (col < minval(this%col_idx(this%row_ptr(row):this%row_ptr(row + 1) - 1))) return

      do i = this%row_ptr(row), this%row_ptr(row + 1) - 1
         if (this%col_idx(i) == col) then
            this%val(i) = val
            return
         end if
      end do

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

! ##################################################
   subroutine addCRS(this, row, col, val)
      class(CRS_), intent(inout) :: this
      integer(int32), intent(in) :: row, col
      real(real64), intent(in) :: val
      integer(int64) :: i, j

      ! update but ignore fill-in
      if (col > maxval(this%col_idx(this%row_ptr(row):this%row_ptr(row + 1) - 1))) return
      if (col < minval(this%col_idx(this%row_ptr(row):this%row_ptr(row + 1) - 1))) return

      do i = this%row_ptr(row), this%row_ptr(row + 1) - 1
         if (this%col_idx(i) == col) then
            this%val(i) = this%val(i) + val
            return
         end if
      end do

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

! ##################################################
   function getCRS(this, row, col) result(val)
      class(CRS_), intent(in) :: this
      integer(int32), intent(in) :: row, col
      real(real64) :: val
      integer(int64) :: i

      ! update but ignore fill-in
      val = 0.0d0
      if (col > maxval(this%col_idx(this%row_ptr(row):this%row_ptr(row + 1) - 1))) return

      if (col < minval(this%col_idx(this%row_ptr(row):this%row_ptr(row + 1) - 1))) return

      do i = this%row_ptr(row), this%row_ptr(row + 1) - 1
         if (this%col_idx(i) == col) then
            val = this%val(i)
            return
         end if
      end do

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

! ##################################################
   logical function is_nonzeroCRS(this, row, col)
      class(CRS_), intent(in) :: this
      integer(int32), intent(in) :: row, col
      real(real64) :: val
      integer(int64) :: i

      is_nonzeroCRS = .false.
      if (row > this%size()) return

      !if (col > maxval(this%col_idx(this%row_ptr(row):this%row_ptr(row+1)-1))  ) return

      !if (col < minval(this%col_idx(this%row_ptr(row):this%row_ptr(row+1)-1))  ) return

      do i = this%row_ptr(row), this%row_ptr(row + 1) - 1
         if (this%col_idx(i) == col) then
            is_nonzeroCRS = .true.
            return
         end if
      end do

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

   function lumpedCRS(this) result(ret)
      class(CRS_), intent(in) :: this
      real(real64), allocatable  :: ret(:)

      ret = this%diag(cell_centered=.true.)

   end function

! ##################################################
   function diagCRS(this, cell_centered) result(diag_vec)
      class(CRS_), intent(in) :: this
      real(real64), allocatable  :: diag_vec(:)
      logical, optional, intent(in) :: cell_centered
      integeR(int64) :: i, j

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

            diag_vec = zeros(this%size())
            do i = 1, this%size() ! for rows
               do j = this%row_ptr(i), this%row_ptr(i + 1) - 1 ! for columns
                  diag_vec(i) = diag_vec(i) + this%val(j)
               end do
            end do
            return
         end if
      end if

      diag_vec = zeros(this%size())
      do i = 1, this%size()
         do j = this%row_ptr(i), this%row_ptr(i + 1) - 1
            if (this%col_idx(j) == i) then
               diag_vec(i) = this%val(j)
               exit
            end if
         end do
      end do

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

   recursive subroutine ILUCRS(this, fill_in_order, RHS, debug)
      class(CRS_), intent(inout) :: this
      integer(int32), intent(in) :: fill_in_order !
      real(real64), optional, intent(inout)   :: RHS(:) ! RHS vector
      real(real64), allocatable  :: diag_vec(:)
      integer(int32), allocatable  :: col_line(:), range_col(:, :)
      type(CCS_) :: ccs
      logical, optional, intent(in) :: debug
      real(real64) :: A_k_j, A_ik
      integer(int32) :: i, k, l, n, m, row, col, col_idx, row_idx, kk, jj
      integer(int64) :: j
      integer(int32), allocatable :: nonzero_k_list(:), nonzero_j_list(:)
      logical :: debug_mode_on = .false.

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

      if (present(RHS)) then
         ! Given: Ax = y
         ! A = LU
         call this%ILU(fill_in_order)

         ! Forward substitution
         do i = 2, this%size()
            do j = this%row_ptr(i), this%row_ptr(i + 1) - 1
               if (this%col_idx(j) < i) then
                  ! execute Forward substitution
                  !print *, i, this%col_idx(j),this%val(j)
                  RHS(i) = RHS(i) - this%val(j)*RHS(this%col_idx(j))
               end if
            end do
         end do

         ! Backward substitution
         diag_vec = this%diag()
         do i = this%size(), 1, -1
            if (diag_vec(i) == 0.0d0) cycle
            RHS(i) = RHS(i)/diag_vec(i)
         end do

         ! Forward substitution
         do i = this%size(), 1, -1
            do j = this%row_ptr(i), this%row_ptr(i + 1) - 1
               if (this%col_idx(j) > i) then
                  ! execute Forward substitution
                  RHS(i) = RHS(i) - this%val(j)/diag_vec(i)*RHS(this%col_idx(j))
               end if
            end do
         end do

         return
      end if

      select case (fill_in_order)
      case default
         print *, "fill-in order", fill_in_order, "is not implemented. Try 0"
         stop
      case (0)
         n = this%size()
         if (debug_mode_on) then
            print *, "[ILU(0)] >> started"
         end if
         diag_vec = this%diag()

         range_col = zeros(n, 2)

            !!$OMP parallel do
         !do row=1,n
         !    range_col(row,1) = minval(this%col_idx(this%row_ptr(row):this%row_ptr(row+1)-1))
         !    range_col(row,2) = maxval(this%col_idx(this%row_ptr(row):this%row_ptr(row+1)-1))
         !enddo
            !!$OMP end parallel do

         do i = 2, n
            if (debug_mode_on) then
               print *, "[ILU(0)] >> U", i, "/", n
            end if
            ! >>>>>>> slow
            nonzero_k_list = this%col_idx(this%row_ptr(i):this%row_ptr(i + 1) - 1)
                !! >> notice!!
                !! each col should be sorted.

            !do k = range_col(i,1) , i-1
            do kk = 1, size(nonzero_k_list)
               k = nonzero_k_list(kk)

               if (k >= i) cycle

               A_ik = this%get(i, k)/diag_vec(k)
               call this%update(row=i, col=k, val=A_ik)

               !$OMP parallel default(shared) private(j)
               !$OMP do
               do jj = 1, size(nonzero_k_list)
                  j = nonzero_k_list(jj)
                  if (j < k + 1) cycle
                  call this%add(row=int(i), col=int(j), val=-A_ik*this%get(int(k), int(j)))
               end do
               !$OMP end do
               !$OMP end parallel

               diag_vec(i) = this%get(int(i), int(i))

            end do
         end do
         ! <<<<<<<<< slow

!            !! ccs
         !ccs = this%to_CCS()

         ! [5,  6,  7]
         ! [2,  8,  9]
         ! [3,  4, 10]

         ! [   5,  6,  7]
         ! [10/5, 20, 23]
         ! [15/5, 50, 67]

!!          ! 高速バージョン

!            do i=2,n
!                if(debug_mode_on)then
!                    print *, "[ILU(0)] >> U",i,"/",n
!                endif
!                ! >>>>>>> slow
!                k = i-1
!                do row_idx = ccs%col_ptr(i-1),ccs%col_ptr(i)-1
!                    if(ccs%row_idx(row_idx) < i ) cycle
!                    if(i==3 ) then
!                        print *,"dbg",ccs%row_idx(row_idx),i-1,ccs%val(row_idx), diag_vec(i-1)
!
!                    endif
!
!                    ccs%val(row_idx) = ccs%val(row_idx)/diag_vec(i-1)
!                    if(ccs%row_idx(row_idx) == i ) then
!                        A_ik = ccs%val(row_idx)
!                    endif
!                enddo
!                if(i==3 ) stop
!
!                do j = this%row_ptr(i),this%row_ptr(i+1)-1
!                    if(this%col_idx(j) < i ) then
!                        cycle
!                    elseif(this%col_idx(j) == i ) then
!                        diag_vec(i) = diag_vec(i) - A_ik*this%get( k,i )
!                        call this%update(i,i,diag_vec(i) )
!                    else
!                        call this%update(row=i,col= this%col_idx(j) ,&
!                            val=this%get(i,this%col_idx(j) ) - A_ik*this%get(k,this%col_idx(j) ) )
!                    endif
!                enddo
!
!            enddo
!            ! [   5,  6,  7]
!            ! [10/5, 20, 23]
!            ! [15/5, 50, 67]
!            call this%load(ccs=ccs,position="L")
!            return

!
      end select

   end subroutine

! ################################################
   subroutine ILU_matvecCRS(this, old_vector, new_vector)
      class(CRS_), intent(in) :: this ! ILU factorlized matrix
      integer(int64) :: i, j
      real(real64), intent(in) :: old_vector(:)
      real(real64), allocatable, intent(inout) :: new_vector(:)
      real(real64), allocatable :: diag_vec(:)
      real(real64) :: new_vector_i
      new_vector = old_vector

      ! Forward substitution
      do i = 2, this%size()
         new_vector_i = new_vector(i)
         !$OMP parallel
         !$OMP do reduction(+:new_vector_i)
         do j = this%row_ptr(i), this%row_ptr(i + 1) - 1
            if (this%col_idx(j) < i) then
               ! execute Forward substitution
               new_vector_i = new_vector_i - this%val(j)*new_vector(this%col_idx(j))
            end if
         end do
         !$OMP end do
         !$OMP end parallel
         new_vector(i) = new_vector_i
      end do
      ! Backward substitution
      diag_vec = this%diag()

      !$OMP parallel
      !$OMP do
      do i = this%size(), 1, -1
         if (diag_vec(i) == 0.0d0) then
            diag_vec(i) = 1.0d0
         end if
         new_vector(i) = new_vector(i)/diag_vec(i)
      end do
      !$OMP end do
      !$OMP end parallel

      ! Forward substitution
      do i = this%size(), 1, -1

         new_vector_i = new_vector(i)
         !$OMP parallel
         !$OMP do reduction(+:new_vector_i)
         do j = this%row_ptr(i), this%row_ptr(i + 1) - 1
            if (this%col_idx(j) > i) then
               ! execute Forward substitution
               new_vector(i) = new_vector(i) - this%val(j)/diag_vec(i)*new_vector(this%col_idx(j))
            end if
         end do
         !$OMP end do
         !$OMP end parallel
         new_vector(i) = new_vector_i

      end do

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

   function to_CCSCRS(this) result(CCS)
      class(CRS_), intent(in) :: this
      type(CCS_) :: CCS
      integer(int64) :: i, j, n, col
      integer(int32), allocatable :: inst_counter(:)

      inst_counter = int(zeros(this%size()))
      CCS%col_ptr = int(zeros(size(this%row_ptr)))
      CCS%row_idx = int(zeros(size(this%col_idx)))
      CCS%val = zeros(size(this%val))

      do i = 1, size(this%col_idx)
         CCS%col_ptr(this%col_idx(i)) = CCS%col_ptr(this%col_idx(i)) + 1
      end do

      ![2,3,3, 3, 2, 0]
      do i = 1, size(CCS%col_ptr) - 1
         CCS%col_ptr(i + 1) = CCS%col_ptr(i + 1) + CCS%col_ptr(i)
      end do
      ! [2,3,3, 3, 2,  0]
      !>[2,5,8,11,13, 13]
      do i = size(CCS%col_ptr), 2, -1
         CCS%col_ptr(i) = CCS%col_ptr(i - 1)
      end do
      ! [2,5,8,11,13, 13]
      !>[2,2,5, 8,11, 13]
      CCS%col_ptr(:) = CCS%col_ptr(:) + 1

      ! [2,2,5, 8,11, 13]
      !>[2,3,6, 9,12, 14]
      CCS%col_ptr(1) = 1
      ! [2,3,6, 9,12, 14]
      !>[1,3,6, 9,12, 14]

      do i = 1, size(this%row_ptr) - 1
         do j = this%row_ptr(i), this%row_ptr(i + 1) - 1
            col = this%col_idx(j)
            CCS%row_idx(CCS%col_ptr(col) + inst_counter(col)) = i
            CCS%val(CCS%col_ptr(col) + inst_counter(col)) = this%val(j)
            inst_counter(col) = inst_counter(col) + 1
         end do
      end do

   end function
!####################################################
   subroutine loadCRS(this, CCS, Position)
      class(CRS_), intent(inout) :: this
      type(CCS_), optional, intent(in) :: CCS
      character(*), optional, intent(in) :: Position
      integer(int32) :: i, j

      if (present(CCS)) then
         if (Position == "L") then
            print *, CCS%col_ptr
            print *, CCS%row_idx
            print *, CCS%val

            do i = 1, size(CCS%col_ptr) - 1
               do j = ccs%col_ptr(i), ccs%col_ptr(i + 1) - 1

                  if (ccs%row_idx(j) <= i) then
                     cycle
                  else
                     print *, ccs%row_idx(j), i
                     call this%update(row=ccs%row_idx(j), col=i, val=ccs%val(j))
                  end if
               end do
            end do
         end if
      end if
   end subroutine

   subroutine eyesCOO(this, n)
      class(COO_) :: this
      integer(int32), intent(in) :: n
      integer(int32) :: i

      call this%init(n)
      do i = 1, n
         call this%set(i, i, 1.0d0)
      end do

   end subroutine

   subroutine poissonCOO(this, n)
      class(COO_) :: this
      integer(int32), intent(in) :: n
      integer(int32) :: i

      call this%init(n)
      call this%set(1, 1, -2.0d0)
      call this%set(1, 1 + 1, 1.0d0)
      do i = 2, n - 1
         call this%set(i, i - 1, 1.0d0)
         call this%set(i, i, -4.0d0)
         call this%set(i, i + 1, 1.0d0)
      end do
      call this%set(n, n - 1, 1.0d0)
      call this%set(n, n, -2.0d0)

   end subroutine

   subroutine LOBPCG_single_CRS(A, B, lambda, X, alpha, tol, debug)
      type(CRS_), intent(in) :: A, B
      real(real64), allocatable, intent(inout) :: X(:)
      real(real64), allocatable, intent(inout) :: lambda
      real(real64), intent(in) :: alpha
      logical, optional, intent(in) :: debug
      logical :: debug_mode

      type(Random_) :: random

      real(real64), allocatable :: r(:), rho
      integer(int32) :: n, i
      integer(int32) :: MAX_ITR = 1000000
      real(real64), intent(in) :: tol

      debug_mode = input(default=.false., option=debug)
      ! BLOPEX IN HYPRE AND PETSC, SIAM Eq. (2.2)
      ! number of eigen value :: m
      ! initialize X and lambda
      n = A%size()
      X = 2.0d0*random%randn(n) - 1.0d0

      lambda = 0.0d0

      !Single-vector version
      do i = 1, MAX_ITR
         rho = dot_product(x, A%matmul(x))/dot_product(x, B%matmul(x))
         r = A%matmul(x) - rho*B%matmul(x)
         x = x - alpha*r
         if (debug_mode) then
            print *, i, norm(r)
         end if
         if (norm(r) < tol) then
            if (debug_mode) then
               print *, "[OK] converged."
            end if
            exit
         else
            cycle
         end if
      end do
      if (i == MAX_ITR) print *, "[ERROR] LOBPCG NOT CONVERGED."
      lambda = rho
   end subroutine

   subroutine LOBPCG_CRS(A, B, lambda, X, m, MAX_ITR, TOL, debug)
      type(CRS_), intent(in) :: A, B
      real(real64), allocatable, intent(out) :: X(:, :)
      real(real64), allocatable, intent(out) :: lambda(:)
      real(real64), intent(in) :: TOL
      integer(int32), intent(in) :: m, MAX_ITR
      logical, optional, intent(in) :: debug
      logical :: debug_mode

      type(Random_) :: random

      real(real64), allocatable :: r(:, :), p(:, :), V(:, :), A_small(:, :), AV(:, :), &
                                   b_small(:, :), lambda_small(:), XAX(:, :), AX(:, :), OR(:, :), BX(:, :), XBX(:, :), &
                                   Bm_small(:, :), BV(:, :)
      real(real64)   :: initial_R
      integer(int32) :: n, i, j, k
      !integer(int32) :: MAX_ITR = 1000000

      ! References:
      ! Knyazev, A. V., Argentati, M. E., Lashuk, I. & Ovtchinnikov, E. E. Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) in hypre and PETSc. SIAM J. Sci. Comput. 29, 2224–2239 (2007).
      ! https://en.wikipedia.org/wiki/LOBPCG
      ! https://webpark1378.sakura.ne.jp/nagai/note_141009_LOBPCG.pdf

      debug_mode = input(default=.false., option=debug)
      ! BLOPEX IN HYPRE AND PETSC, SIAM Eq. (2.2)
      ! number of eigen value :: m

      ! debug

      ! initialize X and lambda
      n = A%size()

      !X = random%randn(n,m)
      X = random%randn(n, m)! + random%randn(n,m) !+ random%randn(n,m)
      do i = 1, m
         X(:, i) = X(:, i)/norm(X(:, i))
      end do

      call GramSchmidt(X, size(X, 1), size(X, 2), X)

      lambda = zeros(m)

      R = zeros(n, m)
      P = zeros(n, m)
      AV = zeros(n, 3*m)
      AX = zeros(n, m)
      A_small = zeros(3*m, 3*m)
      V = zeros(n, 3*m)

      ! m x m 固有値問題

      do i = 1, m
         AX(:, i) = A%matmul(X(:, i))
      end do

      BX = zeros(size(AX, 1), size(AX, 2))
      do i = 1, m
         BX(:, i) = B%matmul(X(:, i))
      end do

      XAX = matmul(transpose(X(:, 1:m)), AX)
      XBX = matmul(transpose(X(:, 1:m)), BX)
      call LAPACK_EIG(XAX, XBX, b_small, lambda_small)
      X = matmul(X, b_small)

      do i = 1, m
         AX(:, i) = A%matmul(X(:, i))
      end do
      BX = zeros(size(AX, 1), size(AX, 2))
      do i = 1, m
         BX(:, i) = B%matmul(X(:, i))
      end do

      do i = 1, m
         R(:, i) = AX(:, i) - BX(:, i)*lambda_small(i)
      end do

      V = zeros(n, 2*m)

      V(:, 1:m) = X(:, :)
      V(:, m + 1:2*m) = R(:, :)

      call GramSchmidt(V, size(V, 1), size(V, 2), V)
      do k = 1, size(X, 2)
         V(:, k) = V(:, k)/norm(V(:, k))
      end do

      X(:, :) = V(:, 1:m)
      R(:, :) = V(:, m + 1:2*m)

      AV = zeros(n, 2*m)
      do j = 1, 2*m  ! n x n, n x 2m
         AV(:, j) = A%matmul(V(:, j))
      end do

      BV = zeros(n, 2*m)
      do j = 1, 2*m  ! n x n, n x 2m
         BV(:, j) = B%matmul(V(:, j))
      end do

      A_small = matmul(transpose(V), AV)
      Bm_small = matmul(transpose(V), BV)

      call LAPACK_EIG(A_small, Bm_small, b_small, lambda_small)

      do j = 1, m
         X(:, j) = matmul(V, b_small(:, j))
      end do
      AX = zeros(n, m)
      do i = 1, m
         AX(:, i) = A%matmul(X(:, i))
      end do

      BX = zeros(size(AX, 1), size(AX, 2))
      do j = 1, m
         BX(:, j) = B%matmul(X(:, j))
      end do

      R = zeros(n, m)
      do j = 1, m
         R(:, j) = AX(:, j) - BX(:, j)*lambda_small(j)
      end do

      OR = zeros(n, 2*m)
      OR(:, 1:m) = 0.0d0
      OR(:, m + 1:2*m) = R(:, :)

      P = zeros(n, m)
      do i = 1, m
         P(:, i) = matmul(OR, b_small(:, i))
      end do

      V = zeros(n, 3*m)
      AV = zeros(n, 3*m)
      BV = zeros(n, 3*m)

      AX = zeros(n, m)
      BX = zeros(n, m)

      do i = 1, MAX_ITR
         V(:, 1:m) = X(:, :)
         V(:, m + 1:2*m) = R(:, :)
         V(:, 2*m + 1:3*m) = P(:, :)

         ! Gram-Scmidtを計算する．
         call GramSchmidt(V, size(V, 1), size(V, 2), V)
         do k = 1, size(X, 2)
            V(:, k) = V(:, k)/norm(V(:, k))
         end do
         do k = size(X, 2) + size(R, 2) + 1, size(V, 2)
            V(:, k) = V(:, k)/norm(V(:, k))
         end do

         X(:, :) = V(:, 1:m)
         R(:, :) = V(:, m + 1:2*m)
         P(:, :) = V(:, 2*m + 1:3*m)

         !$OMP parallel do
         do j = 1, 3*m
            AV(:, j) = A%matmul(V(:, j))
         end do
         !$OMP end parallel do

         !$OMP parallel do
         do j = 1, 3*m
            BV(:, j) = B%matmul(V(:, j))
         end do
         !$OMP end parallel do

         A_small = matmul(transpose(V), AV)
         Bm_small = matmul(transpose(V), BV)

         ! LAPACKの固有値求めるやつを呼ぶ
         call LAPACK_EIG(A_small, Bm_small, b_small, lambda_small)

         X = matmul(V, b_small(:, 1:m))

         !$OMP parallel do
         do j = 1, m
            AX(:, j) = A%matmul(X(:, j))
         end do
         !$OMP end parallel do

         !$OMP parallel do
         do j = 1, m
            BX(:, j) = B%matmul(X(:, j))
         end do
         !$OMP end parallel do

         !$OMP parallel do
         do j = 1, m
            R(:, j) = AX(:, j) - BX(:, j)*lambda_small(j)
         end do
         !$OMP end parallel do

         V(:, 1:m) = 0.0d0

         ! matmul:: (n,3*m) x (n,3*m)
         P = matmul(V, b_small(:, 1:m))

         ! Detect convergence
         if (i == 1) then
            initial_R = maxval(abs(R))
         end if

         if (debug_mode) then
            print *, i, maxval(abs(R))/initial_R
         end if
         if (maxval(abs(R))/initial_R < tol) then
            if (debug_mode) then
               print *, "[OK] converged."
            end if
            lambda = lambda_small(1:m)

            exit
         else
            cycle
         end if
      end do

      if (i == MAX_ITR) print *, "[ERROR] LOBPCG NOT CONVERGED."

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

! #################################################
   subroutine GramSchmidt(mat_v, m, n, mat_v_out)
      integer, intent(in)::m, n
      real(real64), intent(in)::mat_v(1:m, 1:n)
      real(real64), intent(inout)::mat_v_out(1:m, 1:n)
      integer::i, j
      real(real64) :: nai
      !real(real64) :: v(1:m),vi(1:m),viold(1:m)
      real(real64), allocatable :: v(:), vi(:)
      real(real64), allocatable :: viold(:)
      real(real64)::norm_v

      mat_v_out = mat_v
      allocate (viold(1:m))
      allocate (v(1:m))
      allocate (vi(1:m))

      do i = 1, n
         viold = mat_v_out(:, i)

         if (dot_product(viold, viold) == 0.0d0) cycle

         do j = 1, i - 1
            nai = dot_product(mat_v_out(1:m, j), viold(1:m))
            vi = viold - nai*mat_v_out(:, j)
            viold = vi
         end do

         norm_v = sqrt(dble(dot_product(viold, viold)))
         mat_v_out(:, j) = viold/norm_v

      end do

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

! #################################################
   subroutine LAPACK_EIG_DENSE(A, B, x, lambda, debug)
      real(real64), intent(in) :: A(:, :), B(:, :)
      real(real64), allocatable, intent(inout):: x(:, :), lambda(:)
      logical, optional, intent(in) :: debug
      logical :: debug_mode
      !>>>>>>>>>>>>>> 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(:)
      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(CRS_) :: crs

      debug_mode = input(default=.false., option=debug)
      !>>>>>>>>>>>>>> INPUT
      N = size(A, 1)
      LDZ = N
      LWORK = 1 + 6*N + 2*N**2
      LIWORK = 3 + 5*N
      !<<<<<<<<<<<<<< INPUT

      !>>>>>>>>>>>>>>  INPUT/OUTPUT
      AP = zeros(N*(N + 1)/2)
      BP = zeros(N*(N + 1)/2)
      ! Upper triangle matrix

      AP = to_UpperTriangle(A)!UpperTriangularMatrix(CRS_val=this%A_CRS_val,CRS_col=this%A_CRS_index_col,&
      !CRS_rowptr=this%A_CRS_index_row)
      BP = to_UpperTriangle(B)!UpperTriangularMatrix(CRS_val=this%B_CRS_val,CRS_col=this%B_CRS_index_col,&
      !CRS_rowptr=this%B_CRS_index_row)
      !<<<<<<<<<<<<<< INPUT/OUTPUT

      !>>>>>>>>>>>>>> INPUT
      N = size(A, 1)
      LDZ = N
      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

      if (debug_mode) then
         print *, ">> Solver :: LAPACK/DSPGVD"
      end if
      call DSPGVD(ITYPE, JOBZ, UPLO, N, AP, BP, W, Z, LDZ, WORK, &
                  LWORK, IWORK, LIWORK, INFO)

      X = Z
      lambda = W

   end subroutine
! #################################################
! #################################################
   subroutine LAPACK_EIG_SPARSE(A, B, x, lambda, debug)
      type(CRS_), intent(in) :: A, B
      real(real64), allocatable :: Ad(:, :), Bd(:, :)
      real(real64), allocatable, intent(inout):: x(:, :), lambda(:)
      logical, optional, intent(in) :: debug
      logical :: debug_mode
      !>>>>>>>>>>>>>> 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(:)
      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(CRS_) :: crs

      Ad = A%to_dense()
      Bd = B%to_dense()
      debug_mode = input(default=.false., option=debug)
      !>>>>>>>>>>>>>> INPUT
      N = size(Ad, 1)
      LDZ = N
      LWORK = 1 + 6*N + 2*N**2
      LIWORK = 3 + 5*N
      !<<<<<<<<<<<<<< INPUT

      !>>>>>>>>>>>>>>  INPUT/OUTPUT
      AP = zeros(N*(N + 1)/2)
      BP = zeros(N*(N + 1)/2)
      ! Upper triangle matrix

      AP = to_UpperTriangle(Ad)!UpperTriangularMatrix(CRS_val=this%A_CRS_val,CRS_col=this%A_CRS_index_col,&
      !CRS_rowptr=this%A_CRS_index_row)
      BP = to_UpperTriangle(Bd)!UpperTriangularMatrix(CRS_val=this%B_CRS_val,CRS_col=this%B_CRS_index_col,&
      !CRS_rowptr=this%B_CRS_index_row)
      !<<<<<<<<<<<<<< INPUT/OUTPUT

      !>>>>>>>>>>>>>> INPUT
      N = size(Ad, 1)
      LDZ = N
      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

      if (debug_mode) then
         print *, ">> Solver :: LAPACK/DSPGVD"
      end if
      call DSPGVD(ITYPE, JOBZ, UPLO, N, AP, BP, W, Z, LDZ, WORK, &
                  LWORK, IWORK, LIWORK, INFO)

      X = Z
      lambda = W

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

   function to_UpperTriangle(A) result(ret)
      real(real64), intent(in) :: A(:, :)
      real(real64), allocatable :: ret(:)
      integer(int32) :: i, j, n, m

      n = size(A, 1)
      ret = zeros(n*(n + 1)/2)

      m = 0
      do i = 1, size(A, 2)
         do j = 1, i
            m = m + 1
            ret(m) = A(j, i)
         end do
      end do

   end function

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

   subroutine BICGSTAB_CRSSparse(this, x, b, debug, tol)
      class(CRS_), intent(inout) :: this
      real(real64), allocatable, intent(inout) :: x(:)
      logical, optional, intent(in) :: debug
      real(real64), optional, intent(in) :: tol
      real(real64) :: er
      real(real64), intent(in) :: b(:)

      if (.not. allocated(x)) then
         x = zeros(size(b))
      end if

      er = dble(1.0e-14)
      if (present(tol)) then
         er = tol
      end if

      call bicgstab_CRS_SparseClass( &
         a=this%val, &
         ptr_i=this%row_ptr, &
         index_j=this%col_idx, &
         x=x, &
         b=b, &
         itrmax=1000000, &
         er=er, &
         relative_er=er, &
         debug=debug &
         )

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

! #####################################################
   subroutine bicgstab_CRS_SparseClass(a, ptr_i, index_j, x, b, itrmax, er, relative_er, debug)
      integer(int64), intent(in) :: ptr_i(:)
      integer(int32), intent(in) :: index_j(:)
      integer(int32), intent(in) :: itrmax
      real(real64), intent(in) :: a(:)
      real(real64), intent(in) :: b(:)
      real(real64), intent(in) :: 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 SpMV_CRS_Sparse(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 SpMV_CRS_Sparse(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 SpMV_CRS_Sparse(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 SpMV_CRS_Sparse(CRS_value, CRS_col, CRS_row_ptr, old_vector, new_vector, precondition)
      real(real64), intent(in)  :: CRS_value(:), Old_vector(:)
      integeR(int32), intent(in):: CRS_col(:)
      integeR(int64), intent(in):: CRS_row_ptr(:)
      type(CRS_), optional, intent(in) :: precondition
      real(real64), allocatable, intent(inout) :: new_vector(:)
      real(real64), allocatable :: precon_old_vector(:)
      integer(int32) :: i, j, n, gid, lid, row, col
      integer(int64) :: CRS_id
      !> x_i = A_ij b_j

      n = size(CRS_row_ptr) - 1
      !if (size(old_vector) /= n) then
      !   print *, "ERROR crs_matvec :: inconsistent size for old_vector"
      !   return
      !end if

      if (.not. allocated(new_vector)) then
         new_vector = zeros(n)
      else
         new_vector(:) = 0.0d0
      end if

      if (present(precondition)) then
         call precondition%ILU_matvec(old_vector=old_vector, new_vector=precon_old_vector)

         !$OMP parallel do default(shared) private(CRS_id,col)
         do row = 1, n
            do CRS_id = CRS_row_ptr(row), CRS_row_ptr(row + 1) - 1
               col = CRS_col(CRS_id)
               !$OMP atomic
               new_vector(row) = new_vector(row) + CRS_value(CRS_id)*precon_old_vector(col)
            end do
         end do
         !$OMP end parallel do
      else

         !$OMP parallel do default(shared) private(CRS_id,col)
         do row = 1, n
            do CRS_id = CRS_row_ptr(row), CRS_row_ptr(row + 1) - 1
               col = CRS_col(CRS_id)
               !$OMP atomic
               new_vector(row) = new_vector(row) + CRS_value(CRS_id)*old_vector(col)
            end do
         end do
         !$OMP end parallel do

      end if

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

   function divide_by_CRS(this, diag_vector) result(ret_crs)
      class(CRS_), intent(in) :: this
      real(real64), intent(in) :: diag_vector(:)
      type(CRS_) :: ret_crs
      integer(int64) :: i, j
      ret_crs = this

      if (size(diag_vector) == 1) then
         ret_crs%val(:) = ret_crs%val(:)/diag_vector(1)
         return
      end if

      if (this%size() /= size(diag_vector)) then
         print *, "ERROR >> divide_by_CRS >> this%size()/=size(diag_vector) "
         stop
      end if

      do i = 1, ret_crs%size()
         do j = ret_crs%row_ptr(i), ret_crs%row_ptr(i + 1) - 1
            ret_crs%val(j) = ret_crs%val(j)/diag_vector(i)
         end do
      end do

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

   subroutine divide_by_vector_CRS(this, diag_vector) 
      class(CRS_), intent(inout) :: this
      real(real64), intent(in) :: diag_vector(:)
      integer(int64) :: i, j

      
      if (size(diag_vector) == 1) then
         this%val(:) = this%val(:)/diag_vector(1)
         return
      end if

      if (this%size() /= size(diag_vector)) then
         print *, "ERROR >> divide_by_CRS >> this%size()/=size(diag_vector) "
         stop
      end if

      do i = 1, this%size()
         do j = this%row_ptr(i), this%row_ptr(i + 1) - 1
            this%val(j) = this%val(j)/diag_vector(i)
         end do
      end do

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

   function mult_by_CRS(this, diag_vector) result(ret_crs)
      class(CRS_), intent(in) :: this
      real(real64), intent(in) :: diag_vector(:)
      type(CRS_) :: ret_crs
      integer(int64) :: i, j
      ret_crs = this

      if (size(diag_vector) == 1) then
         ret_crs%val(:) = ret_crs%val(:)/diag_vector(1)
         return
      end if

      if (this%size() /= size(diag_vector)) then
         print *, "ERROR >> divide_by_CRS >> this%size()/=size(diag_vector) "
         stop
      end if

      do i = 1, ret_crs%size()
         do j = ret_crs%row_ptr(i), ret_crs%row_ptr(i + 1) - 1
            ret_crs%val(j) = ret_crs%val(j)*diag_vector(i)
         end do
      end do

   end function
! ######################################################
   function matmul_CRS_vec(A, b) result(vec)
      type(CRS_), intent(in) :: A
      real(real64), intent(in) :: b(:)
      real(real64), allocatable :: vec(:)

      vec = A%matmul(b)

   end function

! ######################################################
   function matmul_CRS_CRS(A, B) result(CRS)
      type(CRS_), intent(in) :: A, B
      type(COO_) :: COO
      type(CRS_) :: CRS
      type(CCS_) :: B_CCS
      real(real64), allocatable :: col_vec(:)
      integer(int32) :: col, row, col_ptr, i, j
      integer(int32), allocatable :: A_col_idx(:), B_row_idx(:)
      real(real64), allocatable :: A_col_val(:), B_row_val(:)

      ! Considerably slow!! DO NOT USE!
      ! Considerably slow!! DO NOT USE!
      ! Considerably slow!! DO NOT USE!
      ! Considerably slow!! DO NOT USE!
      ! Considerably slow!! DO NOT USE!
      ! Considerably slow!! DO NOT USE!
      ! Considerably slow!! DO NOT USE!
      ! Considerably slow!! DO NOT USE!
      ! Considerably slow!! DO NOT USE!
      ! Considerably slow!! DO NOT USE!
      ! Considerably slow!! DO NOT USE!

      ! matrix multiplication
      ! C_ij = A_ik B_kj
      call COO%init(A%size())
      B_CCS = B%to_CCS()

!    ! avoiding fill-in (inaccurate)
!    if(A%size() > B%size() )then
!        CRS = A
!        CRS%val(:) = 0.0d0
!    else
!        CRS = B
!        CRS%val(:) = 0.0d0
!    endif
!
!    do row=1,size(CRS%row_ptr)-1
!        do col_ptr=CRS%row_ptr(row),CRS%row_ptr(row+1)-1
!            col = CRS%col_idx(col_ptr)
!            A_col_idx = A%col_idx( A%row_ptr(row):A%row_ptr(row+1)-1 )
!            B_row_idx = B_CCS%row_idx( B_CCS%col_ptr(col):B_CCS%col_ptr(col+1)-1 )
!            A_col_val = A%val( A%row_ptr(row):A%row_ptr(row+1)-1 )
!            B_row_val = B_CCS%val( B_CCS%col_ptr(col):B_CCS%col_ptr(col+1)-1 )
!            if(minval(A_col_idx) > maxval(B_row_idx))cycle
!            if(maxval(A_col_idx) < minval(B_row_idx))cycle
!            do i=1,size(A_col_idx)
!                do j=1, size(B_row_idx)
!                    if(A_col_idx(i)==B_row_idx(j) )then
!                        CRS%val(col) = CRS%val(col) + A_col_val(i)*B_row_val(j)
!                    endif
!                enddo
!            enddo
!        enddo
!    enddo

      do col = 1, size(B_CCS%col_ptr) - 1
         col_vec = A%matmul(B_CCS%get_column(col))
        !!$OMP parallel do
         do row = 1, size(col_vec)
            if (col_vec(row) /= 0.0d0) then
                !!$OMP atomic
               call COO%set(row, col, col_vec(row))
            end if
         end do
        !!$OMP end parallel do
      end do
      CRS = COO%to_CRS()

   end function
! ###################################################
   function get_column_CCS(this, col) result(ret)
      class(CCS_), intent(in) :: this
      integer(int32), intent(in) :: col
      real(real64), allocatable :: ret(:)
      integer(int32) :: row, n

      n = size(this%col_ptr, 1) - 1
      ret = zeros(n)

      do row = this%col_ptr(col), this%col_ptr(col + 1) - 1
         ret(this%row_idx(row)) = this%val(row)
      end do

   end function

! #####################################################
   subroutine randomCOO(this, n, percent)
      class(COO_), intent(inout) :: this
      integer(int32), intent(in) :: n
      real(real32), optional, intent(in) :: percent
      real(real32) :: fill_percent
      integer(int32) :: m, i, row, col
      real(real64) :: val
      type(Random_) :: random

      fill_percent = input(default=10.0, option=percent)
      m = maxval([int(dble(n)*dble(n)*dble(percent)/100.0d0), 1])
      call this%init(n)
      do i = 1, m
         row = int(n*random%random()) + 1
         col = int(n*random%random()) + 1
         val = random%gauss(mu=0.0d0, sigma=1.0d0)
         call this%set(row, col, val)
         call this%set(col, row, val)
      end do

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

   subroutine eyesCRS(this, n)
      class(CRS_), intent(inout) :: This
      type(COO_) :: coo
      type(CRS_) :: buf
      integer(int32), intent(in) :: n
      integer(int32) :: i

      call coo%init(n)
      do i = 1, n
         call coo%set(i, i, 1.0d0)
      end do
      buf = coo%to_CRS()
      this%row_ptr = buf%row_ptr
      this%col_idx = buf%col_idx
      this%val = buf%val
   end subroutine
! ###################################################


! ###################################################
subroutine poissonCRS(this, n)
   class(CRS_),intent(inout) :: this
   type(COO_) :: coo
   type(CRS_) :: buf
   integer(int32),intent(in) :: n

   call coo%poisson(n)
   buf = coo%to_crs()
   this%row_ptr = buf%row_ptr
   this%col_idx = buf%col_idx
   this%val = buf%val

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


! ###################################################
   function tensor_exponential_crs(this, itr_tol, tol, x, dt, fix_idx, fix_val) result(expA_v)
      class(CRS_), intent(in) :: this
      real(real64), intent(in) :: x(:)
      real(real64), optional, intent(in) :: dt
      real(real64), allocatable :: expA_v(:), increA_v(:), bhat(:)

      real(real64), intent(in) :: tol
      integer(int32), intent(in)::itr_tol
!    real(real64),allocatable::increA(:,:)
      real(real64)   :: increment, NN, t
      integer(int32) :: i, j, n

      real(real64), optional, intent(in) :: fix_val(:)
      integer(int32), optional, intent(in)::fix_idx(:)
      type(CRS_) :: Amatrix
!    if(.not. allocated(expA) )allocate(expA(size(A,1),size(A,2) ))
!    allocate(increA(size(A,1),size(A,2) ))
!    if(size(A,1)/=size(A,2)) stop "ERROR:tensor exp is not a square matrix"

!    expA_v = x
!    do n=1,size(expA,1)
!        expA(n,n)=1.0d0
!    enddo
      t = input(default=1.0d0, option=dt)
      ! exp(At) = I + At + 1/2*(At)^2 + 1/6*(At)^3 + ....
      ! exp(At)v = Iv + (At)v + 1/2*(At)^2 v + 1/6*(At)^3 v+ ....
      ! exp(At)v = v + (At)v + 1/2*(At)^2 v + 1/6*(At)^3 v+ ....
!    expA_v = x
!    do i=1,itr_tol
!        increA_v = x
!        do j=1,i
!            increA_v = t/dble(j)*this%matmul(increA_v)
!        enddo
!        expA_v = expA_v + increA_v
!        if(dot_product(increA_v,increA_v) < tol )exit
!
!    enddo
!    return

      ! i==1,2
      ! 0-th order term
      if (present(fix_idx)) then
         Amatrix = this
         bhat = zeros(this%size())
         call Amatrix%fix(idx=fix_idx, RHS=bhat, val=fix_val)

         increA_v = x
         expA_v = zeros(size(x))
         expA_v = expA_v + increA_v
         ! 1-st order term
         increA_v = Amatrix%matmul(t*increA_v) - bhat*dt
         expA_v = expA_v + increA_v

         do i = 2, itr_tol
            increA_v = 1.0d0/(i)*Amatrix%matmul(t*increA_v) !- bhat*dt
            expA_v = expA_v + increA_v

            if (dot_product(increA_v, increA_v) < tol) exit
         end do
      else

         increA_v = x
         expA_v = zeros(size(x))
         expA_v = expA_v + increA_v
         ! 1-st order term
         increA_v = this%matmul(t*increA_v)
         expA_v = expA_v + increA_v

         do i = 2, itr_tol
            increA_v = 1.0d0/(i)*this%matmul(t*increA_v)
            expA_v = expA_v + increA_v

            if (dot_product(increA_v, increA_v) < tol) exit
         end do
      end if
   end function

! ###################################################
   function tensor_exp_sqrt_crs(this, v, tol, coeff, itrmax, binomial) result(exp_sqrtA_v)
      class(CRS_), intent(in) :: this
      real(real64), intent(in) :: v(:)
      real(real64), allocatable :: dv(:), exp_sqrtA_v(:)
      logical, optional, intent(in) :: binomial
      complex(real64), optional, intent(in) :: coeff
      integer(int32) :: i
      integer(int32), intent(in) :: itrmax
      real(real64), intent(in) :: tol
      complex(real64) :: coeffi

      if (present(coeff)) then
         coeffi = coeff
      else
         coeffi = 1.0d0
      end if

      dv = v
      exp_sqrtA_v = zeros(size(v))
      exp_sqrtA_v = exp_sqrtA_v + dv

      ! 1-st order term
      dv = coeffi*this%tensor_sqrt(v=dv, tol=tol, itrmax=itrmax, binomial=binomial)

      exp_sqrtA_v = exp_sqrtA_v + dv

      do i = 2, itrmax
         if (i == 1) then
            cycle
         end if

         dv = 1.0d0/(i)*coeffi*this%tensor_sqrt(v=dv, tol=tol, itrmax=itrmax, binomial=binomial)
         exp_sqrtA_v = exp_sqrtA_v + dv

         if (dot_product(dv, dv) < tol) exit
      end do

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

! ###################################################
   function tensor_sqrt_crs(this, v, tol, itrmax, binomial, r) result(sqrtA_v)
      class(CRS_), intent(in) :: this
      real(real64), intent(in) :: v(:)
      real(real64), allocatable :: dv(:), sqrtA_v(:), logA_v(:)
      integer(int32) :: k
      integer(int32), intent(in) :: itrmax
      real(real64), intent(in) :: tol
      real(real64) :: coeff, nCr
      logical, optional, intenT(in) :: binomial
      real(real64), optional, intent(in) :: r

      if (present(binomial)) then
         if (binomial) then

            !  二項係数によるsqrt(A)の近似
            if (present(r)) then
               coeff = 1.0d0/maxval(this%val*r)
            else
               coeff = 1.0d0/maxval(this%val*5)
            end if

            sqrtA_v = 0.0d0*v
            dv = v
            do k = 0, itrmax
               nCr = gamma(0.50d0 + 1.0d0)/gamma(k + 1.0d0)/gamma(0.50d0 - k + 1)
               sqrtA_v = sqrtA_v + nCr*dv
               dv = (-dv + coeff*this%matmul(dv))
               if (norm(dv)/size(dv) < tol) exit
            end do
            sqrtA_v = sqrtA_v/sqrt(coeff)
            return
         end if
      end if

      dv = v
      sqrtA_v = zeros(size(v))
      sqrtA_v = sqrtA_v + dv
      ! 1-st order term
      dv = 0.50d0*this%tensor_log(v=dv, tol=tol, itrmax=itrmax)
      sqrtA_v = sqrtA_v + dv

      do k = 2, itrmax
         if (k == 1) then
            cycle
         end if

         dv = 1.0d0/dble(k)*0.50d0*this%tensor_log(v=dv, tol=tol, itrmax=itrmax)
         sqrtA_v = sqrtA_v + dv

         if (dot_product(dv, dv) < tol) exit
      end do

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

! ###################################################
   function tensor_log_crs_modified_Sparse(this, v, tol, itrmax, r) result(sqrtA_v)
      class(CRS_), intent(in) :: this
      real(real64), intent(in) :: v(:)
      real(real64), allocatable :: dv(:), sqrtA_v(:)
      integer(int32) :: k, i
      integer(int32), intent(in) :: itrmax
      real(real64), intent(in) :: tol
      real(real64), optional, intent(in) :: r

      !real(real64) :: c
      real(real64), allocatable :: c(:), lambda_mat(:, :), X(:, :), lambda(:)
      real(real64) :: c_val

      type(CRS_) :: I_CRS

!    call I_CRS%eyes(this%size() )
!    call LAPACK_EIG(A=this,B=I_CRS,x=x,lambda=lambda)
!    lambda_mat = zeros(this%size(),this%size())
!    do i=1,this%size()
!        lambda_mat(i,i) = log(lambda(i))
!    enddo
!
!    ! shor algorithm
!    sqrtA_v = matmul(matmul( x, matmul(lambda_mat,transpose(x)) ),v )
!    return

      !c = maxval(this%val)*1.0d0
      !c = this%diag()*10.0d0 ! テイラー展開中心をdiagで決める．
      if (present(r)) then
         c_val = r
      else
         c_val = maxval(abs(this%val))
      end if
      !c = ones(this%size() )
      !c = 1.0d0
      ! k = 0
      ! Martin, J. F. P. On the exponential representation of solutions of linear differential equations. J. Differ. Equ. 4, 257–279 (1968).
      ! ln(this) = Ln(this) = ln(exp(C) X ) - C
      sqrtA_v = -(c_val)*v

      do k = 1, itrmax
         if (k == 1) then
            !dv = k/c_val*(-c_val*v + this%matmul(v)  )
            dv = 1.0d0/k*(exp(c_val)*this%matmul(v) - v)
         else
            dv = (-1.0d0)*dble(k)/dble(k + 1)*(exp(c_val)*this%matmul(dv) - dv)
         end if
         sqrtA_v = sqrtA_v + dv

         if (norm(abs(dv)) < tol) return
      end do

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

! ###################################################
   function tensor_log_crs(this, v, tol, itrmax, r) result(sqrtA_v)

      ! wikipedia's difinition
      ! it is accurate only when "this" is a diagonally dominant matrix
      class(CRS_), intent(in) :: this
      real(real64), intent(in) :: v(:)
      real(real64), allocatable :: dv(:), sqrtA_v(:)
      integer(int32) :: k
      integer(int32), intent(in) :: itrmax
      real(real64), intent(in) :: tol
      !real(real64) :: c
      real(real64) :: c
      real(real64), optional, intent(in) :: r
      !real(real64) :: c
      !c = maxval(this%diag())
      !c = maxval(this%val)*1.0d0
      c = maxval(abs(this%diag()))*2.0d0 ! テイラー展開中心をdiagで決める．
      !c = 1.0d0
      ! k = 0
      if (present(r)) then
         c = r
      end if

      sqrtA_v = log(c)*v

      do k = 1, itrmax
         if (k == 1) then
            dv = (1.0d0)/k*(v - 1.0d0/c*this%matmul(v))
         else
            dv = (1.0d0)*dble(k)/dble(k + 1)*(dv - 1.0d0/c*this%matmul(dv))
         end if
         sqrtA_v = sqrtA_v - dv

         if (norm(abs(dv)) < tol) return
      end do

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

! #####################################################
   subroutine fixCRS(this, idx, val, RHS, only_row)
      class(CRS_), intent(inout) :: this
      integer(int32), intent(in) :: idx(:)
      real(real64), optional, intent(inout) :: RHS(:)
      real(real64), intent(in) :: val(:)
      integer(int64) :: i, j, k, id

      logical, optional, intent(in) :: only_row

      if (present(only_row)) then
         if (only_row) then

            do i = 1, size(idx)
               id = idx(i)
               do j = this%row_ptr(id), this%row_ptr(id + 1) - 1
                  if (this%col_idx(j) == id) then
                     this%val(j) = 1.0d0
                  else
                     this%val(j) = 0.0d0
                  end if
               end do
            end do
            return
         end if
      end if

      do i = 1, size(idx)
         id = idx(i)
         do j = 1, size(this%row_ptr) - 1
            if (j == id) then
               do k = this%row_ptr(j), this%row_ptr(j + 1) - 1
                  if (this%col_idx(k) == id) then
                     this%val(k) = 0.0d0
                  else
                     this%val(k) = 0.0d0
                  end if
               end do
            else
               do k = this%row_ptr(j), this%row_ptr(j + 1) - 1
                  if (this%col_idx(k) == id) then
                     RHS(j) = RHS(j) - this%val(k)*val(i)
                     this%val(k) = 0.0d0
                  end if
               end do
            end if
         end do
      end do

      do i = 1, size(idx)
         id = idx(i)
         do j = 1, size(this%row_ptr) - 1
            if (j == id) then
               do k = this%row_ptr(j), this%row_ptr(j + 1) - 1
                  if (this%col_idx(k) == id) then
                     this%val(k) = 1.0d0
                  end if
               end do
            end if
         end do
      end do

      do i = 1, size(idx)
         RHS(id) = val(i)
      end do
   end subroutine
! #####################################################

! #####################################################
   subroutine fix_zeroCRS(this, idx)
      class(CRS_), intent(inout) :: this
      integer(int32), intent(in) :: idx(:)

      integer(int64) :: i, j, k, id
      do i = 1, size(idx)
         id = idx(i)
         do j = 1, size(this%row_ptr) - 1
            if (j == id) then
               do k = this%row_ptr(j), this%row_ptr(j + 1) - 1
                  if (this%col_idx(k) == id) then
                     this%val(k) = 0.0d0
                  else
                     this%val(k) = 0.0d0
                  end if
               end do
            else
               do k = this%row_ptr(j), this%row_ptr(j + 1) - 1
                  if (this%col_idx(k) == id) then

                     this%val(k) = 0.0d0
                  end if
               end do
            end if
         end do
      end do

      do i = 1, size(idx)
         id = idx(i)
         do j = 1, size(this%row_ptr) - 1
            if (j == id) then
               do k = this%row_ptr(j), this%row_ptr(j + 1) - 1
                  if (this%col_idx(k) == id) then
                     this%val(k) = 1.0d0
                  end if
               end do
            end if
         end do
      end do

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

! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> complex >>>>>>>>>>>>>>>>>>>>>>>

! ###################################################
   function tensor_exponential_complex64_crs(this, itr_tol, tol, x, dt, fix_idx, fix_val) result(expA_v)
      class(CRS_), intent(in) :: this
      complex(real64), intent(in) :: x(:)
      real(real64), optional, intent(in) :: dt
      complex(real64), allocatable :: expA_v(:), increA_v(:), bhat(:)

      real(real64), intent(in) :: tol
      integer(int32), intent(in)::itr_tol
!    complex(real64),allocatable::increA(:,:)
      real(real64)   :: increment, NN, t
      integer(int32) :: i, j, n

      complex(real64), optional, intent(in) :: fix_val(:)
      integer(int32), optional, intent(in)::fix_idx(:)
      type(CRS_) :: Amatrix
!    if(.not. allocated(expA) )allocate(expA(size(A,1),size(A,2) ))
!    allocate(increA(size(A,1),size(A,2) ))
!    if(size(A,1)/=size(A,2)) stop "ERROR:tensor exp is not a square matrix"

!    expA_v = x
!    do n=1,size(expA,1)
!        expA(n,n)=1.0d0
!    enddo
      t = input(default=1.0d0, option=dt)
      ! exp(At) = I + At + 1/2*(At)^2 + 1/6*(At)^3 + ....
      ! exp(At)v = Iv + (At)v + 1/2*(At)^2 v + 1/6*(At)^3 v+ ....
      ! exp(At)v = v + (At)v + 1/2*(At)^2 v + 1/6*(At)^3 v+ ....
!    expA_v = x
!    do i=1,itr_tol
!        increA_v = x
!        do j=1,i
!            increA_v = t/dble(j)*this%matmul(increA_v)
!        enddo
!        expA_v = expA_v + increA_v
!        if(dot_product(increA_v,increA_v) < tol )exit
!
!    enddo
!    return

      ! i==1,2
      ! 0-th order term
      if (present(fix_idx)) then
         Amatrix = this
         bhat = zeros(this%size())
         call Amatrix%fix(idx=fix_idx, RHS=bhat, val=fix_val)

         increA_v = x
         expA_v = zeros(size(x))
         expA_v = expA_v + increA_v
         ! 1-st order term
         increA_v = Amatrix%matmul(t*increA_v) - bhat*dt
         expA_v = expA_v + increA_v

         do i = 2, itr_tol
            if (i == 1) then
               cycle
            end if

            increA_v = 1.0d0/(i)*Amatrix%matmul(t*increA_v) !- bhat*dt
            expA_v = expA_v + increA_v

            if (abs(dot_product(increA_v, increA_v)) < tol) exit
         end do
      else

         increA_v = x
         expA_v = zeros(size(x))
         expA_v = expA_v + increA_v
         ! 1-st order term
         increA_v = this%matmul(t*increA_v)
         expA_v = expA_v + increA_v

         do i = 2, itr_tol
            if (i == 1) then

               cycle
            end if

            increA_v = 1.0d0/(i)*this%matmul(t*increA_v)
            expA_v = expA_v + increA_v

            if (abs(dot_product(increA_v, increA_v)) < tol) exit
         end do
      end if
   end function

! ###################################################
   function tensor_exp_sqrt_complex64_crs(this, v, tol, itrmax, coeff, fix_idx, binomial, r) result(exp_sqrtA_v)
      class(CRS_), intent(in) :: this
      complex(real64), intent(in) :: v(:)
      complex(real64), optional, intent(in) :: coeff
      complex(real64), allocatable :: dv(:), exp_sqrtA_v(:), bhat(:), buf(:)
      integer(int32) :: i, k
      integer(int32), intent(in) :: itrmax
      real(real64), intent(in) :: tol
      logical, optional, intent(in) :: binomial
      real(real64), optional, intent(in) :: r

      !real(real64),optional,intent(in) :: fix_val(:) ! only 0-fix is available
      integer(int32), optional, intent(in)::fix_idx(:)
      complex(real64) :: coeffi

      type(CRS_) :: Amatrix
      type(Math_) :: math

      if (present(coeff)) then
         coeffi = coeff
      else
         coeffi = 1.0d0
      end if

      if (present(fix_idx)) then
         !Amatrix = this
         !bhat = zeros(this%size() )
         !call Amatrix%fix(idx=fix_idx,RHS=bhat,val=fix_val)

         !increA_v = x
         !expA_v = zeros(size(x) )
         !expA_v = expA_v + increA_v
         ! 1-st order term
         !increA_v = Amatrix%matmul(t*increA_v) - bhat*dt
         !expA_v = expA_v + increA_v

         !do i=2,itr_tol
         !    if(i==1)then
         !        cycle
         !    endif
         !    increA_v = 1.0d0/(i)*Amatrix%matmul(t*increA_v) !- bhat*dt
         !    expA_v = expA_v + increA_v
         !
         !    if(abs(dot_product(increA_v,increA_v)) < tol )exit
         !enddo

         Amatrix = this

         dv = zeros(size(v) - size(fix_idx))

         k = 0
         do i = 1, size(v)
            if (i.in.fix_idx) then
               cycle
            else
               k = k + 1
               dv(k) = v(i)
            end if
         end do

         !dv(fix_idx(:) ) = fix_val(:)!/exp(1.0d0)

         !call Amatrix%fix(idx=fix_idx,RHS=bhat,val=fix_val+0.0d0*math%i)

         !call Amatrix%fix(idx=fix_idx,RHS=bhat,val=0.0d0*zeros(size(fix_idx) ) +0.0d0*math%i)

         call Amatrix%remove(idx=fix_idx) ! remove row and columns

         !call Amatrix%fix(idx=fix_idx,val=fix_val+0.0d0*math%i,only_row=.true.)
         !call Amatrix%fix(idx=fix_idx,RHS=dv,val=fix_val+0.0d0*math%i)
         !dv = dv + bhat
         !<test>
         !dv(fix_idx) = fix_val

         ! ignore row and columns
         exp_sqrtA_v = zeros(size(v) - size(fix_idx))
         exp_sqrtA_v = exp_sqrtA_v + dv
         !exp_sqrtA_v(fix_idx(:) ) = fix_val(:)

         ! 1-st order term
         dv = coeffi*Amatrix%tensor_sqrt(v=dv, tol=tol, itrmax=itrmax, binomial=binomial, r=r)
         exp_sqrtA_v = exp_sqrtA_v + dv
         !exp_sqrtA_v(fix_idx) = fix_val

         do i = 2, itrmax

            dv = 1.0d0/(i)*coeffi*Amatrix%tensor_sqrt(v=dv, tol=tol, itrmax=itrmax, binomial=binomial, r=r)
            !dv(fix_idx) = 0.0d0
            exp_sqrtA_v = exp_sqrtA_v + dv
            !exp_sqrtA_v(fix_idx) = fix_val
            if (abs(dot_product(dv, dv)) < tol) exit
         end do
         !exp_sqrtA_v(fix_idx(:) ) = fix_val(:)
         buf = exp_sqrtA_v
         exp_sqrtA_v = 0.0d0*v

         k = 0
         do i = 1, size(v)
            if (i.in.fix_idx) then
               exp_sqrtA_v(i) = 0.0d0
            else
               k = k + 1
               exp_sqrtA_v(i) = buf(k)
            end if
         end do

      else

         dv = v
         exp_sqrtA_v = zeros(size(v))
         exp_sqrtA_v = exp_sqrtA_v + dv

         ! 1-st order term
         dv = coeffi*this%tensor_sqrt(v=dv, tol=tol, itrmax=itrmax, binomial=binomial, r=r)
         exp_sqrtA_v = exp_sqrtA_v + dv

         do i = 2, itrmax

            dv = 1.0d0/(i)*coeffi*this%tensor_sqrt(v=dv, tol=tol, itrmax=itrmax, binomial=binomial, r=r)
            exp_sqrtA_v = exp_sqrtA_v + dv

            if (abs(dot_product(dv, dv)) < tol) exit
         end do
      end if

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

! ###################################################
   function tensor_sqrt_complex64_crs(this, v, tol, itrmax, binomial, r) result(sqrtA_v)
      class(CRS_), intent(in) :: this
      complex(real64), intent(in) :: v(:)
      complex(real64), allocatable :: dv(:), sqrtA_v(:), logA_v(:)
      integer(int32) :: k
      integer(int32), intent(in) :: itrmax
      real(real64), intent(in) :: tol
      real(real64) :: coeff, nCr
      logical, optional, intenT(in) :: binomial
      real(real64), optional, intent(in) :: r
      if (present(binomial)) then
         if (binomial) then

            !  二項係数によるsqrt(A)の近似
            if (present(r)) then
               coeff = 1.0d0/maxval(this%val*r)
            else
               coeff = 1.0d0/maxval(this%val*5)
            end if

            sqrtA_v = 0.0d0*v
            dv = v
            do k = 0, itrmax
               nCr = gamma(0.50d0 + 1.0d0)/gamma(k + 1.0d0)/gamma(0.50d0 - k + 1)
               sqrtA_v = sqrtA_v + nCr*dv
               dv = (-dv + coeff*this%matmul(dv))
               if (norm(dble(dv))/size(dv) < tol) exit
            end do
            sqrtA_v = sqrtA_v/sqrt(coeff)
            return
         end if
      end if

      dv = v
      sqrtA_v = zeros(size(v))
      sqrtA_v = sqrtA_v + dv
      ! 1-st order term
      dv = 0.50d0*this%tensor_log(v=dv, tol=tol, itrmax=itrmax)
      sqrtA_v = sqrtA_v + dv

      do k = 2, itrmax
         if (k == 1) then

            cycle
         end if

         dv = 1.0d0/dble(k)*0.50d0*this%tensor_log(v=dv, tol=tol, itrmax=itrmax)
         sqrtA_v = sqrtA_v + dv

         if (abs(dot_product(dv, dv)) < tol) exit
      end do

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

! ###################################################
   function tensor_log_complex64_crs(this, v, tol, itrmax) result(sqrtA_v)
      class(CRS_), intent(in) :: this
      complex(real64), intent(in) :: v(:)
      complex(real64), allocatable :: dv(:), sqrtA_v(:)
      integer(int32) :: k
      integer(int32), intent(in) :: itrmax
      real(real64), intent(in) :: tol
      !complex(real64) :: c
      real(real64), allocatable :: c(:)

      !c = maxval(this%val)*1.
      c = this%diag() ! テイラー展開中心をdiagで決める．
      !c = maxval(this%diag())*3.0d0*ones(size(v)) ! テイラー展開中心をdiagで決める．
      !c = maxval(this%diag(cell_centered=.true.))*ones(size(v)) ! テイラー展開中心をdiagで決める．

      !c = 1.0d0
      ! k = 0
      sqrtA_v = zeros(size(v))
      sqrtA_v = log(c)*v

      do k = 1, itrmax
         if (k == 1) then
            dv = k/c*(-c*v + this%matmul(v))
         else
            dv = (-1.0d0)/c*dble(k)/dble(k + 1)*(-c*dv + this%matmul(dv))
         end if
         sqrtA_v = sqrtA_v + dv

         if (norm(abs(dv)) < tol) return
      end do

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

! ###################################################
   function tensor_sinc_sqrt_complex64_crs(this, v, tol, itrmax, coeff, debug, fix_idx) result(retA_v)
      class(CRS_), intent(in) :: this
      complex(real64), intent(in) :: v(:)
      complex(real64), optional, intent(in) :: coeff
      complex(real64), allocatable :: dv(:), retA_v(:)
      integer(int32), optional, intent(in) :: fix_idx(:)
      integer(int32) :: k
      integer(int32), intent(in) :: itrmax
      logical, optional, intent(in) :: debug
      real(real64), intent(in) :: tol
      type(Math_) :: math
      complex(real64) :: a, b

      if (present(coeff)) then
         a = coeff
      else
         a = 1.0d0
      end if
      ! sinc(sqrt(A) ) = \sum_{n=0}^{\infty}
      !c = 1.0d0
      ! k = 0

      !k=0
      k = 0
      retA_v = v

      ! zero-fixed Dirichlet boundary
      if (present(fix_idx)) then
         retA_v(fix_idx) = 0
      end if

      ! k=1
      k = 1
      dv = a*a*this%matmul(v)

      ! zero-fixed Dirichlet boundary
      if (present(fix_idx)) then
         dv(fix_idx) = 0
      end if

      retA_v = retA_v + ((-1.0d0)**k)*(sqrt(math%pi)*2.0d0**(-1 - 2*k)) &
               /gamma(1.0d0 + k)/gamma(1.5d0 + k)*dv

      do k = 2, itrmax

         dv = a*a*this%matmul(dv)

         ! zero-fixed Dirichlet boundary
         if (present(fix_idx)) then
            dv(fix_idx) = 0
         end if
         b = ((-1.0d0)**k)*(sqrt(math%pi)*2.0d0**(-1 - 2*k)) &
             /gamma(1.0d0 + k)/gamma(1.5d0 + k)
         retA_v = retA_v + b*dv

         if (norm(abs(b*dv)) < tol) return

         if (present(debug)) then
            if (debug) then
               print *, k, norm(abs(b*dv))
            end if
         end if

      end do

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

! ###################################################
   function tensor_sinc_sqrt_real64_crs(this, v, tol, itrmax, coeff, debug, fix_idx) result(retA_v)
      class(CRS_), intent(in) :: this
      real(real64), intent(in) :: v(:)
      real(real64), optional, intent(in) :: coeff
      real(real64), allocatable :: dv(:), retA_v(:)
      integer(int32), optional, intent(in) :: fix_idx(:)
      integer(int32) :: k
      integer(int32), intent(in) :: itrmax
      logical, optional, intent(in) :: debug
      real(real64), intent(in) :: tol
      type(Math_) :: math
      real(real64) :: a, b

      if (present(coeff)) then
         a = coeff
      else
         a = 1.0d0
      end if
      ! sinc(sqrt(A) ) = \sum_{n=0}^{\infty}
      !c = 1.0d0
      ! k = 0

      !k=0
      k = 0
      retA_v = v

      ! zero-fixed Dirichlet boundary
      if (present(fix_idx)) then
         retA_v(fix_idx) = 0
      end if

      ! k=1
      k = 1
      dv = a*a*this%matmul(v)

      ! zero-fixed Dirichlet boundary
      if (present(fix_idx)) then
         dv(fix_idx) = 0
      end if

      retA_v = retA_v + ((-1.0d0)**k)*(sqrt(math%pi)*2.0d0**(-1 - 2*k)) &
               /gamma(1.0d0 + k)/gamma(1.5d0 + k)*dv

      do k = 2, itrmax

         dv = a*a*this%matmul(dv)

         ! zero-fixed Dirichlet boundary
         if (present(fix_idx)) then
            dv(fix_idx) = 0
         end if
         b = ((-1.0d0)**k)*(sqrt(math%pi)*2.0d0**(-1 - 2*k)) &
             /gamma(1.0d0 + k)/gamma(1.5d0 + k)
         retA_v = retA_v + b*dv

         if (norm(abs(b*dv)) < tol) return

         if (present(debug)) then
            if (debug) then
               print *, k, norm(abs(b*dv))
            end if
         end if

      end do

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

   function tensor_wave_kernel_complex_64_crs(this, u0, v0, tol, itrmax, h, t, fix_idx, debug) result(u)
      class(CRS_), intent(in) :: this
      complex(real64), intent(in) :: u0(:), v0(:)
      complex(real64), allocatable:: Adu(:), Adv(:), u(:)
      complex(real64) :: cos_coeff, sinc_coeff
      real(real64), intent(in) :: h, t
      logical, optional, intent(in) :: debug

      integer(int32), optional, intent(in) :: itrmax
      real(real64), optional, intent(in) :: tol
      integer(int32), optional, intent(in) :: fix_idx(:)

      integer(int32) :: itr_max = 100
      real(real64)   :: itr_tol = dble(1.0e-16)
      logical :: debug_mode
      type(Math_) :: math
      integer(int32) :: n

      if (present(itrmax)) then
         itr_max = itrmax
      end if

      if (present(tol)) then
         itr_tol = tol
      end if

      ! a + 2 h M^{-1} v + M^{-1} K u = 0
      ! u(t) = exp(-ht)( cos(t*sqrt(M^{-1} K - h^2 I)  ) u
      !      + t*sinc( t*sqrt(M^{-1} K - h^2 I) ) v

      u = exp(-h*t)*this%tensor_cos_sqrt(v=u0, tol=itr_tol, itrmax=itr_max, &
                                         coeff=t + 0*math%i, debug=debug, fix_idx=fix_idx) &
          + exp(-h*t)*t*this%tensor_sinc_sqrt(v=v0, tol=itr_tol, itrmax=itr_max, &
                                              coeff=t + 0*math%i, debug=debug, fix_idx=fix_idx)

   end function

   function tensor_wave_kernel_real_64_crs(this, u0, v0, tol, itrmax, h, t, fix_idx, debug) result(u)
      class(CRS_), intent(inout) :: this
      real(real64), intent(in) :: u0(:), v0(:)
      real(real64), intent(in) :: h, t
      logical, optional, intent(in) :: debug

      integer(int32), optional, intent(in) :: itrmax
      real(real64), optional, intent(in) :: tol
      integer(int32), optional, intent(in) :: fix_idx(:)

      real(real64), allocatable:: Adu(:), Adv(:), u(:)
      real(real64) :: cos_coeff, sinc_coeff

      integer(int32) :: itr_max = 100
      real(real64)   :: itr_tol = dble(1.0e-16)
      logical :: debug_mode
      type(Math_) :: math
      integer(int32) :: n

      if (present(itrmax)) then
         itr_max = itrmax
      end if

      if (present(tol)) then
         itr_tol = tol
      end if

      if (present(fix_idx)) then
         call this%fix(fix_idx)
      end if
      ! a + 2 h M^{-1} v + M^{-1} K u = 0
      ! u(t) = exp(-ht)( cos(t*sqrt(M^{-1} K - h^2 I)  ) u
      !      + t*sinc( t*sqrt(M^{-1} K - h^2 I) ) v

      u = exp(-h*t)*this%tensor_cos_sqrt(v=u0, tol=itr_tol, itrmax=itr_max, &
                                         coeff=t, debug=debug, fix_idx=fix_idx) &
          + exp(-h*t)*t*this%tensor_sinc_sqrt(v=v0, tol=itr_tol, itrmax=itr_max, &
                                              coeff=t, debug=debug, fix_idx=fix_idx)

   end function

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

   function tensor_wave_kernel_LPF_real_64_crs(this, u0, v0, tol, itrmax, h, t, fix_idx, debug, cutoff_frequency) result(u)
      class(CRS_), intent(in) :: this
      real(real64), intent(in) :: u0(:), v0(:)
      real(real64), allocatable:: Adu(:), Adv(:), u(:)
      real(real64) :: cos_coeff, sinc_coeff
      real(real64), intent(in) :: h, t, cutoff_frequency
      logical, optional, intent(in) :: debug

      integer(int32), optional, intent(in) :: itrmax
      real(real64), optional, intent(in) :: tol
      integer(int32), optional, intent(in) :: fix_idx(:)

      integer(int32) :: itr_max = 100
      real(real64)   :: itr_tol = dble(1.0e-16)
      logical :: debug_mode
      type(Math_) :: math
      integer(int32) :: n

      if (present(itrmax)) then
         itr_max = itrmax
      end if

      if (present(tol)) then
         itr_tol = tol
      end if

      ! a + 2 h M^{-1} v + M^{-1} K u = 0
      ! u(t) = exp(-ht)( cos(t*sqrt(M^{-1} K - h^2 I)  ) u
      !      + t*sinc( t*sqrt(M^{-1} K - h^2 I) ) v

      u = exp(-h*t)*this%tensor_cos_sqrt_LPF(v=u0, tol=itr_tol, itrmax=itr_max, &
                                             coeff=t, debug=debug, fix_idx=fix_idx, cutoff_frequency=cutoff_frequency) &
          + exp(-h*t)*this%tensor_t_sinc_sqrt_LPF(v=v0, tol=itr_tol, itrmax=itr_max, &
                                                  coeff=t, debug=debug, fix_idx=fix_idx, cutoff_frequency=cutoff_frequency)

   end function

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

   function tensor_wave_kernel_RHS_real_64_crs(this, RHS, t, tol, itrmax, fix_idx, debug) result(u)
      class(CRS_), intent(in) :: this
      real(real64), intent(in) :: RHS(:)
      real(real64), allocatable:: u(:), du(:)
      real(real64) :: cos_coeff, sinc_coeff
      real(real64), intent(in) :: t
      logical, optional, intent(in) :: debug

      integer(int32), optional, intent(in) :: itrmax
      real(real64), optional, intent(in) :: tol
      integer(int32), optional, intent(in) :: fix_idx(:)

      integer(int32) :: itr_max = 100
      real(real64)   :: itr_tol = dble(1.0e-16)
      logical :: debug_mode
      type(Math_) :: math
      integer(int32) :: n

      if (present(itrmax)) then
         itr_max = itrmax
      end if

      if (present(tol)) then
         itr_tol = tol
      end if

      ! Force-induced displacement
      ! RHS is constant for interval [0, t]

      !u = -1.0d0/2.0d0*t*t*RHS
      u = 0.0d0*RHS
      du = -t*t/2.0d0*RHS
      u = u + du
      do n = 1, itr_max
         du = -t*t/dble(2*n + 1)/dble(2*n + 2)*this%matmul(du)
         u = u + du
         if (present(debug)) then
            if (debug) then
               print *, n, norm(du)
            end if
         end if
         if (norm(du) < itr_tol) exit
      end do

   end function

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

   function tensor_wave_kernel_RHS_LPF_real_64_crs(this, RHS, t, tol, itrmax, fix_idx, debug, cutoff_frequency, weights) result(u)
      class(CRS_), intent(in) :: this
      real(real64), intent(in) :: RHS(:), cutoff_frequency
      real(real64), allocatable:: u(:), du(:)
      real(real64) :: cos_coeff, sinc_coeff
      real(real64), intent(in) :: t
      logical, optional, intent(in) :: debug

      integer(int32), optional, intent(in) :: itrmax
      real(real64), optional, intent(in) :: tol
      integer(int32), optional, intent(in) :: fix_idx(:)
      real(real64), optional, intent(in) :: weights(:)

      integer(int32) :: itr_max = 100
      real(real64)   :: itr_tol = dble(1.0e-16)
      real(real64)   :: ddt, hat_t
      logical :: debug_mode
      type(Math_) :: math
      integer(int32) :: n, m
      real(real64), allocatable :: w(:)

      integer(int32) :: lbd, ubd
      

      if (present(weights)) then
         ubd = size(weights) - (size(weights) + 1)/2
         lbd = -(size(weights) - (size(weights) + 1)/2)
         allocate (w(lbd:ubd))
         w(lbd:ubd) = weights(1:)
      else
         allocate (w(-1:1))
         w(-1) = 0.250d0
         w(0) = 0.500d0
         w(1) = 0.250d0
      end if

      if (present(itrmax)) then
         itr_max = itrmax
      end if

      if (present(tol)) then
         itr_tol = tol
      end if

      ! Force-induced displacement
      ! RHS is constant for interval [0, t]

      !ddt = 1.0d0/4.0d0/cutoff_frequency
      ddt = 1.0d0/cutoff_frequency/2.0d0/math%pi*acos(sqrt(2.0d0) - 1.0d0)

      
      !u = -1.0d0/2.0d0*t*t*RHS

      u = 0.0d0*RHS

      hat_t = 0.0d0
      !print *,"dbg",lbound(w, 1), ubound(w, 1) 
      do m = lbound(w, 1), ubound(w, 1)
         hat_t = hat_t + w(m)*(t + m*ddt)**2
      end do
      !hat_t = 0.250d0*(t-ddt)**2 &
      !    + 0.50d0*(t)**2 &
      !    + 0.250d0*(t+ddt)**2

      du = -1.0d0/2.0d0*RHS
      u = u + hat_t*du


      

      do n = 1, itr_max
         
         hat_t = 0.0d0
         do m = lbound(w, 1), ubound(w, 1)
            hat_t = hat_t + w(m)*(t + m*ddt)**(2*n + 2)
         end do

         !hat_t = 0.250d0*(t-ddt)**(2*n+2) &
         !     + 0.50d0*(t)**(2*n+2)&
         !     + 0.250d0*(t+ddt)**(2*n+2)

         du = -1.0d0/dble(2*n + 1)/dble(2*n + 2)*this%matmul(du)
         u = u + hat_t*du
         
         if (present(debug)) then
            if (debug) then
               print *, n, norm(hat_t*du)
            end if
         end if

         if (norm(hat_t*du) < itr_tol) exit
      end do

   end function

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

   function tensor_wave_kernel_RHS_complex_64_crs(this, RHS, t, tol, itrmax, fix_idx, debug) result(u)
      class(CRS_), intent(in) :: this
      complex(real64), intent(in) :: RHS(:)
      complex(real64), allocatable:: u(:), u_n(:), du(:)

      real(real64), intent(in) :: t
      logical, optional, intent(in) :: debug

      integer(int32), optional, intent(in) :: itrmax
      real(real64), optional, intent(in) :: tol
      integer(int32), optional, intent(in) :: fix_idx(:)

      integer(int32) :: itr_max = 100
      real(real64)   :: itr_tol = dble(1.0e-16)
      logical :: debug_mode
      type(Math_) :: math
      integer(int32) :: n

      if (present(itrmax)) then
         itr_max = itrmax
      end if

      if (present(tol)) then
         itr_tol = tol
      end if

      ! Force-induced displacement
      ! RHS is constant for interval [0, t]
      u = 0.0d0*RHS
      du = -t*t/2.0d0*RHS
      u = u + du
      do n = 1, itr_max
         du = -t*t/dble(2*n + 1)/dble(2*n + 2)*this%matmul(du)
         u = u + du
         if (present(debug)) then
            if (debug) then
               print *, n, sqrt(dble(dot_product(du, du)))
            end if
         end if
         if (sqrt(dble(dot_product(du, du))) < itr_tol) exit
      end do

   end function

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

! ###################################################
   function tensor_d1_wave_kernel_complex64_crs(this, u, v, tol, itrmax, coeff) result(retA_v)
      class(CRS_), intent(in) :: this
      complex(real64), intent(in) :: u(:), v(:)
      complex(real64), optional, intent(in) :: coeff
      complex(real64), allocatable :: dv(:), retA_v(:), du(:)
      integer(int32) :: k
      integer(int32), intent(in) :: itrmax
      real(real64), intent(in) :: tol
      type(Math_) :: math
      complex(real64) :: a

      ! first derivative of wave kernel function
      ! v - t u w^2 - 1/2 t^2 (v w^2) + 1/6 t^3 u w^4 + 1/24 t^4 v w^4 - 1/120 t^5 (u w^6) + O(t^6)
      !(テイラー級数)
      ! sinc(sqrt(A) ) = \sum_{n=0}^{\infty}
      !c = 1.0d0
      ! k = 0

      if (present(coeff)) then
         a = coeff
      else
         a = 1.0d0
      end if

      !1,2
      du = this%matmul(u)
      retA_v = v - a*du

      !3
      dv = this%matmul(v)
      retA_v = retA_v - 1.0d0/2.0d0*a*a*dv

      !4
      du = this%matmul(du)
      retA_v = retA_v + 1.0d0/6.0d0*a*a*a*du

      !5
      dv = this%matmul(dv)
      retA_v = retA_v + 1.0d0/24.0d0*a*a*a*a*dv

      !6
      du = this%matmul(du)
      retA_v = retA_v - 1.0d0/120.0d0*a*a*a*a*a*du

      ! 6th-order apploximation

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

! ###################################################
   function tensor_cos_sqrt_complex64_crs(this, v, tol, itrmax, coeff, debug, fix_idx) result(retA_v)
      class(CRS_), intent(in) :: this
      complex(real64), intent(in) :: v(:)
      complex(real64), allocatable :: dv(:), retA_v(:)
      complex(real64), optional, intent(in) :: coeff
      integer(int32), optional, intent(in) :: fix_idx(:)
      logical, optional, intent(in) :: debug
      integer(int32) :: k, row, CRS_id
      integer(int32), intent(in) :: itrmax
      real(real64), intent(in) :: tol
      type(Math_) :: math
      complex(real64) :: a

      if (present(coeff)) then
         a = coeff
      else
         a = 1.0d0
      end if
      ! sinc(sqrt(A) ) = \sum_{n=0}^{\infty}
      !c = 1.0d0
      ! k = 0

      !k=0
      k = 0
      retA_v = v
      ! zero-fixed Dirichlet boundary
      if (present(fix_idx)) then
         retA_v(fix_idx) = 0
      end if

      ! a=t

      ! k=1
      k = 1
      dv = this%matmul(v)
      dv = -0.50d0*a*a*dv

      ! zero-fixed Dirichlet boundary
      if (present(fix_idx)) then
         dv(fix_idx) = 0
      end if
      retA_v = retA_v + dv

      !

      do k = 2, itrmax
         ! k = k + 1
         dv = (a*a)*this%matmul(dv)/(2*k)/(2*k - 1)*(-1.0d0)

         ! zero-fixed Dirichlet boundary
         if (present(fix_idx)) then
            dv(fix_idx) = 0
         end if
         retA_v = retA_v + dv

         if (norm(abs(dv)) < tol) return

         if (present(debug)) then
            if (debug) then
               print *, k, norm(abs(dv))
            end if
         end if

      end do

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

! #####################################################
   function wave_kernel_hanning_coefficient(dt, num_sample, t_power) result(ret)
      real(real64), intent(in)::dt
      integer, intent(in) :: num_sample, t_power
      integer(int32) :: m, k
      real(real64) :: ret, t, period_T
      type(Math_) :: math

      k = num_sample/2
      ret = 0.0d0
      period_T = 1.0d0/dt*2
      do m = -k, k
         t = -num_sample*dt + (m - 1)*dt
         ret = ret + (0.50d0 - 0.50d0*cos(2.0d0*math%pi*(t)/period_T)) &
               *(m*dt)**(t_power)
      end do
   end function
! #####################################################

! ###################################################
   function tensor_cos_sqrt_real64_crs(this, v, tol, itrmax, coeff, debug, fix_idx) result(retA_v)
      class(CRS_), intent(in) :: this
      real(real64), intent(in) :: v(:)
      real(real64), allocatable :: dv(:), retA_v(:)
      real(real64), optional, intent(in) :: coeff
      integer(int32), optional, intent(in) :: fix_idx(:)
      logical, optional, intent(in) :: debug
      integer(int32) :: k, row, CRS_id
      integer(int32), intent(in) :: itrmax
      real(real64), intent(in) :: tol
      type(Math_) :: math
      real(real64) :: a

      if (present(coeff)) then
         a = coeff
      else
         a = 1.0d0
      end if
      ! sinc(sqrt(A) ) = \sum_{n=0}^{\infty}
      !c = 1.0d0
      ! k = 0

      !k=0
      k = 0
      retA_v = v
      ! zero-fixed Dirichlet boundary
      if (present(fix_idx)) then
         retA_v(fix_idx) = 0
      end if

      if (present(debug)) then
         if (debug) then
            print *, k, norm(abs(v))

         end if
      end if

      ! k=1
      k = 1
      dv = this%matmul(v)
      dv = -0.50d0*a*a*dv

      ! zero-fixed Dirichlet boundary
      if (present(fix_idx)) then
         dv(fix_idx) = 0
      end if
      retA_v = retA_v + dv

      if (present(debug)) then
         if (debug) then
            print *, k, norm(abs(dv))
         end if
      end if

      do k = 2, itrmax
         ! k = k + 1
         dv = a*a*this%matmul(dv)/(2*k)/(2*k - 1)*(-1.0d0)

         ! zero-fixed Dirichlet boundary
         if (present(fix_idx)) then
            dv(fix_idx) = 0
         end if
         retA_v = retA_v + dv

         if (norm(abs(dv)) < tol) return

         if (present(debug)) then
            if (debug) then
               print *, k, norm(abs(dv))
            end if
         end if

      end do

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

! ###################################################
   function tensor_cos_sqrt_LPF_real64_crs(this, v, tol, itrmax, coeff, debug, fix_idx, cutoff_frequency) result(retA_v)
      class(CRS_), intent(in) :: this
      real(real64), intent(in) :: v(:)
      real(real64), allocatable :: dv(:), retA_v(:)
      real(real64), optional, intent(in) :: coeff
      integer(int32), optional, intent(in) :: fix_idx(:)

      logical, optional, intent(in) :: debug
      integer(int32) :: k, row, CRS_id, num_hanning_sample
      integer(int32), intent(in) :: itrmax
      real(real64), intent(in) :: tol, cutoff_frequency
      type(Math_) :: math
      real(real64) :: a, dt, b, ddt

      !ddt = 1.0d0/(cutoff_frequency/4.0d0)
      !ddt = 1.0d0/cutoff_frequency/math%pi*acos( (sqrt(2.0d0)-1.0d0)/2.0d0 )
      ddt = 1.0d0/cutoff_frequency/2.0d0/math%pi*acos(sqrt(2.0d0) - 1.0d0)
      if (present(coeff)) then
         a = coeff
      else
         a = 1.0d0
      end if
      ! sinc(sqrt(A) ) = \sum_{n=0}^{\infty}
      !c = 1.0d0
      ! k = 0
      dt = a
      dv = v
      !k=0
      k = 0
      retA_v = v
      ! zero-fixed Dirichlet boundary
      if (present(fix_idx)) then
         retA_v(fix_idx) = 0
      end if

      do k = 1, itrmax
         ! k = k + 1

         !dv = a*a*this%matmul(dv)/(2*k)/(2*k-1)*(-1.0d0)
         dv = this%matmul(dv)

         ! zero-fixed Dirichlet boundary
         if (present(fix_idx)) then
            dv(fix_idx) = 0
         end if
         a = ((-1.0d0)**k)/gamma(2.0d0*k + 1.0d0)
         b = 0.250d0*((dt - ddt)**(2*k)) + 0.50d0*((dt)**(2*k)) + 0.250d0*((dt + ddt)**(2*k))
         retA_v = retA_v + a*b*dv

         if (norm(abs(a*b*dv)) < tol) return

         if (present(debug)) then
            if (debug) then
               print *, k, norm(abs(a*b*dv))
            end if
         end if

      end do

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

! ###################################################
   function tensor_t_sinc_sqrt_LPF_real64_crs(this, v, tol, itrmax, coeff, debug, fix_idx, cutoff_frequency) &
      result(retA_v)
      class(CRS_), intent(in) :: this
      real(real64), intent(in) :: v(:)
      real(real64), optional, intent(in) :: coeff
      real(real64), allocatable :: dv(:), retA_v(:)
      integer(int32), optional, intent(in) :: fix_idx(:)
      integer(int32) :: k, num_hanning_sample
      integer(int32), intent(in) :: itrmax
      logical, optional, intent(in) :: debug
      real(real64), intent(in) :: tol, cutoff_frequency
      type(Math_) :: math
      real(real64) :: a, b, dt, ddt
      !num_hanning_sample = input(default=16,option=window_size)
      !ddt = 1.0d0/(cutoff_frequency/4.0d0)
      !ddt = 1.0d0/cutoff_frequency/math%pi*acos( (sqrt(2.0d0)-1.0d0)/2.0d0 )
      ddt = 1.0d0/cutoff_frequency/2.0d0/math%pi*acos(sqrt(2.0d0) - 1.0d0)
      if (present(coeff)) then
         a = coeff
      else
         a = 1.0d0
      end if
      ! sinc(sqrt(A) ) = \sum_{n=0}^{\infty}
      !c = 1.0d0
      ! k = 0
      dt = a

      !k=0
      k = 0
      retA_v = v
      dv = v

      if (present(fix_idx)) then
         dv(fix_idx) = 0
      end if

      retA_v = retA_v + dt*dv

      do k = 1, itrmax

         a = ((-1.0d0)**k)*(sqrt(math%pi)*2.0d0**(-1 - 2*k)) &
             /gamma(1.0d0 + k)/gamma(1.5d0 + k)
         b = 0.250d0*((dt - ddt)**(2*k + 1)) + 0.50d0*((dt)**(2*k + 1)) &
             + 0.250d0*((dt + ddt)**(2*k + 1))

         dv = this%matmul(dv)

         ! zero-fixed Dirichlet boundary
         if (present(fix_idx)) then
            dv(fix_idx) = 0
         end if

         retA_v = retA_v + a*b*dv

         if (norm(abs(a*b*dv)) < tol) return

         if (present(debug)) then
            if (debug) then
               print *, k, norm(abs(a*b*dv))
            end if
         end if

      end do

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

! #####################################################
   subroutine fix_complex64_CRS(this, idx, val, RHS, only_row)
      class(CRS_), intent(inout) :: this
      integer(int32), intent(in) :: idx(:)
      complex(real64), optional, intent(inout) :: RHS(:)
      complex(real64), intent(in) :: val(:)
      integer(int64) :: i, j, k, id

      logical, optional, intent(in) :: only_row

      if (present(only_row)) then
         if (only_row) then

            do i = 1, size(idx)
               id = idx(i)
               do j = this%row_ptr(id), this%row_ptr(id + 1) - 1
                  if (this%col_idx(j) == id) then
                     this%val(j) = 1.0d0
                  else
                     this%val(j) = 0.0d0
                  end if
               end do
            end do
            return
         end if
      end if

      do i = 1, size(idx)
         id = idx(i)
         do j = 1, size(this%row_ptr) - 1
            if (j == id) then
               do k = this%row_ptr(j), this%row_ptr(j + 1) - 1
                  if (this%col_idx(k) == id) then
                     this%val(k) = 0.0d0
                  else
                     this%val(k) = 0.0d0
                  end if
               end do
            else
               do k = this%row_ptr(j), this%row_ptr(j + 1) - 1
                  if (this%col_idx(k) == id) then
                     RHS(j) = RHS(j) - this%val(k)*val(i)
                     this%val(k) = 0.0d0
                  end if
               end do
            end if
         end do
      end do

      do i = 1, size(idx)
         id = idx(i)
         do j = 1, size(this%row_ptr) - 1
            if (j == id) then
               do k = this%row_ptr(j), this%row_ptr(j + 1) - 1
                  if (this%col_idx(k) == id) then
                     this%val(k) = 1.0d0
                  end if
               end do
            end if
         end do
      end do

      do i = 1, size(idx)
         RHS(id) = val(i)
      end do
   end subroutine
! #####################################################

   subroutine removeCRS(this, idx)
      class(CRS_), intent(inout) :: this
      integer(int32), optional, intent(in) :: idx(:)
      integer(int64), allocatable :: num_col(:), col_idx(:), only_k(:), &
                                     copy_idx(:), row_ptr(:)
      real(real64), allocatable :: val(:)
      integer(int64) :: row, col_id, count_id, n, k, i, j

      if (.not. present(idx)) then
         if (allocated(this%col_idx)) then
            deallocate (this%col_idx)
         end if
         if (allocated(this%row_ptr)) then
            deallocate (this%row_ptr)
         end if
         if (allocated(this%val)) then
            deallocate (this%val)
         end if

         return
      end if

      copy_idx = idx

      ! remove row and column listed in idx(:)
      do i = 1, size(idx)
         k = maxval(copy_idx)

         num_col = int(zeros(this%size()))

         ! zerofill col_idx if row==k or col==k
         do row = 1, size(this%row_ptr) - 1
            num_col(row) = this%row_ptr(row + 1) - this%row_ptr(row)
            do col_id = this%row_ptr(row), this%row_ptr(row + 1) - 1
               if (row == k .or. this%col_idx(col_id) == k) then
                  this%col_idx(col_id) = 0
                  num_col(row) = num_col(row) - 1
               end if
            end do
         end do
         !call print(num_col)
         ! create new row_ptr
         row_ptr = int(zeros(this%size() - 1 + 1))
         row_ptr(1) = 1
         n = 0
         do row = 1, size(this%row_ptr) - 1
            if (row == k) then
               cycle
            else
               n = n + 1
               row_ptr(n + 1) = row_ptr(n) + num_col(row)
            end if
         end do

         ! remove zero
         col_idx = int(zeros(sum(num_col)))
         val = zeros(sum(num_col))
         n = 0
         do row = 1, k - 1
            do col_id = this%row_ptr(row), this%row_ptr(row + 1) - 1
               if (this%col_idx(col_id) /= 0) then

                  n = n + 1
                  col_idx(n) = this%col_idx(col_id)
                  val(n) = this%val(col_id)
               end if
            end do
         end do

         do row = k + 1, size(this%row_ptr) - 1
            do col_id = this%row_ptr(row), this%row_ptr(row + 1) - 1
               if (this%col_idx(col_id) /= 0) then
                  n = n + 1
                  col_idx(n) = this%col_idx(col_id)
                  val(n) = this%val(col_id)
               end if
            end do
         end do

         this%col_idx = col_idx
         do j = 1, size(this%col_idx)
            if (this%col_idx(j) >= k) then
               this%col_idx(j) = this%col_idx(j) - 1
            end if
         end do
         this%val = val
         this%row_ptr = row_ptr
         if (size(copy_idx) == 1) then
            exit
         else
            copy_idx = removeif(copy_idx, equal_to=k)
            cycle
         end if

      end do

   end subroutine

   function sum_sum(a, b, a_params, b_params, i_plus_j, ir, jr) result(ret)
      interface ! coefficient #1
         function a(k, params) result(ret)
            use iso_fortran_env
            integer(int32), intent(in) :: k
            complex(real64), intent(in) :: params(:)
            complex(real64) :: ret

         end function
      end interface

      interface ! coefficient #2
         function b(k, params) result(ret)
            use iso_fortran_env
            integer(int32), intent(in) :: k
            complex(real64), intent(in) :: params(:)
            complex(real64) :: ret

         end function
      end interface
      integer(int32), intent(in) :: i_plus_j ! i+j
      complex(real64), intent(in) :: a_params(:), b_params(:)
      integer(int32), intent(in) :: ir(1:2) ! range of i
      integer(int32), intent(in) :: jr(1:2) ! range of j
      integer(int32) :: i, j
      complex(real64) :: ret

      ret = 0.0d0

      do i = ir(1), i_plus_j
         j = i_plus_j - i
         if (j < jr(1) .or. jr(2) < j) cycle
         ret = ret + a(i, a_params)*b(j, b_params)
      end do

   end function

! ###################################################
   function tensor_cos_sqrt_cos_sqrt_complex64_crs(this, v, tol, itrmax, coeff_1, coeff_2) result(retA_v)
      class(CRS_), intent(in) :: this
      complex(real64), intent(in) :: v(:)
      complex(real64), allocatable :: dv(:), retA_v(:)
      complex(real64), intent(in) :: coeff_1, coeff_2
      integer(int32) :: k
      integer(int32), intent(in) :: itrmax
      real(real64), intent(in) :: tol
      type(Math_) :: math
      complex(real64) :: coeff

      do k = 0, itrmax
         ! k = k + 1
         coeff = sum_sum( &
                 a=cos_sqrt_Taylor, &
                 b=cos_sqrt_Taylor, &
                 a_params=[coeff_1], &
                 b_params=[coeff_2], &
                 i_plus_j=k, &
                 ir=[0, itrmax], &
                 jr=[0, itrmax] &
                 )
         if (k == 0) then
            dv = v
            retA_v = coeff*dv
         else
            dv = this%matmul(dv)
            retA_v = retA_v + coeff*dv
         end if

         if (norm(abs(dv)) < tol) return
      end do

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

! ###################################################
   function tensor_sinc_sqrt_cos_sqrt_complex64_crs(this, v, tol, itrmax, coeff_1, coeff_2, debug) result(retA_v)
      class(CRS_), intent(in) :: this
      complex(real64), intent(in) :: v(:)
      complex(real64), allocatable :: dv(:), retA_v(:)
      complex(real64), intent(in) :: coeff_1, coeff_2
      logical, optional, intent(in) :: debug
      integer(int32) :: k

      integer(int32), intent(in) :: itrmax
      real(real64), intent(in) :: tol
      type(Math_) :: math
      complex(real64) :: coeff

      do k = 0, itrmax
         ! k = k + 1
         coeff = sum_sum( &
                 a=sinc_sqrt_Taylor, &
                 b=cos_sqrt_Taylor, &
                 a_params=[coeff_1], &
                 b_params=[coeff_2], &
                 i_plus_j=k, &
                 ir=[0, itrmax], &
                 jr=[0, itrmax] &
                 )
         if (k == 0) then
            dv = v
            retA_v = coeff*dv
         else
            dv = this%matmul(dv)
            retA_v = retA_v + coeff*dv
         end if

         if (norm(abs(dv)) < tol) return
         if (present(debug)) then
            if (debug) then
               print *, k, norm(abs(dv))
            end if
         end if

      end do

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

! ###################################################
   function tensor_cos_sqrt_sinc_sqrt_complex64_crs(this, v, tol, itrmax, coeff_1, coeff_2) result(retA_v)
      class(CRS_), intent(in) :: this
      complex(real64), intent(in) :: v(:)
      complex(real64), allocatable :: dv(:), retA_v(:)
      complex(real64), intent(in) :: coeff_1, coeff_2
      integer(int32) :: k
      integer(int32), intent(in) :: itrmax
      real(real64), intent(in) :: tol
      type(Math_) :: math
      complex(real64) :: coeff

      do k = 0, itrmax
         ! k = k + 1
         coeff = sum_sum( &
                 a=cos_sqrt_Taylor, &
                 b=sinc_sqrt_Taylor, &
                 a_params=[coeff_1], &
                 b_params=[coeff_2], &
                 i_plus_j=k, &
                 ir=[0, itrmax], &
                 jr=[0, itrmax] &
                 )
         if (k == 0) then
            dv = v
            retA_v = coeff*dv
         else
            dv = this%matmul(dv)
            retA_v = retA_v + coeff*dv
         end if

         if (norm(abs(dv)) < tol) return
      end do

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

! ###################################################
   function tensor_sinc_sqrt_sinc_sqrt_complex64_crs(this, v, tol, itrmax, coeff_1, coeff_2) result(retA_v)
      class(CRS_), intent(in) :: this
      complex(real64), intent(in) :: v(:)
      complex(real64), allocatable :: dv(:), retA_v(:)
      complex(real64), intent(in) :: coeff_1, coeff_2
      integer(int32) :: k
      integer(int32), intent(in) :: itrmax
      real(real64), intent(in) :: tol
      type(Math_) :: math
      complex(real64) :: coeff

      do k = 0, itrmax
         ! k = k + 1
         coeff = sum_sum( &
                 a=sinc_sqrt_Taylor, &
                 b=sinc_sqrt_Taylor, &
                 a_params=[coeff_1], &
                 b_params=[coeff_2], &
                 i_plus_j=k, &
                 ir=[0, itrmax], &
                 jr=[0, itrmax] &
                 )
         if (k == 0) then
            dv = v
            retA_v = coeff*dv
         else
            dv = this%matmul(dv)
            retA_v = retA_v + coeff*dv
         end if

         if (norm(abs(dv)) < tol) return
      end do

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

   function cos_sqrt_Taylor(k, params) result(ret)
      integer(int32), intent(in) :: k
      complex(real64), intent(in) :: params(:)
      complex(real64) :: ret
      ! Taylor expanision coefficient for cos(a*sqrt(x) ) around x=0
      ret = ((-1.0d0)**k)*(params(1)**(2*k))/gamma(2.0d0*k + 1.0d0)

   end function

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

   function sinc_sqrt_Taylor(k, params) result(ret)
      integer(int32), intent(in) :: k
      complex(real64), intent(in) :: params(:)
      complex(real64) :: ret
      type(Math_) :: math
      ! a = params(1)
      ! Taylor expanision coefficient for sinc(a*sqrt(x) ) around x=0
      ret = dble(2.0d0**(-1 - 2*k))*dble((-1.0d0)**k)*(params(1)**(2*k))*sqrt(math%pi) &
            /Gamma(1.0d0 + k*1.0d0)/Gamma(1.50d0 + k*1.0d0)
      !ret = dble(2.0d0**(-1-2*k))*dble((-1.0d0)**k)!*(params(1)**(2*k) )

   end function

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

   function sinc_complex64(x) result(ret)
      complex(real64), intent(in) :: x
      complex(real64) :: ret
      if (abs(x) == 0.0d0) then
         ret = 1.0d0
      else
         ret = sin(x)/x
      end if
   end function

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

   function sinc_real64(x) result(ret)
      real(real64), intent(in) :: x
      real(real64) :: ret
      if (abs(x) == 0.0d0) then
         ret = 1.0d0
      else
         ret = sin(x)/x
      end if
   end function
! #####################################################

   function sinc_real64_vec(x) result(ret)
      real(real64), intent(in) :: x(:)
      real(real64), allocatable :: ret(:)
      integer(int32) :: i

      allocate (ret(size(x)))

      do i = 1, size(x)
         ret(i) = sinc(x(i))
      end do
   end function

   function to_diag_vector_to_CRS(diag_vec) result(ret)
      real(real64), intent(in) :: diag_vec(:)
      type(CRS_) :: ret

      call ret%eyes(size(diag_vec))
      ret%val = diag_vec

   end function

! ###################################################
   subroutine to_real32_CRS(this)
      class(CRS_), intent(inout) :: this

      this%val_real32 = real(this%val)
      deallocate (this%val)
      this%dtype = real32

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

!! ###################################################
!function divide_by_diagonal_vector_CRS(this,diag_vector) result(ret)
!    type(CRS_),intent(in) :: this
!    type(CRS_) :: ret
!    real(real64),intent(in)  :: diag_vector(:)
!    integer(int32) :: i
!
!    ret = this
!    do i=1,ret%size()
!        ret%val(ret%row_ptr(i):ret%row_ptr(i+1)-1) = &
!            ret%val(ret%row_ptr(i):ret%row_ptr(i+1)-1)/diag_vector(i)
!    enddo
!end function
!! ###################################################
!
!function divide_by_scalar_CRS(this,scalar_value) result(ret)
!    type(CRS_),intent(in) :: this
!    real(real64),intent(in)  :: scalar_value
!    type(CRS_) :: ret
!
!    ret = this
!    ret%val = ret%val/scalar_value
!end function
!! ###################################################

! ###################################################
   function to_COO_from_DenseMatrix(dense_matrix) result(ret)
      type(COO_) :: ret
      real(real64), intent(in) :: dense_matrix(:, :)
      integer(int32) :: DOF, i, j

      DOF = size(dense_matrix, 1)
      call ret%init(DOF)
      do i = 1, size(dense_matrix, 1)
         do j = 1, size(dense_matrix, 2)
            call ret%add(i, j, dense_matrix(i, j))
         end do
      end do

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

! ###################################################
   function to_COO_from_ArrayObject(arrayobject) result(ret)
      type(COO_) :: ret
      type(Array_), intent(in) :: arrayobject
      real(real64), allocatable :: dense_matrix(:, :)
      integer(int32) :: DOF, i, j

      dense_matrix = arrayobject
      DOF = size(dense_matrix, 1)
      call ret%init(DOF)
      do i = 1, size(dense_matrix, 1)
         do j = 1, size(dense_matrix, 1)
            if (dense_matrix(i, j) == 0.0d0) cycle
            call ret%add(i, j, dense_matrix(i, j))
         end do
      end do

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

! ###################################################
   function to_CRS_from_DenseMatrix(dense_matrix) result(ret)
      type(COO_) :: coo
      type(CRS_) :: ret
      real(real64), intent(in) :: dense_matrix(:, :)

      coo = to_COO(dense_matrix)
      ret = coo%to_CRS()

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

! ###################################################
   function to_CRS_from_ArrayObject(arrayobject) result(ret)
      type(COO_) :: coo
      type(CRS_) :: ret
      type(Array_), intent(in) :: arrayobject
      real(real64), allocatable :: dense_matrix(:, :)

      dense_matrix = arrayobject
      coo = to_COO(dense_matrix)
      ret = coo%to_CRS()

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

! ###################################################
   function imag_ICOO(n) result(ret)
      integer(int32), intent(in) :: n
      type(COO_) :: ret
      integer(int32) :: i
      type(Math_) :: math

      call ret%init(n)
      if (n == 1) then
         call ret%add(1, 1, math%i)
         return
      end if

      if (mod(n, 2) == 0) then
         ! even
         do i = 1, n/2
            call ret%add(i, n - i + 1, -1.0d0)
         end do
         do i = n/2 + 1, n
            call ret%add(i, n - i + 1, 1.0d0)
         end do
      else
         ! odd
         ! even
         do i = 1, (n - 1)/2
            call ret%add(i, n, -1.0d0)
         end do
         call ret%add((n - 1)/2 + 1, n, math%i)
         do i = (n - 1)/2 + 2, n
            call ret%add(i, n, 1.0d0)
         end do
      end if

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

! ###################################################
   function speyeCOO(n) result(ret)
      integeR(int32), intent(in) :: n
      type(COO_) :: ret

      call ret%eyes(n)

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

! ###################################################
   subroutine fillSymmetric_CRS(this)
      class(CRS_),intent(in) :: this
      ! fill-in to symmetric parts
      ! under implementation!


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


! ###################################################
   function allocated_CRS(this) result(ret)
      class(CRS_),intent(in) :: this
      logical :: ret

      ret = allocated(this%val)
      ret = ret .and. allocated(this%row_ptr)
      ret = ret .and. allocated(this%col_idx)

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

! ###################################################
! >> BCRS procedures
! ###################################################

! ###################################################
!   subroutine assemble_CRS(this,crs_matrices,matrices_shape)
!      class(CRS_),intent(inout) :: this
!      type(CRS_),intent(in) :: crs_matrices(:)
!      integer(int32),intent(in) :: matrices_shape(1:2)
!
!
!   end subroutine
! ###################################################

   subroutine setBCRS(this,block_position,CRS) 
      class(BCRS_),intent(inout) :: this
      integer(int32) :: block_position(1:2)
      type(CRS_),intent(in) :: CRS
      type(BCRS_) :: BCRS_buffer

      if(.not.allocated(this%CRS))then
         allocate(this%CRS(block_position(1),block_position(2)))
      endif
      if(&
         block_position(1)>size(this%CRS,1) &
            .or. &
         block_position(2)>size(this%CRS,2) )then
         BCRS_buffer = this
         deallocate(this%CRS)
         allocate(this%CRS( max(block_position(1),size(BCRS_buffer%CRS,1)),&
            max(block_position(2),size(BCRS_buffer%CRS,2))))
         this%CRS(1:size(BCRS_buffer%CRS,1),1:size(BCRS_buffer%CRS,2)) = &
            BCRS_buffer%CRS(1:size(BCRS_buffer%CRS,1),1:size(BCRS_buffer%CRS,2))
      endif
      
      this%CRS(block_position(1),block_position(2)) = CRS

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


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

   subroutine addBCRS(this,block_position,CRS) 
      class(BCRS_),intent(inout) :: this
      integer(int32) :: block_position(1:2)
      type(CRS_),intent(in) :: CRS
      type(BCRS_) :: BCRS_buffer
      
      if(.not.allocated(this%CRS))then
         allocate(this%CRS(block_position(1),block_position(2)))
      endif
      if(&
         block_position(1)>size(this%CRS,1) &
            .or. &
         block_position(2)>size(this%CRS,2) )then
         BCRS_buffer = this
         deallocate(this%CRS)
         allocate(this%CRS( max(block_position(1),size(BCRS_buffer%CRS,1)),&
            max(block_position(2),size(BCRS_buffer%CRS,2))))
         this%CRS(1:size(BCRS_buffer%CRS,1),1:size(BCRS_buffer%CRS,2)) = &
            BCRS_buffer%CRS(1:size(BCRS_buffer%CRS,1),1:size(BCRS_buffer%CRS,2))
      endif
      
      this%CRS(block_position(1),block_position(2)) = &
         this%CRS(block_position(1),block_position(2)) + CRS

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


! ###################################################
   subroutine showShapeBCRS(this)
      class(BCRS_),intent(in) :: this
      integer(int32) :: i,j
      character(:),allocatable :: msg      
      integer(int32) :: this_shape(1:2)

      print *, "BCRS :: shape >> "
      do i=1,size(this%CRS,1)
         do j = 1, size(this%CRS,2)
            if(allocated(this%CRS(i,j)))then
               this_shape = this%CRS(i,j)%shape()
               msg = "[" + str(this_shape(1))+","+str(this_shape(2)) + "]"
            else
               msg = "[-0-]"
            endif
            write (*,fmt='(a)', advance='no') msg
         enddo
         write (*,fmt='(a)', advance='yes') ""
      enddo
      print *, "BCRS :: shape << "
   end subroutine
! ###################################################

! ###################################################
   function matmulBCRS(this,vec) result(ret)
      class(BCRS_),intent(in) :: this
      real(real64),intent(in) :: vec(:)
      real(real64),allocatable :: ret(:)
      integer(int32),allocatable :: non_zero_idx(:,:)
      integer(int32) :: i,j,row_r(1:2),col_r(1:2),bcrs_shape(1:2),non_zero_counter,row,col
      
      bcrs_shape = this%shape()
      ret = zeros(bcrs_shape(1))

      ! Pick-up non-zero values
      non_zero_counter = 0
      do i=1,size(this%CRS,1)
         do j=1,size(this%CRS,2)
            if(.not.allocated(this%CRS(i,j))) then
               cycle
            else
               non_zero_counter = non_zero_counter + 1
            endif
         enddo
      enddo
      allocate(non_zero_idx(non_zero_counter,6))
      non_zero_counter = 0
      do i=1,size(this%CRS,1)
         do j=1,size(this%CRS,2)
            if(.not.allocated(this%CRS(i,j))) then
               cycle
            else
               non_zero_counter = non_zero_counter + 1
               non_zero_idx(non_zero_counter,1) = i ! row of block
               non_zero_idx(non_zero_counter,2) = j ! column of block
               non_zero_idx(non_zero_counter,3:4) = this%row_range(i,j) ! global row idx
               non_zero_idx(non_zero_counter,5:6) = this%col_range(i,j) ! global column idx
            endif
         enddo
      enddo
            
      !$OMP parallel do private(row,col,row_r,col_r) reduction(+:ret)
      do i=1,non_zero_counter
         row   = non_zero_idx(i,1)
         col   = non_zero_idx(i,2)
         row_r = non_zero_idx(i,3:4)
         col_r = non_zero_idx(i,5:6)
         ret(row_r(1):row_r(2)) = ret(row_r(1):row_r(2))&
             + this%CRS(row,col)%matmul(vec(col_r(1):col_r(2)))
      enddo
      !$OMP end parallel do

   end function

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


! ###################################################
   function row_range_BCRS(this,box_row,box_col) result(ret)
      class(BCRS_),intent(in) :: this
      integer(int32),intent(in) :: box_row,box_col
      integer(int32) :: ret(1:2),this_shape(1:2)
      integer(int32),allocatable :: row_len(:)
      integer(int32) :: i,j

      ret(1:2) = [1,0]

      do i=1,box_row-1 
         ! for each box_row
         row_len = int(zeros(size(this%CRS,1)))
         do j=1,size(this%CRS,2)
            row_len(j) = 1 .of. this%CRS(i,j)%shape()
         enddo
         
         ret(1) = ret(1) + maxval(row_len)
      enddo

      row_len = int(zeros(size(this%CRS,1)))
      do j=1,size(this%CRS,2)
         row_len(j) = 1 .of. this%CRS(box_row,j)%shape()
      enddo
      
      ret(2) = ret(1) + maxval(row_len) -1
      

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


! ###################################################
   function col_range_BCRS(this,box_row,box_col) result(ret)
      class(BCRS_),intent(in) :: this
      integer(int32),intent(in) :: box_row,box_col
      integer(int32) :: ret(1:2),this_shape(1:2)
      integer(int32),allocatable :: col_len(:)
      integer(int32) :: i,j

      ret(1:2) = [1,0]

      do i=1,box_col-1 
         ! for each box_col
         col_len = int(zeros(size(this%CRS,2)))
         do j=1,size(this%CRS,2)
            col_len(j) = 2 .of. this%CRS(j,i)%shape()
         enddo
         
         ret(1) = ret(1) + maxval(col_len)
      enddo

      col_len = int(zeros(size(this%CRS,2)))
      do j=1,size(this%CRS,2)
         col_len(j) = 2 .of. this%CRS(j,box_col)%shape()
      enddo
      
      ret(2) = ret(1) + maxval(col_len) -1
      

      !this_shape = this%CRS(box_row,box_col)%shape()
      !if (this_shape(2)==0)then
      !   ret(2) = ret(1) 
      !else
      !   ret(2) = ret(1) + this_shape(2) -1
      !endif

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

function shapeBCRS(this) result(ret)
   class(BCRS_),intent(in) :: this
   integer(int32) :: ret(1:2)

   ret(1) = maxval(this%row_range(size(this%CRS,1),size(this%CRS,2)))
   ret(2) = maxval(this%col_range(size(this%CRS,1),size(this%CRS,2)))

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

! ###################################################
function to_dense_BCRS(this) result(ret)
   class(BCRS_),intent(in)  :: this
   real(real64),allocatable :: ret(:,:)
   integer(int32) :: n_row,n_col,i,j,bcrs_row,bcrs_col,row_1,row_2,col_1,col_2

   ret = zeros(1 .of. this%shape(), 2 .of. this%shape())


   
   do bcrs_row=1,size(this%crs,1)
      do bcrs_col=1,size(this%crs,2)
         if(.not.allocated(this%crs(bcrs_row,bcrs_col))) cycle
         row_1 = 1 .of. this%row_range(bcrs_row,bcrs_col)
         row_2 = 1 .of. this%row_range(bcrs_row,bcrs_col) + (1 .of. this%crs(bcrs_row,bcrs_col)%shape()) -1
         col_1 = 1 .of. this%col_range(bcrs_row,bcrs_col)
         col_2 = 1 .of. this%col_range(bcrs_row,bcrs_col) + (2 .of. this%crs(bcrs_row,bcrs_col)%shape()) -1 

         

         ret( &
            row_1: row_2 ,col_1: col_2 ) = this%crs(bcrs_row,bcrs_col)%to_dense()
      enddo
   enddo


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


! ###################################################
function expBCRS(this,vec,max_itr,fix_idx) result(b)
   ! tensor exponential 

   ! exp(x) = \sum_{0}^{\infty} \cfrac{1}{n!}x^{n}
   ! let a_{k} = x^{k}, b_{k} = \sum_{0}^{k} \cfrac{1}{n!}x^{n}
   ! b_{0} = 1
   ! a_{0} = 1
   ! a_{k+1} = \cfrac{1}{k} x*a_{k} 
   ! b_{k+1} = b_{k} + a_{k+1}
   class(BCRS_),intent(in) :: this
   real(real64),intent(in) :: vec(:)
   integer(int32),optional,intent(in) :: max_itr,fix_idx(:)
   integer(int32) :: k,n,itr_max,this_shape(1:2)
   real(real64)   :: tol
   real(real64),allocatable :: a(:),b(:)
   itr_max = input(default=100,option=max_itr)

   a = vec(:)
   b = vec(:)

   do k=1,itr_max
      
      a = 1.0d0/dble(k)*this%matmul(a)
      if(present(fix_idx))then
         if(size(fix_idx)>=1)then
            a(fix_idx(:))=0.0d0
            b(fix_idx(:))=0.0d0
         endif
      endif
      b = b + a
   enddo

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

end module SparseClass
