!>
!> @file sparse_mod.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>
!>
MODULE sparse
!
!    SPARSE: Implement sparse matrix using dynamic linked lists
!            as matrix rows.
!
!    T.M. Tran, CRPP-EPFL
!    October 2010
!
  IMPLICIT NONE
!
  TYPE elt
     INTEGER            :: index=0
     DOUBLE PRECISION   :: val=0.0d0
     TYPE(elt), POINTER :: next => NULL()
  END TYPE elt
!
  TYPE zelt
     INTEGER             :: index=0
     DOUBLE COMPLEX      :: val=(0.0d0, 0.0d0)
     TYPE(zelt), POINTER :: next => NULL()
  END TYPE zelt
!
  TYPE sprow
     INTEGER            :: nnz=0           ! Number of non zeros in a row
     TYPE(elt), POINTER :: row0 => NULL()  ! Points to head of a (sparse) row
  END TYPE sprow
!
  TYPE zsprow
     INTEGER             :: nnz=0          ! Number of non zeros in a row
     TYPE(zelt), POINTER :: row0 => NULL() ! Points to head of a (sparse) row
  END TYPE zsprow
!
  TYPE spmat
     INTEGER :: rank
     TYPE(sprow), POINTER :: row(:) => NULL()
  END TYPE spmat
!
  TYPE zspmat
     INTEGER :: rank
     TYPE(zsprow), POINTER :: row(:) => NULL()
  END TYPE zspmat
!
  INTERFACE init
     MODULE PROCEDURE init_spmat, init_zspmat
 END INTERFACE init
!
  INTERFACE updtmat
     MODULE PROCEDURE updt_sp, updt_zsp, updt_spmat, updt_zspmat
  END INTERFACE updtmat
!
  INTERFACE putele
     MODULE PROCEDURE putele_sp, putele_zsp, putele_spmat, putele_zspmat
  END INTERFACE putele
!
  INTERFACE getele
     MODULE PROCEDURE getele_sp, getele_zsp, getele_spmat, getele_zspmat
  END INTERFACE getele
!
  INTERFACE putrow
     MODULE PROCEDURE putrow_csr,  putrow_full,  putrow_spmat, &
          &           putrow_zcsr, putrow_zfull, putrow_zspmat
  END INTERFACE putrow
!
  INTERFACE getrow
     MODULE PROCEDURE getrow_csr,  getrow_full,  getrow_spmat, &
          &           getrow_zcsr, getrow_zfull, getrow_zspmat
  END INTERFACE getrow
!
  INTERFACE putcol
     MODULE PROCEDURE putcol_spmat, putcol_zspmat
  END INTERFACE putcol
!
  INTERFACE getcol
     MODULE PROCEDURE getcol_spmat, getcol_zspmat
  END INTERFACE getcol
!
  INTERFACE get_count
     MODULE PROCEDURE get_count_sp,  get_count_spmat, &
          &           get_count_zsp, get_count_zspmat
  END INTERFACE get_count
!
  INTERFACE destroy
     MODULE PROCEDURE destroy_spmat,  destroy_row,  destroy_node, &
          &           destroy_zspmat, destroy_zrow, destroy_znode
 END INTERFACE destroy
!
CONTAINS
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE init_spmat(n, mat, istart, iend)
!
!   Initial an empty  sparse matrix
!
    INTEGER, INTENT(in)           :: n
    INTEGER, INTENT(in), OPTIONAL :: istart, iend
    TYPE(spmat)                   :: mat
!
    mat%rank = n
    IF(ASSOCIATED(mat%row)) DEALLOCATE(mat%row)
    IF(PRESENT(istart)) THEN
       ALLOCATE(mat%row(istart:iend))
    ELSE
       ALLOCATE(mat%row(n))
    END IF
!
  END SUBROUTINE init_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE init_zspmat(n, mat, istart, iend)
!
!   Initial an empty  sparse matrix
!
    INTEGER, INTENT(in)           :: n
    INTEGER, INTENT(in), OPTIONAL :: istart, iend
    TYPE(zspmat)                  :: mat
