!>
!> @file mumps_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 .
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran
!>
MODULE mumps_bsplines
!
! MUMPS_BSPLINES: Simple interface to the sparse direct solver MUMPS
!
! T.M. Tran, CRPP-EPFL
! June 2011
!
USE sparse
IMPLICIT NONE
INCLUDE 'dmumps_struc.h'
INCLUDE 'zmumps_struc.h'
!
TYPE mumps_mat
INTEGER :: rank, nnz
INTEGER :: nterms, kmat
INTEGER :: istart, iend
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()
DOUBLE PRECISION, POINTER :: val(:) => NULL()
TYPE(dmumps_struc) :: mumps_par
END TYPE mumps_mat
!
TYPE zmumps_mat
INTEGER :: rank, nnz
INTEGER :: nterms, kmat
INTEGER :: istart, iend
INTEGER :: nnz_start, nnz_end, nnz_loc
LOGICAL :: nlsym
LOGICAL :: nlherm
LOGICAL :: nlpos
LOGICAL :: nlforce_zero
TYPE(zspmat), POINTER :: mat => NULL()
INTEGER, POINTER :: cols(:) => NULL()
INTEGER, POINTER :: irow(:) => NULL()
INTEGER, POINTER :: perm(:) => NULL()
DOUBLE COMPLEX, POINTER :: val(:) => NULL()
TYPE(zmumps_struc) :: mumps_par
END TYPE zmumps_mat
!
INTERFACE init
MODULE PROCEDURE init_mumps_mat, init_zmumps_mat
END INTERFACE init
!
INTERFACE clear_mat
MODULE PROCEDURE clear_mumps_mat, clear_zmumps_mat
END INTERFACE clear_mat
!
INTERFACE updtmat
MODULE PROCEDURE updt_mumps_mat, updt_zmumps_mat
END INTERFACE updtmat
!
INTERFACE putele
MODULE PROCEDURE putele_mumps_mat, putele_zmumps_mat
END INTERFACE putele
!
INTERFACE getele
MODULE PROCEDURE getele_mumps_mat, getele_zmumps_mat
END INTERFACE getele
!
INTERFACE putrow
MODULE PROCEDURE putrow_mumps_mat, putrow_zmumps_mat
END INTERFACE putrow
!
INTERFACE getrow
MODULE PROCEDURE getrow_mumps_mat, getrow_zmumps_mat
END INTERFACE getrow
!
INTERFACE putcol
MODULE PROCEDURE putcol_mumps_mat, putcol_zmumps_mat
END INTERFACE putcol
!
INTERFACE getcol
MODULE PROCEDURE getcol_mumps_mat, getcol_zmumps_mat
END INTERFACE getcol
!
INTERFACE get_count
MODULE PROCEDURE get_count_mumps_mat, get_count_zmumps_mat
END INTERFACE get_count
!
INTERFACE to_mat
MODULE PROCEDURE to_mumps_mat, to_zmumps_mat
END INTERFACE to_mat
!
INTERFACE reord_mat
MODULE PROCEDURE reord_mumps_mat, reord_zmumps_mat
END INTERFACE reord_mat
!
INTERFACE numfact
MODULE PROCEDURE numfact_mumps_mat, numfact_zmumps_mat
END INTERFACE numfact
!
INTERFACE factor
MODULE PROCEDURE factor_mumps_mat, factor_zmumps_mat
END INTERFACE factor
!
INTERFACE bsolve
MODULE PROCEDURE bsolve_mumps_mat1, bsolve_mumps_matn, &
& bsolve_zmumps_mat1, bsolve_zmumps_matn
END INTERFACE bsolve
!
INTERFACE vmx
MODULE PROCEDURE vmx_mumps_mat, vmx_mumps_matn, &
& vmx_zmumps_mat, vmx_zmumps_matn
END INTERFACE vmx
!
INTERFACE destroy
MODULE PROCEDURE destroy_mumps_mat, destroy_zmumps_mat
END INTERFACE destroy
!
INTERFACE putmat
MODULE PROCEDURE put_mumps_mat, put_zmumps_mat
END INTERFACE putmat
!
INTERFACE getmat
MODULE PROCEDURE get_mumps_mat, get_zmumps_mat
END INTERFACE getmat
!
INTERFACE mcopy
MODULE PROCEDURE mcopy_mumps_mat, mcopy_zmumps_mat
END INTERFACE mcopy
!
INTERFACE maddto
MODULE PROCEDURE maddto_mumps_mat, maddto_zmumps_mat
END INTERFACE maddto
!
INTERFACE psum_mat
MODULE PROCEDURE psum_mumps_mat, psum_zmumps_mat
END INTERFACE psum_mat
!
INTERFACE p2p_mat
MODULE PROCEDURE p2p_mumps_mat, p2p_zmumps_mat
END INTERFACE p2p_mat
!
CONTAINS
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE init_mumps_mat(n, nterms, mat, kmat, nlsym, nlpos, &
& nlforce_zero, comm_in)
!
! Initialize an empty sparse mumps matrix
!
USE pputils2
INCLUDE 'mpif.h'
INTEGER, INTENT(in) :: n, nterms
TYPE(mumps_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
!
comm = MPI_COMM_SELF ! Default is serial!
IF(PRESENT(comm_in)) comm = comm_in
!
mat%rank = n
mat%nterms = nterms
mat%nnz = 0
mat%nlsym = .FALSE.
mat%nlpos = .TRUE.
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
!
IF(ASSOCIATED(mat%mat)) DEALLOCATE(mat%mat)
ALLOCATE(mat%mat)
CALL init(n, mat%mat, mat%istart, mat%iend)
!
! Initialize a MUMPS instance
!
mat%mumps_par%N = n
mat%mumps_par%NZ = 0
mat%mumps_par%COMM = comm
mat%mumps_par%PAR = 1 ! Host involved in calculations
IF(mat%nlsym) THEN
IF(mat%nlpos) THEN
mat%mumps_par%SYM = 1 ! symmetric, positive definite
ELSE
mat%mumps_par%SYM = 2 ! symmetric, indefinite
END IF
ELSE
mat%mumps_par%SYM = 0 ! unsymmetric
END IF
!
mat%mumps_par%JOB = -1 ! Init phase
CALL dmumps(mat%mumps_par)
!
! Nullify MUMPS pointers for distributed matrix
!
NULLIFY(mat%mumps_par%A_loc)
NULLIFY(mat%mumps_par%IRN_loc)
NULLIFY(mat%mumps_par%JCN_loc)
NULLIFY(mat%mumps_par%RHS)
!
END SUBROUTINE init_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE init_zmumps_mat(n, nterms, mat, kmat, nlsym, nlherm, &
& nlpos, nlforce_zero, comm_in)
!
! Initialize an empty sparse mumps matrix
!
USE pputils2
INCLUDE 'mpif.h'
INTEGER, INTENT(in) :: n, nterms
TYPE(zmumps_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
!
comm = MPI_COMM_SELF ! Default is serial!
IF(PRESENT(comm_in)) comm = comm_in
!
mat%rank = n
mat%nterms = nterms
mat%nnz = 0
mat%nlsym = .FALSE.
mat%nlherm = .FALSE.
mat%nlpos = .TRUE.
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
!
IF(ASSOCIATED(mat%mat)) DEALLOCATE(mat%mat)
ALLOCATE(mat%mat)
CALL init(n, mat%mat, mat%istart, mat%iend)
!
! Initialize a MUMPS instance
!
mat%mumps_par%N = n
mat%mumps_par%NZ = 0
mat%mumps_par%COMM = comm
mat%mumps_par%PAR = 1 ! Host involved in calculations
mat%mumps_par%SYM = 0 ! General unsymmetric
IF(mat%nlsym) THEN
IF(mat%nlpos) THEN
mat%mumps_par%SYM = 1 ! symmetric, positive definite
ELSE
mat%mumps_par%SYM = 2 ! symmetric, indefinite
END IF
END IF
!
mat%mumps_par%JOB = -1 ! Init phase
CALL zmumps(mat%mumps_par)
!
! WARNING: SYM=1 is currently (version 4.10.0) is treated as SYM=2.
! The Hermitian case is not implemented yet!
!
! Nullify MUMPS pointers for distributed matrix
!
NULLIFY(mat%mumps_par%A_loc)
NULLIFY(mat%mumps_par%IRN_loc)
NULLIFY(mat%mumps_par%JCN_loc)
NULLIFY(mat%mumps_par%RHS)
!
END SUBROUTINE init_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE clear_mumps_mat(mat)
!
! Clear matrix, keeping its sparse structure unchanged
!
TYPE(mumps_mat) :: mat
!
mat%val = 0.0d0
END SUBROUTINE clear_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE clear_zmumps_mat(mat)
!
! Clear matrix, keeping its sparse structure unchanged
!
TYPE(zmumps_mat) :: mat
!
mat%val = (0.0d0, 0.0d0)
END SUBROUTINE clear_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE updt_mumps_mat(mat, i, j, val)
!
! Update element Aij of mumps matrix
!
TYPE(mumps_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) - mat%nnz_start + 1
e = mat%irow(i+1) - mat%nnz_start
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 mumps_mod ***'
END IF
END IF
END SUBROUTINE updt_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE updt_zmumps_mat(mat, i, j, val)
!
! Update element Aij of mumps matrix
!
TYPE(zmumps_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(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) - mat%nnz_start + 1
e = mat%irow(i+1) - mat%nnz_start
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 mumps_mod ***'
END IF
END IF
END SUBROUTINE updt_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putele_mumps_mat(mat, i, j, val)
!
! Put element (i,j) of matrix
!
TYPE(mumps_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) - mat%nnz_start + 1
e = mat%irow(iput+1) - mat%nnz_start
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 mumps_mod ***'
END IF
END IF
END IF
END SUBROUTINE putele_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putele_zmumps_mat(mat, i, j, val)
!
! Put element (i,j) of matrix
!
TYPE(zmumps_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%nlherm .OR. mat%nlsym) 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) - mat%nnz_start + 1
e = mat%irow(iput+1) - mat%nnz_start
k = isearch(mat%cols(s:e), jput)
IF( k.GE.0 ) THEN
mat%val(s+k) = valput
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 mumps_mod ***'
END IF
END IF
END IF
END SUBROUTINE putele_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getele_mumps_mat(mat, i, j, val)
!
! Get element (i,j) of sparse matrix
!
TYPE(mumps_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) - mat%nnz_start + 1
e = mat%irow(iget+1) - mat%nnz_start
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_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getele_zmumps_mat(mat, i, j, val)
!
! Get element (i,j) of sparse matrix
!
TYPE(zmumps_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) - mat%nnz_start + 1
e = mat%irow(iget+1) - mat%nnz_start
k = isearch(mat%cols(s:e), jget)
IF( k.GE.0 ) THEN
valget =mat%val(s+k)
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_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putrow_mumps_mat(amat, i, arr)
!
! Put a row into sparse matrix
!
TYPE(mumps_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_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putrow_zmumps_mat(amat, i, arr)
!
! Put a row into sparse matrix
!
TYPE(zmumps_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_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getrow_mumps_mat(amat, i, arr)
!
! Get a row from sparse matrix
!
TYPE(mumps_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_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getrow_zmumps_mat(amat, i, arr)
!
! Get a row from sparse matrix
!
TYPE(zmumps_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_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putcol_mumps_mat(amat, j, arr)
!
! Put a column into sparse matrix
!
TYPE(mumps_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_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putcol_zmumps_mat(amat, j, arr)
!
! Put a column into sparse matrix
!
TYPE(zmumps_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_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getcol_mumps_mat(amat, j, arr)
!
! Get a column from sparse matrix
!
TYPE(mumps_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_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getcol_zmumps_mat(amat, j, arr)
!
! Get a column from sparse matrix
!
TYPE(zmumps_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_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
INTEGER FUNCTION get_count_mumps_mat(mat, nnz)
!
! Number of non-zeros in sparse matrix
!
TYPE(mumps_mat) :: mat
INTEGER, INTENT(out), OPTIONAL :: nnz(:)
INTEGER :: i
!
IF(ASSOCIATED(mat%mat)) THEN
get_count_mumps_mat = get_count(mat%mat, nnz)
ELSE
get_count_mumps_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_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
INTEGER FUNCTION get_count_zmumps_mat(mat, nnz)
!
! Number of non-zeros in sparse matrix
!
TYPE(zmumps_mat) :: mat
INTEGER, INTENT(out), OPTIONAL :: nnz(:)
INTEGER :: i
!
IF(ASSOCIATED(mat%mat)) THEN
get_count_zmumps_mat = get_count(mat%mat, nnz)
ELSE
get_count_zmumps_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_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE to_mumps_mat(mat, nlkeep)
!
! Convert linked list spmat to mumps matrice structure
!
INCLUDE 'mpif.h'
TYPE(mumps_mat) :: mat
LOGICAL, INTENT(in), OPTIONAL :: nlkeep
INTEGER :: i, nnz, rank, s, e
INTEGER :: comm, ierr, nnz_loc
LOGICAL :: nlclean
!
nlclean = .TRUE.
IF(PRESENT(nlkeep)) THEN
nlclean = .NOT. nlkeep
END IF
!
comm = mat%mumps_par%COMM
mat%mumps_par%ICNTL(18) = 3 ! Distributed assembled matrix
!
! Allocate the Mumps matrix structure
! CSR format: (cols, irow, val) or (JCN, irow, A)
! COO format: (IRN, JCN, A) or (IRN, cols, val)
!
rank = mat%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)
!
mat%mumps_par%N = rank
mat%mumps_par%NZ_loc = nnz_loc
!
IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
IF(ASSOCIATED(mat%perm)) DEALLOCATE(mat%perm)
IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
ALLOCATE(mat%val(nnz_loc))
ALLOCATE(mat%cols(nnz_loc))
ALLOCATE(mat%irow(mat%istart:mat%iend+1))
!
IF(ASSOCIATED(mat%mumps_par%IRN_loc)) DEALLOCATE(mat%mumps_par%IRN_loc)
ALLOCATE(mat%mumps_par%IRN_loc(nnz_loc))
mat%mumps_par%A_loc => mat%val
mat%mumps_par%JCN_loc => mat%cols
!
! Fill Mumps structure and deallocate the sparse rows
!
mat%irow(mat%istart) = mat%nnz_start
DO i=mat%istart,mat%iend
mat%irow(i+1) = mat%irow(i) + mat%mat%row(i)%nnz
s = mat%irow(i) - mat%nnz_start + 1
e = mat%irow(i+1) - mat%nnz_start
CALL getrow(mat%mat%row(i), mat%val(s:e), mat%cols(s:e))
mat%mumps_par%IRN_loc(s:e) = i
IF(nlclean) CALL destroy(mat%mat%row(i))
END DO
IF(nlclean) DEALLOCATE(mat%mat)
END SUBROUTINE to_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE to_zmumps_mat(mat, nlkeep)
!
! Convert linked list spmat to mumps matrice structure
!
INCLUDE 'mpif.h'
TYPE(zmumps_mat) :: mat
LOGICAL, INTENT(in), OPTIONAL :: nlkeep
INTEGER :: i, nnz, rank, s, e
INTEGER :: comm, ierr, nnz_loc
LOGICAL :: nlclean
!
nlclean = .TRUE.
IF(PRESENT(nlkeep)) THEN
nlclean = .NOT. nlkeep
END IF
!
comm = mat%mumps_par%COMM
mat%mumps_par%ICNTL(18) = 3 ! Distributed assembled matrix
!
! Allocate the Mumps matrix structure
! CSR format: (cols, irow, val) or (JCN, irow, A)
! COO format: (IRN, JCN, A) or (IRN, cols, val)
!
rank = mat%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)
!
mat%mumps_par%N = rank
mat%mumps_par%NZ_loc = nnz_loc
!
IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
IF(ASSOCIATED(mat%perm)) DEALLOCATE(mat%perm)
IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
ALLOCATE(mat%val(nnz_loc))
ALLOCATE(mat%cols(nnz_loc))
ALLOCATE(mat%irow(mat%istart:mat%iend+1))
!
IF(ASSOCIATED(mat%mumps_par%IRN_loc)) DEALLOCATE(mat%mumps_par%IRN_loc)
ALLOCATE(mat%mumps_par%IRN_loc(nnz_loc))
mat%mumps_par%A_loc => mat%val
mat%mumps_par%JCN_loc => mat%cols
!
! Fill Mumps structure and deallocate the sparse rows
!
mat%irow(mat%istart) = mat%nnz_start
DO i=mat%istart,mat%iend
mat%irow(i+1) = mat%irow(i) + mat%mat%row(i)%nnz
s = mat%irow(i) - mat%nnz_start + 1
e = mat%irow(i+1) - mat%nnz_start
CALL getrow(mat%mat%row(i), mat%val(s:e), mat%cols(s:e))
mat%mumps_par%IRN_loc(s:e) = i
IF(nlclean) CALL destroy(mat%mat%row(i))
END DO
IF(nlclean) DEALLOCATE(mat%mat)
END SUBROUTINE to_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE reord_mumps_mat(mat, nlmetis, debug)
!
! Reordering and symbolic factorization
!
TYPE(mumps_mat) :: mat
LOGICAL, OPTIONAL, INTENT(in) :: nlmetis
LOGICAL, OPTIONAL, INTENT(in) :: debug
!
! Verbose messages
!
mat%mumps_par%ICNTL(3) = 0
IF(PRESENT(debug)) THEN
IF(debug) mat%mumps_par%ICNTL(3) = 6
END IF
!
! Ordering
!
mat%mumps_par%ICNTL(7) = 7 ! Automatic choice
IF(PRESENT(nlmetis)) THEN
IF(nlmetis) mat%mumps_par%ICNTL(7) = 5 ! use METIS nested dissection
END IF
!
mat%mumps_par%JOB = 1
CALL dmumps(mat%mumps_par)
mat%perm => mat%mumps_par%SYM_PERM
IF(mat%mumps_par%INFOG(1).NE.0) THEN
WRITE(*,'(a,2i12)') 'REORD: Reordering failed with error', &
& mat%mumps_par%INFOG(1:2)
STOP
END IF
END SUBROUTINE reord_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE reord_zmumps_mat(mat, nlmetis, debug)
!
! Reordering and symbolic factorization
!
TYPE(zmumps_mat) :: mat
LOGICAL, OPTIONAL, INTENT(in) :: nlmetis
LOGICAL, OPTIONAL, INTENT(in) :: debug
!
! Verbose messages
!
mat%mumps_par%ICNTL(3) = 0
IF(PRESENT(debug)) THEN
IF(debug) mat%mumps_par%ICNTL(3) = 6
END IF
!
! Ordering
!
mat%mumps_par%ICNTL(7) = 7 ! Automatic choice
IF(PRESENT(nlmetis)) THEN
IF(nlmetis) mat%mumps_par%ICNTL(7) = 5 ! use METIS nested dissection
END IF
!
mat%mumps_par%JOB = 1
CALL zmumps(mat%mumps_par)
IF(mat%mumps_par%INFOG(1).NE.0) THEN
WRITE(*,'(a,2i12)') 'REORD: Reordering failed with error', &
& mat%mumps_par%INFOG(1:2)
STOP
END IF
END SUBROUTINE reord_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE numfact_mumps_mat(mat, debug)
!
! Numerical factorization
!
TYPE(mumps_mat) :: mat
LOGICAL, OPTIONAL, INTENT(in) :: debug
!
! Verbose messages
!
mat%mumps_par%ICNTL(3) = 0
IF(PRESENT(debug)) THEN
IF(debug) mat%mumps_par%ICNTL(3) = 6
END IF
!
mat%mumps_par%JOB = 2
CALL dmumps(mat%mumps_par)
IF(mat%mumps_par%INFOG(1).NE.0) THEN
WRITE(*,'(a,2i12)') 'FACTOR: Reordering failed with error', &
& mat%mumps_par%INFOG(1:2)
STOP
END IF
!
END SUBROUTINE numfact_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE numfact_zmumps_mat(mat, debug)
!
! Numerical factorization
!
TYPE(zmumps_mat) :: mat
LOGICAL, OPTIONAL, INTENT(in) :: debug
!
! Verbose messages
!
mat%mumps_par%ICNTL(3) = 0
IF(PRESENT(debug)) THEN
IF(debug) mat%mumps_par%ICNTL(3) = 6
END IF
!
mat%mumps_par%JOB = 2
CALL zmumps(mat%mumps_par)
IF(mat%mumps_par%INFOG(1).NE.0) THEN
WRITE(*,'(a,2i12)') 'FACTOR: Reordering failed with error', &
& mat%mumps_par%INFOG(1:2)
STOP
END IF
!
END SUBROUTINE numfact_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE factor_mumps_mat(mat, nlreord, nlmetis, debug)
!
! Factor (create +reorder + factor) a mumps_mat matrix
!
TYPE(mumps_mat) :: mat
LOGICAL, OPTIONAL, INTENT(in) :: nlreord
LOGICAL, OPTIONAL, INTENT(in) :: nlmetis
LOGICAL, OPTIONAL, INTENT(in) :: debug
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, nlmetis, debug)
END IF
!----------------------------------------------------------------------
! 3.0 Numerical factorization
!
CALL numfact(mat, debug)
END SUBROUTINE factor_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE factor_zmumps_mat(mat, nlreord, nlmetis, debug)
!
! Factor (create +reorder + factor) a mumps_mat matrix
!
TYPE(zmumps_mat) :: mat
LOGICAL, OPTIONAL, INTENT(in) :: nlreord
LOGICAL, OPTIONAL, INTENT(in) :: nlmetis
LOGICAL, OPTIONAL, INTENT(in) :: debug
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, nlmetis, debug)
END IF
!----------------------------------------------------------------------
! 3.0 Numerical factorization
!
CALL numfact(mat, debug)
END SUBROUTINE factor_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE bsolve_mumps_mat1(mat, rhs, sol, nref, debug)
!
! Backsolve, using Mumps
!
INCLUDE 'mpif.h'
TYPE(mumps_mat) :: mat
DOUBLE PRECISION, INTENT(inout) :: rhs(:)
DOUBLE PRECISION, OPTIONAL, INTENT(inout) :: sol(:)
INTEGER, OPTIONAL, INTENT(in) :: nref
LOGICAL, OPTIONAL, INTENT(in) :: debug
!
INTEGER :: nrank, ierr
!
nrank = SIZE(rhs,1)
!
! Verbose messages
!
mat%mumps_par%ICNTL(3) = 0
IF(PRESENT(debug)) THEN
IF(debug) mat%mumps_par%ICNTL(3) = 6
END IF
!
IF(mat%mumps_par%MYID .EQ. 0) THEN
mat%mumps_par%NRHS = 1
mat%mumps_par%LRHS = nrank
mat%mumps_par%ICNTL(10) = 0 ! Max numbers of iterative refinement steps
IF(PRESENT(nref)) mat%mumps_par%ICNTL(10) = nref
!
ALLOCATE(mat%mumps_par%RHS(nrank))
mat%mumps_par%RHS = rhs
END IF
!
mat%mumps_par%JOB = 3
CALL dmumps(mat%mumps_par)
!
! The solution will be broadcasted to everyone
!
IF(PRESENT(sol)) THEN
IF(mat%mumps_par%MYID .EQ. 0) sol=mat%mumps_par%RHS
CALL mpi_bcast(sol, nrank, MPI_DOUBLE_PRECISION, &
& 0, mat%mumps_par%COMM, ierr)
ELSE
IF(mat%mumps_par%MYID .EQ. 0) rhs=mat%mumps_par%RHS
CALL mpi_bcast(rhs, nrank, MPI_DOUBLE_PRECISION, &
& 0, mat%mumps_par%COMM, ierr)
END IF
!
IF(mat%mumps_par%MYID .EQ. 0) THEN
DEALLOCATE(mat%mumps_par%RHS)
IF(mat%mumps_par%INFOG(1).NE.0) THEN
WRITE(*,'(a,2i12)') 'BSOLVE: Reordering failed with error', &
& mat%mumps_par%INFOG(1:2)
STOP
END IF
END IF
END SUBROUTINE bsolve_mumps_mat1
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE bsolve_zmumps_mat1(mat, rhs, sol, nref, debug)
!
! Backsolve, using Mumps
!
INCLUDE 'mpif.h'
TYPE(zmumps_mat) :: mat
DOUBLE COMPLEX, INTENT(inout) :: rhs(:)
DOUBLE COMPLEX, OPTIONAL, INTENT(inout) :: sol(:)
INTEGER, OPTIONAL, INTENT(in) :: nref
LOGICAL, OPTIONAL, INTENT(in) :: debug
!
INTEGER :: nrank, ierr
!
nrank = SIZE(rhs,1)
!
! Verbose messages
!
mat%mumps_par%ICNTL(3) = 0
IF(PRESENT(debug)) THEN
IF(debug) mat%mumps_par%ICNTL(3) = 6
END IF
!
IF(mat%mumps_par%MYID .EQ. 0) THEN
mat%mumps_par%NRHS = 1
mat%mumps_par%LRHS = nrank
mat%mumps_par%ICNTL(10) = 0 ! Max numbers of iterative refinement steps
IF(PRESENT(nref)) mat%mumps_par%ICNTL(10) = nref
!
ALLOCATE(mat%mumps_par%RHS(nrank))
mat%mumps_par%RHS = rhs
END IF
!
mat%mumps_par%JOB = 3
CALL zmumps(mat%mumps_par)
!
! The solution will be broadcasted to everyone
!
IF(PRESENT(sol)) THEN
IF(mat%mumps_par%MYID .EQ. 0) sol=mat%mumps_par%RHS
CALL mpi_bcast(sol, SIZE(rhs), MPI_DOUBLE_COMPLEX, &
& 0, mat%mumps_par%COMM, ierr)
ELSE
IF(mat%mumps_par%MYID .EQ. 0) rhs=mat%mumps_par%RHS
CALL mpi_bcast(rhs, SIZE(rhs), MPI_DOUBLE_COMPLEX, &
& 0, mat%mumps_par%COMM, ierr)
END IF
!
IF(mat%mumps_par%MYID .EQ. 0) THEN
DEALLOCATE(mat%mumps_par%RHS)
IF(mat%mumps_par%INFOG(1).NE.0) THEN
WRITE(*,'(a,2i12)') 'BSOLVE: Reordering failed with error', &
& mat%mumps_par%INFOG(1:2)
STOP
END IF
END IF
END SUBROUTINE bsolve_zmumps_mat1
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE bsolve_mumps_matn(mat, rhs, sol, nref, debug)
!
! Backsolve, using Mumps
!
INCLUDE 'mpif.h'
TYPE(mumps_mat) :: mat
DOUBLE PRECISION, INTENT(inout) :: rhs(:,:)
DOUBLE PRECISION, OPTIONAL, INTENT(inout) :: sol(:,:)
INTEGER, OPTIONAL, INTENT(in) :: nref
LOGICAL, OPTIONAL, INTENT(in) :: debug
!
INTEGER :: nrank, nrhs, ierr
!
nrank = SIZE(rhs,1)
nrhs = SIZE(rhs,2)
!
! Verbose messages
!
mat%mumps_par%ICNTL(3) = 0
IF(PRESENT(debug)) THEN
IF(debug) mat%mumps_par%ICNTL(3) = 6
END IF
!
IF(mat%mumps_par%MYID .EQ. 0) THEN
mat%mumps_par%NRHS = nrhs
mat%mumps_par%LRHS = nrank
mat%mumps_par%ICNTL(10) = 0 ! Max numbers of iterative refinement steps
IF(PRESENT(nref)) mat%mumps_par%ICNTL(10) = nref
!
ALLOCATE(mat%mumps_par%RHS(nrhs*nrank))
mat%mumps_par%RHS = RESHAPE(rhs, SHAPE(mat%mumps_par%RHS))
END IF
!
mat%mumps_par%JOB = 3
CALL dmumps(mat%mumps_par)
!
! The solution will be broadcasted to everyone
!
IF(PRESENT(sol)) THEN
IF(mat%mumps_par%MYID .EQ. 0) sol=RESHAPE(mat%mumps_par%RHS, SHAPE(sol))
CALL mpi_bcast(sol, PRODUCT(SHAPE(rhs)), MPI_DOUBLE_PRECISION, &
& 0, mat%mumps_par%COMM, ierr)
ELSE
IF(mat%mumps_par%MYID .EQ. 0) rhs=RESHAPE(mat%mumps_par%RHS, SHAPE(rhs))
CALL mpi_bcast(rhs, PRODUCT(SHAPE(rhs)), MPI_DOUBLE_PRECISION, &
& 0, mat%mumps_par%COMM, ierr)
END IF
!
IF(mat%mumps_par%MYID .EQ. 0) THEN
DEALLOCATE(mat%mumps_par%RHS)
IF(mat%mumps_par%INFOG(1).NE.0) THEN
WRITE(*,'(a,2i12)') 'BSOLVE: Reordering failed with error', &
& mat%mumps_par%INFOG(1:2)
STOP
END IF
END IF
END SUBROUTINE bsolve_mumps_matn
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE bsolve_zmumps_matn(mat, rhs, sol, nref, debug)
!
! Backsolve, using Mumps
!
INCLUDE 'mpif.h'
TYPE(zmumps_mat) :: mat
DOUBLE COMPLEX, INTENT(inout) :: rhs(:,:)
DOUBLE COMPLEX, OPTIONAL, INTENT(inout) :: sol(:,:)
INTEGER, OPTIONAL, INTENT(in) :: nref
LOGICAL, OPTIONAL, INTENT(in) :: debug
!
INTEGER :: nrank, nrhs, ierr
!
nrank = SIZE(rhs,1)
nrhs = SIZE(rhs,2)
!
! Verbose messages
!
mat%mumps_par%ICNTL(3) = 0
IF(PRESENT(debug)) THEN
IF(debug) mat%mumps_par%ICNTL(3) = 6
END IF
!
IF(mat%mumps_par%MYID .EQ. 0) THEN
mat%mumps_par%NRHS = nrhs
mat%mumps_par%LRHS = nrank
mat%mumps_par%ICNTL(10) = 0 ! Max numbers of iterative refinement steps
IF(PRESENT(nref)) mat%mumps_par%ICNTL(10) = nref
!
ALLOCATE(mat%mumps_par%RHS(nrhs*nrank))
mat%mumps_par%RHS = RESHAPE(rhs, SHAPE(mat%mumps_par%RHS))
END IF
!
mat%mumps_par%JOB = 3
CALL zmumps(mat%mumps_par)
!
! The solution will be broadcasted to everyone
!
IF(PRESENT(sol)) THEN
IF(mat%mumps_par%MYID .EQ. 0) sol=RESHAPE(mat%mumps_par%RHS, SHAPE(sol))
CALL mpi_bcast(sol, PRODUCT(SHAPE(rhs)), MPI_DOUBLE_COMPLEX, &
& 0, mat%mumps_par%COMM, ierr)
ELSE
IF(mat%mumps_par%MYID .EQ. 0) rhs=RESHAPE(mat%mumps_par%RHS, SHAPE(rhs))
CALL mpi_bcast(rhs, PRODUCT(SHAPE(rhs)), MPI_DOUBLE_COMPLEX, &
& 0, mat%mumps_par%COMM, ierr)
END IF
!
IF(mat%mumps_par%MYID .EQ. 0) THEN
DEALLOCATE(mat%mumps_par%RHS)
IF(mat%mumps_par%INFOG(1).NE.0) THEN
WRITE(*,'(a,2i12)') 'BSOLVE: Reordering failed with error', &
& mat%mumps_par%INFOG(1:2)
STOP
END IF
END IF
END SUBROUTINE bsolve_zmumps_matn
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION vmx_mumps_mat(mat, xarr) RESULT(yarr)
!
! Return product mat*x
!
TYPE(mumps_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_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION vmx_zmumps_mat(mat, xarr) RESULT(yarr)
!
! Return product mat*x
!
TYPE(zmumps_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
!
n = mat%rank
!
#ifdef MKL
IF(mat%nlsym) THEN
matdescra = 'sun'
ELSE IF(mat%nlherm) THEN
matdescra = 'hun'
ELSE
matdescra = 'g'
END IF
CALL mkl_zcsrmv('N', 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
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
ELSE IF(mat%nlherm) THEN
DO j=mat%irow(i)+1, mat%irow(i+1)-1
yarr(mat%cols(j)) = yarr(mat%cols(j)) &
& + CONJG(mat%val(j))*xarr(i)
END DO
END IF
END DO
#endif
!
END FUNCTION vmx_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION vmx_mumps_matn(mat, xarr) RESULT(yarr)
!
! Return product mat*x
!
TYPE(mumps_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_mumps_matn
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION vmx_zmumps_matn(mat, xarr) RESULT(yarr)
!
! Return product mat*x
!
TYPE(zmumps_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
!
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
!
CALL mkl_zcsrmm('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,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
ELSE IF(mat%nlherm) THEN
DO j=mat%irow(i)+1, mat%irow(i+1)-1
yarr(mat%cols(j),:) = yarr(mat%cols(j),:) &
& + CONJG(mat%val(j))*xarr(i,:)
END DO
END IF
END DO
#endif
!
END FUNCTION vmx_zmumps_matn
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE destroy_mumps_mat(mat)
!
! Deallocate the sparse matrix mat
!
TYPE(mumps_mat) :: mat
!
IF(ASSOCIATED(mat%mat)) THEN
CALL destroy(mat%mat)
DEALLOCATE(mat%mat)
END IF
!
IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
!
IF(ASSOCIATED(mat%mumps_par%IRN_loc)) DEALLOCATE(mat%mumps_par%IRN_loc)
mat%mumps_par%JOB = -2
CALL dmumps(mat%mumps_par)
IF(mat%mumps_par%INFOG(1).NE.0) THEN
WRITE(*,'(a,2i12)') 'DESTROY: Reordering failed with error', &
& mat%mumps_par%INFOG(1:2)
STOP
END IF
END SUBROUTINE destroy_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE destroy_zmumps_mat(mat)
!
! Deallocate the sparse matrix mat
!
TYPE(zmumps_mat) :: mat
!
IF(ASSOCIATED(mat%mat)) THEN
CALL destroy(mat%mat)
DEALLOCATE(mat%mat)
END IF
!
IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
!
IF(ASSOCIATED(mat%mumps_par%IRN_loc)) DEALLOCATE(mat%mumps_par%IRN_loc)
mat%mumps_par%JOB = -2
CALL zmumps(mat%mumps_par)
IF(mat%mumps_par%INFOG(1).NE.0) THEN
WRITE(*,'(a,2i12)') 'DESTROY: Reordering failed with error', &
& mat%mumps_par%INFOG(1:2)
STOP
END IF
END SUBROUTINE destroy_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE put_mumps_mat(fid, label, mat, str)
!
! Write matrix to hdf5 file
!
USE futils
!
INTEGER, INTENT(in) :: fid
CHARACTER(len=*), INTENT(in) :: label
TYPE(mumps_mat) :: mat
CHARACTER(len=*), OPTIONAL, INTENT(in) :: str
CHARACTER(len=128) :: mumps_grp
!
IF(PRESENT(str)) THEN
CALL creatg(fid, label, str)
ELSE
CALL creatg(fid, label)
END IF
IF(ASSOCIATED(mat%mat)) THEN
CALL to_mat(mat)
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)//'/perm', mat%perm)
CALL putarr(fid, TRIM(label)//'/val', mat%val)
!
mumps_grp = TRIM(label)//'/mumps_par'
CALL creatg(fid, mumps_grp)
CALL attach(fid, mumps_grp, 'PAR', mat%mumps_par%PAR)
CALL attach(fid, mumps_grp, 'SYM', mat%mumps_par%SYM)
CALL putarr(fid, TRIM(mumps_grp)//'/IRN', mat%mumps_par%IRN)
!
END SUBROUTINE put_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE put_zmumps_mat(fid, label, mat, str)
!
! Write matrix to hdf5 file
!
USE futils
!
INTEGER, INTENT(in) :: fid
CHARACTER(len=*), INTENT(in) :: label
TYPE(zmumps_mat) :: mat
CHARACTER(len=*), OPTIONAL, INTENT(in) :: str
CHARACTER(len=128) :: mumps_grp
!
IF(PRESENT(str)) THEN
CALL creatg(fid, label, str)
ELSE
CALL creatg(fid, label)
END IF
IF(ASSOCIATED(mat%mat)) THEN
CALL to_mat(mat)
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)
!
mumps_grp = TRIM(label)//'/mumps_par'
CALL creatg(fid, mumps_grp)
CALL attach(fid, mumps_grp, 'PAR', mat%mumps_par%PAR)
CALL attach(fid, mumps_grp, 'SYM', mat%mumps_par%SYM)
CALL putarr(fid, TRIM(mumps_grp)//'/IRN', mat%mumps_par%IRN)
!
END SUBROUTINE put_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE get_mumps_mat(fid, label, mat)
!
! Read matrix from hdf5 file
!
USE futils
!
INTEGER, INTENT(in) :: fid
CHARACTER(len=*), INTENT(in) :: label
TYPE(mumps_mat) :: mat
CHARACTER(len=128) :: mumps_grp
!
CALL getatt(fid, label, 'RANK', mat%rank)
CALL getatt(fid, label, 'NNZ', mat%nnz)
CALL getatt(fid, label, 'NLSYM', mat%nlsym)
CALL getarr(fid, TRIM(label)//'/irow', mat%irow)
CALL getarr(fid, TRIM(label)//'/cols', mat%cols)
CALL getarr(fid, TRIM(label)//'/perm', mat%perm)
CALL getarr(fid, TRIM(label)//'/val', mat%val)
!
mumps_grp = TRIM(label)//'/mumps_par'
CALL getatt(fid, mumps_grp, 'PAR', mat%mumps_par%PAR)
CALL getatt(fid, mumps_grp, 'SYM', mat%mumps_par%SYM)
CALL getarr(fid, TRIM(mumps_grp)//'/IRN', mat%mumps_par%IRN)
!
END SUBROUTINE get_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE get_zmumps_mat(fid, label, mat)
!
! Read matrix from hdf5 file
!
USE futils
!
INTEGER, INTENT(in) :: fid
CHARACTER(len=*), INTENT(in) :: label
TYPE(zmumps_mat) :: mat
CHARACTER(len=128) :: mumps_grp
!
CALL getatt(fid, label, 'RANK', mat%rank)
CALL getatt(fid, label, 'NNZ', mat%nnz)
CALL getatt(fid, label, 'NLSYM', mat%nlsym)
CALL getarr(fid, TRIM(label)//'/irow', mat%irow)
CALL getarr(fid, TRIM(label)//'/cols', mat%cols)
CALL getarr(fid, TRIM(label)//'/perm', mat%perm)
CALL getarr(fid, TRIM(label)//'/val', mat%val)
!
mumps_grp = TRIM(label)//'/mumps_par'
CALL getatt(fid, mumps_grp, 'PAR', mat%mumps_par%PAR)
CALL getatt(fid, mumps_grp, 'SYM', mat%mumps_par%SYM)
CALL getarr(fid, TRIM(mumps_grp)//'/IRN', mat%mumps_par%IRN)
!
END SUBROUTINE get_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE mcopy_mumps_mat(mata, matb)
!
! Matrix copy: B = A (assume that B is already initialize)
!
TYPE(mumps_mat) :: mata, matb
INTEGER :: n, nnz, nnz_loc
!
IF(ASSOCIATED(matb%mat)) THEN ! Sparse linled list not needed
CALL destroy(matb%mat)
DEALLOCATE(matb%mat)
END IF
!
n = mata%rank
nnz = mata%nnz
nnz_loc = mata%nnz_loc
matb%nnz = nnz
matb%nnz_loc = nnz_loc
matb%nnz_start = mata%nnz_start
matb%nnz_end = mata%nnz_end
matb%istart = mata%istart
matb%iend = mata%iend
!
matb%mumps_par%NZ_loc = mata%mumps_par%NZ_loc
!
IF(ASSOCIATED(matb%val)) DEALLOCATE(matb%val)
IF(ASSOCIATED(matb%cols)) DEALLOCATE(matb%cols)
IF(ASSOCIATED(matb%irow)) DEALLOCATE(matb%irow)
ALLOCATE(matb%val(nnz_loc)); matb%val = mata%val
ALLOCATE(matb%cols(nnz_loc)); matb%cols = mata%cols
ALLOCATE(matb%irow(matb%istart:matb%iend+1)); matb%irow = mata%irow
!
ALLOCATE(matb%mumps_par%IRN_loc(nnz_loc))
matb%mumps_par%IRN_loc = mata%mumps_par%IRN_loc
matb%mumps_par%A_loc => matb%val
matb%mumps_par%JCN_loc => matb%cols
END SUBROUTINE mcopy_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE mcopy_zmumps_mat(mata, matb)
!
! Matrix copy: B = A (assume that B is already initialize)
!
TYPE(zmumps_mat) :: mata, matb
INTEGER :: n, nnz, nnz_loc
!
IF(ASSOCIATED(matb%mat)) THEN ! Sparse linled list not needed
CALL destroy(matb%mat)
DEALLOCATE(matb%mat)
END IF
!
n = mata%rank
nnz = mata%nnz
nnz_loc = mata%nnz_loc
matb%nnz = nnz
matb%nnz_loc = nnz_loc
matb%nnz_start = mata%nnz_start
matb%nnz_end = mata%nnz_end
matb%istart = mata%istart
matb%iend = mata%iend
!
matb%mumps_par%NZ_loc = mata%mumps_par%NZ_loc
!
IF(ASSOCIATED(matb%val)) DEALLOCATE(matb%val)
IF(ASSOCIATED(matb%cols)) DEALLOCATE(matb%cols)
IF(ASSOCIATED(matb%irow)) DEALLOCATE(matb%irow)
ALLOCATE(matb%val(nnz_loc)); matb%val = mata%val
ALLOCATE(matb%cols(nnz_loc)); matb%cols = mata%cols
ALLOCATE(matb%irow(matb%istart:matb%iend+1)); matb%irow = mata%irow
!
ALLOCATE(matb%mumps_par%IRN_loc(nnz_loc))
matb%mumps_par%IRN_loc = mata%mumps_par%IRN_loc
matb%mumps_par%A_loc => matb%val
matb%mumps_par%JCN_loc => matb%cols
END SUBROUTINE mcopy_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE maddto_mumps_mat(mata, alpha, matb)
!
! A <- A + alpha*B
!
TYPE(mumps_mat) :: mata, matb
DOUBLE PRECISION :: alpha
!
mata%val = mata%val + alpha*matb%val
END SUBROUTINE maddto_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE maddto_zmumps_mat(mata, alpha, matb)
!
! A <- A + alpha*B
!
TYPE(zmumps_mat) :: mata, matb
DOUBLE COMPLEX :: alpha
!
mata%val = mata%val + alpha*matb%val
END SUBROUTINE maddto_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE psum_mumps_mat(mat, comm)
!
! Parallel sum of sparse matrices
!
INCLUDE "mpif.h"
!
TYPE(mumps_mat) :: mat
INCLUDE 'psum_mat.tpl'
END SUBROUTINE psum_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE psum_zmumps_mat(mat, comm)
!
! Parallel sum of sparse matrices
!
INCLUDE "mpif.h"
!
TYPE(zmumps_mat) :: mat
INCLUDE 'psum_mat.tpl'
END SUBROUTINE psum_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE p2p_mumps_mat(mat, dest, extyp, op, comm)
!
! Point-to-point combine sparse matrix between 2 processes
!
INCLUDE "mpif.h"
!
TYPE(mumps_mat) :: mat
DOUBLE PRECISION, ALLOCATABLE :: val(:)
INTEGER :: mpi_type=MPI_DOUBLE_PRECISION
!
INCLUDE "p2p_mat.tpl"
END SUBROUTINE p2p_mumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE p2p_zmumps_mat(mat, dest, extyp, op, comm)
!
! Point-to-point combine sparse matrix between 2 processes
!
INCLUDE "mpif.h"
!
TYPE(zmumps_mat) :: mat
DOUBLE COMPLEX, ALLOCATABLE :: val(:)
INTEGER :: mpi_type=MPI_DOUBLE_COMPLEX
!
INCLUDE "p2p_mat.tpl"
END SUBROUTINE p2p_zmumps_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
END MODULE mumps_bsplines