!>
!> @file pwsmp_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 .
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran
!>
MODULE pwsmp_bsplines
!
! PWSMP_BSPLINES: Simple interface to the parallel sparse direct
! solver PWSMP.
!
! T.M. Tran, CRPP-EPFL
! December 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
INTEGER :: comm
INTEGER :: istart, iend, rank_loc
INTEGER :: nnz_start, nnz_end, nnz_loc
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
INTEGER :: comm
INTEGER :: istart, iend, rank_loc
INTEGER :: nnz_start, nnz_end, nnz_loc
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, comm_in)
!
! Initialize an empty sparse wsmp matrix
!
USE pputils2
INCLUDE 'mpif.h'
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, OPTIONAL, INTENT(in) :: comm_in
!
INTEGER :: comm, nloc
INTEGER :: info
INTEGER :: idummy = 0
DOUBLE PRECISION :: dummy = 0.0d0
!
comm = MPI_COMM_WORLD
IF(PRESENT(comm_in)) comm = comm_in
mat%comm = comm
!
! Store away (valid) current matrix id
!
IF(current_matid .GE. 0) THEN
CALL wstoremat(current_matid, info)
IF(info.NE.0) THEN
WRITE(*,'(a,i12)') '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
!
! Matrix partition
!
CALL dist1d(comm, 1, n, mat%istart, nloc)
mat%iend = mat%istart + nloc - 1
mat%rank_loc = nloc
!
IF(ASSOCIATED(mat%mat)) DEALLOCATE(mat%mat)
ALLOCATE(mat%mat)
CALL init(n, mat%mat, mat%istart, mat%iend)
!
! Fill 'iparm' and 'dparm' with default values
!
mat%p%iparm(1:3) = 0
IF(mat%nlsym) THEN
CALL pwssmp(mat%rank_loc, mat%irow, mat%cols, mat%val, mat%diag, &
& mat%perm, mat%invp, dummy, mat%rank_loc, 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
CALL pwgsmp(mat%rank_loc, 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,i12)') 'INIT: Initialization failed with error', &
& mat%p%iparm(64)
STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
END IF
!
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, comm_in)
!
! Initialize an empty sparse wsmp matrix
!
USE pputils2
INCLUDE 'mpif.h'
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, OPTIONAL, INTENT(in) :: comm_in
!
INTEGER :: comm, nloc
INTEGER :: info
INTEGER :: idummy = 0
DOUBLE COMPLEX :: dummy = 0.0d0
!
comm = MPI_COMM_WORLD
IF(PRESENT(comm_in)) comm = comm_in
mat%comm = comm
!
! Store away (valid) current matrix id
!
IF(current_matid .GE. 0) THEN
CALL wstoremat(current_matid, info)
IF(info.NE.0) THEN
WRITE(*,'(a,i12)') '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
!
! Matrix partition
!
CALL dist1d(comm, 1, n, mat%istart, nloc)
mat%iend = mat%istart + nloc - 1
mat%rank_loc = nloc
!
IF(ASSOCIATED(mat%mat)) DEALLOCATE(mat%mat)
ALLOCATE(mat%mat)
CALL init(n, mat%mat, mat%istart, mat%iend)
!
! Fill 'iparm' and 'dparm' with default values
!
mat%p%iparm(1:3) = 0
IF(mat%nlherm .OR. mat%nlsym) THEN
CALL pzssmp(mat%rank_loc, mat%irow, mat%cols, mat%val, mat%diag, &
& mat%perm, mat%invp, dummy, mat%rank_loc, 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 pzgsmp(mat%rank_loc, 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,i12)') '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,i12)') '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,i12)') '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,i12)') '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,i12)') '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(i.LT.mat%istart .OR. i.GT.mat%iend) THEN
WRITE(*,'(a,2i6)') 'UPDT: i, j out of range ', i, j
WRITE(*,'(a,2i6)') ' istart, iend ', mat%istart, mat%iend
STOP '*** Abnormal EXIT in MODULE mumps_mod ***'
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
!
! Do nothing if outside
IF(iput.LT.mat%istart .OR. iput.GT.mat%iend) RETURN
!
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
!
! Do nothing if outside
IF(iput.LT.mat%istart .OR. iput.GT.mat%iend) RETURN
!
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
!
val = 0.0d0 ! Assume zero val if outside
IF( iget.LT.mat%istart .OR. iget.GT.mat%iend ) RETURN
!
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
!
val = (0.0d0, 0.0d0) ! Assume zero val if outside
IF( iget.LT.mat%istart .OR. iget.GT.mat%iend ) RETURN
!
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=amat%istart,amat%iend
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=amat%istart,amat%iend
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=amat%istart,amat%iend
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=amat%istart,amat%iend
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=mat%istart,mat%iend
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=mat%istart,mat%iend
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
!
INCLUDE 'mpif.h'
TYPE(wsmp_mat) :: mat
LOGICAL, INTENT(in), OPTIONAL :: nlkeep
INTEGER :: i, nnz, rank, s, e
INTEGER :: comm, ierr, nnz_loc, rank_loc
LOGICAL :: nlclean
!
nlclean = .TRUE.
IF(PRESENT(nlkeep)) THEN
nlclean = .NOT. nlkeep
END IF
!
comm = mat%comm
!
! Allocate the WSMP matrix structure
!
rank = mat%rank
rank_loc = mat%rank
!
nnz_loc = get_count(mat)
mat%nnz_start = 0
CALL mpi_exscan(nnz_loc, mat%nnz_start, 1, MPI_INTEGER, MPI_SUM, comm, ierr)
mat%nnz_start = mat%nnz_start + 1
mat%nnz_end = mat%nnz_start + nnz_loc - 1
mat%nnz_loc = nnz_loc
CALL mpi_allreduce(nnz_loc, mat%nnz, 1, MPI_INTEGER, MPI_SUM, comm, ierr)
!
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 LOCAL irow, cols and val
IF(mat%nlsym) THEN
ALLOCATE(mat%perm(rank))
ALLOCATE(mat%invp(rank))
END IF
ALLOCATE(mat%val(nnz_loc))
ALLOCATE(mat%cols(nnz_loc))
ALLOCATE(mat%irow(mat%istart:mat%iend+1))
!
! Fill WSMP structure and deallocate the sparse rows
!
mat%irow(mat%istart) = 1
DO i=mat%istart,mat%iend
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
!
INCLUDE 'mpif.h'
TYPE(zwsmp_mat) :: mat
LOGICAL, INTENT(in), OPTIONAL :: nlkeep
INTEGER :: i, nnz, rank, s, e
INTEGER :: comm, ierr, nnz_loc, rank_loc
LOGICAL :: nlclean
!
nlclean = .TRUE.
IF(PRESENT(nlkeep)) THEN
nlclean = .NOT. nlkeep
END IF
!
comm = mat%comm
!
! Allocate the WSMP matrix structure
!
rank = mat%rank
rank_loc = mat%rank
!
nnz_loc = get_count(mat)
mat%nnz_start = 0
CALL mpi_exscan(nnz_loc, mat%nnz_start, 1, MPI_INTEGER, MPI_SUM, comm, ierr)
mat%nnz_start = mat%nnz_start + 1
mat%nnz_end = mat%nnz_start + nnz_loc - 1
mat%nnz_loc = nnz_loc
CALL mpi_allreduce(nnz_loc, mat%nnz, 1, MPI_INTEGER, MPI_SUM, comm, ierr)
!
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 LOCAL irow, cols and val
IF(mat%nlsym) THEN
ALLOCATE(mat%perm(rank))
ALLOCATE(mat%invp(rank))
END IF
ALLOCATE(mat%val(nnz_loc))
ALLOCATE(mat%cols(nnz_loc))
ALLOCATE(mat%irow(mat%istart:mat%iend+1))
!
! Fill WSMP structure and deallocate the sparse rows
!
mat%irow(mat%istart) = 1
DO i=mat%istart,mat%iend
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 pwssmp(mat%rank_loc, mat%irow, mat%cols, mat%val, mat%diag, &
& mat%perm, mat%invp, dummy, mat%rank_loc, 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 pwgsmp(mat%rank_loc, mat%irow, mat%cols, mat%val, dummy, mat%rank_loc, &
& mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
END IF
IF(mat%p%iparm(64).NE.0) THEN
WRITE(*,'(a,i12)') '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_loc, mat%irow, mat%cols, mat%val, mat%diag, &
& mat%perm, mat%invp, dummy, mat%rank_loc, 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_loc, mat%irow, mat%cols, mat%val, dummy, mat%rank_loc, &
& mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
END IF
IF(mat%p%iparm(64).NE.0) THEN
WRITE(*,'(a,i12)') '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 pwssmp(mat%rank_loc, mat%irow, mat%cols, mat%val, mat%diag, &
& mat%perm, mat%invp, dummy, mat%rank_loc, 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 pwgsmp(mat%rank_loc, mat%irow, mat%cols, mat%val, dummy, mat%rank_loc, &
& mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
END IF
IF(mat%p%iparm(64).NE.0) THEN
WRITE(*,'(a,i12)') '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_loc, mat%irow, mat%cols, mat%val, mat%diag, &
& mat%perm, mat%invp, dummy, mat%rank_loc, 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_loc, mat%irow, mat%cols, mat%val, dummy, mat%rank_loc, &
& mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
END IF
IF(mat%p%iparm(64).NE.0) THEN
WRITE(*,'(a,i12)') '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
!
INCLUDE 'mpif.h'
TYPE(wsmp_mat) :: mat
DOUBLE PRECISION :: rhs(:)
DOUBLE PRECISION, OPTIONAL :: sol(:)
INTEGER, OPTIONAL :: nref
!
DOUBLE PRECISION :: sol_loc(mat%rank_loc)
INTEGER :: nloc, me, nprocs, ierr, i
INTEGER, ALLOCATABLE :: nlocs(:), displs(:)
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
!
! Extract local rhs from global rhs
!
sol_loc = rhs(mat%istart:mat%iend)
!
IF(mat%nlsym) THEN
CALL pwssmp(mat%rank_loc, mat%irow, mat%cols, mat%val, mat%diag, &
& mat%perm, mat%invp, sol_loc, mat%rank_loc, mat%nrhs, &
& mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
& mat%p%dparm)
ELSE
CALL pwgsmp(mat%rank_loc, mat%irow, mat%cols, mat%val, sol_loc, &
& mat%rank_loc, mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
END IF
!
IF(mat%p%iparm(64).NE.0) THEN
WRITE(*,'(a,i12)') 'BSOLVE: Failed with error', mat%p%iparm(64)
STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
END IF
!
! Allgatherv local sol
!
CALL mpi_comm_rank(mat%comm, me, ierr)
CALL mpi_comm_size(mat%comm, nprocs, ierr)
!
ALLOCATE(displs(0:nprocs))
ALLOCATE(nlocs(0:nprocs-1))
!
nloc = mat%rank_loc
CALL mpi_allgather(nloc, 1, MPI_INTEGER, nlocs, 1, MPI_INTEGER, &
& mat%comm, ierr)
!
displs(0) = 0
DO i=0,nprocs-1
displs(i+1) = displs(i)+nlocs(i)
END DO
!
IF(PRESENT(sol)) THEN
CALL mpi_allgatherv(sol_loc, nloc, MPI_DOUBLE_PRECISION, &
& sol, nlocs, displs, MPI_DOUBLE_PRECISION, &
& mat%comm, ierr)
ELSE
CALL mpi_allgatherv(sol_loc, nloc, MPI_DOUBLE_PRECISION, &
& rhs, nlocs, displs, MPI_DOUBLE_PRECISION, &
& mat%comm, ierr)
END IF
!
DEALLOCATE(nlocs)
DEALLOCATE(displs)
!
END SUBROUTINE bsolve_wsmp_mat1
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE bsolve_zwsmp_mat1(mat, rhs, sol, nref)
!
! Backsolve, using Wsmp
!
INCLUDE 'mpif.h'
TYPE(zwsmp_mat) :: mat
DOUBLE COMPLEX :: rhs(:)
DOUBLE COMPLEX, OPTIONAL :: sol(:)
INTEGER, OPTIONAL :: nref
!
DOUBLE COMPLEX :: sol_loc(mat%rank_loc)
INTEGER :: nloc, me, nprocs, ierr, i
INTEGER, ALLOCATABLE :: nlocs(:), displs(:)
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
!
! Extract local rhs from global rhs
!
sol_loc = rhs(mat%istart:mat%iend)
!
IF(mat%nlsym .OR. mat%nlherm) THEN
CALL pzssmp(mat%rank_loc, mat%irow, mat%cols, mat%val, mat%diag, &
& mat%perm, mat%invp, sol_loc, mat%rank_loc, mat%nrhs, &
& mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
& mat%p%dparm)
ELSE
CALL pzgsmp(mat%rank_loc, mat%irow, mat%cols, mat%val, sol_loc, &
& mat%rank_loc, mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
END IF
!
IF(mat%p%iparm(64).NE.0) THEN
WRITE(*,'(a,i12)') 'BSOLVE: Failed with error', mat%p%iparm(64)
STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
END IF
!
! Allgatherv local sol
!
CALL mpi_comm_rank(mat%comm, me, ierr)
CALL mpi_comm_size(mat%comm, nprocs, ierr)
!
ALLOCATE(displs(0:nprocs))
ALLOCATE(nlocs(0:nprocs-1))
!
nloc = mat%rank_loc
CALL mpi_allgather(nloc, 1, MPI_INTEGER, nlocs, 1, MPI_INTEGER, &
& mat%comm, ierr)
!
displs(0) = 0
DO i=0,nprocs-1
displs(i+1) = displs(i)+nlocs(i)
END DO
!
IF(PRESENT(sol)) THEN
CALL mpi_allgatherv(sol_loc, nloc, MPI_DOUBLE_COMPLEX, &
& sol, nlocs, displs, MPI_DOUBLE_COMPLEX, &
& mat%comm, ierr)
ELSE
CALL mpi_allgatherv(sol_loc, nloc, MPI_DOUBLE_COMPLEX, &
& rhs, nlocs, displs, MPI_DOUBLE_COMPLEX, &
& mat%comm, ierr)
END IF
!
DEALLOCATE(nlocs)
DEALLOCATE(displs)
!
END SUBROUTINE bsolve_zwsmp_mat1
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE bsolve_wsmp_matn(mat, rhs, sol, nref)
!
! Backsolve, using Wsmp, multiple RHS
!
INCLUDE 'mpif.h'
TYPE(wsmp_mat) :: mat
DOUBLE PRECISION :: rhs(:,:)
DOUBLE PRECISION, OPTIONAL :: sol(:,:)
INTEGER, OPTIONAL :: nref
!
DOUBLE PRECISION :: sol_loc(mat%rank_loc,SIZE(rhs,2))
INTEGER :: nloc, me, nprocs, ierr, i
INTEGER, ALLOCATABLE :: nlocs(:), displs(:)
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)
!
! Extract local rhs from global rhs
!
sol_loc(:,:) = rhs(mat%istart:mat%iend,:)
!
IF(mat%nlsym) THEN
CALL pwssmp(mat%rank_loc, mat%irow, mat%cols, mat%val, mat%diag, &
& mat%perm, mat%invp, sol_loc, mat%rank_loc, mat%nrhs, &
& mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
& mat%p%dparm)
ELSE
CALL pwgsmp(mat%rank_loc, mat%irow, mat%cols, mat%val, sol_loc, &
& mat%rank_loc, mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
END IF
!
IF(mat%p%iparm(64).NE.0) THEN
WRITE(*,'(a,i12)') 'BSOLVE: Failed with error', mat%p%iparm(64)
STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
END IF
!
! Allgatherv local sol
!
CALL mpi_comm_rank(mat%comm, me, ierr)
CALL mpi_comm_size(mat%comm, nprocs, ierr)
!
ALLOCATE(displs(0:nprocs))
ALLOCATE(nlocs(0:nprocs-1))
!
nloc = mat%rank_loc
CALL mpi_allgather(nloc, 1, MPI_INTEGER, nlocs, 1, MPI_INTEGER, &
& mat%comm, ierr)
!
displs(0) = 0
DO i=0,nprocs-1
displs(i+1) = displs(i)+nlocs(i)
END DO
!
DO i=1,mat%nrhs
IF(PRESENT(sol)) THEN
CALL mpi_allgatherv(sol_loc(1,i), nloc, MPI_DOUBLE_PRECISION, &
& sol(1,i), nlocs, displs, MPI_DOUBLE_PRECISION, &
& mat%comm, ierr)
ELSE
CALL mpi_allgatherv(sol_loc(1,i), nloc, MPI_DOUBLE_PRECISION, &
& rhs(1,i), nlocs, displs, MPI_DOUBLE_PRECISION, &
& mat%comm, ierr)
END IF
END DO
!
DEALLOCATE(nlocs)
DEALLOCATE(displs)
!
END SUBROUTINE bsolve_wsmp_matn
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE bsolve_zwsmp_matn(mat, rhs, sol, nref)
!
! Backsolve, using Wsmp, multiple RHS
!
INCLUDE 'mpif.h'
TYPE(zwsmp_mat) :: mat
DOUBLE COMPLEX :: rhs(:,:)
DOUBLE COMPLEX, OPTIONAL :: sol(:,:)
INTEGER, OPTIONAL :: nref
!
DOUBLE COMPLEX :: sol_loc(mat%rank_loc,SIZE(rhs,2))
INTEGER :: nloc, me, nprocs, ierr, i
INTEGER, ALLOCATABLE :: nlocs(:), displs(:)
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)
!
! Extract local rhs from global rhs
!
sol_loc(:,:) = rhs(mat%istart:mat%iend,:)
!
IF(mat%nlsym) THEN
CALL pzssmp(mat%rank_loc, mat%irow, mat%cols, mat%val, mat%diag, &
& mat%perm, mat%invp, sol_loc, mat%rank_loc, mat%nrhs, &
& mat%aux, SIZE(mat%aux), mat%mrp, mat%p%iparm, &
& mat%p%dparm)
ELSE
CALL pzgsmp(mat%rank_loc, mat%irow, mat%cols, mat%val, sol_loc, &
& mat%rank_loc, mat%nrhs, dummy, mat%p%iparm, mat%p%dparm)
END IF
!
IF(mat%p%iparm(64).NE.0) THEN
WRITE(*,'(a,i12)') 'BSOLVE: Failed with error', mat%p%iparm(64)
STOP '*** Abnormal EXIT in MODULE wsmp_mod ***'
END IF
!
! Allgatherv local sol
!
CALL mpi_comm_rank(mat%comm, me, ierr)
CALL mpi_comm_size(mat%comm, nprocs, ierr)
!
ALLOCATE(displs(0:nprocs))
ALLOCATE(nlocs(0:nprocs-1))
!
nloc = mat%rank_loc
CALL mpi_allgather(nloc, 1, MPI_INTEGER, nlocs, 1, MPI_INTEGER, &
& mat%comm, ierr)
!
displs(0) = 0
DO i=0,nprocs-1
displs(i+1) = displs(i)+nlocs(i)
END DO
!
DO i=1,mat%nrhs
IF(PRESENT(sol)) THEN
CALL mpi_allgatherv(sol_loc(1,i), nloc, MPI_DOUBLE_COMPLEX, &
& sol(1,i), nlocs, displs, MPI_DOUBLE_COMPLEX, &
& mat%comm, ierr)
ELSE
CALL mpi_allgatherv(sol_loc(1,i), nloc, MPI_DOUBLE_COMPLEX, &
& rhs(1,i), nlocs, displs, MPI_DOUBLE_COMPLEX, &
& mat%comm, ierr)
END IF
END DO
!
DEALLOCATE(nlocs)
DEALLOCATE(displs)
!
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 pwsmp_bsplines