!
    mat%rank = n
    IF(ASSOCIATED(mat%row)) DEALLOCATE(mat%row)
    IF(PRESENT(istart)) THEN
       ALLOCATE(mat%row(istart:iend))
    ELSE
       ALLOCATE(mat%row(n))
    END IF
!
  END SUBROUTINE init_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE updt_sp(arow, j, val)
!
!   Update element j of row arow or insert it in an increasing "index"
!
    TYPE(sprow), TARGET          :: arow
    INTEGER, INTENT(in)          :: j
    DOUBLE PRECISION, INTENT(in) :: val
!
    TYPE(elt), TARGET :: pre_root
    TYPE(elt), POINTER :: t, p
!
    pre_root%next => arow%row0 ! pre_root is linked to the head of the list.
    t => pre_root
    DO WHILE( ASSOCIATED(t%next) )
       p => t%next
       IF( p%index .EQ. j ) THEN
          p%val = p%val+val
          RETURN
       END IF
       IF( p%index .GT. j ) EXIT
       t => t%next
    END DO
    ALLOCATE(p)
    p = elt(j, val, t%next)
    t%next => p
!
    arow%nnz = arow%nnz+1
    arow%row0 => pre_root%next ! In case the head is altered
  END SUBROUTINE updt_sp
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE updt_zsp(arow, j, val)
!
!   Update element j of row arow or insert it in an increasing "index"
!
    TYPE(zsprow), TARGET       :: arow
    INTEGER, INTENT(in)        :: j
    DOUBLE COMPLEX, INTENT(in) :: val
!
    TYPE(zelt), TARGET :: pre_root
    TYPE(zelt), POINTER :: t, p
!
    pre_root%next => arow%row0 ! pre_root is linked to the head of the list.
    t => pre_root
    DO WHILE( ASSOCIATED(t%next) )
       p => t%next
       IF( p%index .EQ. j ) THEN
          p%val = p%val+val
          RETURN
       END IF
       IF( p%index .GT. j ) EXIT
       t => t%next
    END DO
    ALLOCATE(p)
    p = zelt(j, val, t%next)
    t%next => p
!
    arow%nnz = arow%nnz+1
    arow%row0 => pre_root%next ! In case the head is altered
  END SUBROUTINE updt_zsp
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE updt_spmat(mat, i, j, val)
!
!   Update element Aij of sparse matrix
!
    TYPE(spmat)                  :: mat
    INTEGER, INTENT(in)          :: i, j
    DOUBLE PRECISION, INTENT(in) :: val
!
    CALL updt_sp(mat%row(i), j, val)
  END SUBROUTINE updt_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE updt_zspmat(mat, i, j, val)
!
!   Update element Aij of sparse matrix
!
    TYPE(zspmat)               :: mat
    INTEGER, INTENT(in)        :: i, j
    DOUBLE COMPLEX, INTENT(in) :: val
!
    CALL updt_zsp(mat%row(i), j, val)
  END SUBROUTINE updt_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getele_sp(arow, j, val, found)
!
!   Get element j from row arow
!
    TYPE(sprow), TARGET            :: arow
    INTEGER, INTENT(in)            :: j
    DOUBLE PRECISION, INTENT(out)  :: val
    LOGICAL, INTENT(out), OPTIONAL :: found
!
    TYPE(elt), POINTER :: t
    INTEGER :: i
!
    val = 0.0d0
    t => arow%row0  ! Start of a row
    DO WHILE( ASSOCIATED(t) )
       IF(t%index .EQ. j) THEN
          val = t%val
          IF(PRESENT(found)) found = .TRUE.
          RETURN
       END IF
       t => t%next
    END DO
    IF(PRESENT(found)) found = .FALSE.
  END SUBROUTINE getele_sp
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getele_zsp(arow, j, val, found)
!
!   Get element j from row arow
!
    TYPE(zsprow), TARGET           :: arow
    INTEGER, INTENT(in)            :: j
    DOUBLE COMPLEX, INTENT(out)    :: val
    LOGICAL, INTENT(out), OPTIONAL :: found
!
    TYPE(zelt), POINTER :: t
    INTEGER :: i
