diff --git a/src/csr_mod.f90 b/src/csr_mod.f90
index 23afb72..ab2f57d 100644
--- a/src/csr_mod.f90
+++ b/src/csr_mod.f90
@@ -1,1255 +1,1255 @@
 !>
 !> @file csr_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/>.
 !>
 !> @author
 !> (in alphabetical order)
 !> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
 !>
 MODULE csr
 !
 !   CSR: Implement CSR (Compressed Sparse Row) matrice
 !
 !    T.M. Tran, CRPP-EPFL
 !    October 2012
 !
   USE sparse
   USE mumps_bsplines
   IMPLICIT NONE
 !
   TYPE, EXTENDS(spmat) :: csr_mat
      INTEGER                   :: mrows, ncols
      INTEGER                   :: nnz = 0            ! Number of non-zeros
      INTEGER                   :: nterms             ! Number of terms in weak form
      LOGICAL                   :: nlforce_zero       ! Keep exixting nodes with zero value if .true.
      INTEGER, POINTER          :: irow(:) => NULL()  ! points to start of row
      INTEGER, POINTER          :: idiag(:) => NULL() ! points to diagonal element
      INTEGER, POINTER          :: cols(:) => NULL()  ! Column indices
      DOUBLE PRECISION, POINTER :: val(:) => NULL()   ! Elelement values
      TYPE(mumps_mat), ALLOCATABLE :: mumps
   END TYPE csr_mat
 !
   TYPE, EXTENDS(zspmat) :: zcsr_mat
      INTEGER                   :: mrows, ncols
      INTEGER                   :: nnz = 0            ! Number of non-zeros
      INTEGER                   :: nterms             ! Number of terms in weak form
      LOGICAL                   :: nlforce_zero       ! Keep exixting nodes with zero value if .true.
      INTEGER, POINTER          :: irow(:) => NULL()  ! points to start of row
      INTEGER, POINTER          :: idiag(:) => NULL() ! points to diagonal element
      INTEGER, POINTER          :: cols(:) => NULL()  ! Column indices
      DOUBLE COMPLEX, POINTER :: val(:) => NULL()   ! Elelement values
 !
      TYPE(zmumps_mat), ALLOCATABLE :: mumps
   END TYPE zcsr_mat
 !
   INTERFACE init
      MODULE PROCEDURE init_csr_mat, init_zcsr_mat
   END INTERFACE init
   INTERFACE clear_mat
      MODULE PROCEDURE clear_csr_mat, clear_zcsr_mat
   END INTERFACE clear_mat
   INTERFACE updtmat
      MODULE PROCEDURE updt_csr_mat, updt_zcsr_mat
   END INTERFACE updtmat
   INTERFACE putele
      MODULE PROCEDURE putele_csr_mat, putele_zcsr_mat
   END INTERFACE putele
   INTERFACE getele
      MODULE PROCEDURE getele_csr_mat, getele_zcsr_mat
   END INTERFACE getele
   INTERFACE putrow
      MODULE PROCEDURE putrow_csr_mat, putrow_zcsr_mat
   END INTERFACE putrow
   INTERFACE getrow
      MODULE PROCEDURE getrow_csr_mat, getrow_zcsr_mat
   END INTERFACE getrow
   INTERFACE getdiag
      MODULE PROCEDURE getdiag_csr_mat, getdiag_zcsr_mat
   END INTERFACE getdiag
   INTERFACE putcol
      MODULE PROCEDURE putcol_csr_mat, putcol_zcsr_mat
   END INTERFACE putcol
   INTERFACE getcol
      MODULE PROCEDURE getcol_csr_mat, getcol_zcsr_mat
   END INTERFACE getcol
   INTERFACE to_mat
      MODULE PROCEDURE to_csr_mat, to_zcsr_mat
   END INTERFACE to_mat
   INTERFACE vmx
      MODULE PROCEDURE vmx_csr_mat,  vmx_csr_matn, vmx_zcsr_mat,  vmx_zcsr_matn
   END INTERFACE vmx
   INTERFACE destroy
      MODULE PROCEDURE destroy_csr_mat, destroy_zcsr_mat
   END INTERFACE destroy
   INTERFACE putmat
      MODULE PROCEDURE put_csr_mat, put_zcsr_mat
   END INTERFACE putmat
 !>>>>>
 !>>>>> CONMAT
 !>>>>
   INTERFACE conmat
      MODULE PROCEDURE conmat_1d_csr, conmat_2d_csr, conmat_1d_zcsr, conmat_2d_zcsr
   END INTERFACE conmat
 !>>>>
 !>>>> MULTIGRID_MOD
 !>>>>
   INTERFACE femat
      MODULE PROCEDURE femat_csr
   END INTERFACE femat
   INTERFACE matnorm
      MODULE PROCEDURE matnorm_csr
   END INTERFACE matnorm
   INTERFACE kron
      MODULE PROCEDURE kron_csr
   END INTERFACE kron
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 CONTAINS
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE init_csr_mat(n, nterms, mat, nlforce_zero, ncols)
 !
 !   Initialize an empty CSR matrix
 !
     INTEGER, INTENT(in)           :: n, nterms
     TYPE(csr_mat)                 :: mat
     LOGICAL, OPTIONAL, INTENT(in) :: nlforce_zero
     INTEGER, OPTIONAL, intent(in) :: ncols
 !
     CALL init(n, mat%spmat)
     mat%mrows = n
     mat%ncols = n
     IF(PRESENT(ncols)) mat%ncols = ncols
     mat%nterms = nterms
     mat%nlforce_zero = .TRUE. ! Do not remove existing nodes with zero val
     IF(PRESENT(nlforce_zero)) mat%nlforce_zero = nlforce_zero
 !
   END SUBROUTINE init_csr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE init_zcsr_mat(n, nterms, mat, nlforce_zero, ncols)
 !
 !   Initialize an empty CSR matrix
 !
     INTEGER, INTENT(in)           :: n, nterms
     TYPE(zcsr_mat)                :: mat
     LOGICAL, OPTIONAL, INTENT(in) :: nlforce_zero
     INTEGER, OPTIONAL, intent(in) :: ncols
 !
     CALL init(n, mat%zspmat)
     mat%mrows = n
     mat%ncols = n
     IF(PRESENT(ncols)) mat%ncols = ncols
     mat%nterms = nterms
     mat%nlforce_zero = .TRUE. ! Do not remove existing nodes with zero val
     IF(PRESENT(nlforce_zero)) mat%nlforce_zero = nlforce_zero
 !
   END SUBROUTINE init_zcsr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE clear_zcsr_mat(mat)
 !
 !   Clear matrix, keeping its sparse structure unchanged
 !
     TYPE(zcsr_mat) :: mat
 !
     mat%val = (0.0d0,0.0d0)
   END SUBROUTINE clear_zcsr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE clear_csr_mat(mat)
 !
 !   Clear matrix, keeping its sparse structure unchanged
 !
     TYPE(csr_mat) :: mat
 !
     mat%val = 0.0d0
   END SUBROUTINE clear_csr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE updt_csr_mat(mat, i, j, val)
 !
 !   Update element Aij of csr  matrix
 !
     TYPE(csr_mat)                :: mat
     INTEGER, INTENT(in)          :: i, j
     DOUBLE PRECISION, INTENT(in) :: val
     INTEGER :: k, s, e
 !
     IF(mat%nnz.EQ.0) THEN   ! Still using linked lists
        CALL updtmat(mat%spmat, i, j, val)
     ELSE
        s = mat%irow(i)
        e = mat%irow(i+1)-1
        k = isearch(mat%cols(s:e), j)
        IF( k.GE.0 ) THEN
           mat%val(s+k) =  mat%val(s+k)+val
        ELSE
           WRITE(*,'(a,2i6)') 'UPDT: i, j out of range ', i, j
           STOP '*** Abnormal EXIT in MODULE csr_mod ***'
        END IF
     END IF
   END SUBROUTINE updt_csr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE updt_zcsr_mat(mat, i, j, val)
 !
 !   Update element Aij of csr  matrix
 !
     TYPE(zcsr_mat)               :: mat
     INTEGER, INTENT(in)          :: i, j
     DOUBLE COMPLEX, INTENT(in)   :: val
     INTEGER :: k, s, e
 !
     IF(mat%nnz.EQ.0) THEN   ! Still using linked lists
        CALL updtmat(mat%zspmat, i, j, val)
     ELSE
        s = mat%irow(i)
        e = mat%irow(i+1)-1
        k = isearch(mat%cols(s:e), j)
        IF( k.GE.0 ) THEN
           mat%val(s+k) =  mat%val(s+k)+val
        ELSE
           WRITE(*,'(a,2i6)') 'UPDT: i, j out of range ', i, j
           STOP '*** Abnormal EXIT in MODULE csr_mod ***'
        END IF
     END IF
   END SUBROUTINE updt_zcsr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE putele_csr_mat(mat, i, j, val)
 !
 !  Put element (i,j) of matrix
 !
     TYPE(csr_mat)                 :: mat
     INTEGER, INTENT(in)          :: i, j
     DOUBLE PRECISION, INTENT(in) :: val
     INTEGER :: iput, jput
     INTEGER :: k, s, e
 !
     iput = i
     jput = j
 !
     IF(mat%nnz.EQ.0) THEN   ! Still using linked lists
        CALL putele(mat%spmat, iput, jput, val, &
             &      nlforce_zero=mat%nlforce_zero)
     ELSE
        s = mat%irow(iput)
        e = mat%irow(iput+1)-1
        k = isearch(mat%cols(s:e), jput)
        IF( k.GE.0 ) THEN
           mat%val(s+k) =  val
        ELSE
           IF(ABS(val) .GT. EPSILON(0.0d0)) THEN   ! Ok for zero val
              WRITE(*,'(a,2i6)') 'PUTELE: i, j out of range ', i, j
              PRINT*, 'val', val
              PRINT*, 'matrix m, n', mat%mrows, mat%ncols
              STOP '*** Abnormal EXIT in MODULE csr_mod ***'
           END IF
        END IF
    END IF
   END SUBROUTINE putele_csr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE putele_zcsr_mat(mat, i, j, val)
 !
 !  Put element (i,j) of matrix
 !
     TYPE(zcsr_mat)                :: mat
     INTEGER, INTENT(in)          :: i, j
     DOUBLE COMPLEX, INTENT(in)  :: val
     INTEGER :: iput, jput
     INTEGER :: k, s, e
 !
     iput = i
     jput = j
 !
     IF(mat%nnz.EQ.0) THEN   ! Still using linked lists
        CALL putele(mat%zspmat, iput, jput, val, &
             &      nlforce_zero=mat%nlforce_zero)
     ELSE
        s = mat%irow(iput)
        e = mat%irow(iput+1)-1
        k = isearch(mat%cols(s:e), jput)
        IF( k.GE.0 ) THEN
           mat%val(s+k) =  val
        ELSE
           IF(ABS(val) .GT. EPSILON(0.0d0)) THEN   ! Ok for zero val
              WRITE(*,'(a,2i6)') 'PUTELE: i, j out of range ', i, j
              PRINT*, 'val', val
              PRINT*, 'matrix m, n', mat%mrows, mat%ncols
              STOP '*** Abnormal EXIT in MODULE csr_mod ***'
           END IF
        END IF
    END IF
  END SUBROUTINE putele_zcsr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE getele_csr_mat(mat, i, j, val)
 !
 !   Get element (i,j) of sparse matrix
 !
     TYPE(csr_mat)                 :: mat
     INTEGER, INTENT(in)           :: i, j
     DOUBLE PRECISION, INTENT(out) :: val
     INTEGER :: iget, jget
     INTEGER :: k, s, e
 !
     iget = i
     jget = j
     IF(mat%nnz.EQ.0) THEN   ! Still using linked lists
        CALL getele(mat%spmat, iget, jget, val)
     ELSE
        s = mat%irow(iget)
        e = mat%irow(iget+1)-1
        k = isearch(mat%cols(s:e), jget)
        IF( k.GE.0 ) THEN
