!>
!> @file wsmp_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 wsmp_bsplines
!
!    WSMP_BSPLINES: Simple interface to the sparse direct solver WSMP.
!
!    T.M. Tran, CRPP-EPFL
!    November 2011
!
  USE sparse
  IMPLICIT NONE
!
  INTEGER, SAVE :: current_matid = -1
  INTEGER, SAVE :: last_matid = -1
!
  TYPE wsmp_param
     INTEGER          :: iparm(64)
     DOUBLE PRECISION :: dparm(64)
  END TYPE wsmp_param
!
  TYPE wsmp_mat
     INTEGER :: matid=-1
     INTEGER :: rank=0, nnz
     INTEGER :: nterms, kmat, nrhs
     LOGICAL :: nlsym
     LOGICAL :: nlpos
     LOGICAL :: nlforce_zero
     TYPE(spmat), POINTER      :: mat => NULL()
     INTEGER, POINTER          :: cols(:) => NULL()
     INTEGER, POINTER          :: irow(:) => NULL()
     INTEGER, POINTER          :: perm(:) => NULL()
     INTEGER, POINTER          :: invp(:) => NULL()
     INTEGER, POINTER          :: mrp(:) => NULL()
     DOUBLE PRECISION, POINTER :: diag(:) => NULL()
     DOUBLE PRECISION, POINTER :: val(:) => NULL()
     DOUBLE PRECISION, POINTER :: aux(:) => NULL()
     TYPE(wsmp_param)       :: p
  END TYPE wsmp_mat
!
  TYPE zwsmp_mat
     INTEGER :: matid=-1
     INTEGER :: rank=0, nnz
     INTEGER :: nterms, kmat, nrhs
     LOGICAL :: nlherm
     LOGICAL :: nlsym
     LOGICAL :: nlpos
     LOGICAL :: nlforce_zero
     TYPE(zspmat), POINTER     :: mat => NULL()
     INTEGER, POINTER          :: cols(:) => NULL()
     INTEGER, POINTER          :: irow(:) => NULL()
     INTEGER, POINTER          :: perm(:) => NULL()
     INTEGER, POINTER          :: invp(:) => NULL()
     INTEGER, POINTER          :: mrp(:) => NULL()
     DOUBLE COMPLEX, POINTER :: diag(:) => NULL()
     DOUBLE COMPLEX, POINTER :: val(:) => NULL()
     DOUBLE COMPLEX, POINTER :: aux(:) => NULL()
     TYPE(wsmp_param)       :: p
  END TYPE zwsmp_mat
!
  INTERFACE init
     MODULE PROCEDURE init_wsmp_mat, init_zwsmp_mat
  END INTERFACE init
!
  INTERFACE check_mat
     MODULE PROCEDURE check_wsmp_mat, check_zwsmp_mat
  END INTERFACE check_mat
!
  INTERFACE clear_mat
     MODULE PROCEDURE clear_wsmp_mat, clear_zwsmp_mat
  END INTERFACE clear_mat
!
  INTERFACE updtmat
     MODULE PROCEDURE updt_wsmp_mat, updt_zwsmp_mat
  END INTERFACE updtmat
!
  INTERFACE putele
     MODULE PROCEDURE putele_wsmp_mat, putele_zwsmp_mat
  END INTERFACE putele
!
  INTERFACE getele
     MODULE PROCEDURE getele_wsmp_mat, getele_zwsmp_mat
  END INTERFACE getele
!
  INTERFACE putrow
     MODULE PROCEDURE putrow_wsmp_mat, putrow_zwsmp_mat
  END INTERFACE putrow
!
  INTERFACE getrow
     MODULE PROCEDURE getrow_wsmp_mat, getrow_zwsmp_mat
  END INTERFACE getrow
!
  INTERFACE putcol
     MODULE PROCEDURE putcol_wsmp_mat, putcol_zwsmp_mat
  END INTERFACE putcol
!
  INTERFACE getcol
     MODULE PROCEDURE getcol_wsmp_mat, getcol_zwsmp_mat
  END INTERFACE getcol
!
  INTERFACE get_count
     MODULE PROCEDURE get_count_wsmp_mat, get_count_zwsmp_mat
  END INTERFACE get_count
!
  INTERFACE to_mat
     MODULE PROCEDURE to_wsmp_mat, to_zwsmp_mat
  END INTERFACE to_mat
!
  INTERFACE reord_mat
     MODULE PROCEDURE reord_wsmp_mat, reord_zwsmp_mat
  END INTERFACE reord_mat
!
  INTERFACE numfact
     MODULE PROCEDURE numfact_wsmp_mat, numfact_zwsmp_mat
  END INTERFACE numfact
!
  INTERFACE factor
     MODULE PROCEDURE factor_wsmp_mat, factor_zwsmp_mat
  END INTERFACE factor
!
  INTERFACE bsolve
     MODULE PROCEDURE bsolve_wsmp_mat1,  bsolve_wsmp_matn, &
          &           bsolve_zwsmp_mat1, bsolve_zwsmp_matn
  END INTERFACE bsolve
!
  INTERFACE vmx
     MODULE PROCEDURE vmx_wsmp_mat,  vmx_wsmp_matn, &
          &           vmx_zwsmp_mat, vmx_zwsmp_matn
  END INTERFACE vmx
!
  INTERFACE destroy
     MODULE PROCEDURE destroy_wsmp_mat, destroy_zwsmp_mat
 END INTERFACE destroy
!
 INTERFACE putmat
    MODULE PROCEDURE put_wsmp_mat, put_zwsmp_mat
 END INTERFACE putmat
!
 INTERFACE getmat
    MODULE PROCEDURE get_wsmp_mat, get_zwsmp_mat
 END INTERFACE getmat
!
 INTERFACE mcopy
    MODULE PROCEDURE mcopy_wsmp_mat,  mcopy_zwsmp_mat
 END INTERFACE mcopy
!
 INTERFACE maddto
    MODULE PROCEDURE maddto_wsmp_mat, maddto_zwsmp_mat
 END INTERFACE maddto
!
 INTERFACE psum_mat
    MODULE PROCEDURE psum_wsmp_mat, psum_zwsmp_mat
 END INTERFACE psum_mat
!
 INTERFACE p2p_mat
    MODULE PROCEDURE p2p_wsmp_mat, p2p_zwsmp_mat
 END INTERFACE p2p_mat
!
CONTAINS
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE init_wsmp_mat(n, nterms, mat, kmat, nlsym, nlpos, &
       &                   nlforce_zero)
!
!   Initialize an empty sparse wsmp matrix
!
    INTEGER, INTENT(in)           :: n, nterms
    TYPE(wsmp_mat)                :: mat
    INTEGER, OPTIONAL, INTENT(in) :: kmat
    LOGICAL, OPTIONAL, INTENT(in) :: nlsym
    LOGICAL, OPTIONAL, INTENT(in) :: nlpos
    LOGICAL, OPTIONAL, INTENT(in) :: nlforce_zero
    INTEGER :: info
    INTEGER :: idummy = 0
    DOUBLE PRECISION :: dummy = 0.0d0
!
!  Store away (valid) current matrix id
!
    IF(current_matid .GE. 0) THEN
       CALL wstoremat(current_matid, info)
       IF(info.NE.0) THEN
          WRITE(*,'(a,i4)') 'INIT: WSTOREMAT failed WITH error', info
          STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
       END IF
    END IF
    last_matid = last_matid+1
    mat%matid = last_matid
    current_matid = mat%matid
!
!  Initialize sparse matrice structure
!
    mat%rank = n
    mat%nterms = nterms
    mat%nnz = 0
    mat%nlsym = .FALSE.
    mat%nlpos = .TRUE.
    mat%nrhs = 1
    mat%nlforce_zero = .TRUE. ! Do not remove existing nodes with zero val
    IF(PRESENT(kmat)) mat%kmat = kmat
    IF(PRESENT(nlsym)) mat%nlsym = nlsym
    IF(PRESENT(nlpos)) mat%nlpos = nlpos
    IF(PRESENT(nlforce_zero)) mat%nlforce_zero = nlforce_zero
    IF(ASSOCIATED(mat%mat)) DEALLOCATE(mat%mat)
    ALLOCATE(mat%mat)
    CALL init(n, mat%mat)