!
    val = 0.0d0
    t => arow%row0  ! Start of a row
    DO WHILE( ASSOCIATED(t) )
       IF(t%index .EQ. j) THEN
          val = t%val
          IF(PRESENT(found)) found = .TRUE.
          RETURN
       END IF
       t => t%next
    END DO
    IF(PRESENT(found)) found = .FALSE.
  END SUBROUTINE getele_zsp
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getele_spmat(mat, i, j, val)
!
!   Get element (i,j) of sparse matrix
!
    TYPE(spmat)                   :: mat
    INTEGER, INTENT(in)           :: i, j
    DOUBLE PRECISION, INTENT(out) :: val
!
    CALL getele(mat%row(i), j, val)
  END SUBROUTINE getele_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getele_zspmat(mat, i, j, val)
!
!   Get element (i,j) of sparse matrix
!
    TYPE(zspmat)                :: mat
    INTEGER, INTENT(in)         :: i, j
    DOUBLE COMPLEX, INTENT(out) :: val
!
    CALL getele(mat%row(i), j, val)
  END SUBROUTINE getele_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putele_sp(arow, j, val, nlforce_zero)
!
!   Put (overwrite) element j of row arow or insert it in an increasing "index"
!
    TYPE(sprow), TARGET           :: arow
    INTEGER, INTENT(in)           :: j
    DOUBLE PRECISION, INTENT(in)  :: val
    LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
    TYPE(elt), TARGET :: pre_root
    TYPE(elt), POINTER :: t, p
    LOGICAL :: rmnode
!
    pre_root%next => arow%row0 ! pre_root is linked to the head of the list.
    t => pre_root
!
!    Remove node which has zero val or not?
!    But never create new node with zero val
!
    rmnode = .TRUE.
    IF(PRESENT(nlforce_zero)) rmnode = .NOT.nlforce_zero
!
    DO WHILE( ASSOCIATED(t%next) )
       p => t%next
       IF( p%index .EQ. j ) THEN
          IF(ABS(val).LE.EPSILON(0.0d0) .AND. rmnode) THEN ! Remove the node for zero val!
             t%next => p%next
             arow%nnz = arow%nnz-1
             arow%row0 => pre_root%next ! In case the head is altered
             DEALLOCATE(p)
          ELSE
             p%val = val
          END IF
          RETURN
       END IF
       IF( p%index .GT. j ) EXIT
       t => t%next
    END DO
!
!    Never create new node with zero val
!
    IF(ABS(val).GT.EPSILON(0.0d0)) THEN
       ALLOCATE(p)
       p = elt(j, val, t%next)
       t%next => p
       arow%nnz = arow%nnz+1
       arow%row0 => pre_root%next
    END IF
  END SUBROUTINE putele_sp
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putele_zsp(arow, j, val, nlforce_zero)
!
!   Put (overwrite) element j of row arow or insert it in an increasing "index"
!
    TYPE(zsprow), TARGET       :: arow
    INTEGER, INTENT(in)        :: j
    DOUBLE COMPLEX, INTENT(in) :: val
    LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
    TYPE(zelt), TARGET :: pre_root
    TYPE(zelt), POINTER :: t, p
    LOGICAL :: rmnode
!
    pre_root%next => arow%row0 ! pre_root is linked to the head of the list.
    t => pre_root
!
!    Remove node which has zero val or not?
!    But never create new node with zero val
!
    rmnode = .TRUE.
    IF(PRESENT(nlforce_zero)) rmnode = .NOT.nlforce_zero
!
    DO WHILE( ASSOCIATED(t%next) )
       p => t%next
       IF( p%index .EQ. j ) THEN
          IF(ABS(val).LE.EPSILON(0.0d0) .AND. rmnode) THEN ! Remove the node for zero val!
             t%next => p%next
             arow%nnz = arow%nnz-1
             arow%row0 => pre_root%next ! In case the head is altered
             DEALLOCATE(p)
          ELSE
             p%val = val
          END IF
          RETURN
       END IF
       IF( p%index .GT. j ) EXIT
       t => t%next
    END DO