-          val =mat%val(s+k) 
+          val =mat%val(s+k)
        ELSE
           val = 0.0d0   ! Assume zero val if not found
        END IF
     END IF
   END SUBROUTINE getele_csr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE getele_zcsr_mat(mat, i, j, val)
 !
 !   Get element (i,j) of sparse matrix
 !
     TYPE(zcsr_mat)                :: mat
     INTEGER, INTENT(in)           :: i, j
     DOUBLE COMPLEX, INTENT(out)   :: val
     INTEGER :: iget, jget
     INTEGER :: k, s, e
 !
     iget = i
     jget = j
     IF(mat%nnz.EQ.0) THEN   ! Still using linked lists
        CALL getele(mat%zspmat, iget, jget, val)
     ELSE
        s = mat%irow(iget)
        e = mat%irow(iget+1)-1
        k = isearch(mat%cols(s:e), jget)
        IF( k.GE.0 ) THEN
-          val =mat%val(s+k) 
+          val =mat%val(s+k)
        ELSE
           val = 0.0d0   ! Assume zero val if not found
        END IF
     END IF
   END SUBROUTINE getele_zcsr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE putrow_csr_mat(amat, i, arr)
 !
 ! Put a row into sparse matrix
 !
     TYPE(csr_mat), INTENT(inout) :: amat
     INTEGER, INTENT(in)          :: i
     DOUBLE PRECISION, INTENT(in) :: arr(:)
     INTEGER :: s, e, j
 !
     IF(amat%nnz.EQ.0) THEN   ! Still using linked lists
        DO j=1,amat%ncols
           CALL putele(amat, i, j, arr(j))
        END DO
     ELSE
        s = amat%irow(i)
        e = amat%irow(i+1)-1
        DO j=s,e
           amat%val(j) = arr(amat%cols(j))
        END DO
     END IF
   END SUBROUTINE putrow_csr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE putrow_zcsr_mat(amat, i, arr)
 !
 ! Put a row into sparse matrix
 !
     TYPE(zcsr_mat), INTENT(inout) :: amat
     INTEGER, INTENT(in)           :: i
     DOUBLE COMPLEX, INTENT(in)    :: arr(:)
     INTEGER :: s, e, j
 !
     IF(amat%nnz.EQ.0) THEN   ! Still using linked lists
        DO j=1,amat%ncols
           CALL putele(amat, i, j, arr(j))
        END DO
     ELSE
        s = amat%irow(i)
        e = amat%irow(i+1)-1
        DO j=s,e
           amat%val(j) = arr(amat%cols(j))
        END DO
     END IF
   END SUBROUTINE putrow_zcsr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE getrow_csr_mat(amat, i, arr)
 !
 !   Get a row from sparse matrix
 !
     TYPE(csr_mat), INTENT(in)     :: amat
     INTEGER, INTENT(in)           :: i
     DOUBLE PRECISION, INTENT(out) :: arr(:)
     INTEGER :: s, e, j
 !
     arr = 0.0d0
     IF(amat%nnz.EQ.0) THEN   ! Still using linked lists
        DO j=1,amat%ncols
           CALL getele(amat, i, j, arr(j))
        END DO
     ELSE
        s = amat%irow(i)
        e = amat%irow(i+1)-1
        DO j=s,e
           arr(amat%cols(j)) = amat%val(j)
        END DO
     END IF
   END SUBROUTINE getrow_csr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE getrow_zcsr_mat(amat, i, arr)
 !
 !   Get a row from sparse matrix
 !
     TYPE(zcsr_mat), INTENT(in)    :: amat
     INTEGER, INTENT(in)           :: i
     DOUBLE COMPLEX, INTENT(out)   :: arr(:)
     INTEGER :: s, e, j
 !
     arr = 0.0d0
     IF(amat%nnz.EQ.0) THEN   ! Still using linked lists
        DO j=1,amat%ncols
           CALL getele(amat, i, j, arr(j))
        END DO
     ELSE
        s = amat%irow(i)
        e = amat%irow(i+1)-1
        DO j=s,e
           arr(amat%cols(j)) = amat%val(j)
        END DO
     END IF
   END SUBROUTINE getrow_zcsr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE getdiag_csr_mat(amat, arr)
 !
 !   Get the diagonal from matrix
 !
     TYPE(csr_mat), INTENT(in)     :: amat
     DOUBLE PRECISION, INTENT(out) :: arr(:)
 !
 !  WARNING: assume that CSR matrix has been converted from linked lists
 !
     arr(:) = amat%val(amat%idiag(:))
   END SUBROUTINE getdiag_csr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE getdiag_zcsr_mat(amat, arr)
 !
 !   Get the diagonal from matrix
 !
     TYPE(zcsr_mat), INTENT(in)    :: amat
     DOUBLE COMPLEX, INTENT(out)   :: arr(:)
 !
 !  WARNING: assume that CSR matrix has been converted from linked lists
 !
     arr(:) = amat%val(amat%idiag(:))
   END SUBROUTINE getdiag_zcsr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE putcol_csr_mat(amat, j, arr)
 !
 ! Put a column into sparse matrix
 !
     TYPE(csr_mat), INTENT(inout)     :: amat
     INTEGER, INTENT(in)              :: j
     DOUBLE PRECISION, INTENT(in)     :: arr(:)
     INTEGER :: i
 !
     DO i=1,amat%mrows
        CALL putele(amat, i, j, arr(i))
     END DO
   END SUBROUTINE putcol_csr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE putcol_zcsr_mat(amat, j, arr)
 !
 ! Put a column into sparse matrix
 !
     TYPE(zcsr_mat), INTENT(inout)    :: amat
     INTEGER, INTENT(in)              :: j
     DOUBLE COMPLEX, INTENT(in)       :: arr(:)
     INTEGER :: i
 !
     DO i=1,amat%mrows
        CALL putele(amat, i, j, arr(i))
     END DO
   END SUBROUTINE putcol_zcsr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE getcol_csr_mat(amat, j, arr)
 !
 !   Get a column from sparse matrix
 !
     TYPE(csr_mat), INTENT(in)     :: amat
     INTEGER, INTENT(in)           :: j
     DOUBLE PRECISION, INTENT(out) :: arr(:)
     INTEGER :: i
 !
     DO i=1,amat%mrows
        CALL getele(amat, i, j, arr(i))
     END DO
   END SUBROUTINE getcol_csr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE getcol_zcsr_mat(amat, j, arr)
 !
 !   Get a column from sparse matrix
 !
     TYPE(zcsr_mat), INTENT(in)    :: amat
     INTEGER, INTENT(in)           :: j
     DOUBLE COMPLEX, INTENT(out)   :: arr(:)
     INTEGER :: i
 !
     DO i=1,amat%mrows
        CALL getele(amat, i, j, arr(i))
     END DO
   END SUBROUTINE getcol_zcsr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE to_csr_mat(mat, nlkeep)
 !
 !   Convert linked list spmat to csr matrice structure
 !
     TYPE(csr_mat)             :: mat
     LOGICAL, INTENT(in), OPTIONAL :: nlkeep
     INTEGER :: nnz_arr(mat%rank)
     INTEGER :: i, nnz, rank, s, e
     LOGICAL :: nlclean
 !
     nlclean = .TRUE.
     IF(PRESENT(nlkeep)) THEN
        nlclean = .NOT. nlkeep
     END IF
 !
 !    Allocate the csr matrix structure
 !
     nnz = get_count(mat%spmat, nnz_arr)
     rank = mat%rank
     mat%nnz = nnz
     IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
     IF(ASSOCIATED(mat%idiag)) DEALLOCATE(mat%idiag)
     IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
     IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
     ALLOCATE(mat%irow(rank+1))
     ALLOCATE(mat%idiag(rank))
     ALLOCATE(mat%cols(nnz))
     ALLOCATE(mat%val(nnz))
 !
 !    Fill csr structure and optionally deallocate the sparse rows
 !
     mat%irow = 1
     DO i=1,rank
        mat%irow(i+1) = mat%irow(i) + nnz_arr(i)
        s = mat%irow(i)
        e = mat%irow(i+1)-1
        CALL getrow(mat%spmat%row(i), mat%val(s:e), mat%cols(s:e))
        mat%idiag(i) = isearch(mat%cols(s:e), i) + s
        IF(nlclean) CALL destroy(mat%spmat%row(i))
     END DO
 !!$!
 !!$!    MUMPS mat for direct solver
 !!$!
 !!$    ALLOCATE(mat%mumps)
 !!$    CALL csr2mumps(mat, mat%mumps)
   END SUBROUTINE to_csr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE to_zcsr_mat(mat, nlkeep)
 !
 !   Convert linked list spmat to csr matrice structure
 !
     TYPE(zcsr_mat)                :: mat
     LOGICAL, INTENT(in), OPTIONAL :: nlkeep
     INTEGER :: nnz_arr(mat%rank)
     INTEGER :: i, nnz, rank, s, e
     LOGICAL :: nlclean
 !
     nlclean = .TRUE.
     IF(PRESENT(nlkeep)) THEN
        nlclean = .NOT. nlkeep
     END IF
 !
 !    Allocate the csr matrix structure
 !
     nnz = get_count(mat%zspmat, nnz_arr)
     rank = mat%rank
     mat%nnz = nnz
     IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
     IF(ASSOCIATED(mat%idiag)) DEALLOCATE(mat%idiag)
     IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
     IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
     ALLOCATE(mat%irow(rank+1))
     ALLOCATE(mat%idiag(rank))
     ALLOCATE(mat%cols(nnz))
     ALLOCATE(mat%val(nnz))
 !
 !    Fill csr structure and optionally deallocate the sparse rows
 !
     mat%irow = 1
     DO i=1,rank
        mat%irow(i+1) = mat%irow(i) + nnz_arr(i)
        s = mat%irow(i)
        e = mat%irow(i+1)-1
        CALL getrow(mat%zspmat%row(i), mat%val(s:e), mat%cols(s:e))
        mat%idiag(i) = isearch(mat%cols(s:e), i) + s
        IF(nlclean) CALL destroy(mat%zspmat%row(i))
     END DO
 !!$!
 !!$!    MUMPS mat for direct solver
 !!$!
 !!$    ALLOCATE(mat%mumps)
 !!$    CALL csr2mumps(mat, mat%mumps)
   END SUBROUTINE to_zcsr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE csr2mumps(mat, mat_mumps)
 !
 !   Fill mumps structure (based on routine to_mumps_mat)
 !
     INCLUDE 'mpif.h'
     TYPE(csr_mat)   :: mat
     TYPE(mumps_mat) :: mat_mumps
 !
     INTEGER :: i, rank, s, e
     INTEGER :: comm, ierr, nnz_loc
 !
     CALL init(mat%rank, mat%nterms, mat_mumps)
 !
     comm = mat_mumps%mumps_par%COMM
     mat_mumps%mumps_par%ICNTL(18) = 3 ! Distributed assembled matrix
 !
 !    Allocate the Mumps matrix structure
 !      CSR format: (cols, irow, val) or (JCN, irow, A)
 !      COO format: (IRN, JCN, A)  or (IRN, cols, val)
 !
     rank = mat_mumps%rank
     nnz_loc = mat%nnz
     mat_mumps%nnz_start = 0
     CALL mpi_exscan(nnz_loc, mat_mumps%nnz_start, 1, MPI_INTEGER, MPI_SUM, comm, ierr)
     mat_mumps%nnz_start = mat_mumps%nnz_start + 1
     mat_mumps%nnz_end = mat_mumps%nnz_start + nnz_loc - 1
     mat_mumps%nnz_loc = nnz_loc
     CALL mpi_allreduce(nnz_loc, mat_mumps%nnz, 1, MPI_INTEGER, MPI_SUM, comm, ierr)
 !
     mat_mumps%mumps_par%N = rank
     mat_mumps%mumps_par%NZ_loc = nnz_loc
 !
     mat_mumps%cols => mat%cols
     mat_mumps%irow => mat%irow
     mat_mumps%val => mat%val
 !
 !     (A,JCN) picked from CSR mat
     mat_mumps%mumps_par%A_loc => mat_mumps%val
     mat_mumps%mumps_par%JCN_loc => mat_mumps%cols
 !
 !    Determine IRN array
     IF(ASSOCIATED(mat_mumps%mumps_par%IRN_loc)) DEALLOCATE(mat_mumps%mumps_par%IRN_loc)
     ALLOCATE(mat_mumps%mumps_par%IRN_loc(nnz_loc))
     DO i=mat_mumps%istart,mat_mumps%iend
        s = mat_mumps%irow(i) - mat_mumps%nnz_start + 1
        e = mat_mumps%irow(i+1) - mat_mumps%nnz_start
        mat_mumps%mumps_par%IRN_loc(s:e) = i
     END DO
     CALL destroy(mat_mumps%mat)
     NULLIFY(mat_mumps%mat)
