!>
!> @file ibcmat.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> spclibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> spclibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
SUBROUTINE ibcmat(mat, ny)
!
!   Apply BC on matrix
!
  USE matrix
  IMPLICIT NONE
  TYPE(gbmat), INTENT(inout) :: mat
  INTEGER, INTENT(in) :: ny
  INTEGER :: kl, ku, nrank, i, j
  DOUBLE PRECISION, ALLOCATABLE :: zsum(:), arr(:)
  INTEGER :: i0, ii
  INTEGER :: i0_arr(ny), col(ny)
!===========================================================================
!              1.0 Prologue
!

  kl = mat%kl
  ku = mat%ku
  nrank = mat%rank
  ALLOCATE(zsum(nrank), arr(nrank))
!
  i0 = nrank - ku
  WRITE(*,'(a,i6)') 'Estimated i0', i0
  DO i=1,ny
     arr = 0.0d0
     col(i) = nrank-ny+i
     CALL getcol(mat, nrank-ny+i, arr)
     DO ii=1,nrank
        i0_arr(i)=ii
        IF(arr(ii) .NE. 0.0d0) EXIT
     END DO
  END DO
!!$  WRITE(*,'(a/(10i6))') 'col', col
!!$  WRITE(*,'(a/(10i6))') 'i0_arr', i0_arr
!
!===========================================================================
!              2.0 Unicity at the axis
!
!   The vertical sum on the NY-th row
!
  zsum = 0.0d0
  DO i=1,ny
     arr = 0.0d0
     CALL getrow(mat, i, arr)
     DO j=1,ny+ku
        zsum(j) = zsum(j) + arr(j)
     END DO
  END DO
  CALL putrow(mat, ny, zsum)
!
!   The horizontal sum on the NY-th column
!
  zsum = 0.0d0
  DO j=1,ny
     arr = 0.0d0
     CALL getcol(mat, j, arr)
     DO i=ny,ny+kl
        zsum(i) = zsum(i) + arr(i)
     END DO
  END DO
  CALL putcol(mat, ny, zsum)
!
!   The away operator
!
  DO j = 1,ny-1
     arr = 0.0d0; arr(j) = 1.0d0
     CALL putcol(mat, j, arr)     
  END DO
!
  DO i = 1,ny-1
     arr = 0.0d0; arr(i) = 1.0d0
     CALL putrow(mat, i, arr)     
  END DO
!  
!===========================================================================
!              3.0 Dirichlet on right boundary
!
  DO j = nrank, nrank-ny+1, -1
     arr = 0.0d0; arr(j) = 1.0d0
     CALL putcol(mat, j, arr)     
  END DO
!
  DO i = nrank, nrank-ny+1, -1
     arr = 0.0d0; arr(i) = 1.0d0
     CALL putrow(mat, i, arr)     
  END DO
!===========================================================================
!              9.0  Epilogue
!
  DEALLOCATE(zsum, arr)
!
END SUBROUTINE ibcmat
!+++
SUBROUTINE ibcrhs(rhs, ny)
!
!   Apply BC on RHS
!
  IMPLICIT NONE
  DOUBLE PRECISION, INTENT(inout) :: rhs(:)
  INTEGER, INTENT(in) :: ny
  INTEGER :: nrank
  DOUBLE PRECISION :: zsum
!===========================================================================
!              1.0 Prologue
!
  nrank = SIZE(rhs,1)
!===========================================================================
!              2.0 Unicity at the axis
!
!   The vertical sum
!
  zsum = SUM(rhs(1:ny))
  rhs(ny) = zsum
  rhs(1:ny-1) = 0.0d0
!===========================================================================
!              3.0 Dirichlet on right boundary
!
  rhs(nrank-ny+1:nrank) = 0.0d0
END SUBROUTINE ibcrhs
!+++
SUBROUTINE ibcrhs3(rhs, ny)
!
!   Apply BC on RHS
!
  IMPLICIT NONE
  DOUBLE PRECISION, INTENT(inout) :: rhs(:,:)
  INTEGER, INTENT(in) :: ny
  INTEGER :: nrank, nz, k
  DOUBLE PRECISION :: zsum
!===========================================================================
!              1.0 Prologue
!
  nrank = SIZE(rhs,1)
  nz = SIZE(rhs,2)
!===========================================================================
!              2.0 Unicity at the axis
!
!   The vertical sum
!
  DO k=1,nz
     zsum = SUM(rhs(1:ny,k))
     rhs(ny,k) = zsum
     rhs(1:ny-1,k) = 0.0d0
  END DO
!===========================================================================
!              3.0 Dirichlet on right boundary
!
  DO k=1,nz
     rhs(nrank-ny+1:nrank,k) = 0.0d0
  END DO
END SUBROUTINE ibcrhs3