!
!    Never create new node with zero val
!
    IF(ABS(val).GT.EPSILON(0.0d0)) THEN
       ALLOCATE(p)
       p = zelt(j, val, t%next)
       t%next => p
       arow%nnz = arow%nnz+1
       arow%row0 => pre_root%next
    END IF
  END SUBROUTINE putele_zsp
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putele_spmat(mat, i, j, val, nlforce_zero)
!
!  Put element (i,j) of matrix
!
    TYPE(spmat)                   :: mat
    INTEGER, INTENT(in)           :: i, j
    DOUBLE PRECISION, INTENT(in)  :: val
    LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
    CALL putele(mat%row(i), j, val, nlforce_zero)
  END SUBROUTINE putele_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putele_zspmat(mat, i, j, val, nlforce_zero)
!
!  Put element (i,j) of matrix
!
    TYPE(zspmat)                  :: mat
    INTEGER, INTENT(in)           :: i, j
    DOUBLE COMPLEX, INTENT(in)    :: val
    LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
    CALL putele(mat%row(i), j, val, nlforce_zero)
  END SUBROUTINE putele_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  INTEGER FUNCTION get_count_sp(arow)
!
!  Number of elements in arow
!
    TYPE(sprow), INTENT(in) :: arow
    TYPE(elt), POINTER :: t
    INTEGER :: i
!
    t => arow%row0  ! Start of a row
    i = 0
    DO WHILE( ASSOCIATED(t) )
       i=i+1
       t => t%next
    END DO
    get_count_sp = i
  END FUNCTION get_count_sp
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  INTEGER FUNCTION get_count_zsp(arow)
!
!  Number of elements in arow
!
    TYPE(zsprow), INTENT(in) :: arow
    TYPE(zelt), POINTER :: t
    INTEGER :: i
!
    t => arow%row0  ! Start of a row
    i = 0
    DO WHILE( ASSOCIATED(t) )
       i=i+1
       t => t%next
    END DO
    get_count_zsp = i
  END FUNCTION get_count_zsp
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  INTEGER FUNCTION get_count_spmat(mat, nnz)
!
!   Number of non-zeros in sparse matrix
!
    TYPE(spmat)                    :: mat
    INTEGER, INTENT(out), OPTIONAL :: nnz(:)
!
    INTEGER :: i, c(LBOUND(mat%row,1):UBOUND(mat%row,1))
    DO i=LBOUND(mat%row,1),UBOUND(mat%row,1)
       c(i) = get_count_sp(mat%row(i))
    END DO
    IF(PRESENT(nnz)) nnz = c
    get_count_spmat = SUM(c)
  END FUNCTION get_count_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  INTEGER FUNCTION get_count_zspmat(mat, nnz)
!
!   Number of non-zeros in sparse matrix
!
    TYPE(zspmat)                   :: mat
    INTEGER, INTENT(out), OPTIONAL :: nnz(:)
!
    INTEGER :: i, c(LBOUND(mat%row,1):UBOUND(mat%row,1))
    DO i=LBOUND(mat%row,1),UBOUND(mat%row,1)
       c(i) = get_count_zsp(mat%row(i))
    END DO
    IF(PRESENT(nnz)) nnz = c
    get_count_zspmat = SUM(c)
  END FUNCTION get_count_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getrow_csr(arow, arr, col, count)
!
!   Get a row from sparse row arow and put it in a CSR format
!
    TYPE(sprow), INTENT(in)        :: arow
    DOUBLE PRECISION, INTENT(out)  :: arr(:)
    INTEGER, INTENT(out)           :: col(:)
    INTEGER, OPTIONAL, INTENT(out) :: count
!
    TYPE(elt), POINTER :: t
    INTEGER :: i
!
    t => arow%row0  ! Start of a row
    i = 0
    DO WHILE( ASSOCIATED(t) )
       i=i+1
       col(i) = t%index
       arr(i) = t%val
       t => t%next
    END DO
    IF(PRESENT(count)) count = i
  END SUBROUTINE getrow_csr
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getrow_zcsr(arow, arr, col, count)
!
!   Get a row from sparse row arow and put it in a CSR format
!
    TYPE(zsprow), INTENT(in)       :: arow
    DOUBLE COMPLEX, INTENT(out)    :: arr(:)
    INTEGER, INTENT(out)           :: col(:)
    INTEGER, OPTIONAL, INTENT(out) :: count