-!    
+!
   END SUBROUTINE csr2mumps
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE zcsr2mumps(mat, mat_mumps)
 !
 !   Fill mumps structure (based on routine to_mumps_mat)
 !
     INCLUDE 'mpif.h'
     TYPE(zcsr_mat)   :: mat
     TYPE(zmumps_mat) :: mat_mumps
 !
     INTEGER :: i, rank, s, e
     INTEGER :: comm, ierr, nnz_loc
 !
     CALL init(mat%rank, mat%nterms, mat_mumps)
 !
     comm = mat_mumps%mumps_par%COMM
     mat_mumps%mumps_par%ICNTL(18) = 3 ! Distributed assembled matrix
 !
 !    Allocate the Mumps matrix structure
 !      CSR format: (cols, irow, val) or (JCN, irow, A)
 !      COO format: (IRN, JCN, A)  or (IRN, cols, val)
 !
     rank = mat_mumps%rank
     nnz_loc = mat%nnz
     mat_mumps%nnz_start = 0
     CALL mpi_exscan(nnz_loc, mat_mumps%nnz_start, 1, MPI_INTEGER, MPI_SUM, comm, ierr)
     mat_mumps%nnz_start = mat_mumps%nnz_start + 1
     mat_mumps%nnz_end = mat_mumps%nnz_start + nnz_loc - 1
     mat_mumps%nnz_loc = nnz_loc
     CALL mpi_allreduce(nnz_loc, mat_mumps%nnz, 1, MPI_INTEGER, MPI_SUM, comm, ierr)
 !
     mat_mumps%mumps_par%N = rank
     mat_mumps%mumps_par%NZ_loc = nnz_loc
 !
     mat_mumps%cols => mat%cols
     mat_mumps%irow => mat%irow
     mat_mumps%val => mat%val
 !
 !     (A,JCN) picked from CSR mat
     mat_mumps%mumps_par%A_loc => mat_mumps%val
     mat_mumps%mumps_par%JCN_loc => mat_mumps%cols
 !
 !    Determine IRN array
     IF(ASSOCIATED(mat_mumps%mumps_par%IRN_loc)) DEALLOCATE(mat_mumps%mumps_par%IRN_loc)
     ALLOCATE(mat_mumps%mumps_par%IRN_loc(nnz_loc))
     DO i=mat_mumps%istart,mat_mumps%iend
        s = mat_mumps%irow(i) - mat_mumps%nnz_start + 1
        e = mat_mumps%irow(i+1) - mat_mumps%nnz_start
        mat_mumps%mumps_par%IRN_loc(s:e) = i
     END DO
     CALL destroy(mat_mumps%mat)
     NULLIFY(mat_mumps%mat)