!
!  Fill 'iparm' and 'dparm' with default values
!
    mat%p%iparm(1:3) = 0
    IF(mat%nlsym) THEN
       CALL wssmp(mat%rank, mat%irow, mat%cols, mat%val, mat%diag, &
            &      mat%perm, mat%invp, dummy, mat%rank, mat%nrhs, &
            &      mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
            &      mat%p%dparm)
       IF(mat%nlpos) THEN
          mat%p%iparm(31) = 0
       ELSE
!!$          mat%p%iparm(31) = 1  ! LDL^T without pivoting
          mat%p%iparm(31) = 2  ! LDL^T with pivoting
       END IF
    ELSE
       CALL wgsmp(mat%rank, mat%irow, mat%cols, mat%val, dummy, mat%rank, &
            &     mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
   END IF
    IF(mat%p%iparm(64).NE.0) THEN
       WRITE(*,'(a,i4)') 'INIT: Initialization failed with error', &
            &             mat%p%iparm(64)
       STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
    END IF
!!$    WRITE(*,'(/a/(10i8))') 'iparm', mat%p%iparm
!!$    WRITE(*,'(/a/(10(1pe8.1)))') 'dparm', mat%p%dparm
!
    CALL setup_wsmp(mat%p%iparm, mat%p%dparm)
!
  CONTAINS
    SUBROUTINE setup_wsmp(iparm, dparm)
      INTEGER          :: iparm(:)
      DOUBLE PRECISION :: dparm(:)
    END SUBROUTINE setup_wsmp
  END SUBROUTINE init_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE init_zwsmp_mat(n, nterms, mat, kmat, nlsym, nlherm, &
       &                   nlpos, nlforce_zero)
!
!   Initialize an empty sparse wsmp matrix
!
    INTEGER, INTENT(in)           :: n, nterms
    TYPE(zwsmp_mat)               :: mat
    INTEGER, OPTIONAL, INTENT(in) :: kmat
    LOGICAL, OPTIONAL, INTENT(in) :: nlsym
    LOGICAL, OPTIONAL, INTENT(in) :: nlherm
    LOGICAL, OPTIONAL, INTENT(in) :: nlpos
    LOGICAL, OPTIONAL, INTENT(in) :: nlforce_zero
    INTEGER :: info
    INTEGER :: idummy = 0
    DOUBLE COMPLEX :: dummy = 0.0d0
!
!  Store away (valid) current matrix id
!
    IF(current_matid .GE. 0) THEN
       CALL wstoremat(current_matid, info)
       IF(info.NE.0) THEN
          WRITE(*,'(a,i4)') 'INIT: WSTOREMAT failed WITH error', info
          STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
       END IF
    END IF
    last_matid = last_matid+1
    mat%matid = last_matid
    current_matid = mat%matid
!
!  Initialize sparse matrice structure
!
    mat%rank = n
    mat%nterms = nterms
    mat%nnz = 0
    mat%nlsym = .FALSE.
    mat%nlherm = .FALSE.
    mat%nlpos = .TRUE.
    mat%nrhs = 1
    mat%nlforce_zero = .TRUE. ! Do not remove existing nodes with zero val
    IF(PRESENT(kmat)) mat%kmat = kmat
    IF(PRESENT(nlsym)) mat%nlsym = nlsym
    IF(PRESENT(nlherm)) mat%nlherm = nlherm
    IF(PRESENT(nlpos)) mat%nlpos = nlpos
    IF(PRESENT(nlforce_zero)) mat%nlforce_zero = nlforce_zero
    IF(ASSOCIATED(mat%mat)) DEALLOCATE(mat%mat)
    ALLOCATE(mat%mat)
    CALL init(n, mat%mat)
!
!  Fill 'iparm' and 'dparm' with default values
!
    mat%p%iparm(1:3) = 0
    IF(mat%nlherm .OR. mat%nlsym) THEN
       CALL zssmp(mat%rank, mat%irow, mat%cols, mat%val, mat%diag, &
            &      mat%perm, mat%invp, dummy, mat%rank, mat%nrhs, &
            &      mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
            &      mat%p%dparm)
       IF(mat%nlherm) THEN
          IF(mat%nlpos) THEN
             mat%p%iparm(31) = 0  ! hermitian, positive definite 
          ELSE
             mat%p%iparm(31) = 2  ! hermitian, no-definite, LDL^T with pivoting
          END IF
       ELSE
          mat%p%iparm(31) = 3     ! non-hermitian, symmetric
       END IF
    ELSE
       CALL zgsmp(mat%rank, mat%irow, mat%cols, mat%val, dummy, mat%rank, &
            &     mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
   END IF
    IF(mat%p%iparm(64).NE.0) THEN
       WRITE(*,'(a,i4)') 'INIT: Initialization failed with error', &
            &             mat%p%iparm(64)
       STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
    END IF
!!$    WRITE(*,'(/a/(10i8))') 'iparm', mat%p%iparm
!!$    WRITE(*,'(/a/(10(1pe8.1)))') 'dparm', mat%p%dparm
!
    CALL setup_wsmp(mat%p%iparm, mat%p%dparm)
!
  CONTAINS
    SUBROUTINE setup_wsmp(iparm, dparm)
      INTEGER          :: iparm(:)
      DOUBLE PRECISION :: dparm(:)
    END SUBROUTINE setup_wsmp
  END SUBROUTINE init_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE check_wsmp_mat(mat)
!
!  Check matrice id and  recall the matrice if not current
!
    TYPE(wsmp_mat) :: mat
    INTEGER :: info
!
    IF(.NOT.mat%nlsym) THEN
       IF( mat%matid.NE.current_matid ) THEN      
          WRITE(*,'(a)') "Processing multi matrices is not possible "// &
               &         "for non-symetric matrices."
          STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
       ELSE
          RETURN
       END IF
    END IF
!
    IF( mat%matid.NE.current_matid ) THEN
       IF(current_matid .GE. 0) THEN
          CALL wstoremat(current_matid, info)
          IF(info.NE.0) THEN
             WRITE(*,'(a,i3,a,i4)') 'Store matrix', current_matid, &
                  &               ' failed with error', info
             STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
          END IF
       END IF
       CALL wrecallmat(mat%matid, info)
       IF(info.NE.0) THEN
          WRITE(*,'(a,i3,a,i4)') 'Recall matrix', mat%matid, &
               &              ' failed with error', info
          STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
       END IF
       current_matid = mat%matid
    END IF
  END SUBROUTINE check_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE check_zwsmp_mat(mat)
!
!  Check matrice id and  recall the matrice if not current
!
    TYPE(zwsmp_mat) :: mat
    INTEGER :: info
!
    IF(.NOT.mat%nlsym .AND. .NOT.mat%nlherm ) THEN
       IF( mat%matid.NE.current_matid ) THEN      
          WRITE(*,'(a)') "Processing multi matrices is not possible "// &
               &         "for non-symetric/non-hermitian matrices."
          STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
       ELSE
          RETURN
       END IF
    END IF
!
    IF( mat%matid.NE.current_matid ) THEN
       IF(current_matid .GE. 0) THEN
          CALL wstoremat(current_matid, info)
          IF(info.NE.0) THEN
             WRITE(*,'(a,i3,a,i4)') 'Store matrix', current_matid, &
                  &               ' failed with error', info
             STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
          END IF
       END IF
       CALL wrecallmat(mat%matid, info)
       IF(info.NE.0) THEN
          WRITE(*,'(a,i3,a,i4)') 'Recall matrix', mat%matid, &
               &              ' failed with error', info
          STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
       END IF
       current_matid = mat%matid
    END IF
  END SUBROUTINE check_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE clear_wsmp_mat(mat)
!
!   Clear matrix, keeping its sparse structure unchanged
!
    TYPE(wsmp_mat) :: mat
!
    mat%val = 0.0d0
  END SUBROUTINE clear_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE clear_zwsmp_mat(mat)
!
!   Clear matrix, keeping its sparse structure unchanged
!
    TYPE(zwsmp_mat) :: mat
!
    mat%val = (0.0d0, 0.0d0)
  END SUBROUTINE clear_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE updt_wsmp_mat(mat, i, j, val)
!
!   Update element Aij of wsmp  matrix
!
    TYPE(wsmp_mat)            :: mat
    INTEGER, INTENT(in)          :: i, j
    DOUBLE PRECISION, INTENT(in) :: val
    INTEGER :: k, s, e
!
    IF(mat%nlsym) THEN   ! Store only upper part for symmetric matrices
       IF(i.GT.j) RETURN
    END IF
!
    IF(ASSOCIATED(mat%mat)) THEN
       CALL updtmat(mat%mat, 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 wsmp_mod ***'
       END IF
    END IF
  END SUBROUTINE updt_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE updt_zwsmp_mat(mat, i, j, val)
!
!   Update element Aij of wsmp  matrix
!
    TYPE(zwsmp_mat)            :: mat
    INTEGER, INTENT(in)        :: i, j
    DOUBLE COMPLEX, INTENT(in) :: val
    INTEGER :: k, s, e
!
    IF(mat%nlherm .OR. mat%nlsym) THEN   ! Store only upper part
       IF(i.GT.j) RETURN
    END IF
!
    IF(ASSOCIATED(mat%mat)) THEN
       CALL updtmat(mat%mat, 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
          IF(mat%nlherm) THEN
             mat%val(s+k) =  mat%val(s+k)+CONJG(val)  ! CSR-UT* = CSC-LT
          ELSE
             mat%val(s+k) =  mat%val(s+k)+val
       END IF
       ELSE
          WRITE(*,'(a,2i6)') 'UPDT: i, j out of range ', i, j
          STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
       END IF
    END IF
  END SUBROUTINE updt_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putele_wsmp_mat(mat, i, j, val)
!
!  Put element (i,j) of matrix
!
    TYPE(wsmp_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%nlsym) THEN
       IF( i.GT.j ) THEN    ! Lower triangular part
          iput = j
          jput = i
       END IF
    END IF
!
    IF(ASSOCIATED(mat%mat)) THEN
       CALL putele(mat%mat, 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
             STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
          END IF
       END IF
   END IF
  END SUBROUTINE putele_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putele_zwsmp_mat(mat, i, j, val)
!
!  Put element (i,j) of matrix
!
    TYPE(zwsmp_mat)            :: mat
    INTEGER, INTENT(in)        :: i, j
    DOUBLE COMPLEX, INTENT(in) :: val
    DOUBLE COMPLEX :: valput
    INTEGER :: iput, jput
    INTEGER :: k, s, e
!
    iput = i
    jput = j
    valput = val
    IF(mat%nlsym .OR. mat%nlherm) THEN
       IF( i.GT.j ) THEN    ! Lower triangular part
          iput = j
          jput = i
          IF(mat%nlherm) THEN
             valput = CONJG(val)
          ELSE
             valput = val
          END IF
       END IF
    END IF
!
    IF(ASSOCIATED(mat%mat)) THEN
       CALL putele(mat%mat, iput, jput, valput, &
            &      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
          IF(mat%nlherm) THEN
             mat%val(s+k) = CONJG(valput)   ! CSR-UT* = CSC-LT
          ELSE
             mat%val(s+k) =  valput
          END IF
       ELSE
          IF(ABS(val) .GT. EPSILON(0.0d0)) THEN   ! Ok for zero val
             WRITE(*,'(a,2i6)') 'PUTELE: i, j out of range ', i, j
             STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
          END IF
       END IF
   END IF
 END SUBROUTINE putele_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getele_wsmp_mat(mat, i, j, val)
!
!   Get element (i,j) of sparse matrix
!
    TYPE(wsmp_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%nlsym) THEN
       IF( i.GT.j ) THEN    ! Lower triangular part
          iget = j
          jget = i
       END IF
    END IF
!
    IF(ASSOCIATED(mat%mat)) THEN
       CALL getele(mat%mat, 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) 
       ELSE
          val = 0.0d0   ! Assume zero val if not found
       END IF
    END IF
  END SUBROUTINE getele_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getele_zwsmp_mat(mat, i, j, val)
!
!   Get element (i,j) of sparse matrix
!
    TYPE(zwsmp_mat)             :: mat
    INTEGER, INTENT(in)         :: i, j
    DOUBLE COMPLEX, INTENT(out) :: val
    DOUBLE COMPLEX :: valget
    INTEGER :: iget, jget
    INTEGER :: k, s, e
!
    iget = i
    jget = j
    IF(mat%nlherm .OR. mat%nlsym) THEN
       IF( i.GT.j ) THEN    ! Lower triangular part
          iget = j
          jget = i
       END IF
    END IF
!
    IF(ASSOCIATED(mat%mat)) THEN
       CALL getele(mat%mat, iget, jget, valget)
    ELSE
       s = mat%irow(iget)
       e = mat%irow(iget+1)-1
       k = isearch(mat%cols(s:e), jget)
       IF( k.GE.0 ) THEN
          IF(mat%nlherm) THEN
             valget = CONJG(mat%val(s+k))   ! CSR-UT* = CSC-LT
          ELSE
             valget = mat%val(s+k) 
          END IF
       ELSE
          valget = (0.0d0,0.0d0)   ! Assume zero val if not found
       END IF
    END IF
    val = valget
    IF( i.GT.j ) THEN
       IF(mat%nlherm) THEN 
          val = CONJG(valget)
       END IF
    END IF
  END SUBROUTINE getele_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putrow_wsmp_mat(amat, i, arr)
!
! Put a row into sparse matrix
!
    TYPE(wsmp_mat), INTENT(inout) :: amat
    INTEGER, INTENT(in)              :: i
    DOUBLE PRECISION, INTENT(in)     :: arr(:)
    INTEGER :: j
!
    DO j=1,amat%rank
       CALL putele(amat, i, j, arr(j))
    END DO
  END SUBROUTINE putrow_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putrow_zwsmp_mat(amat, i, arr)
!
! Put a row into sparse matrix
!
    TYPE(zwsmp_mat), INTENT(inout) :: amat
    INTEGER, INTENT(in)            :: i
    DOUBLE COMPLEX, INTENT(in)     :: arr(:)
    INTEGER :: j
!
    DO j=1,amat%rank
       CALL putele(amat, i, j, arr(j))
    END DO
  END SUBROUTINE putrow_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getrow_wsmp_mat(amat, i, arr)
!
!   Get a row from sparse matrix
!
    TYPE(wsmp_mat), INTENT(in)  :: amat
    INTEGER, INTENT(in)            :: i
    DOUBLE PRECISION, INTENT(out)  :: arr(:)
    INTEGER :: j
!
    DO j=1,amat%rank
       CALL getele(amat, i, j, arr(j))
    END DO
  END SUBROUTINE getrow_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getrow_zwsmp_mat(amat, i, arr)
!
!   Get a row from sparse matrix
!
    TYPE(zwsmp_mat), INTENT(in)  :: amat
    INTEGER, INTENT(in)          :: i
    DOUBLE COMPLEX, INTENT(out)  :: arr(:)
    INTEGER :: j
!
    DO j=1,amat%rank
       CALL getele(amat, i, j, arr(j))
    END DO
  END SUBROUTINE getrow_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putcol_wsmp_mat(amat, j, arr)
!
! Put a column into sparse matrix
!
    TYPE(wsmp_mat), INTENT(inout) :: amat
    INTEGER, INTENT(in)              :: j
    DOUBLE PRECISION, INTENT(in)     :: arr(:)
    INTEGER :: i
!
    DO i=1,amat%rank
       CALL putele(amat, i, j, arr(i))
    END DO
  END SUBROUTINE putcol_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE putcol_zwsmp_mat(amat, j, arr)
!
! Put a column into sparse matrix
!
    TYPE(zwsmp_mat), INTENT(inout) :: amat
    INTEGER, INTENT(in)            :: j
    DOUBLE COMPLEX, INTENT(in)     :: arr(:)
    INTEGER :: i
!
    DO i=1,amat%rank
       CALL putele(amat, i, j, arr(i))
    END DO
  END SUBROUTINE putcol_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getcol_wsmp_mat(amat, j, arr)
!
!   Get a column from sparse matrix
!
    TYPE(wsmp_mat), INTENT(in)  :: amat
    INTEGER, INTENT(in)            :: j
    DOUBLE PRECISION, INTENT(out)  :: arr(:)
    INTEGER :: i
!
    DO i=1,amat%rank
       CALL getele(amat, i, j, arr(i))
    END DO
  END SUBROUTINE getcol_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE getcol_zwsmp_mat(amat, j, arr)
!
!   Get a column from sparse matrix
!
    TYPE(zwsmp_mat), INTENT(in)  :: amat
    INTEGER, INTENT(in)          :: j
    DOUBLE COMPLEX, INTENT(out)  :: arr(:)
    INTEGER :: i
!
    DO i=1,amat%rank
       CALL getele(amat, i, j, arr(i))
    END DO
  END SUBROUTINE getcol_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  INTEGER FUNCTION get_count_wsmp_mat(mat, nnz)
!
!   Number of non-zeros in sparse matrix
!
    TYPE(wsmp_mat)                 :: mat
    INTEGER, INTENT(out), OPTIONAL :: nnz(:)
    INTEGER :: i
!
    IF(ASSOCIATED(mat%mat)) THEN
       get_count_wsmp_mat = get_count(mat%mat, nnz)
    ELSE
       get_count_wsmp_mat = mat%nnz
       IF(PRESENT(nnz)) THEN
          DO i=1,mat%rank
             nnz(i) = mat%irow(i+1)-mat%irow(i)
          END DO
       END IF
    END IF
  END FUNCTION get_count_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  INTEGER FUNCTION get_count_zwsmp_mat(mat, nnz)
!
!   Number of non-zeros in sparse matrix
!
    TYPE(zwsmp_mat)                :: mat
    INTEGER, INTENT(out), OPTIONAL :: nnz(:)
    INTEGER :: i
!
    IF(ASSOCIATED(mat%mat)) THEN
       get_count_zwsmp_mat = get_count(mat%mat, nnz)
    ELSE
       get_count_zwsmp_mat = mat%nnz
       IF(PRESENT(nnz)) THEN
          DO i=1,mat%rank
             nnz(i) = mat%irow(i+1)-mat%irow(i)
          END DO
       END IF
    END IF
  END FUNCTION get_count_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE to_wsmp_mat(mat, nlkeep)
!
!   Convert linked list spmat to wsmp matrice structure
!
    TYPE(wsmp_mat) :: mat
    LOGICAL, INTENT(in), OPTIONAL :: nlkeep
    INTEGER :: i, nnz, rank, s, e
    LOGICAL :: nlclean
!
    nlclean = .TRUE.
    IF(PRESENT(nlkeep)) THEN
       nlclean = .NOT. nlkeep
    END IF
!
!    Allocate the WSMP matrix structure
!
    nnz = get_count(mat)
    rank = mat%rank
    mat%nnz = nnz
    IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
    IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
    IF(ASSOCIATED(mat%perm)) DEALLOCATE(mat%perm)
    IF(ASSOCIATED(mat%invp)) DEALLOCATE(mat%invp)
    IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
    ALLOCATE(mat%val(nnz))
    IF(mat%nlsym) THEN
       ALLOCATE(mat%perm(rank))
       ALLOCATE(mat%invp(rank))
    END IF
    ALLOCATE(mat%irow(rank+1))
    ALLOCATE(mat%cols(nnz))
!
!    Fill WSMP structure and deallocate the sparse rows
!
    mat%irow = 1
    DO i=1,rank
       mat%irow(i+1) = mat%irow(i) + mat%mat%row(i)%nnz
       s = mat%irow(i)
       e = mat%irow(i+1)-1
       CALL getrow(mat%mat%row(i), mat%val(s:e), mat%cols(s:e))
       IF(nlclean) CALL destroy(mat%mat%row(i))
    END DO
    IF(nlclean) DEALLOCATE(mat%mat)
  END SUBROUTINE to_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE to_zwsmp_mat(mat, nlkeep)
!
!   Convert linked list spmat to wsmp matrice structure
!
    TYPE(zwsmp_mat) :: mat
    LOGICAL, INTENT(in), OPTIONAL :: nlkeep
    INTEGER :: i, nnz, rank, s, e
    LOGICAL :: nlclean
!
    nlclean = .TRUE.
    IF(PRESENT(nlkeep)) THEN
       nlclean = .NOT. nlkeep
    END IF
!
!    Allocate the WSMP matrix structure
!
    nnz = get_count(mat)
    rank = mat%rank
    mat%nnz = nnz
    IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
    IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
    IF(ASSOCIATED(mat%perm)) DEALLOCATE(mat%perm)
    IF(ASSOCIATED(mat%invp)) DEALLOCATE(mat%invp)
    IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
    ALLOCATE(mat%val(nnz))
    IF(mat%nlsym) THEN
       ALLOCATE(mat%perm(rank))
       ALLOCATE(mat%invp(rank))
    END IF
    ALLOCATE(mat%irow(rank+1))
    ALLOCATE(mat%cols(nnz))
!
!    Fill WSMP structure and deallocate the sparse rows
!
    mat%irow = 1
    DO i=1,rank
       mat%irow(i+1) = mat%irow(i) + mat%mat%row(i)%nnz
       s = mat%irow(i)
       e = mat%irow(i+1)-1
       CALL getrow(mat%mat%row(i), mat%val(s:e), mat%cols(s:e))
       IF(nlclean) CALL destroy(mat%mat%row(i))
    END DO
    IF(mat%nlherm) THEN
       mat%val(:) = CONJG(mat%val(:))  ! CSR-UT* = CSC-LT
    END IF
    IF(nlclean) DEALLOCATE(mat%mat)
  END SUBROUTINE to_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE reord_wsmp_mat(mat)
!
!   Reordering and symbolic factorization
!
    TYPE(wsmp_mat) :: mat
    DOUBLE PRECISION :: dummy
!
!    Recall the matrice if not current
!
    CALL check_mat(mat)
!
    IF(mat%nlsym) THEN
       mat%p%iparm(2) = 1  ! Ordering
       mat%p%iparm(3) = 2  ! Symbolic factorization
       CALL wssmp(mat%rank, mat%irow, mat%cols, mat%val, mat%diag, &
            &      mat%perm, mat%invp, dummy, mat%rank, mat%nrhs, &
            &      mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
            &      mat%p%dparm)
    ELSE
       mat%p%iparm(2) = 1  ! Analysis and reordering
       mat%p%iparm(3) = 1
       CALL wgsmp(mat%rank, mat%irow, mat%cols, mat%val, dummy, mat%rank, &
            &     mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
    END IF
    IF(mat%p%iparm(64).NE.0) THEN
       WRITE(*,'(a,i4)') 'REORD: Reordering failed with error', mat%p%iparm(64)
       STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
    END IF
  END SUBROUTINE reord_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE reord_zwsmp_mat(mat)
!
!   Reordering and symbolic factorization
!
    TYPE(zwsmp_mat)             :: mat
    DOUBLE COMPLEX :: dummy
!
!    Recall the matrice if not current
!
    CALL check_mat(mat)
!
    IF(mat%nlsym .OR. mat%nlherm) THEN
       mat%p%iparm(2) = 1  ! Ordering
       mat%p%iparm(3) = 2  ! Symbolic factorization
       CALL zssmp(mat%rank, mat%irow, mat%cols, mat%val, mat%diag, &
            &      mat%perm, mat%invp, dummy, mat%rank, mat%nrhs, &
            &      mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
            &      mat%p%dparm)
!!$       WRITE(*,'(a,i3/(10i8))') 'REORD: matrice', mat%matid, mat%perm
    ELSE
       mat%p%iparm(2) = 1  ! Analysis and reordering
       mat%p%iparm(3) = 1
       CALL zgsmp(mat%rank, mat%irow, mat%cols, mat%val, dummy, mat%rank, &
            &     mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
    END IF
    IF(mat%p%iparm(64).NE.0) THEN
       WRITE(*,'(a,i4)') 'REORD: Reordering failed with error', mat%p%iparm(64)
       STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
    END IF
  END SUBROUTINE reord_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE numfact_wsmp_mat(mat)
!
!   Numerical factorization
!
    TYPE(wsmp_mat) :: mat
    DOUBLE PRECISION :: dummy
!
!    Recall the matrice if not current
!
    CALL check_mat(mat)
!
    IF(mat%nlsym) THEN
       mat%p%iparm(2) = 3   ! Numerical factorization
       mat%p%iparm(3) = 3
       CALL wssmp(mat%rank, mat%irow, mat%cols, mat%val, mat%diag, &
            &      mat%perm, mat%invp, dummy, mat%rank, mat%nrhs, &
            &      mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
            &      mat%p%dparm)
    ELSE
       mat%p%iparm(2) = 2  ! Factorization
       mat%p%iparm(3) = 2
       CALL wgsmp(mat%rank, mat%irow, mat%cols, mat%val, dummy, mat%rank, &
            &     mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
    END IF
    IF(mat%p%iparm(64).NE.0) THEN
       WRITE(*,'(a,i4)') 'NUMFACT: Num. Factor failed with error', mat%p%iparm(64)
       STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
    END IF
  END SUBROUTINE numfact_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE numfact_zwsmp_mat(mat)
!
!   Numerical factorization
!
    TYPE(zwsmp_mat) :: mat
    DOUBLE COMPLEX :: dummy
!
!    Recall the matrice if not current
!
    CALL check_mat(mat)
!
    IF(mat%nlsym .OR. mat%nlherm) THEN
       mat%p%iparm(2) = 3   ! Numerical factorization
       mat%p%iparm(3) = 3
       CALL zssmp(mat%rank, mat%irow, mat%cols, mat%val, mat%diag, &
            &      mat%perm, mat%invp, dummy, mat%rank, mat%nrhs, &
            &      mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
            &      mat%p%dparm)
    ELSE
       mat%p%iparm(2) = 2  ! Factorization
       mat%p%iparm(3) = 2
       CALL zgsmp(mat%rank, mat%irow, mat%cols, mat%val, dummy, mat%rank, &
            &     mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
    END IF
    IF(mat%p%iparm(64).NE.0) THEN
       WRITE(*,'(a,i4)') 'NUMFACT: Num. Factor failed with error', mat%p%iparm(64)
       STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
    END IF
  END SUBROUTINE numfact_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE factor_wsmp_mat(mat, nlreord)
!
!  Factor (create  +reorder + factor) a wsmp_mat matrix
!
    TYPE(wsmp_mat)                :: mat
    LOGICAL, OPTIONAL, INTENT(in) :: nlreord
    LOGICAL :: mlreord
!----------------------------------------------------------------------
!               1.0  Creation from the sparse rows
!
    IF(ASSOCIATED(mat%mat)) THEN
       CALL to_mat(mat)
    END IF
!----------------------------------------------------------------------
!               2.0  Reordering and symbolic factorization phase
!
    mlreord = .TRUE.
    IF(PRESENT(nlreord)) mlreord = nlreord
    IF(mlreord) THEN
       CALL reord_mat(mat)
    END IF
!----------------------------------------------------------------------
!               3.0  Numerical factorization
!
    CALL numfact(mat)
  END SUBROUTINE factor_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE factor_zwsmp_mat(mat, nlreord)
!
!  Factor (create  +reorder + factor) a wsmp_mat matrix
!
    TYPE(zwsmp_mat)               :: mat
    LOGICAL, OPTIONAL, INTENT(in) :: nlreord
    LOGICAL :: mlreord
!----------------------------------------------------------------------
!               1.0  Creation from the sparse rows
!
    IF(ASSOCIATED(mat%mat)) THEN
       CALL to_mat(mat)
    END IF
!----------------------------------------------------------------------
!               2.0  Reordering and symbolic factorization phase
!
    mlreord = .TRUE.
    IF(PRESENT(nlreord)) mlreord = nlreord
    IF(mlreord) THEN
       CALL reord_mat(mat)
    END IF
!----------------------------------------------------------------------
!               3.0  Numerical factorization
!
    CALL numfact(mat)
  END SUBROUTINE factor_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE bsolve_wsmp_mat1(mat, rhs, sol, nref)
!
!   Backsolve, using Wsmp
!
    TYPE(wsmp_mat)             :: mat
    DOUBLE PRECISION           :: rhs(:)
    DOUBLE PRECISION, OPTIONAL :: sol(:)
    INTEGER, OPTIONAL          :: nref
    DOUBLE PRECISION :: dummy
!
!    Recall the matrice if not current
!
    CALL check_mat(mat)
!
    IF(mat%nlsym) THEN
       mat%p%iparm(2) = 4  ! Back substitution
       mat%p%iparm(3) = 4
    ELSE
       mat%p%iparm(2) = 3  ! Back substitution
       mat%p%iparm(3) = 3
    END IF
    mat%p%iparm(6) = 0  ! Max numbers of iterative refinement steps
    IF(PRESENT(nref)) THEN
       IF(mat%nlsym) THEN
          mat%p%iparm(3) = 5
       ELSE
          mat%p%iparm(3) = 4
       END IF
       mat%p%iparm(6) = nref
    END IF
    mat%nrhs = 1
    IF(PRESENT(sol)) THEN
       sol = rhs
       IF(mat%nlsym) THEN
          CALL wssmp(mat%rank, mat%irow, mat%cols, mat%val, mat%diag, &
               &      mat%perm, mat%invp, sol, mat%rank, mat%nrhs, &
               &      mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
               &      mat%p%dparm)
       ELSE
          CALL wgsmp(mat%rank, mat%irow, mat%cols, mat%val, sol, mat%rank, &
               &     mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
       END IF
    ELSE
       IF(mat%nlsym) THEN
          CALL wssmp(mat%rank, mat%irow, mat%cols, mat%val, mat%diag, &
               &      mat%perm, mat%invp, rhs, mat%rank, mat%nrhs, &
               &      mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
               &      mat%p%dparm)
       ELSE
          CALL wgsmp(mat%rank, mat%irow, mat%cols, mat%val, rhs, mat%rank, &
               &     mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
       END IF
    END IF
    IF(mat%p%iparm(64).NE.0) THEN
       WRITE(*,'(a,i4)') 'BSOLVE: Failed with error', mat%p%iparm(64)
       STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
    END IF
!
   END SUBROUTINE bsolve_wsmp_mat1
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE bsolve_zwsmp_mat1(mat, rhs, sol, nref)
!
!   Backsolve, using Wsmp
!
    TYPE(zwsmp_mat)          :: mat
    DOUBLE COMPLEX           :: rhs(:)
    DOUBLE COMPLEX, OPTIONAL :: sol(:)
    INTEGER, OPTIONAL        :: nref
    DOUBLE COMPLEX :: dummy
!
!    Recall the matrice if not current
!
    CALL check_mat(mat)
!
    IF(mat%nlsym .OR. mat%nlherm) THEN
       mat%p%iparm(2) = 4  ! Back substitution
       mat%p%iparm(3) = 4
    ELSE
       mat%p%iparm(2) = 3  ! Back substitution
       mat%p%iparm(3) = 3
    END IF
    mat%p%iparm(6) = 0  ! Max numbers of iterative refinement steps
    IF(PRESENT(nref)) THEN
       IF(mat%nlsym .OR. mat%nlherm) THEN
          mat%p%iparm(3) = 5
       ELSE
          mat%p%iparm(3) = 4
       END IF
       mat%p%iparm(6) = nref
    END IF
    mat%nrhs = 1
    IF(PRESENT(sol)) THEN
       sol = rhs
       IF(mat%nlsym .OR. mat%nlherm) THEN
          CALL zssmp(mat%rank, mat%irow, mat%cols, mat%val, mat%diag, &
               &      mat%perm, mat%invp, sol, mat%rank, mat%nrhs, &
               &      mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
               &      mat%p%dparm)
       ELSE
          CALL zgsmp(mat%rank, mat%irow, mat%cols, mat%val, sol, mat%rank, &
               &     mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
       END IF
    ELSE
       IF(mat%nlsym .OR. mat%nlherm) THEN
          CALL zssmp(mat%rank, mat%irow, mat%cols, mat%val, mat%diag, &
               &      mat%perm, mat%invp, rhs, mat%rank, mat%nrhs, &
               &      mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
               &      mat%p%dparm)
       ELSE
          CALL zgsmp(mat%rank, mat%irow, mat%cols, mat%val, rhs, mat%rank, &
               &     mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
       END IF
    END IF
    IF(mat%p%iparm(64).NE.0) THEN
       WRITE(*,'(a,i4)') 'BSOLVE: Failed with error', mat%p%iparm(64)
       STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
    END IF
!
  END SUBROUTINE bsolve_zwsmp_mat1
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE bsolve_wsmp_matn(mat, rhs, sol, nref)
!
!   Backsolve, using Wsmp, multiple RHS
!
    TYPE(wsmp_mat)             :: mat
    DOUBLE PRECISION           :: rhs(:,:)
    DOUBLE PRECISION, OPTIONAL :: sol(:,:)
    INTEGER, OPTIONAL          :: nref
    DOUBLE PRECISION :: dummy
!
!    Recall the matrice if not current
!
    CALL check_mat(mat)
!
    IF(mat%nlsym) THEN
       mat%p%iparm(2) = 4  ! Back substitution
       mat%p%iparm(3) = 4
    ELSE
       mat%p%iparm(2) = 3  ! Back substitution
       mat%p%iparm(3) = 3
    END IF
    mat%p%iparm(6) = 0  ! Max numbers of iterative refinement steps
    IF(PRESENT(nref)) THEN
       IF(mat%nlsym) THEN
          mat%p%iparm(3) = 5
       ELSE
          mat%p%iparm(3) = 4
       END IF
       mat%p%iparm(6) = nref
    END IF
    mat%nrhs = SIZE(rhs,2)
    IF(PRESENT(sol)) THEN
       sol = rhs
       IF(mat%nlsym) THEN
          CALL wssmp(mat%rank, mat%irow, mat%cols, mat%val, mat%diag, &
               &      mat%perm, mat%invp, sol, mat%rank, mat%nrhs, &
               &      mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
               &      mat%p%dparm)
       ELSE
          CALL wgsmp(mat%rank, mat%irow, mat%cols, mat%val, sol, mat%rank, &
               &     mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
       END IF
    ELSE
       IF(mat%nlsym) THEN
          CALL wssmp(mat%rank, mat%irow, mat%cols, mat%val, mat%diag, &
               &      mat%perm, mat%invp, rhs, mat%rank, mat%nrhs, &
               &      mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
               &      mat%p%dparm)
       ELSE
          CALL wgsmp(mat%rank, mat%irow, mat%cols, mat%val, rhs, mat%rank, &
               &     mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
       END IF
    END IF
    IF(mat%p%iparm(64).NE.0) THEN
       WRITE(*,'(a,i4)') 'BSOLVE: Failed with error', mat%p%iparm(64)
       STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
    END IF
!
  END SUBROUTINE bsolve_wsmp_matn
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE bsolve_zwsmp_matn(mat, rhs, sol, nref)
!
!   Backsolve, using Wsmp, multiple RHS
!
    TYPE(zwsmp_mat)          :: mat
    DOUBLE COMPLEX           :: rhs(:,:)
    DOUBLE COMPLEX, OPTIONAL :: sol(:,:)
    INTEGER, OPTIONAL        :: nref
    DOUBLE COMPLEX :: dummy
!
!    Recall the matrice if not current
!
    CALL check_mat(mat)
!
    IF(mat%nlsym .or. mat%nlherm) THEN
       mat%p%iparm(2) = 4  ! Back substitution
       mat%p%iparm(3) = 4
    ELSE
       mat%p%iparm(2) = 3  ! Back substitution
       mat%p%iparm(3) = 3
    END IF
    mat%p%iparm(6) = 0  ! Max numbers of iterative refinement steps
    IF(PRESENT(nref)) THEN
       IF(mat%nlsym .OR. mat%nlherm) THEN
          mat%p%iparm(3) = 5
       ELSE
          mat%p%iparm(3) = 4
       END IF
       mat%p%iparm(6) = nref
    END IF
    mat%nrhs = SIZE(rhs,2)
    IF(PRESENT(sol)) THEN
       sol = rhs
       IF(mat%nlsym .OR. mat%nlherm) THEN
          CALL zssmp(mat%rank, mat%irow, mat%cols, mat%val, mat%diag, &
               &      mat%perm, mat%invp, sol, mat%rank, mat%nrhs, &
               &      mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
               &      mat%p%dparm)
       ELSE
          CALL zgsmp(mat%rank, mat%irow, mat%cols, mat%val, sol, mat%rank, &
               &     mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
       END IF
    ELSE
       IF(mat%nlsym .OR. mat%nlherm) THEN
          CALL zssmp(mat%rank, mat%irow, mat%cols, mat%val, mat%diag, &
               &      mat%perm, mat%invp, rhs, mat%rank, mat%nrhs, &
               &      mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
               &      mat%p%dparm)
       ELSE
          CALL zgsmp(mat%rank, mat%irow, mat%cols, mat%val, rhs, mat%rank, &
               &     mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
       END IF
    END IF
    IF(mat%p%iparm(64).NE.0) THEN
       WRITE(*,'(a,i4)') 'BSOLVE: Failed with error', mat%p%iparm(64)
       STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
    END IF
!
  END SUBROUTINE bsolve_zwsmp_matn
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  FUNCTION vmx_wsmp_mat(mat, xarr) RESULT(yarr)
!
!   Return product mat*x
!
    TYPE(wsmp_mat)                :: mat
    DOUBLE PRECISION, INTENT(in)  :: xarr(:)
    DOUBLE PRECISION              :: yarr(SIZE(xarr))
    DOUBLE PRECISION :: alpha=1.0d0, beta=0.0d0
    CHARACTER(len=6) :: matdescra
    INTEGER :: n, i, j
!
    n = mat%rank
!
#ifdef MKL
    IF(mat%nlsym) THEN
       matdescra = 'sun'
    ELSE
       matdescra = 'g'
    END IF
!
    CALL mkl_dcsrmv('N', n, n, alpha, matdescra, mat%val, &
         &          mat%cols, mat%irow(1), mat%irow(2), xarr, &
         &          beta, yarr)
#else
    yarr = 0.0d0
    DO i=1,n
       DO j=mat%irow(i), mat%irow(i+1)-1
          yarr(i) = yarr(i) + mat%val(j)*xarr(mat%cols(j))
       END DO
       IF(mat%nlsym) THEN
          DO j=mat%irow(i)+1, mat%irow(i+1)-1
             yarr(mat%cols(j)) = yarr(mat%cols(j)) &
                  &              + mat%val(j)*xarr(i)
          END DO
       END IF
    END DO
#endif
!
  END FUNCTION vmx_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  FUNCTION vmx_zwsmp_mat(mat, xarr) RESULT(yarr)
!
!   Return product mat*x
!
    TYPE(zwsmp_mat)             :: mat
    DOUBLE COMPLEX, INTENT(in)  :: xarr(:)
    DOUBLE COMPLEX              :: yarr(SIZE(xarr))
    DOUBLE COMPLEX :: alpha=(1.0d0,0.0d0), beta=(0.0d0,0.0d0)
    INTEGER :: n, i, j
    CHARACTER(len=6) :: matdescra
    CHARACTER(len=1) :: transa
!
    n = mat%rank
!
#ifdef MKL
    IF(mat%nlsym) THEN
       matdescra = 'sun'
    ELSE IF(mat%nlherm) THEN
       matdescra = 'hun'
    ELSE
       matdescra = 'g'
    END IF
    transa='N'
    IF(mat%nlherm) THEN
       transa='T'   ! Transpose(CSR-UT*) = CSC-LT* = CSR-UT
    END IF
    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,0.0d0)
    DO i=1,n
       IF(mat%nlherm) THEN ! CSR-UT = (CSR-UT*)*
          DO j=mat%irow(i), mat%irow(i+1)-1
             yarr(i) = yarr(i) + CONJG(mat%val(j))*xarr(mat%cols(j))
          END DO
       ELSE
          DO j=mat%irow(i), mat%irow(i+1)-1
             yarr(i) = yarr(i) + mat%val(j)*xarr(mat%cols(j))
          END DO
       END IF
       IF(mat%nlsym) THEN
          DO j=mat%irow(i)+1, mat%irow(i+1)-1
             yarr(mat%cols(j)) = yarr(mat%cols(j)) &
                  &              + mat%val(j)*xarr(i)
          END DO
       ELSE IF(mat%nlherm) THEN ! CSR-UT = (CSR-UT*)*
          DO j=mat%irow(i)+1, mat%irow(i+1)-1
             yarr(mat%cols(j)) = yarr(mat%cols(j)) &
                  &              + mat%val(j)*xarr(i)
          END DO
       END IF
    END DO
#endif
!
  END FUNCTION vmx_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  FUNCTION vmx_wsmp_matn(mat, xarr) RESULT(yarr)
!
!   Return product mat*x
!
    TYPE(wsmp_mat)                :: mat
    DOUBLE PRECISION, INTENT(in)  :: xarr(:,:)
    DOUBLE PRECISION              :: yarr(SIZE(xarr,1),SIZE(xarr,2))
!
    DOUBLE PRECISION :: alpha=1.0d0, beta=0.0d0
    INTEGER :: n, nrhs, i, j
    CHARACTER(len=6) :: matdescra
!
    n = mat%rank
    nrhs = SIZE(xarr,2)
!
#ifdef MKL
    IF(mat%nlsym) THEN
       matdescra = 'sun'
    ELSE
       matdescra = 'g'
    END IF
!
    CALL mkl_dcsrmm('N', 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
          yarr(i,:) = yarr(i,:) + mat%val(j)*xarr(mat%cols(j),:)
       END DO
       IF(mat%nlsym) THEN
          DO j=mat%irow(i)+1, mat%irow(i+1)-1
             yarr(mat%cols(j),:) = yarr(mat%cols(j),:) &
                  &                + mat%val(j)*xarr(i,:)
          END DO
       END IF
    END DO
#endif
!
  END FUNCTION vmx_wsmp_matn
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  FUNCTION vmx_zwsmp_matn(mat, xarr) RESULT(yarr)
!
!   Return product mat*x
!
    TYPE(zwsmp_mat)             :: mat
    DOUBLE COMPLEX, INTENT(in)  :: xarr(:,:)
    DOUBLE COMPLEX              :: yarr(SIZE(xarr,1),SIZE(xarr,2))
!
    DOUBLE COMPLEX :: alpha=(1.0d0,0.0d0), beta=(0.0d0,0.0d0)
    INTEGER :: n, nrhs, i, j
    CHARACTER(len=6) :: matdescra
    CHARACTER(len=1) :: transa
!
    n = mat%rank
    nrhs = SIZE(xarr,2)
!
#ifdef MKL
    IF(mat%nlsym) THEN
       matdescra = 'sun'
    ELSE IF(mat%nlherm) THEN
       matdescra = 'hun'
    ELSE
       matdescra = 'g'
    END IF
    transa='N'
    IF(mat%nlherm) THEN
       transa='T' ! Transpose(CSR-UT*) = CSC-LT* = CSR-UT
    END IF
!
    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
       IF(mat%nlherm) THEN ! CSR-UT = (CSR-UT*)*
          DO j=mat%irow(i), mat%irow(i+1)-1
             yarr(i,:) = yarr(i,:) + CONJG(mat%val(j))*xarr(mat%cols(j),:)
          END DO
       ELSE
          DO j=mat%irow(i), mat%irow(i+1)-1
             yarr(i,:) = yarr(i,:) + mat%val(j)*xarr(mat%cols(j),:)
          END DO
       END IF
       IF(mat%nlsym) THEN
          DO j=mat%irow(i)+1, mat%irow(i+1)-1
             yarr(mat%cols(j),:) = yarr(mat%cols(j),:) &
                  &                + mat%val(j)*xarr(i,:)
          END DO
       ELSE IF(mat%nlherm) THEN ! CSR-UT = (CSR-UT*)*
          DO j=mat%irow(i)+1, mat%irow(i+1)-1
             yarr(mat%cols(j),:) = yarr(mat%cols(j),:) &
                  &               + mat%val(j)*xarr(i,:)
          END DO
       END IF
    END DO
#endif
!
  END FUNCTION vmx_zwsmp_matn
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE destroy_wsmp_mat(mat)
!
!  Deallocate the sparse matrix mat
!
    TYPE(wsmp_mat) :: mat
!
    IF(ASSOCIATED(mat%mat)) THEN
       CALL destroy(mat%mat)
       DEALLOCATE(mat%mat)
    END IF
!
!  Release memory for factors for symmetric matrix
    IF(mat%nlsym) THEN
       CALL check_mat(mat)
       CALL wsffree
    END IF
!
    IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
    IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
    IF(ASSOCIATED(mat%perm)) DEALLOCATE(mat%perm)
    IF(ASSOCIATED(mat%invp)) DEALLOCATE(mat%invp)
    IF(ASSOCIATED(mat%mrp)) DEALLOCATE(mat%mrp)
    IF(ASSOCIATED(mat%diag)) DEALLOCATE(mat%diag)
    IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
    IF(ASSOCIATED(mat%aux)) DEALLOCATE(mat%aux)
  END SUBROUTINE destroy_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE destroy_zwsmp_mat(mat)
!
!  Deallocate the sparse matrix mat
!
    TYPE(zwsmp_mat) :: mat
!
    IF(ASSOCIATED(mat%mat)) THEN
       CALL destroy(mat%mat)
       DEALLOCATE(mat%mat)
    END IF
!
!  Release memory for factors for symmetric/hermitian matrix
    IF(mat%nlsym .OR. mat%nlherm) THEN
       CALL check_mat(mat)
       CALL wsffree
    END IF
!
    IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
    IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
    IF(ASSOCIATED(mat%perm)) DEALLOCATE(mat%perm)
    IF(ASSOCIATED(mat%invp)) DEALLOCATE(mat%invp)
    IF(ASSOCIATED(mat%mrp)) DEALLOCATE(mat%mrp)
    IF(ASSOCIATED(mat%diag)) DEALLOCATE(mat%diag)
    IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
    IF(ASSOCIATED(mat%aux)) DEALLOCATE(mat%aux)
  END SUBROUTINE destroy_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE put_wsmp_mat(fid, label, mat, str)
!
!   Write matrix to hdf5 file
!
    USE futils
!
    INTEGER, INTENT(in)                    :: fid
    CHARACTER(len=*), INTENT(in)           :: label
    TYPE(wsmp_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 attach(fid, label, 'NLSYM', mat%nlsym)
    CALL attach(fid, label, 'NLPOS', mat%nlpos)
    CALL putarr(fid, TRIM(label)//'/irow', mat%irow)
    CALL putarr(fid, TRIM(label)//'/cols', mat%cols)
    CALL putarr(fid, TRIM(label)//'/val', mat%val)
!
    CALL creatg(fid, TRIM(label)//'/p')
    CALL putarr(fid, TRIM(label)//'/p/iparm', mat%p%iparm)
    CALL putarr(fid, TRIM(label)//'/p/dparm', mat%p%dparm)
  END SUBROUTINE put_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE put_zwsmp_mat(fid, label, mat, str)
!
!   Write matrix to hdf5 file
!
    USE futils
!
    INTEGER, INTENT(in)                    :: fid
    CHARACTER(len=*), INTENT(in)           :: label
    TYPE(zwsmp_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 attach(fid, label, 'NLSYM', mat%nlsym)
    CALL attach(fid, label, 'NLPOS', mat%nlpos)
    CALL attach(fid, label, 'NLHERM', mat%nlherm)
    CALL putarr(fid, TRIM(label)//'/irow', mat%irow)
    CALL putarr(fid, TRIM(label)//'/cols', mat%cols)
    CALL putarr(fid, TRIM(label)//'/val', mat%val)
!
    CALL creatg(fid, TRIM(label)//'/p')
    CALL putarr(fid, TRIM(label)//'/p/iparm', mat%p%iparm)
    CALL putarr(fid, TRIM(label)//'/p/dparm', mat%p%dparm)
  END SUBROUTINE put_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE get_wsmp_mat(fid, label, mat)
!
!   Read matrix from hdf5 file
!
    USE futils
!
    INTEGER, INTENT(in)          :: fid
    CHARACTER(len=*), INTENT(in) :: label
    TYPE(wsmp_mat)               :: mat
!
    CALL getatt(fid, label, 'RANK', mat%rank)
    CALL getatt(fid, label, 'NNZ',  mat%nnz)
    CALL getatt(fid, label, 'NLSYM', mat%nlsym)
    CALL getatt(fid, label, 'NLPOS', mat%nlpos)
    CALL getarr(fid, TRIM(label)//'/irow', mat%irow)
    CALL getarr(fid, TRIM(label)//'/cols', mat%cols)
    IF(mat%nlsym) THEN
       CALL getarr(fid, TRIM(label)//'/perm', mat%perm)
       CALL getarr(fid, TRIM(label)//'/invp', mat%invp)
    END IF
    CALL getarr(fid, TRIM(label)//'/val', mat%val)
!
    CALL getarr(fid, TRIM(label)//'/p/iparm', mat%p%iparm)
    CALL getarr(fid, TRIM(label)//'/p/dparm', mat%p%dparm)
  END SUBROUTINE get_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE get_zwsmp_mat(fid, label, mat)
!
!   Read matrix from hdf5 file
!
    USE futils
!
    INTEGER, INTENT(in)          :: fid
    CHARACTER(len=*), INTENT(in) :: label
    TYPE(zwsmp_mat)              :: mat
!
    CALL getatt(fid, label, 'RANK', mat%rank)
    CALL getatt(fid, label, 'NNZ',  mat%nnz)
    CALL getatt(fid, label, 'NLSYM', mat%nlsym)
    CALL getatt(fid, label, 'NLPOS', mat%nlpos)
    CALL getatt(fid, label, 'NLHERM', mat%nlherm)
    CALL getarr(fid, TRIM(label)//'/irow', mat%irow)
    CALL getarr(fid, TRIM(label)//'/cols', mat%cols)
    IF(mat%nlsym) THEN
       CALL getarr(fid, TRIM(label)//'/perm', mat%perm)
       CALL getarr(fid, TRIM(label)//'/invp', mat%invp)
    END IF
    CALL getarr(fid, TRIM(label)//'/val', mat%val)
!
    CALL getarr(fid, TRIM(label)//'/p/iparm', mat%p%iparm)
    CALL getarr(fid, TRIM(label)//'/p/dparm', mat%p%dparm)
  END SUBROUTINE get_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE mcopy_wsmp_mat(mata, matb)
!
!   Matrix copy: B = A
!
    TYPE(wsmp_mat) :: mata, matb
    INTEGER :: n, nnz
!
! Assume that matb was already initialized by init_wsmp_mat.
    IF(matb%rank.LE.0) THEN
       WRITE(*,'(a)') 'MCOPY: Mat B is not initialized with INIT'
       STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
    END IF
!
    IF(ASSOCIATED(matb%mat)) THEN 
       CALL destroy(matb%mat)
       DEALLOCATE(matb%mat)
    END IF
!
    n = mata%rank
    nnz = mata%nnz
    matb%rank = n
    matb%nnz = nnz
    matb%nlsym = mata%nlsym
    matb%nlpos = mata%nlpos
    matb%nlforce_zero = mata%nlforce_zero
!
    IF(ASSOCIATED(matb%val)) DEALLOCATE(matb%val)
    IF(ASSOCIATED(matb%cols)) DEALLOCATE(matb%cols)
    IF(ASSOCIATED(matb%irow)) DEALLOCATE(matb%irow)
    IF(ASSOCIATED(matb%perm)) DEALLOCATE(matb%perm)
    IF(ASSOCIATED(matb%invp)) DEALLOCATE(matb%invp)
    ALLOCATE(matb%val(nnz)); matb%val = mata%val
    ALLOCATE(matb%cols(nnz)); matb%cols = mata%cols
    ALLOCATE(matb%irow(n+1)); matb%irow = mata%irow
    ALLOCATE(matb%perm(n))
    IF(matb%nlsym) THEN
       ALLOCATE(matb%perm(n))
       ALLOCATE(matb%invp(n))
    END IF
  END SUBROUTINE mcopy_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE mcopy_zwsmp_mat(mata, matb)
!
!   Matrix copy: B = A
!
    TYPE(zwsmp_mat) :: mata, matb
    INTEGER :: n, nnz
!
! Assume that matb was already initialized by init_wsmp_mat.
    IF(matb%rank.LE.0) THEN
       WRITE(*,'(a)') 'MCOPY: Mat B is not initialized with INIT'
       STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
    END IF
!
    IF(ASSOCIATED(matb%mat)) THEN 
       CALL destroy(matb%mat)
       DEALLOCATE(matb%mat)
    END IF
!
    n = mata%rank
    nnz = mata%nnz
    matb%rank = n
    matb%nnz = nnz
    matb%nlsym = mata%nlsym
    matb%nlherm = mata%nlherm
    matb%nlpos = mata%nlpos
    matb%nlforce_zero = mata%nlforce_zero
!
    IF(ASSOCIATED(matb%val)) DEALLOCATE(matb%val)
    IF(ASSOCIATED(matb%cols)) DEALLOCATE(matb%cols)
    IF(ASSOCIATED(matb%irow)) DEALLOCATE(matb%irow)
    IF(ASSOCIATED(matb%perm)) DEALLOCATE(matb%perm)
    IF(ASSOCIATED(matb%invp)) DEALLOCATE(matb%invp)
    ALLOCATE(matb%val(nnz)); matb%val = mata%val
    ALLOCATE(matb%cols(nnz)); matb%cols = mata%cols
    ALLOCATE(matb%irow(n+1)); matb%irow = mata%irow
    ALLOCATE(matb%perm(n))
    IF(matb%nlsym) THEN
       ALLOCATE(matb%perm(n))
       ALLOCATE(matb%invp(n))
    END IF
  END SUBROUTINE mcopy_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE maddto_wsmp_mat(mata, alpha, matb)
!
!   A <- A + alpha*B
!
    TYPE(wsmp_mat)   :: mata, matb
    DOUBLE PRECISION :: alpha
!
    mata%val = mata%val + alpha*matb%val
  END SUBROUTINE maddto_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE maddto_zwsmp_mat(mata, alpha, matb)
!
!   A <- A + alpha*B
!
    TYPE(zwsmp_mat) :: mata, matb
    DOUBLE COMPLEX  :: alpha
!
    mata%val = mata%val + alpha*matb%val
  END SUBROUTINE maddto_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE psum_wsmp_mat(mat, comm)
!
!   Parallel sum of sparse matrices
!
    INCLUDE "mpif.h"
!
    TYPE(wsmp_mat)   :: mat
    INCLUDE 'psum_mat.tpl'
  END SUBROUTINE psum_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE psum_zwsmp_mat(mat, comm)
!
!   Parallel sum of sparse matrices
!
    INCLUDE "mpif.h"
!
    TYPE(zwsmp_mat)   :: mat
    INCLUDE 'psum_mat.tpl'
  END SUBROUTINE psum_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE p2p_wsmp_mat(mat, dest, extyp, op, comm)
!
!   Point-to-point combine sparse matrix between 2 processes
!
    INCLUDE "mpif.h"
!
    TYPE(wsmp_mat)            :: mat
    DOUBLE PRECISION, ALLOCATABLE :: val(:)
    INTEGER :: mpi_type=MPI_DOUBLE_PRECISION
!
    INCLUDE "p2p_mat.tpl"
  END SUBROUTINE p2p_wsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  SUBROUTINE p2p_zwsmp_mat(mat, dest, extyp, op, comm)
!
!   Point-to-point combine sparse matrix between 2 processes
!
    INCLUDE "mpif.h"
!
    TYPE(zwsmp_mat)           :: mat
    DOUBLE COMPLEX, ALLOCATABLE :: val(:)
    INTEGER :: mpi_type=MPI_DOUBLE_COMPLEX
!
    INCLUDE "p2p_mat.tpl"
  END SUBROUTINE p2p_zwsmp_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
END MODULE wsmp_bsplines