!
    TYPE(zelt), POINTER :: t
    INTEGER :: i
!
    t => arow%row0  ! Start of a row
    i = 0
    DO WHILE( ASSOCIATED(t) )
       i=i+1
       col(i) = t%index
       arr(i) = t%val
       t => t%next
    END DO
    IF(PRESENT(count)) count = i
  END SUBROUTINE getrow_zcsr
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getrow_full(arow, arr, count)
!
!   Get a row from sparse row arow and put it in an full row
!
    TYPE(sprow), INTENT(in)        :: arow
    DOUBLE PRECISION, INTENT(out)  :: arr(:)
    INTEGER, OPTIONAL, INTENT(out) :: count
!
    TYPE(elt), POINTER :: t
    INTEGER :: n, i, j
!
    n = SIZE(arr)
    arr = 0.0d0
    t => arow%row0  ! Start of a row
    i = 0
    DO WHILE( ASSOCIATED(t) )
       i=i+1
       j = t%index
       IF(j.LE.n) THEN
          arr(j) = t%val
          t => t%next
       ELSE
          WRITE(*,'(a)') 'GETROW_FULL: size of input ARR too small!'
          STOP
       END IF
    END DO
    IF(PRESENT(count)) count = i
  END SUBROUTINE getrow_full
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getrow_zfull(arow, arr, count)
!
!   Get a row from sparse row arow and put it in an full row
!
    TYPE(zsprow), INTENT(in)       :: arow
    DOUBLE COMPLEX, INTENT(out)    :: arr(:)
    INTEGER, OPTIONAL, INTENT(out) :: count
!
    TYPE(zelt), POINTER :: t
    INTEGER :: n, i, j
!
    n = SIZE(arr)
    arr = 0.0d0
    t => arow%row0  ! Start of a row
    i = 0
    DO WHILE( ASSOCIATED(t) )
       i=i+1
       j = t%index
       IF(j.LE.n) THEN
          arr(j) = t%val
          t => t%next
       ELSE
          WRITE(*,'(a)') 'GETROW_FULL: size of input ARR too small!'
          STOP
       END IF
    END DO
    IF(PRESENT(count)) count = i
  END SUBROUTINE getrow_zfull
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getrow_spmat(mat, i, arr, col)
!
!   Get a row from sparse matrix
!
    TYPE(spmat), INTENT(in)        :: mat
    INTEGER, INTENT(in)            :: i
    DOUBLE PRECISION, INTENT(out)  :: arr(:)
    INTEGER, INTENT(out), OPTIONAL :: col(:)
!
    IF(PRESENT(col)) THEN ! The output row is defined by (col, arr)
       CALL getrow_csr(mat%row(i), arr, col)
    ELSE
       CALL getrow_full(mat%row(i), arr)
    END IF
  END SUBROUTINE getrow_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getrow_zspmat(mat, i, arr, col)
!
!   Get a row from sparse matrix
!
    TYPE(zspmat), INTENT(in)       :: mat
    INTEGER, INTENT(in)            :: i
    DOUBLE COMPLEX, INTENT(out)    :: arr(:)
    INTEGER, INTENT(out), OPTIONAL :: col(:)
!
    IF(PRESENT(col)) THEN ! The output row is defined by (col, arr)
       CALL getrow_zcsr(mat%row(i), arr, col)
    ELSE
       CALL getrow_zfull(mat%row(i), arr)
    END IF
  END SUBROUTINE getrow_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putrow_csr(arow, arr, col, nlforce_zero)
!
!   Put a row from sparse row arow and put it in a CSR format
!
    TYPE(sprow), INTENT(inout)    :: arow
    DOUBLE PRECISION, INTENT(in)  :: arr(:)
    INTEGER, INTENT(in)           :: col(:)
    LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
    INTEGER :: n, i