-!    
+!
   END SUBROUTINE zcsr2mumps
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE destroy_csr_mat(mat)
 !
 !   Deallocate csr mat
 !
     TYPE(csr_mat) :: mat
 !
     CALL destroy(mat%spmat)
     IF(mat%nnz.GT.0) THEN
        DEALLOCATE(mat%irow)
        DEALLOCATE(mat%idiag)
        DEALLOCATE(mat%cols)
        DEALLOCATE(mat%val)
     END IF
     mat%nnz = 0
   END SUBROUTINE destroy_csr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE destroy_zcsr_mat(mat)
 !
 !   Deallocate csr mat
 !
     TYPE(zcsr_mat) :: mat
 !
     CALL destroy(mat%zspmat)
     IF(mat%nnz.GT.0) THEN
        DEALLOCATE(mat%irow)
        DEALLOCATE(mat%idiag)
        DEALLOCATE(mat%cols)
        DEALLOCATE(mat%val)
     END IF
     mat%nnz = 0
   END SUBROUTINE destroy_zcsr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   FUNCTION vmx_csr_mat(mat, xarr, transa_in) RESULT(yarr)
 !
 !   Return product mat*x
 !
     TYPE(csr_mat)                          :: mat
     DOUBLE PRECISION, INTENT(in)           :: xarr(:)
     CHARACTER(len=1), OPTIONAL, INTENT(in) :: transa_in
     DOUBLE PRECISION                       :: yarr(SIZE(xarr))
 !
     DOUBLE PRECISION :: alpha=1.0d0, beta=0.0d0
     CHARACTER(len=1) :: transa
     CHARACTER(len=6) :: matdescra
     INTEGER :: n, i, j
 !
     n = mat%rank
     transa = 'N'
     IF(PRESENT(transa_in)) transa = transa_in
 !
 #ifdef MKL
     matdescra = 'g'
     CALL mkl_dcsrmv(transa, n, n, alpha, matdescra, mat%val, &
          &          mat%cols, mat%irow(1), mat%irow(2), xarr, &
          &          beta, yarr)
 #else
     yarr = 0.0d0
     DO i=1,mat%rank
        DO j=mat%irow(i), mat%irow(i+1)-1
           IF(transa .EQ. 'N') THEN
              yarr(i) = yarr(i) + mat%val(j)*xarr(mat%cols(j))
           ELSE
              yarr(mat%cols(j)) = yarr(mat%cols(j)) + mat%val(j)*xarr(i)
           END IF
        END DO
     END DO
 #endif
 !
   END FUNCTION vmx_csr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   FUNCTION vmx_zcsr_mat(mat, xarr, transa_in) RESULT(yarr)
 !
 !   Return product mat*x
 !
     TYPE(zcsr_mat)                         :: mat
     DOUBLE COMPLEX, INTENT(in)             :: xarr(:)
     CHARACTER(len=1), OPTIONAL, INTENT(in) :: transa_in
     DOUBLE COMPLEX                         :: yarr(SIZE(xarr))
 !
     DOUBLE COMPLEX :: alpha=(1.0d0,0.0d0), beta=(0.0d0,0.0d0)
     CHARACTER(len=1) :: transa
     CHARACTER(len=6) :: matdescra
     INTEGER :: n, i, j
 !
     n = mat%rank
     transa = 'N'
     IF(PRESENT(transa_in)) transa = transa_in
 !
 #ifdef MKL
     matdescra = 'g'
     CALL mkl_zcsrmv(transa, n, n, alpha, matdescra, mat%val, &
          &          mat%cols, mat%irow(1), mat%irow(2), xarr, &
          &          beta, yarr)
 #else
     yarr = 0.0d0
     DO i=1,mat%rank
        DO j=mat%irow(i), mat%irow(i+1)-1
           IF(transa .EQ. 'N') THEN
              yarr(i) = yarr(i) + mat%val(j)*xarr(mat%cols(j))
           ELSE
              yarr(mat%cols(j)) = yarr(mat%cols(j)) + mat%val(j)*xarr(i)
           END IF
        END DO
     END DO
 #endif
 !
   END FUNCTION vmx_zcsr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   FUNCTION vmx_csr_matn(mat, xarr, transa_in) RESULT(yarr)
 !
 !   Return product mat*x
 !
     TYPE(csr_mat)                          :: mat
     DOUBLE PRECISION, INTENT(in)           :: xarr(:,:)
     CHARACTER(len=1), OPTIONAL, INTENT(in) :: transa_in
     DOUBLE PRECISION                       :: yarr(SIZE(xarr,1),SIZE(xarr,2))
 !
     DOUBLE PRECISION :: alpha=1.0d0, beta=0.0d0
     CHARACTER(len=1) :: transa
     INTEGER :: n, nrhs, i, j
     CHARACTER(len=6) :: matdescra
 !
     n = mat%rank
     nrhs = SIZE(xarr,2)
     transa = 'N'
     IF(PRESENT(transa_in)) transa = transa_in
 !
 #ifdef MKL
     matdescra = 'g'
     CALL mkl_dcsrmm(transa, n, nrhs, n, alpha, matdescra, mat%val,&
          &           mat%cols, mat%irow(1), mat%irow(2), xarr, &
          &           n, beta, yarr, n)
 #else
     yarr = 0.0d0
     DO i=1,n
        DO j=mat%irow(i), mat%irow(i+1)-1
           IF(transa .EQ. 'N') THEN
              yarr(i,:) = yarr(i,:) + mat%val(j)*xarr(mat%cols(j),:)
           ELSE
              yarr(mat%cols(j),:) = yarr(mat%cols(j),:) + mat%val(j)*xarr(i,:)
           END IF
        END DO
     END DO
 #endif
 !
   END FUNCTION vmx_csr_matn
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   FUNCTION vmx_zcsr_matn(mat, xarr, transa_in) RESULT(yarr)
 !
 !   Return product mat*x
 !
     TYPE(zcsr_mat)                         :: mat
     DOUBLE COMPLEX, INTENT(in)             :: xarr(:,:)
     CHARACTER(len=1), OPTIONAL, INTENT(in) :: transa_in
     DOUBLE COMPLEX                          :: yarr(SIZE(xarr,1),SIZE(xarr,2))
 !
     DOUBLE COMPLEX   :: alpha=(1.0d0,0.0d0), beta=(0.0d0,0.0d0)
     CHARACTER(len=1) :: transa
     INTEGER :: n, nrhs, i, j
     CHARACTER(len=6) :: matdescra
 !
     n = mat%rank
     nrhs = SIZE(xarr,2)
     transa = 'N'
     IF(PRESENT(transa_in)) transa = transa_in
 !
 #ifdef MKL
     matdescra = 'g'
     CALL mkl_zcsrmm(transa, n, nrhs, n, alpha, matdescra, mat%val,&
          &           mat%cols, mat%irow(1), mat%irow(2), xarr, &
          &           n, beta, yarr, n)
 #else
     yarr = (0.0d0,0.0d0)
     DO i=1,n
        DO j=mat%irow(i), mat%irow(i+1)-1
           IF(transa .EQ. 'N') THEN
              yarr(i,:) = yarr(i,:) + mat%val(j)*xarr(mat%cols(j),:)
           ELSE
              yarr(mat%cols(j),:) = yarr(mat%cols(j),:) + mat%val(j)*xarr(i,:)
           END IF
        END DO
     END DO
 #endif
 !
   END FUNCTION vmx_zcsr_matn
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE conmat_1d_csr(spl, mat, coefeq, maxder)
 !
 !   Construction of FE matrix mat for 1D differential operator
 !   using spline spl
 !
     USE bsplines
     TYPE(csr_mat)              :: mat
     TYPE(spline1d), INTENT(in) :: spl
 !
