Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/src/pardiso_mod.f90
8a373078a934feature/refactoring

Last Change
Raw File

pardiso_mod.f90
View Options

!>
!> @file pardiso_mod.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
MODULE pardiso_bsplines
!
! PARDISO_BSPLINES: Simple interface to the sparse direct solver PARDISO
! (MKL version).
!
! T.M. Tran, CRPP-EPFL
! November 2010
!
USE sparse
IMPLICIT NONE
!
TYPE pardiso_param
INTEGER :: error, mtype, msglvl, phase, maxfct, mnum, nrhs
INTEGER :: iparm(64)
INTEGER*8 :: pt(64)
END TYPE pardiso_param
!
TYPE pardiso_mat
INTEGER :: rank, nnz
INTEGER :: nterms, kmat
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(pardiso_param) :: p
END TYPE pardiso_mat
!
TYPE zpardiso_mat
INTEGER :: rank, nnz
INTEGER :: nterms, kmat
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(pardiso_param) :: p
END TYPE zpardiso_mat
!
INTERFACE init
MODULE PROCEDURE init_pardiso_mat, init_zpardiso_mat
END INTERFACE init
!
INTERFACE clear_mat
MODULE PROCEDURE clear_pardiso_mat, clear_zpardiso_mat
END INTERFACE clear_mat
!
INTERFACE updtmat
MODULE PROCEDURE updt_pardiso_mat, updt_zpardiso_mat
END INTERFACE updtmat
!
INTERFACE putele
MODULE PROCEDURE putele_pardiso_mat, putele_zpardiso_mat
END INTERFACE putele
!
INTERFACE getele
MODULE PROCEDURE getele_pardiso_mat, getele_zpardiso_mat
END INTERFACE getele
!
INTERFACE putrow
MODULE PROCEDURE putrow_pardiso_mat, putrow_zpardiso_mat
END INTERFACE putrow
!
INTERFACE getrow
MODULE PROCEDURE getrow_pardiso_mat, getrow_zpardiso_mat
END INTERFACE getrow
!
INTERFACE putcol
MODULE PROCEDURE putcol_pardiso_mat, putcol_zpardiso_mat
END INTERFACE putcol
!
INTERFACE getcol
MODULE PROCEDURE getcol_pardiso_mat, getcol_zpardiso_mat
END INTERFACE getcol
!
INTERFACE get_count
MODULE PROCEDURE get_count_pardiso_mat, get_count_zpardiso_mat
END INTERFACE get_count
!
INTERFACE to_mat
MODULE PROCEDURE to_pardiso_mat, to_zpardiso_mat
END INTERFACE to_mat
!
INTERFACE reord_mat
MODULE PROCEDURE reord_pardiso_mat, reord_zpardiso_mat
END INTERFACE reord_mat
!
INTERFACE numfact
MODULE PROCEDURE numfact_pardiso_mat, numfact_zpardiso_mat
END INTERFACE numfact
!
INTERFACE factor
MODULE PROCEDURE factor_pardiso_mat, factor_zpardiso_mat
END INTERFACE factor
!
INTERFACE bsolve
MODULE PROCEDURE bsolve_pardiso_mat1, bsolve_pardiso_matn, &
& bsolve_zpardiso_mat1, bsolve_zpardiso_matn
END INTERFACE bsolve
!
INTERFACE vmx
MODULE PROCEDURE vmx_pardiso_mat, vmx_pardiso_matn, &
& vmx_zpardiso_mat, vmx_zpardiso_matn
END INTERFACE vmx
!
INTERFACE destroy
MODULE PROCEDURE destroy_pardiso_mat, destroy_zpardiso_mat
END INTERFACE destroy
!
INTERFACE putmat
MODULE PROCEDURE put_pardiso_mat, put_zpardiso_mat
END INTERFACE putmat
!
INTERFACE getmat
MODULE PROCEDURE get_pardiso_mat, get_zpardiso_mat
END INTERFACE getmat
!
INTERFACE mcopy
MODULE PROCEDURE mcopy_pardiso_mat, mcopy_zpardiso_mat
END INTERFACE mcopy
!
INTERFACE maddto
MODULE PROCEDURE maddto_pardiso_mat, maddto_zpardiso_mat
END INTERFACE maddto
!
INTERFACE psum_mat
MODULE PROCEDURE psum_pardiso_mat, psum_zpardiso_mat
END INTERFACE psum_mat
!
INTERFACE p2p_mat
MODULE PROCEDURE p2p_pardiso_mat, p2p_zpardiso_mat
END INTERFACE p2p_mat
!
CONTAINS
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE init_pardiso_mat(n, nterms, mat, kmat, nlsym, nlpos, &
& nlforce_zero)
!
! Initialize an empty sparse pardiso matrix
!
INTEGER, INTENT(in) :: n, nterms
TYPE(pardiso_mat) :: mat
INTEGER, OPTIONAL, INTENT(in) :: kmat
LOGICAL, OPTIONAL, INTENT(in) :: nlsym
LOGICAL, OPTIONAL, INTENT(in) :: nlpos
LOGICAL, OPTIONAL, INTENT(in) :: nlforce_zero
!
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
!
IF(ASSOCIATED(mat%mat)) DEALLOCATE(mat%mat)
ALLOCATE(mat%mat)
CALL init(n, mat%mat)
!
mat%p%iparm = 0
CALL setup_pardiso(mat%p%iparm)
mat%p%maxfct = 1 ! Max number of factorizations
mat%p%mnum = 1 ! Actual matrix, shoild be 1<= num <= maxfct
mat%p%error = 0 ! initialize error flag
mat%p%msglvl = 1 ! print statistical information (0: no stat)
IF(mat%nlsym) THEN
IF(mat%nlpos) THEN
mat%p%mtype = 2 ! symmetric, positive definite
ELSE
mat%p%mtype = -2 ! symmetric, indefinite
END IF
ELSE
mat%p%mtype = 11 ! unsymmetric
END IF
mat%p%nrhs = 1 ! number of RHSs
mat%p%pt(1:64) = 0 ! Initialize Pardiso address pointer (handle)
!
CONTAINS
SUBROUTINE setup_pardiso(iparm)
INTEGER :: iparm(:)
iparm(1) = 1 ! no solver default
!!$ iparm(2) = 2 ! fill-in reordering from METIS
iparm(2) = 0 ! Minimum degree fill-in reordering
iparm(3) = 1 ! numbers of processors
iparm(4) = 0 ! no iterative-direct algorithm
iparm(5) = 0 ! no user fill-in reducing permutation
iparm(6) = 0 ! =0 solution on the first n compoments of x
iparm(7) = 0 ! not in use
iparm(8) = 9 ! numbers of iterative refinement steps
iparm(9) = 0 ! not in use
iparm(10) = 13 ! perturbe the pivot elements with 1E-13
iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
iparm(12) = 0 ! not in use
iparm(13) = 0 ! not in use
iparm(14) = 0 ! Output: number of perturbed pivots
iparm(15) = 0 ! not in use
iparm(16) = 0 ! not in use
iparm(17) = 0 ! not in use
iparm(18) = -1 ! Output: number of nonzeros in the factor LU
iparm(19) = -1 ! Output: Mflops for LU factorization
iparm(20) = 0 ! Output: Numbers of CG Iterations
END SUBROUTINE setup_pardiso
END SUBROUTINE init_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE init_zpardiso_mat(n, nterms, mat, kmat, nlsym, nlherm, &
& nlpos, nlforce_zero)
!
! Initialize an empty sparse pardiso matrix
!
INTEGER, INTENT(in) :: n, nterms
TYPE(zpardiso_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
!
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(nlherm)) mat%nlherm = nlherm
IF(PRESENT(nlpos)) mat%nlpos = nlpos
IF(PRESENT(nlforce_zero)) mat%nlforce_zero = nlforce_zero
!
IF(ASSOCIATED(mat%mat)) DEALLOCATE(mat%mat)
ALLOCATE(mat%mat)
CALL init(n, mat%mat)
!
mat%p%iparm = 0
CALL setup_pardiso(mat%p%iparm)
mat%p%maxfct = 1 ! Max number of factorizations
mat%p%mnum = 1 ! Actual matrix, shoild be 1<= num <= maxfct
mat%p%error = 0 ! initialize error flag
mat%p%msglvl = 1 ! print statistical information (0: no stat)
IF(mat%nlherm) THEN
IF(mat%nlpos) THEN
mat%p%mtype = 4 ! hermitian, positive definite
ELSE
mat%p%mtype = -4 ! hermitian, indefinite
END IF
ELSE IF(mat%nlsym) THEN
mat%p%mtype = 6 ! symmetric
ELSE
mat%p%mtype = 13 ! unsymmetric
END IF
mat%p%nrhs = 1 ! number of RHSs
mat%p%pt(1:64) = 0 ! Initialize Pardiso address pointer (handle)
!
CONTAINS
SUBROUTINE setup_pardiso(iparm)
INTEGER :: iparm(:)
iparm(1) = 1 ! no solver default
!!$ iparm(2) = 2 ! fill-in reordering from METIS
iparm(2) = 0 ! Minimum degree fill-in reordering
iparm(3) = 1 ! numbers of processors
iparm(4) = 0 ! no iterative-direct algorithm
iparm(5) = 0 ! no user fill-in reducing permutation
iparm(6) = 0 ! =0 solution on the first n compoments of x
iparm(7) = 0 ! not in use
iparm(8) = 9 ! numbers of iterative refinement steps
iparm(9) = 0 ! not in use
iparm(10) = 13 ! perturbe the pivot elements with 1E-13
iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
iparm(12) = 0 ! not in use
iparm(13) = 0 ! not in use
iparm(14) = 0 ! Output: number of perturbed pivots
iparm(15) = 0 ! not in use
iparm(16) = 0 ! not in use
iparm(17) = 0 ! not in use
iparm(18) = -1 ! Output: number of nonzeros in the factor LU
iparm(19) = -1 ! Output: Mflops for LU factorization
iparm(20) = 0 ! Output: Numbers of CG Iterations
END SUBROUTINE setup_pardiso
END SUBROUTINE init_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE clear_pardiso_mat(mat)
!
! Clear matrix, keeping its sparse structure unchanged
!
TYPE(pardiso_mat) :: mat
!
mat%val = 0.0d0
mat%perm = 0
END SUBROUTINE clear_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE clear_zpardiso_mat(mat)
!
! Clear matrix, keeping its sparse structure unchanged
!
TYPE(zpardiso_mat) :: mat
!
mat%val = (0.0d0, 0.0d0)
mat%perm = 0
END SUBROUTINE clear_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE updt_pardiso_mat(mat, i, j, val)
!
! Update element Aij of pardiso matrix
!
TYPE(pardiso_mat) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE PRECISION, INTENT(in) :: val
INTEGER :: k, s, e
!
IF(mat%nlsym) THEN ! Store only upper part for symmetric matrices
IF(i.GT.j) RETURN
END IF
!
IF(ASSOCIATED(mat%mat)) THEN
CALL updtmat(mat%mat, i, j, val)
ELSE
s = mat%irow(i)
e = mat%irow(i+1)-1
k = isearch(mat%cols(s:e), j)
IF( k.GE.0 ) THEN
mat%val(s+k) = mat%val(s+k)+val
ELSE
WRITE(*,'(a,2i6)') 'UPDT: i, j out of range ', i, j
STOP '*** Abnormal EXIT in MODULE pardiso_mod ***'
END IF
END IF
END SUBROUTINE updt_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE updt_zpardiso_mat(mat, i, j, val)
!
! Update element Aij of pardiso matrix
!
TYPE(zpardiso_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
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 pardiso_mod ***'
END IF
END IF
END SUBROUTINE updt_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putele_pardiso_mat(mat, i, j, val)
!
! Put element (i,j) of matrix
!
TYPE(pardiso_mat) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE PRECISION, INTENT(in) :: val
INTEGER :: iput, jput
INTEGER :: k, s, e
!
iput = i
jput = j
IF(mat%nlsym) THEN
IF( i.GT.j ) THEN ! Lower triangular part
iput = j
jput = i
END IF
END IF
!
IF(ASSOCIATED(mat%mat)) THEN
CALL putele(mat%mat, iput, jput, val, &
& nlforce_zero=mat%nlforce_zero)
ELSE
s = mat%irow(iput)
e = mat%irow(iput+1)-1
k = isearch(mat%cols(s:e), jput)
IF( k.GE.0 ) THEN
mat%val(s+k) = val
ELSE
IF(ABS(val) .GT. EPSILON(0.0d0)) THEN ! Ok for zero val
WRITE(*,'(a,2i6)') 'PUTELE: i, j out of range ', i, j
STOP '*** Abnormal EXIT in MODULE pardiso_mod ***'
END IF
END IF
END IF
END SUBROUTINE putele_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putele_zpardiso_mat(mat, i, j, val)
!
! Put element (i,j) of matrix
!
TYPE(zpardiso_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
!
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
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 pardiso_mod ***'
END IF
END IF
END IF
END SUBROUTINE putele_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getele_pardiso_mat(mat, i, j, val)
!
! Get element (i,j) of sparse matrix
!
TYPE(pardiso_mat) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE PRECISION, INTENT(out) :: val
INTEGER :: iget, jget
INTEGER :: k, s, e
!
iget = i
jget = j
IF(mat%nlsym) THEN
IF( i.GT.j ) THEN ! Lower triangular part
iget = j
jget = i
END IF
END IF
!
IF(ASSOCIATED(mat%mat)) THEN
CALL getele(mat%mat, iget, jget, val)
ELSE
s = mat%irow(iget)
e = mat%irow(iget+1)-1
k = isearch(mat%cols(s:e), jget)
IF( k.GE.0 ) THEN
val =mat%val(s+k)
ELSE
val = 0.0d0 ! Assume zero val if not found
END IF
END IF
END SUBROUTINE getele_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getele_zpardiso_mat(mat, i, j, val)
!
! Get element (i,j) of sparse matrix
!
TYPE(zpardiso_mat) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE COMPLEX, INTENT(out) :: val
DOUBLE COMPLEX :: valget
INTEGER :: iget, jget
INTEGER :: k, s, e
!
iget = i
jget = j
IF(mat%nlherm .OR. mat%nlsym) THEN
IF( i.GT.j ) THEN ! Lower triangular part
iget = j
jget = i
END IF
END IF
!
IF(ASSOCIATED(mat%mat)) THEN
CALL getele(mat%mat, iget, jget, valget)
ELSE
s = mat%irow(iget)
e = mat%irow(iget+1)-1
k = isearch(mat%cols(s:e), jget)
IF( k.GE.0 ) THEN
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_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putrow_pardiso_mat(amat, i, arr)
!
! Put a row into sparse matrix
!
TYPE(pardiso_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_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putrow_zpardiso_mat(amat, i, arr)
!
! Put a row into sparse matrix
!
TYPE(zpardiso_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_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getrow_pardiso_mat(amat, i, arr)
!
! Get a row from sparse matrix
!
TYPE(pardiso_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_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getrow_zpardiso_mat(amat, i, arr)
!
! Get a row from sparse matrix
!
TYPE(zpardiso_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_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putcol_pardiso_mat(amat, j, arr)
!
! Put a column into sparse matrix
!
TYPE(pardiso_mat), INTENT(inout) :: amat
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, INTENT(in) :: arr(:)
INTEGER :: i
!
DO i=1,amat%rank
CALL putele(amat, i, j, arr(i))
END DO
END SUBROUTINE putcol_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putcol_zpardiso_mat(amat, j, arr)
!
! Put a column into sparse matrix
!
TYPE(zpardiso_mat), INTENT(inout) :: amat
INTEGER, INTENT(in) :: j
DOUBLE COMPLEX, INTENT(in) :: arr(:)
INTEGER :: i
!
DO i=1,amat%rank
CALL putele(amat, i, j, arr(i))
END DO
END SUBROUTINE putcol_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getcol_pardiso_mat(amat, j, arr)
!
! Get a column from sparse matrix
!
TYPE(pardiso_mat), INTENT(in) :: amat
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, INTENT(out) :: arr(:)
INTEGER :: i
!
DO i=1,amat%rank
CALL getele(amat, i, j, arr(i))
END DO
END SUBROUTINE getcol_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getcol_zpardiso_mat(amat, j, arr)
!
! Get a column from sparse matrix
!
TYPE(zpardiso_mat), INTENT(in) :: amat
INTEGER, INTENT(in) :: j
DOUBLE COMPLEX, INTENT(out) :: arr(:)
INTEGER :: i
!
DO i=1,amat%rank
CALL getele(amat, i, j, arr(i))
END DO
END SUBROUTINE getcol_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
INTEGER FUNCTION get_count_pardiso_mat(mat, nnz)
!
! Number of non-zeros in sparse matrix
!
TYPE(pardiso_mat) :: mat
INTEGER, INTENT(out), OPTIONAL :: nnz(:)
INTEGER :: i
!
IF(ASSOCIATED(mat%mat)) THEN
get_count_pardiso_mat = get_count(mat%mat, nnz)
ELSE
get_count_pardiso_mat = mat%nnz
IF(PRESENT(nnz)) THEN
DO i=1,mat%rank
nnz(i) = mat%irow(i+1)-mat%irow(i)
END DO
END IF
END IF
END FUNCTION get_count_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
INTEGER FUNCTION get_count_zpardiso_mat(mat, nnz)
!
! Number of non-zeros in sparse matrix
!
TYPE(zpardiso_mat) :: mat
INTEGER, INTENT(out), OPTIONAL :: nnz(:)
INTEGER :: i
!
IF(ASSOCIATED(mat%mat)) THEN
get_count_zpardiso_mat = get_count(mat%mat, nnz)
ELSE
get_count_zpardiso_mat = mat%nnz
IF(PRESENT(nnz)) THEN
DO i=1,mat%rank
nnz(i) = mat%irow(i+1)-mat%irow(i)
END DO
END IF
END IF
END FUNCTION get_count_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE to_pardiso_mat(mat, nlkeep)
!
! Convert linked list spmat to pardiso matrice structure
!
TYPE(pardiso_mat) :: mat
LOGICAL, INTENT(in), OPTIONAL :: nlkeep
INTEGER :: i, nnz, rank, s, e
LOGICAL :: nlclean
!
nlclean = .TRUE.
IF(PRESENT(nlkeep)) THEN
nlclean = .NOT. nlkeep
END IF
!
! Allocate the Pardiso matrix structure
!
nnz = get_count(mat)
rank = mat%rank
mat%nnz = nnz
IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
IF(ASSOCIATED(mat%perm)) DEALLOCATE(mat%perm)
IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
ALLOCATE(mat%val(nnz))
ALLOCATE(mat%perm(rank))
ALLOCATE(mat%irow(rank+1))
ALLOCATE(mat%cols(nnz))
!
! Fill Pardiso structure and optionnaly deallocate the sparse rows
!
mat%irow = 1
DO i=1,rank
mat%irow(i+1) = mat%irow(i) + mat%mat%row(i)%nnz
s = mat%irow(i)
e = mat%irow(i+1)-1
CALL getrow(mat%mat%row(i), mat%val(s:e), mat%cols(s:e))
IF(nlclean) CALL destroy(mat%mat%row(i))
END DO
IF(nlclean) DEALLOCATE(mat%mat)
END SUBROUTINE to_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE to_zpardiso_mat(mat, nlkeep)
!
! Convert linked list spmat to pardiso matrice structure
!
TYPE(zpardiso_mat) :: mat
LOGICAL, INTENT(in), OPTIONAL :: nlkeep
INTEGER :: i, nnz, rank, s, e
LOGICAL :: nlclean
!
nlclean = .TRUE.
IF(PRESENT(nlkeep)) THEN
nlclean = .NOT. nlkeep
END IF
!
! Allocate the Pardiso matrix structure
!
nnz = get_count(mat)
rank = mat%rank
mat%nnz = nnz
IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
IF(ASSOCIATED(mat%perm)) DEALLOCATE(mat%perm)
IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
ALLOCATE(mat%val(nnz))
ALLOCATE(mat%perm(rank))
ALLOCATE(mat%irow(rank+1))
ALLOCATE(mat%cols(nnz))
!
! Fill Pardiso structure and deallocate the sparse rows
!
mat%irow = 1
DO i=1,rank
mat%irow(i+1) = mat%irow(i) + mat%mat%row(i)%nnz
s = mat%irow(i)
e = mat%irow(i+1)-1
CALL getrow(mat%mat%row(i), mat%val(s:e), mat%cols(s:e))
IF(nlclean) CALL destroy(mat%mat%row(i))
END DO
IF(nlclean) DEALLOCATE(mat%mat)
END SUBROUTINE to_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE reord_pardiso_mat(mat, nlmetis, debug)
!
! Reordering and symbolic factorization
!
TYPE(pardiso_mat) :: mat
LOGICAL, OPTIONAL, INTENT(in) :: nlmetis
LOGICAL, OPTIONAL, INTENT(in) :: debug
DOUBLE PRECISION :: dummy
!
mat%p%iparm(2) = 0 ! use minimum degree algorithm
IF(PRESENT(nlmetis)) THEN
IF(nlmetis) mat%p%iparm(2) = 2 ! use METIS nested dissection
END IF
mat%p%iparm(5)= 2 ! return the permutation vector in mat%perm
mat%p%phase = 11 ! Reordering and symbolic factorization
mat%p%msglvl = 0
IF(PRESENT(debug)) THEN
IF(debug) mat%p%msglvl = 1
END IF
CALL pardiso(mat%p%pt, mat%p%maxfct, mat%p%mnum, mat%p%mtype, &
& mat%p%phase, mat%rank, mat%val, mat%irow, mat%cols, &
& mat%perm, mat%p%nrhs, mat%p%iparm, mat%p%msglvl, &
& dummy, dummy, mat%p%error)
IF(mat%p%error.NE.0) THEN
WRITE(*,'(a,i4)') 'REORD: Reordering failed with error', mat%p%error
END IF
END SUBROUTINE reord_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE reord_zpardiso_mat(mat, nlmetis, debug)
!
! Reordering and symbolic factorization
!
TYPE(zpardiso_mat) :: mat
LOGICAL, OPTIONAL, INTENT(in) :: nlmetis
LOGICAL, OPTIONAL, INTENT(in) :: debug
DOUBLE COMPLEX :: dummy
!
mat%p%iparm(2) = 0 ! use minimum degree algorithm
IF(PRESENT(nlmetis)) THEN
IF(nlmetis) mat%p%iparm(2) = 2 ! use METIS nested dissection
END IF
mat%p%iparm(5)= 2 ! return the permutation vector in mat%perm
mat%p%phase = 11 ! Reordering and symbolic factorization
mat%p%msglvl = 0
IF(PRESENT(debug)) THEN
IF(debug) mat%p%msglvl = 1
END IF
CALL pardiso(mat%p%pt, mat%p%maxfct, mat%p%mnum, mat%p%mtype, &
& mat%p%phase, mat%rank, mat%val, mat%irow, mat%cols, &
& mat%perm, mat%p%nrhs, mat%p%iparm, mat%p%msglvl, &
& dummy, dummy, mat%p%error)
IF(mat%p%error.NE.0) THEN
WRITE(*,'(a,i4)') 'REORD: Reordering failed with error', mat%p%error
END IF
END SUBROUTINE reord_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE numfact_pardiso_mat(mat, debug)
!
! Numerical factorization
!
TYPE(pardiso_mat) :: mat
LOGICAL, OPTIONAL, INTENT(in) :: debug
DOUBLE PRECISION :: dummy
!
mat%p%phase = 22 ! Numerical factorization
mat%p%msglvl = 0
IF(PRESENT(debug)) THEN
IF(debug) mat%p%msglvl = 1
END IF
CALL pardiso(mat%p%pt, mat%p%maxfct, mat%p%mnum, mat%p%mtype, &
& mat%p%phase, mat%rank, mat%val, mat%irow, mat%cols, &
& mat%perm, mat%p%nrhs, mat%p%iparm, mat%p%msglvl, &
& dummy, dummy, mat%p%error)
IF(mat%p%error.NE.0) THEN
WRITE(*,'(a,i4)') 'FACTOR: Factorization failed with error', mat%p%error
END IF
END SUBROUTINE numfact_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE numfact_zpardiso_mat(mat, debug)
!
! Numerical factorization
!
TYPE(zpardiso_mat) :: mat
LOGICAL, OPTIONAL, INTENT(in) :: debug
DOUBLE COMPLEX :: dummy
!
mat%p%phase = 22 ! Numerical factorization
mat%p%msglvl = 0
IF(PRESENT(debug)) THEN
IF(debug) mat%p%msglvl = 1
END IF
CALL pardiso(mat%p%pt, mat%p%maxfct, mat%p%mnum, mat%p%mtype, &
& mat%p%phase, mat%rank, mat%val, mat%irow, mat%cols, &
& mat%perm, mat%p%nrhs, mat%p%iparm, mat%p%msglvl, &
& dummy, dummy, mat%p%error)
IF(mat%p%error.NE.0) THEN
WRITE(*,'(a,i4)') 'FACTOR: Factorization failed with error', mat%p%error
END IF
END SUBROUTINE numfact_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE factor_pardiso_mat(mat, nlreord, nlmetis, debug)
!
! Factor (create +reorder + factor) a pardiso_mat matrix
!
TYPE(pardiso_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_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE factor_zpardiso_mat(mat, nlreord, nlmetis, debug)
!
! Factor (create +reorder + factor) a pardiso_mat matrix
!
TYPE(zpardiso_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_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE bsolve_pardiso_mat1(mat, rhs, sol, nref, debug)
!
! Backsolve, using Pardiso
!
TYPE(pardiso_mat) :: mat
DOUBLE PRECISION :: rhs(:)
DOUBLE PRECISION, OPTIONAL :: sol(:)
INTEGER, OPTIONAL :: nref
LOGICAL, OPTIONAL, INTENT(in) :: debug
DOUBLE PRECISION :: dummy(SIZE(rhs))
!
mat%p%iparm(8) = 0 ! Max numbers of iterative refinement steps
IF(PRESENT(nref)) mat%p%iparm(8) = nref
mat%p%phase = 33 ! Backsolve
mat%p%nrhs = 1
mat%p%msglvl = 0
IF(PRESENT(debug)) THEN
IF(debug) THEN
mat%p%msglvl = 1
END IF
END IF
IF(PRESENT(sol)) THEN
mat%p%iparm(6) = 0
CALL pardiso(mat%p%pt, mat%p%maxfct, mat%p%mnum, mat%p%mtype, &
& mat%p%phase, mat%rank, mat%val, mat%irow, mat%cols, &
& mat%perm, mat%p%nrhs, mat%p%iparm, mat%p%msglvl, &
& rhs, sol, mat%p%error)
ELSE
mat%p%iparm(6) = 1
CALL pardiso(mat%p%pt, mat%p%maxfct, mat%p%mnum, mat%p%mtype, &
& mat%p%phase, mat%rank, mat%val, mat%irow, mat%cols, &
& mat%perm, mat%p%nrhs, mat%p%iparm, mat%p%msglvl, &
& rhs, dummy, mat%p%error)
END IF
IF(mat%p%error.NE.0) THEN
WRITE(*,'(a,i4)') 'BSOLVE: Failed with error', mat%p%error
END IF
END SUBROUTINE bsolve_pardiso_mat1
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE bsolve_zpardiso_mat1(mat, rhs, sol, nref, debug)
!
! Backsolve, using Pardiso
!
TYPE(zpardiso_mat) :: mat
DOUBLE COMPLEX :: rhs(:)
DOUBLE COMPLEX, OPTIONAL :: sol(:)
INTEGER, OPTIONAL :: nref
LOGICAL, OPTIONAL, INTENT(in) :: debug
DOUBLE COMPLEX :: dummy(SIZE(rhs))
!
mat%p%iparm(8) = 0 ! Max numbers of iterative refinement steps
IF(PRESENT(nref)) mat%p%iparm(8) = nref
mat%p%phase = 33 ! Backsolve
mat%p%nrhs = 1
mat%p%msglvl = 0
IF(PRESENT(debug)) THEN
IF(debug) THEN
mat%p%msglvl = 1
END IF
END IF
IF(PRESENT(sol)) THEN
mat%p%iparm(6) = 0
CALL pardiso(mat%p%pt, mat%p%maxfct, mat%p%mnum, mat%p%mtype, &
& mat%p%phase, mat%rank, mat%val, mat%irow, mat%cols, &
& mat%perm, mat%p%nrhs, mat%p%iparm, mat%p%msglvl, &
& rhs, sol, mat%p%error)
ELSE
mat%p%iparm(6) = 1
CALL pardiso(mat%p%pt, mat%p%maxfct, mat%p%mnum, mat%p%mtype, &
& mat%p%phase, mat%rank, mat%val, mat%irow, mat%cols, &
& mat%perm, mat%p%nrhs, mat%p%iparm, mat%p%msglvl, &
& rhs, dummy, mat%p%error)
END IF
IF(mat%p%error.NE.0) THEN
WRITE(*,'(a,i4)') 'BSOLVE: Failed with error', mat%p%error
END IF
END SUBROUTINE bsolve_zpardiso_mat1
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE bsolve_pardiso_matn(mat, rhs, sol, nref, debug)
!
! Backsolve, using Pardiso, multiple RHS
!
TYPE(pardiso_mat) :: mat
DOUBLE PRECISION :: rhs(:,:)
DOUBLE PRECISION, OPTIONAL :: sol(:,:)
INTEGER, OPTIONAL :: nref
LOGICAL, OPTIONAL, INTENT(in) :: debug
DOUBLE PRECISION :: dummy(SIZE(rhs,1),SIZE(rhs,2))
!
mat%p%iparm(8) = 0 ! Max numbers of iterative refinement steps
IF(PRESENT(nref)) mat%p%iparm(8) = nref
mat%p%phase = 33 ! Backsolve
mat%p%nrhs = SIZE(rhs,2)
mat%p%msglvl = 0
IF(PRESENT(debug)) THEN
IF(debug) THEN
mat%p%msglvl = 1
END IF
END IF
IF(PRESENT(sol)) THEN
mat%p%iparm(6) = 0
CALL pardiso(mat%p%pt, mat%p%maxfct, mat%p%mnum, mat%p%mtype, &
& mat%p%phase, mat%rank, mat%val, mat%irow, mat%cols, &
& mat%perm, mat%p%nrhs, mat%p%iparm, mat%p%msglvl, &
& rhs, sol, mat%p%error)
ELSE
mat%p%iparm(6) = 1
CALL pardiso(mat%p%pt, mat%p%maxfct, mat%p%mnum, mat%p%mtype, &
& mat%p%phase, mat%rank, mat%val, mat%irow, mat%cols, &
& mat%perm, mat%p%nrhs, mat%p%iparm, mat%p%msglvl, &
& rhs, dummy, mat%p%error)
END IF
IF(mat%p%error.NE.0) THEN
WRITE(*,'(a,i4)') 'BSOLVE: Failed with error', mat%p%error
END IF
END SUBROUTINE bsolve_pardiso_matn
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE bsolve_zpardiso_matn(mat, rhs, sol, nref, debug)
!
! Backsolve, using Pardiso, multiple RHS
!
TYPE(zpardiso_mat) :: mat
DOUBLE COMPLEX :: rhs(:,:)
DOUBLE COMPLEX, OPTIONAL :: sol(:,:)
INTEGER, OPTIONAL :: nref
LOGICAL, OPTIONAL, INTENT(in) :: debug
DOUBLE COMPLEX :: dummy(SIZE(rhs,1),SIZE(rhs,2))
!
mat%p%iparm(8) = 0 ! Max numbers of iterative refinement steps
IF(PRESENT(nref)) mat%p%iparm(8) = nref
mat%p%phase = 33 ! Backsolve
mat%p%nrhs = SIZE(rhs,2)
mat%p%msglvl = 0
IF(PRESENT(debug)) THEN
IF(debug) THEN
mat%p%msglvl = 1
END IF
END IF
IF(PRESENT(sol)) THEN
mat%p%iparm(6) = 0
CALL pardiso(mat%p%pt, mat%p%maxfct, mat%p%mnum, mat%p%mtype, &
& mat%p%phase, mat%rank, mat%val, mat%irow, mat%cols, &
& mat%perm, mat%p%nrhs, mat%p%iparm, mat%p%msglvl, &
& rhs, sol, mat%p%error)
ELSE
mat%p%iparm(6) = 1
CALL pardiso(mat%p%pt, mat%p%maxfct, mat%p%mnum, mat%p%mtype, &
& mat%p%phase, mat%rank, mat%val, mat%irow, mat%cols, &
& mat%perm, mat%p%nrhs, mat%p%iparm, mat%p%msglvl, &
& rhs, dummy, mat%p%error)
END IF
IF(mat%p%error.NE.0) THEN
WRITE(*,'(a,i4)') 'BSOLVE: Failed with error', mat%p%error
END IF
END SUBROUTINE bsolve_zpardiso_matn
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION vmx_pardiso_mat(mat, xarr) RESULT(yarr)
!
! Return product mat*x
!
TYPE(pardiso_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_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION vmx_zpardiso_mat(mat, xarr) RESULT(yarr)
!
! Return product mat*x
!
TYPE(zpardiso_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_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION vmx_pardiso_matn(mat, xarr) RESULT(yarr)
!
! Return product mat*x
!
TYPE(pardiso_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_pardiso_matn
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION vmx_zpardiso_matn(mat, xarr) RESULT(yarr)
!
! Return product mat*x
!
TYPE(zpardiso_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_zpardiso_matn
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE destroy_pardiso_mat(mat)
!
! Deallocate the sparse matrix mat
!
TYPE(pardiso_mat) :: mat
DOUBLE PRECISION :: dummy
!
IF(ASSOCIATED(mat%mat)) THEN
CALL destroy(mat%mat)
DEALLOCATE(mat%mat)
END IF
!
IF(mat%p%phase .GT. 0) THEN
mat%p%phase = 0 ! Release memory for factors
CALL pardiso(mat%p%pt, mat%p%maxfct, mat%p%mnum, mat%p%mtype, &
& mat%p%phase, mat%rank, mat%val, mat%irow, mat%cols, &
& mat%perm, mat%p%nrhs, mat%p%iparm, mat%p%msglvl, &
& dummy, dummy, mat%p%error)
IF(mat%p%error.NE.0) THEN
WRITE(*,'(a,i4)') 'DESTROY: Mem release failed with error', mat%p%error
END IF
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%val)) DEALLOCATE(mat%val)
END SUBROUTINE destroy_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE destroy_zpardiso_mat(mat)
!
! Deallocate the sparse matrix mat
!
TYPE(zpardiso_mat) :: mat
DOUBLE COMPLEX :: dummy
!
IF(ASSOCIATED(mat%mat)) THEN
CALL destroy(mat%mat)
DEALLOCATE(mat%mat)
END IF
!
IF(mat%p%phase .GT. 0) THEN
mat%p%phase = 0 ! Release memory for factors
CALL pardiso(mat%p%pt, mat%p%maxfct, mat%p%mnum, mat%p%mtype, &
& mat%p%phase, mat%rank, mat%val, mat%irow, mat%cols, &
& mat%perm, mat%p%nrhs, mat%p%iparm, mat%p%msglvl, &
& dummy, dummy, mat%p%error)
IF(mat%p%error.NE.0) THEN
WRITE(*,'(a,i4)') 'DESTROY: Mem release failed with error', mat%p%error
END IF
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%val)) DEALLOCATE(mat%val)
END SUBROUTINE destroy_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE put_pardiso_mat(fid, label, mat, str)
!
! Write matrix to hdf5 file
!
USE futils
!
INTEGER, INTENT(in) :: fid
CHARACTER(len=*), INTENT(in) :: label
TYPE(pardiso_mat) :: mat
CHARACTER(len=*), OPTIONAL, INTENT(in) :: str
!
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)
!
CALL creatg(fid, TRIM(label)//'/p')
CALL attach(fid, TRIM(label)//'/p', 'error', mat%p%error)
CALL attach(fid, TRIM(label)//'/p', 'mtype', mat%p%mtype)
CALL attach(fid, TRIM(label)//'/p', 'phase', mat%p%phase)
CALL attach(fid, TRIM(label)//'/p', 'msglv', mat%p%msglvl)
CALL attach(fid, TRIM(label)//'/p', 'maxfct', mat%p%maxfct)
CALL attach(fid, TRIM(label)//'/p', 'mnum', mat%p%mnum)
CALL attach(fid, TRIM(label)//'/p', 'nrhs', mat%p%nrhs)
CALL putarr(fid, TRIM(label)//'/p/iparm', mat%p%iparm)
END SUBROUTINE put_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE put_zpardiso_mat(fid, label, mat, str)
!
! Write matrix to hdf5 file
!
USE futils
!
INTEGER, INTENT(in) :: fid
CHARACTER(len=*), INTENT(in) :: label
TYPE(zpardiso_mat) :: mat
CHARACTER(len=*), OPTIONAL, INTENT(in) :: str
!
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 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)//'/perm', mat%perm)
CALL putarr(fid, TRIM(label)//'/val', mat%val)
!
CALL creatg(fid, TRIM(label)//'/p')
CALL attach(fid, TRIM(label)//'/p', 'error', mat%p%error)
CALL attach(fid, TRIM(label)//'/p', 'mtype', mat%p%mtype)
CALL attach(fid, TRIM(label)//'/p', 'phase', mat%p%phase)
CALL attach(fid, TRIM(label)//'/p', 'msglv', mat%p%msglvl)
CALL attach(fid, TRIM(label)//'/p', 'maxfct', mat%p%maxfct)
CALL attach(fid, TRIM(label)//'/p', 'mnum', mat%p%mnum)
CALL attach(fid, TRIM(label)//'/p', 'nrhs', mat%p%nrhs)
CALL putarr(fid, TRIM(label)//'/p/iparm', mat%p%iparm)
END SUBROUTINE put_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE get_pardiso_mat(fid, label, mat)
!
! Read matrix from hdf5 file
!
USE futils
!
INTEGER, INTENT(in) :: fid
CHARACTER(len=*), INTENT(in) :: label
TYPE(pardiso_mat) :: mat
!
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)
!
CALL getatt(fid, TRIM(label)//'/p', 'error', mat%p%error)
CALL getatt(fid, TRIM(label)//'/p', 'mtype', mat%p%mtype)
CALL getatt(fid, TRIM(label)//'/p', 'phase', mat%p%phase)
CALL getatt(fid, TRIM(label)//'/p', 'msglv', mat%p%msglvl)
CALL getatt(fid, TRIM(label)//'/p', 'maxfct', mat%p%maxfct)
CALL getatt(fid, TRIM(label)//'/p', 'mnum', mat%p%mnum)
CALL getatt(fid, TRIM(label)//'/p', 'nrhs', mat%p%nrhs)
CALL getarr(fid, TRIM(label)//'/p/iparm', mat%p%iparm)
END SUBROUTINE get_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE get_zpardiso_mat(fid, label, mat)
!
! Read matrix from hdf5 file
!
USE futils
!
INTEGER, INTENT(in) :: fid
CHARACTER(len=*), INTENT(in) :: label
TYPE(zpardiso_mat) :: mat
!
CALL getatt(fid, label, 'RANK', mat%rank)
CALL getatt(fid, label, 'NNZ', mat%nnz)
CALL getatt(fid, label, 'SYM', mat%nlsym)
CALL getatt(fid, label, 'HERM', mat%nlherm)
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)
!
CALL getatt(fid, TRIM(label)//'/p', 'error', mat%p%error)
CALL getatt(fid, TRIM(label)//'/p', 'mtype', mat%p%mtype)
CALL getatt(fid, TRIM(label)//'/p', 'phase', mat%p%phase)
CALL getatt(fid, TRIM(label)//'/p', 'msglv', mat%p%msglvl)
CALL getatt(fid, TRIM(label)//'/p', 'maxfct', mat%p%maxfct)
CALL getatt(fid, TRIM(label)//'/p', 'mnum', mat%p%mnum)
CALL getatt(fid, TRIM(label)//'/p', 'nrhs', mat%p%nrhs)
CALL getarr(fid, TRIM(label)//'/p/iparm', mat%p%iparm)
END SUBROUTINE get_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE mcopy_pardiso_mat(mata, matb)
!
! Matrix copy: B = A
!
TYPE(pardiso_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 pardiso_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)
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))
END SUBROUTINE mcopy_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE mcopy_zpardiso_mat(mata, matb)
!
! Matrix copy: B = A
!
TYPE(zpardiso_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 pardiso_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)
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))
END SUBROUTINE mcopy_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE maddto_pardiso_mat(mata, alpha, matb)
!
! A <- A + alpha*B
!
TYPE(pardiso_mat) :: mata, matb
DOUBLE PRECISION :: alpha
!
mata%val = mata%val + alpha*matb%val
END SUBROUTINE maddto_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE maddto_zpardiso_mat(mata, alpha, matb)
!
! A <- A + alpha*B
!
TYPE(zpardiso_mat) :: mata, matb
DOUBLE COMPLEX :: alpha
!
mata%val = mata%val + alpha*matb%val
END SUBROUTINE maddto_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE psum_pardiso_mat(mat, comm)
!
! Parallel sum of sparse matrices
!
INCLUDE "mpif.h"
!
TYPE(pardiso_mat) :: mat
INCLUDE 'psum_mat.tpl'
END SUBROUTINE psum_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE psum_zpardiso_mat(mat, comm)
!
! Parallel sum of sparse matrices
!
INCLUDE "mpif.h"
!
TYPE(zpardiso_mat) :: mat
INCLUDE 'psum_mat.tpl'
END SUBROUTINE psum_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE p2p_pardiso_mat(mat, dest, extyp, op, comm)
!
! Point-to-point combine sparse matrix between 2 processes
!
INCLUDE "mpif.h"
!
TYPE(pardiso_mat) :: mat
DOUBLE PRECISION, ALLOCATABLE :: val(:)
INTEGER :: mpi_type=MPI_DOUBLE_PRECISION
!
INCLUDE "p2p_mat.tpl"
END SUBROUTINE p2p_pardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE p2p_zpardiso_mat(mat, dest, extyp, op, comm)
!
! Point-to-point combine sparse matrix between 2 processes
!
INCLUDE "mpif.h"
!
TYPE(zpardiso_mat) :: mat
DOUBLE COMPLEX, ALLOCATABLE :: val(:)
INTEGER :: mpi_type=MPI_DOUBLE_COMPLEX
!
INCLUDE "p2p_mat.tpl"
END SUBROUTINE p2p_zpardiso_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
END MODULE pardiso_bsplines
c4science · Help