!
    n=SIZE(arr)
    DO i=1,n
       CALL putele(arow, col(i), arr(i), nlforce_zero)
    END DO
  END SUBROUTINE putrow_csr
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putrow_zcsr(arow, arr, col, nlforce_zero)
!
!   Put a row from sparse row arow and put it in a CSR format
!
    TYPE(zsprow), INTENT(inout)   :: arow
    DOUBLE COMPLEX, INTENT(in)    :: arr(:)
    INTEGER, INTENT(in)           :: col(:)
    LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
    INTEGER :: n, i
!
    n=SIZE(arr)
    DO i=1,n
       CALL putele(arow, col(i), arr(i), nlforce_zero)
    END DO
  END SUBROUTINE putrow_zcsr
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putrow_full(arow, arr, nlforce_zero)
!
!   Put a row from sparse row arow and put it in a full row
!
    TYPE(sprow), INTENT(inout)    :: arow
    DOUBLE PRECISION, INTENT(in)  :: arr(:)
    LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
    INTEGER :: n, i
!
    n=SIZE(arr)
    DO i=1,n
       CALL putele(arow, i, arr(i), nlforce_zero)
    END DO
  END SUBROUTINE putrow_full
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putrow_zfull(arow, arr, nlforce_zero)
!
!   Put a row from sparse row arow and put it in a full row
!
    TYPE(zsprow), INTENT(inout)    :: arow
    DOUBLE COMPLEX, INTENT(in)  :: arr(:)
    LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
    INTEGER :: n, i
!
    n=SIZE(arr)
    DO i=1,n
       CALL putele(arow, i, arr(i), nlforce_zero)
    END DO
  END SUBROUTINE putrow_zfull
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putrow_spmat(mat, i, arr, col, nlforce_zero)
!
!   Put a row to matrix
!
    TYPE(spmat)                   :: mat
    INTEGER, intent(in)           :: i
    DOUBLE PRECISION, INTENT(in)  :: arr(:)
    INTEGER, INTENT(in), OPTIONAL :: col(:)
    LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
    IF(PRESENT(col)) THEN ! The input row is defined by (col, arr)
       CALL putrow_csr(mat%row(i), arr, col, nlforce_zero)
    ELSE
       CALL putrow_full(mat%row(i), arr, nlforce_zero)
    END IF
  END SUBROUTINE putrow_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putrow_zspmat(mat, i, arr, col, nlforce_zero)
!
!   Put a row to matrix
!
    TYPE(zspmat)                  :: mat
    INTEGER, intent(in)           :: i
    DOUBLE COMPLEX, INTENT(in)    :: arr(:)
    INTEGER, INTENT(in), OPTIONAL :: col(:)
    LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
    IF(PRESENT(col)) THEN ! The input row is defined by (col, arr)
       CALL putrow_zcsr(mat%row(i), arr, col, nlforce_zero)
    ELSE
       CALL putrow_zfull(mat%row(i), arr, nlforce_zero)
    END IF
  END SUBROUTINE putrow_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putcol_spmat(mat, j, arr, nlforce_zero)
!
!   Put a column to mtarix
!
    TYPE(spmat)                   :: mat
    INTEGER, INTENT(in)           :: j
    DOUBLE PRECISION, INTENT(in)  :: arr(:)
    LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
    INTEGER :: i
    DO i=1,mat%rank
       CALL putele(mat, i, j, arr(i), nlforce_zero)
    END DO
  END SUBROUTINE putcol_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putcol_zspmat(mat, j, arr, nlforce_zero)
!
!   Put a column to mtarix
!
    TYPE(zspmat)                  :: mat
    INTEGER, INTENT(in)           :: j
    DOUBLE COMPLEX, INTENT(in)    :: arr(:)
    LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
    INTEGER :: i
    DO i=1,mat%rank
       CALL putele(mat, i, j, arr(i), nlforce_zero)
    END DO
  END SUBROUTINE putcol_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getcol_spmat(mat, j, arr)
!
!   Get column j of matrix
!
    TYPE(spmat)                   :: mat
    INTEGER, INTENT(in)           :: j
    DOUBLE PRECISION, INTENT(out) :: arr(:)
    INTEGER :: i
    DO i=1,mat%rank
       CALL getele(mat, i, j, arr(i))
    END DO
  END SUBROUTINE getcol_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getcol_zspmat(mat, j, arr)