-    INCLUDE '../../bsplines/src/conmat_1d.tpl'
+    INCLUDE 'conmat_1d.tpl'
   END SUBROUTINE conmat_1d_csr
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE conmat_1d_zcsr(spl, mat, coefeq, maxder)
 !
 !   Construction of FE matrix mat for 1D differential operator
 !   using spline spl
 !
     USE bsplines
     TYPE(zcsr_mat)              :: mat
     TYPE(spline1d), INTENT(in) :: spl
 !
-    INCLUDE '../../bsplines/src/zconmat_1d.tpl'
+    INCLUDE 'zconmat_1d.tpl'
   END SUBROUTINE conmat_1d_zcsr
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE conmat_2d_csr(spl, mat, coefeq, maxder, nat_order)
 !
 !   Construction of FE matrix mat for 2D differential operator
 !   using spline spl
 !
     USE bsplines
     TYPE(spline2d), INTENT(in) :: spl
     TYPE(csr_mat)              :: mat
 !
-    INCLUDE '../../bsplines/src/conmat.tpl'
+    INCLUDE 'conmat.tpl'
   END SUBROUTINE conmat_2d_csr
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE conmat_2d_zcsr(spl, mat, coefeq, maxder, nat_order)
 !
 !   Construction of FE matrix mat for 2D differential operator
 !   using spline spl
 !
     USE bsplines
     TYPE(spline2d), INTENT(in) :: spl
     TYPE(zcsr_mat)             :: mat
 !
-    INCLUDE '../../bsplines/src/zconmat.tpl'
+    INCLUDE 'zconmat.tpl'
   END SUBROUTINE conmat_2d_zcsr
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE femat_csr(spl, mat, coefeq, nterms)
 !
 !   Compute fe matrix
 !
     USE bsplines
     TYPE(spline1d), INTENT(in)   :: spl
     TYPE(csr_mat), INTENT(inout) :: mat
     INTEGER, INTENT(in)          :: nterms
     INTERFACE
        SUBROUTINE coefeq(x, idt, idw, c)
          DOUBLE PRECISION, INTENT(in) :: x
          INTEGER, INTENT(out) :: idt(:), idw(:)
          DOUBLE PRECISION, INTENT(out) :: c(:)
        END SUBROUTINE coefeq
     END INTERFACE
 !
     INTEGER :: nrank, nx, nidbas
 !
     CALL get_dim(spl, nrank, nx, nidbas)
     IF(spl%period) nrank = nx
     IF(mat%nnz.EQ.0) THEN
        WRITE(*,'(a,i0,a)') 'FEMAT: Initialize mat with ', &
             &    nterms, ' terms ...'
        CALL init(nrank, nterms, mat)
     END IF
     CALL conmat(spl, mat, coefeq)
   END SUBROUTINE femat_csr
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE put_csr_mat(fid, label, mat, str)
 !
 !   Write matrix to hdf5 file
 !
     USE futils
 !
     INTEGER, INTENT(in)                    :: fid
     CHARACTER(len=*), INTENT(in)           :: label
     TYPE(csr_mat)                          :: mat
     CHARACTER(len=*), OPTIONAL, INTENT(in) :: str
 !
     IF(PRESENT(str)) THEN
        CALL creatg(fid, label, str)
     ELSE
        CALL creatg(fid, label)
     END IF
     CALL attach(fid, label, 'RANK', mat%rank)
     CALL attach(fid, label, 'NNZ',  mat%nnz)
     CALL putarr(fid, TRIM(label)//'/irow', mat%irow)
     CALL putarr(fid, TRIM(label)//'/cols', mat%cols)
     CALL putarr(fid, TRIM(label)//'/idiag', mat%idiag)
     CALL putarr(fid, TRIM(label)//'/val', mat%val)
   END SUBROUTINE put_csr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE put_zcsr_mat(fid, label, mat, str)
 !
 !   Write matrix to hdf5 file
 !
     USE futils
 !
     INTEGER, INTENT(in)                    :: fid
     CHARACTER(len=*), INTENT(in)           :: label
     TYPE(zcsr_mat)                         :: mat
     CHARACTER(len=*), OPTIONAL, INTENT(in) :: str
 !
     IF(PRESENT(str)) THEN
        CALL creatg(fid, label, str)
     ELSE
        CALL creatg(fid, label)
     END IF
     CALL attach(fid, label, 'RANK', mat%rank)
     CALL attach(fid, label, 'NNZ',  mat%nnz)
     CALL putarr(fid, TRIM(label)//'/irow', mat%irow)
     CALL putarr(fid, TRIM(label)//'/cols', mat%cols)
     CALL putarr(fid, TRIM(label)//'/idiag', mat%idiag)
     CALL putarr(fid, TRIM(label)//'/val', mat%val)
   END SUBROUTINE put_zcsr_mat
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   DOUBLE PRECISION FUNCTION matnorm_csr(mat, p)
 !
 !   Compute matrix norm
 !
     TYPE(csr_mat), INTENT(in)              :: mat
     CHARACTER(len=*), OPTIONAL, INTENT(in) :: p
 !
     CHARACTER(len=4) :: norm_type
     INTEGER          :: i, j
     DOUBLE PRECISION :: temp(mat%rank)
 !
     norm_type = 'fro'
     IF(PRESENT(p)) norm_type = p
 !
     SELECT CASE (norm_type)
     CASE ('inf')
        DO i=1,mat%rank
           temp(i) = SUM(ABS(mat%val(mat%irow(i):mat%irow(i+1)-1)))
        END DO
        matnorm_csr = MAXVAL(temp)
     CASE ('1')
        temp = 0.0d0
        DO i=1,mat%rank
           DO j=mat%irow(i), mat%irow(i+1)-1
              temp(mat%cols(j)) = temp(mat%cols(j)) + ABS(mat%val(j))
           END DO
        END DO
        matnorm_csr = MAXVAL(temp)
     CASE('fro')
        matnorm_csr = SQRT(DOT_PRODUCT(mat%val, mat%val))
     END SELECT
   END FUNCTION matnorm_csr
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE full_to_csr(fullmat, mat)
 !
 !   Convert full rectangular matrix to csr mat
 !
     DOUBLE PRECISION, INTENT(inout) :: fullmat(:,:)
     TYPE(csr_mat), INTENT(out)      :: mat
 !
     INTEGER :: m, n, nnz
     INTEGER :: i, j, k
 !
     m = SIZE(fullmat,1)
     n = SIZE(fullmat,2)
     CALL init(m, 0, mat, ncols=n)
 !
 !   Determine nnz of matrix
 !
     IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
     IF(ASSOCIATED(mat%idiag)) DEALLOCATE(mat%idiag)
     IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
     IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
 !
     ALLOCATE(mat%irow(m+1))
     ALLOCATE(mat%idiag(m))
 !
 ! Clear matrix small elements of fullmat
     WHERE( ABS(fullmat) < 1.d-8) fullmat=0.0d0
 !
     mat%irow(1) = 1
     nnz = 0
     DO i=1,m
        DO j=1,n
-          IF(ABS(fullmat(i,j)).GT.EPSILON(0.0d0)) THEN 
+          IF(ABS(fullmat(i,j)).GT.EPSILON(0.0d0)) THEN
              nnz = nnz+1
              IF(m.EQ.n .AND. i.EQ.j) THEN ! Only for square matrix
                 mat%idiag(i) = nnz
              END IF
           END IF
        END DO
        mat%irow(i+1) = nnz+1
     END DO
 !
 !    Allocate and fill the csr matrix structure
 !
     mat%nnz = nnz
     ALLOCATE(mat%cols(nnz))
     ALLOCATE(mat%val(nnz))
     k=0
     DO i=1,m
        DO j=1,n
-          IF(ABS(fullmat(i,j)).GT.EPSILON(0.0d0)) THEN 
+          IF(ABS(fullmat(i,j)).GT.EPSILON(0.0d0)) THEN
              k=k+1
              mat%cols(k) = j
              mat%val(k) = fullmat(i,j)
           END IF
        END DO
     END DO
   END SUBROUTINE full_to_csr
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE full_to_zcsr(fullmat, mat)
 !
 !   Convert full rectangular matrix to csr mat
 !
     DOUBLE COMPLEX, INTENT(inout) :: fullmat(:,:)
     TYPE(zcsr_mat), INTENT(out)   :: mat
 !
     INTEGER :: m, n, nnz
     INTEGER :: i, j, k
 !
     m = SIZE(fullmat,1)
     n = SIZE(fullmat,2)
     CALL init(m, 0, mat, ncols=n)
 !
 !   Determine nnz of matrix
 !
     IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
     IF(ASSOCIATED(mat%idiag)) DEALLOCATE(mat%idiag)
     IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
     IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
 !
     ALLOCATE(mat%irow(m+1))
     ALLOCATE(mat%idiag(m))
 !
 ! Clear matrix small elements of fullmat
     WHERE( ABS(fullmat) < 1.d-8) fullmat=(0.0d0,0.0d0)
 !
     mat%irow(1) = 1
     nnz = 0
     DO i=1,m
        DO j=1,n
-          IF(ABS(fullmat(i,j)).GT.EPSILON(0.0d0)) THEN 
+          IF(ABS(fullmat(i,j)).GT.EPSILON(0.0d0)) THEN
              nnz = nnz+1
              IF(m.EQ.n .AND. i.EQ.j) THEN ! Only for square matrix
                 mat%idiag(i) = nnz
              END IF
           END IF
        END DO
        mat%irow(i+1) = nnz+1
     END DO
 !
 !    Allocate and fill the csr matrix structure
 !
     mat%nnz = nnz
     ALLOCATE(mat%cols(nnz))
     ALLOCATE(mat%val(nnz))
     k=0
     DO i=1,m
        DO j=1,n
-          IF(ABS(fullmat(i,j)).GT.EPSILON(0.0d0)) THEN 
+          IF(ABS(fullmat(i,j)).GT.EPSILON(0.0d0)) THEN
              k=k+1
              mat%cols(k) = j
              mat%val(k) = fullmat(i,j)
           END IF
        END DO
     END DO
   END SUBROUTINE full_to_zcsr
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE check_dom(mat)
 !
 !  Check whether mat is strict diagonal dominabce.
-! 
+!
     TYPE(csr_mat), INTENT(in) :: mat
     DOUBLE PRECISION :: arow(mat%rank), asum(mat%rank)
     INTEGER :: n, i, j1, j2, jdiag
 !
     n = mat%rank
     DO i=1,n
        j1 = mat%irow(i)
        jdiag = mat%idiag(i)
        j2 = mat%irow(i+1)-1
        asum(i) = SUM(ABS(mat%val(j1:j2))) / ABS(mat%val(jdiag)) - 1.0d0
     END DO
     WRITE(*,'(/a,1pe12.3)') 'Max of sum of off-diag', MAXVAL(asum)
 !
   END SUBROUTINE check_dom
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE kron_csr(mata, matb, matc)
 !
 !   Kronecker product of 2 CSR matrices
 !
     USE sparse, ONLY : isearch
     TYPE(csr_mat), INTENT(in)  :: mata, matb
     TYPE(csr_mat), INTENT(out) :: matc
 !
     INTEGER :: m1, n1, nnz1, m2, n2, nnz2, m, n, nnz
     INTEGER :: i,i1,i2,j1,s,s1,s2,e,e1,e2,k,nc2
 !
     m1 = mata%mrows
     n1 = mata%ncols
     nnz1 = mata%nnz
     m2 = matb%mrows
     n2 = matb%ncols
     nnz2 = matb%nnz
     m = m1*m2
     n = n1*n2
     nnz = nnz1*nnz2
 !
     CALL init(m, 0, matc, ncols=n)
     matc%nnz = nnz
     IF(ASSOCIATED(matc%irow)) DEALLOCATE(matc%irow)
     IF(ASSOCIATED(matc%idiag)) DEALLOCATE(matc%idiag)
     IF(ASSOCIATED(matc%cols)) DEALLOCATE(matc%cols)
     IF(ASSOCIATED(matc%val)) DEALLOCATE(matc%val)
     ALLOCATE(matc%irow(m+1))
     IF(m.EQ.n) THEN
        ALLOCATE(matc%idiag(m))     ! Only for square matrices
     END IF
     ALLOCATE(matc%cols(nnz))
     ALLOCATE(matc%val(nnz))
 !
     k = 0
     matc%irow(1) = 1
     DO i1=1,m1
        s1=mata%irow(i1)
        e1=mata%irow(i1+1)-1
        DO i2=1,m2
           s2=matb%irow(i2)
           e2=matb%irow(i2+1)-1
           nc2=e2-s2+1
           DO j1=s1,e1
              matc%val(k+1:k+nc2) = mata%val(j1)*matb%val(s2:e2)
              matc%cols(k+1:k+nc2) = (mata%cols(j1)-1)*n2 + matb%cols(s2:e2)
              k = k+nc2
              matc%irow((i1-1)*m2+i2+1) = k+1  ! Points to next row
           END DO
        END DO
     END DO
 !
 !  Search the diagonals
     DO i=1,matc%mrows
        s = matc%irow(i)
        e = matc%irow(i+1)-1
        matc%idiag(i) = isearch(matc%cols(s:e),i) + s
     END DO
   END SUBROUTINE kron_csr
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 END MODULE csr