!
!   Get column j of matrix
!
    TYPE(zspmat)                :: mat
    INTEGER, INTENT(in)         :: j
    DOUBLE COMPLEX, INTENT(out) :: arr(:)
    INTEGER :: i
    DO i=1,mat%rank
       CALL getele(mat, i, j, arr(i))
    END DO
  END SUBROUTINE getcol_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE destroy_row(arow)
!
!  Deallocate the sparse row
!
    TYPE(sprow), INTENT(inout) :: arow
!
    IF(ASSOCIATED(arow%row0)) CALL destroy_node(arow%row0)
    arow%nnz = get_count(arow)
  END SUBROUTINE destroy_row
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE destroy_zrow(arow)
!
!  Deallocate the sparse row
!
    TYPE(zsprow), INTENT(inout) :: arow
!
    IF(ASSOCIATED(arow%row0)) CALL destroy_znode(arow%row0)
    arow%nnz = get_count(arow)
  END SUBROUTINE destroy_zrow
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  RECURSIVE SUBROUTINE destroy_node(p)
!
!  Deallocate recursively the linked list
!
    TYPE(elt), POINTER :: p
!
    IF(ASSOCIATED(p%next)) CALL destroy_node(p%next)
    DEALLOCATE(p)
  END SUBROUTINE destroy_node
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  RECURSIVE SUBROUTINE destroy_znode(p)
!
!  Deallocate recursively the linked list
!
    TYPE(zelt), POINTER :: p
!
    IF(ASSOCIATED(p%next)) CALL destroy_znode(p%next)
    DEALLOCATE(p)
  END SUBROUTINE destroy_znode
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE destroy_spmat(mat)
!
!  Deallocate the sparse matrix
!
    TYPE(spmat) :: mat
    INTEGER :: n, i
!
    n = mat%rank
    DO i=LBOUND(mat%row,1),UBOUND(mat%row,1)
       CALL destroy(mat%row(i))
    END DO
    DEALLOCATE(mat%row)
  END SUBROUTINE destroy_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE destroy_zspmat(mat)
!
!  Deallocate the sparse matrix
!
    TYPE(zspmat) :: mat
    INTEGER :: n, i
!
    n = mat%rank
    DO i=LBOUND(mat%row,1),UBOUND(mat%row,1)
       CALL destroy(mat%row(i))
    END DO
    DEALLOCATE(mat%row)
  END SUBROUTINE destroy_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  INTEGER FUNCTION isearch(karr, k)
!
!   Sequential search an ordered table of integers
!
    INTEGER, INTENT(in) :: karr(0:)
    INTEGER, INTENT(in) :: k
    INTEGER :: n
!
    n=SIZE(karr)
    isearch = -1   ! Failure
    IF( k.GT.karr(n-1)) RETURN
!
    isearch=0
    DO 
       IF( k.LE.karr(isearch)) EXIT
       isearch = isearch+1
    END DO
    IF( k.NE.karr(isearch)) isearch = -1  ! Failure
  END FUNCTION isearch
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  INTEGER FUNCTION isearch_bin(karr, k)
!
!   Binary search an ordered table of integers
!
    INTEGER, INTENT(in) :: karr(0:)
    INTEGER, INTENT(in) :: k
    INTEGER :: n
    INTEGER :: l, u
!
    n=SIZE(karr)
    isearch_bin = -1   ! Failure
    IF( k.LT.karr(0) .OR. k.GT.karr(n-1)) RETURN
!
    l=0; u=n-1
    DO WHILE(l.LE.u)
       isearch_bin = (l+u)/2
       IF(k.EQ.karr(isearch_bin)) THEN 
          RETURN
       ELSE IF(k.LT.karr(isearch_bin)) THEN
          u = isearch_bin-1
       ELSE
          l = isearch_bin+1
       END IF
    END DO
    isearch_bin = -1   ! Failure
  END FUNCTION isearch_bin
!
END MODULE sparse