Homec4science
Diffusion SPClibs (master) Browse

spclibs/bsplines/cds_mod.f90
ca0028807c2emaster

Last Change
Raw File

cds_mod.f90
View Options

!>
!> @file cds_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 cds
!
! CDSMAT: Implement sparse matrix using Compressed
! Diagonal Storage.
!
! T.M. Tran, CRPP-EPFL
! November 2010
!
USE mumps_bsplines
IMPLICIT NONE
!
TYPE cds_mat ! Compressed Diagonal Storage
INTEGER :: rank
INTEGER :: kl, ku, ndiags
INTEGER :: nterms, kmat
INTEGER :: ny
INTEGER, DIMENSION(:), POINTER :: dists => NULL()
DOUBLE PRECISION, DIMENSION(:), POINTER :: rowv => NULL()
DOUBLE PRECISION, DIMENSION(:), POINTER :: colh => NULL()
DOUBLE PRECISION, DIMENSION(:), POINTER :: bal => NULL()
DOUBLE PRECISION, DIMENSION(:,:), POINTER :: val => NULL()
TYPE(mumps_mat), ALLOCATABLE :: mumps
END TYPE cds_mat
!
!--------------------------------------------------------------------------------
INTERFACE init
MODULE PROCEDURE init_cds_mat
END INTERFACE init
INTERFACE clear_mat
MODULE PROCEDURE clear_cds_mat
END INTERFACE clear_mat
INTERFACE destroy
MODULE PROCEDURE destroy_cds_mat
END INTERFACE destroy
INTERFACE updtmat
MODULE PROCEDURE updt_cds
END INTERFACE updtmat
INTERFACE getele
MODULE PROCEDURE getele_cds
END INTERFACE getele
INTERFACE putele
MODULE PROCEDURE putele_cds
END INTERFACE putele
INTERFACE getcol
MODULE PROCEDURE getcol_cds
END INTERFACE getcol
INTERFACE getrow
MODULE PROCEDURE getrow_cds
END INTERFACE getrow
INTERFACE putcol
MODULE PROCEDURE putcol_cds
END INTERFACE putcol
INTERFACE putrow
MODULE PROCEDURE putrow_cds
END INTERFACE putrow
INTERFACE getdiag
MODULE PROCEDURE getdiag_cds
END INTERFACE getdiag
INTERFACE vmx
MODULE PROCEDURE vmx_cds, vmxn_cds
END INTERFACE
INTERFACE putmat
MODULE PROCEDURE putmat_cds
END INTERFACE
INTERFACE getmat
MODULE PROCEDURE getmat_cds
END INTERFACE
INTERFACE flops
MODULE PROCEDURE flops_cds
END INTERFACE flops
INTERFACE matnorm
MODULE PROCEDURE matnorm_cds
END INTERFACE matnorm
!
CONTAINS
!===========================================================================
SUBROUTINE init_cds_mat(rank, dists, nterms, mat, bw0, kmat)
!
! Initialize a CDS matrix obtained for a 2d FE discretization
! using Splines of orders p(1) and p(2).
! Number first the 2nd (periodic) dimension.
!
INTEGER, INTENT(in) :: rank
INTEGER, ALLOCATABLE, INTENT(in) :: dists(:)
INTEGER, INTENT(in) :: nterms
TYPE(cds_mat) :: mat
INTEGER, OPTIONAL, INTENT(in) :: bw0, kmat
!
INTEGER :: kl, ku
!
mat%rank = rank
mat%nterms = nterms
mat%ny = 0 ! Used for unicity condition in cyl. geometry.
IF(PRESENT(kmat)) mat%kmat = kmat
!
kl = -LBOUND(dists,1)
ku = UBOUND(dists,1)
!
mat%kl = kl
mat%ku = ku
mat%ndiags = ku + kl + 1
IF(ASSOCIATED(mat%dists)) DEALLOCATE(mat%dists)
ALLOCATE(mat%dists(-kl:ku))
mat%dists = dists
!
IF(ASSOCIATED(mat%rowv)) DEALLOCATE(mat%rowv)
IF(ASSOCIATED(mat%colh)) DEALLOCATE(mat%colh)
IF(PRESENT(bw0)) THEN
ALLOCATE(mat%rowv(bw0), mat%colh(bw0))
mat%rowv = 0.0d0
mat%colh = 0.0d0
ELSE
ALLOCATE(mat%rowv(0))
ALLOCATE(mat%colh(0))
END IF
!
IF( ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
ALLOCATE(mat%val(rank, -kl:ku))
mat%val = 0.0d0
!
IF(ASSOCIATED(mat%bal)) DEALLOCATE(mat%bal)
ALLOCATE(mat%bal(rank))
mat%bal = 0.0d0
END SUBROUTINE init_cds_mat
!===========================================================================
SUBROUTINE clear_cds_mat(mat)
!
! Clear matrix, keeping its sparse structure unchanged
!
TYPE(cds_mat) :: mat
!
mat%val = 0.0d0
END SUBROUTINE clear_cds_mat
!===========================================================================
SUBROUTINE destroy_cds_mat(mat)
!
! Deallocate pointers in mat
!
TYPE(cds_mat) :: mat
IF( ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
IF( ASSOCIATED(mat%dists)) DEALLOCATE(mat%dists)
IF( ASSOCIATED(mat%rowv)) DEALLOCATE(mat%rowv)
IF( ASSOCIATED(mat%colh)) DEALLOCATE(mat%colh)
END SUBROUTINE destroy_cds_mat
!===========================================================================
SUBROUTINE updt_cds(mat, i, j, val)
!
! Update element Aij into sparse CDS matrix
!
TYPE(cds_mat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE PRECISION, INTENT(in) :: val
INTEGER :: d, k
!
d = j-i
DO k = -mat%kl, mat%ku
IF( d .EQ. mat%dists(k) ) THEN
mat%val(i,k) = mat%val(i,k)+val
RETURN
END IF
END DO
WRITE(*,'(a,2i6)') 'UPDT: i, j out of range ', i, j
WRITE(*,'(a/(10i6))') 'Valid distances', mat%dists
STOP '*** Abnormal EXIT in MODULE matrix ***'
END SUBROUTINE updt_cds
!===========================================================================
SUBROUTINE getele_cds(mat, i, j, val)
!
! Get element Aij of sparse CDS matrix
!
TYPE(cds_mat), INTENT(in) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE PRECISION, INTENT(out) :: val
INTEGER :: d, k
!
d = j-i
DO k = -mat%kl, mat%ku
IF( d .EQ. mat%dists(k) ) THEN
val = mat%val(i,k)
RETURN
END IF
END DO
WRITE(*,'(a,2i6)') 'GETELE: i, j out of range ', i, j
WRITE(*,'(a/(10i6))') 'Valid distances', mat%dists
STOP '*** Abnormal EXIT in MODULE matrix ***'
END SUBROUTINE getele_cds
!===========================================================================
SUBROUTINE putele_cds(mat, i, j, val)
!
! Update element Aij into sparse CDS matrix
!
TYPE(cds_mat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE PRECISION, INTENT(in) :: val
INTEGER :: d, k
!
d = j-i
DO k = -mat%kl, mat%ku
IF( d .EQ. mat%dists(k) ) THEN
mat%val(i,k) = val
RETURN
END IF
END DO
WRITE(*,'(a,2i6)') 'UPDT: i, j out of range ', i, j
WRITE(*,'(a/(10i6))') 'Valid distances', mat%dists
STOP '*** Abnormal EXIT in MODULE matrix ***'
END SUBROUTINE putele_cds
!===========================================================================
SUBROUTINE getcol_cds(mat, j, arr)
!
! Get a column from matrix
!
TYPE(cds_mat), INTENT(in) :: mat
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, INTENT(out) :: arr(:)
INTEGER :: n,i, k
!
n = mat%rank
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'GETCOL: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
arr(1:n) = 0.0d0
DO k=-mat%kl, mat%ku
i = j-mat%dists(k)
IF( i.GE.1 .AND. i.LE.n ) arr(i) = mat%val(i,k)
END DO
END SUBROUTINE getcol_cds
!===========================================================================
SUBROUTINE getrow_cds(mat, i, arr)
!
! Get a row from matrix
!
TYPE(cds_mat), INTENT(in) :: mat
INTEGER, INTENT(in) :: i
DOUBLE PRECISION, INTENT(out) :: arr(:)
INTEGER :: n, j, k
!
n = mat%rank
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'GETROW: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
arr(1:n) = 0.0d0
DO k=-mat%kl, mat%ku
j = i+mat%dists(k)
IF( j.GE.1 .AND. j.LE.n ) arr(j) = mat%val(i,k)
END DO
END SUBROUTINE getrow_cds
!===========================================================================
SUBROUTINE putcol_cds(mat, j, arr)
!
! Put a column to matrix
!
TYPE(cds_mat) :: mat
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, INTENT(in) :: arr(:)
INTEGER :: n,i, k
!
n = mat%rank
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'PUTCOL: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
DO k=-mat%kl, mat%ku
i = j-mat%dists(k)
IF( i.GE. 1 .AND. i.LE.n ) mat%val(i,k) = arr(i)
END DO
END SUBROUTINE putcol_cds
!===========================================================================
SUBROUTINE putrow_cds(mat, i, arr)
!
! Put a row from matrix
!
TYPE(cds_mat) :: mat
INTEGER, INTENT(in) :: i
DOUBLE PRECISION, INTENT(in) :: arr(:)
INTEGER :: n, j, k
!
n = mat%rank
IF( SIZE(arr) .LT. n ) THEN
WRITE(*,*) 'PUTROW: size of arr too small'
STOP '*** Abnormal EXIT in MODULE matrix ***'
END IF
DO k=-mat%kl, mat%ku
j = i+mat%dists(k)
IF( j.GE.1 .AND. j.LE.n ) mat%val(i,k) = arr(j)
END DO
END SUBROUTINE putrow_cds
!===========================================================================
SUBROUTINE getdiag_cds(mat, d)
!
! Returns diagonal of matrix
!
TYPE(cds_mat) :: mat
DOUBLE PRECISION :: d(:)
INTEGER :: ny
!
d(:) = mat%val(:,0)
!
! The extra row and column implied by periodic BC
!!$ ny = mat%ny
!!$ IF( ny .NE. 0 ) THEN
!!$ d(ny) = mat%rowv(ny) + mat%colh(ny)
!!$ END IF
!!$ WRITE(*,'(a/(8(1pe12.3)))') 'd', d
END SUBROUTINE getdiag_cds
!===========================================================================
FUNCTION vmx_cds(mat, xarr)
!
! Return product mat*x
!
TYPE(cds_mat) :: mat
DOUBLE PRECISION, INTENT(in) :: xarr(:)
DOUBLE PRECISION :: vmx_cds(SIZE(xarr))
!
DOUBLE PRECISION :: alpha=1.0d0, beta=1.0d0
INTEGER :: m, bw0, ny, k, d, i, i1, i2
!
m = mat%rank
bw0 = SIZE(mat%rowv)
ny = mat%ny
vmx_cds = 0.0d0
!
IF( ny .NE. 0 ) THEN ! Contributions from unicity BC
vmx_cds(ny:bw0) = mat%colh(ny:bw0)*xarr(ny)
vmx_cds(ny) = vmx_cds(ny) + DOT_PRODUCT(mat%rowv(ny:bw0), xarr(ny:bw0))
END IF
!
#ifdef MKL
CALL mkl_ddiamv('n', m, m, alpha, 'g', mat%val, m, mat%dists, &
& mat%ndiags, xarr, beta, vmx_cds)
#else
DO k=-mat%kl,mat%ku
d = mat%dists(k)
i1 = MAX(1,1-d)
i2 = MIN(mat%rank,mat%rank-d)
DO i=i1,i2
vmx_cds(i) = vmx_cds(i) + mat%val(i,k)*xarr(i+d)
END DO
END DO
#endif
END FUNCTION vmx_cds
!===========================================================================
FUNCTION vmxn_cds(mat, xarr)
!
! Return product mat*x
!
TYPE(cds_mat) :: mat
DOUBLE PRECISION, INTENT(in) :: xarr(:,:)
DOUBLE PRECISION :: vmxn_cds(SIZE(xarr,1),SIZE(xarr,2))
!
DOUBLE PRECISION :: alpha=1.0d0, beta=1.0d0
INTEGER :: m, nrhs, bw0, ny, k, d, i, j, i1, i2
!
m = mat%rank
nrhs = SIZE(xarr,2)
bw0 = SIZE(mat%rowv)
ny = mat%ny
vmxn_cds = 0.0d0
!
IF( ny .NE. 0 ) THEN ! Contributions from unicity BC
DO j=1,nrhs
vmxn_cds(ny:bw0,j) = mat%colh(ny:bw0)*xarr(ny,j)
vmxn_cds(ny,j) = vmxn_cds(ny,j) + &
& DOT_PRODUCT(mat%rowv(ny:bw0), xarr(ny:bw0,j))
END DO
END IF
!
#ifdef MKL
CALL mkl_ddiamm('n', m, nrhs, m, alpha, 'g', mat%val, m, &
& mat%dists, mat%ndiags, xarr, m, beta, vmxn_cds, m)
#else
DO k=-mat%kl,mat%ku
d = mat%dists(k)
i1 = MAX(1,1-d)
i2 = MIN(mat%rank,mat%rank-d)
DO j=1,nrhs
DO i=i1,i2
vmxn_cds(i,j) = vmxn_cds(i,j) + mat%val(i,k)*xarr(i+d,j)
END DO
END DO
END DO
#endif
END FUNCTION vmxn_cds
!===========================================================================
SUBROUTINE getmat_cds(fid, label, mat)
!
! Read in CDS matrix from hdf5 file
!
USE futils
!
INTEGER, INTENT(in) :: fid
CHARACTER(len=*), INTENT(in) :: label
TYPE(cds_mat) :: mat
!
CALL getatt(fid, label, 'RANK', mat%rank)
CALL getatt(fid, label, 'KL', mat%kl)
CALL getatt(fid, label, 'KU', mat%ku)
CALL getatt(fid, label, 'NDIAGS', mat%ndiags)
CALL getatt(fid, label, 'NTERMS', mat%nterms)
CALL getatt(fid, label, 'KMAT', mat%kmat)
CALL getatt(fid, label, 'NY', mat%ny)
IF( ASSOCIATED(mat%dists) ) THEN
CALL getarr(fid, TRIM(label)//'/dists', mat%dists)
END IF
IF(ASSOCIATED(mat%bal)) THEN
CALL getarr(fid, TRIM(label)//'/bal', mat%bal)
END IF
CALL getarr(fid, TRIM(label)//'/vals', mat%val)
IF(ASSOCIATED(mat%rowv)) THEN
CALL getarr(fid, TRIM(label)//'/rowv', mat%rowv)
CALL getarr(fid, TRIM(label)//'/colh', mat%colh)
END IF
END SUBROUTINE getmat_cds
!===========================================================================
SUBROUTINE putmat_cds(fid, label, mat, str)
!
! Write CDS matrix in hdf5 file
!
USE futils
!
INTEGER, INTENT(in) :: fid
CHARACTER(len=*), INTENT(in) :: label
TYPE(cds_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, 'KL', mat%kl)
CALL attach(fid, label, 'KU', mat%ku)
CALL attach(fid, label, 'NDIAGS', mat%ndiags)
CALL attach(fid, label, 'NTERMS', mat%nterms)
CALL attach(fid, label, 'KMAT', mat%kmat)
CALL attach(fid, label, 'NY', mat%ny)
IF( ASSOCIATED(mat%dists) ) THEN
CALL putarr(fid, TRIM(label)//'/dists', mat%dists)
END IF
IF(ASSOCIATED(mat%bal)) THEN
CALL putarr(fid, TRIM(label)//'/bal', mat%bal)
END IF
CALL putarr(fid, TRIM(label)//'/vals', mat%val)
IF(ASSOCIATED(mat%rowv)) THEN
CALL putarr(fid, TRIM(label)//'/rowv', mat%rowv)
CALL putarr(fid, TRIM(label)//'/colh', mat%colh)
END IF
END SUBROUTINE putmat_cds
!===========================================================================
FUNCTION flops_cds(mat, xarr, ny)
!
! Return FLOPS in product mat*x
!
TYPE(cds_mat) :: mat
DOUBLE PRECISION, INTENT(in) :: xarr(:)
DOUBLE PRECISION :: flops_cds
INTEGER, OPTIONAL, INTENT(in) :: ny
!
INTEGER :: k, d, i, i1, i2
!
flops_cds = 0.0d0
IF( PRESENT(ny) ) THEN ! Contributions from unicity BC
flops_cds = 4.0d0*(SIZE(mat%rowv)-ny+1)
END IF
DO k=-mat%kl,mat%ku
d = mat%dists(k)
i1 = MAX(1,1-d)
i2 = MIN(mat%rank,mat%rank-d)
flops_cds = flops_cds + 2.0d0*(i2-i1+1)
END DO
END FUNCTION flops_cds
!===========================================================================
SUBROUTINE cds2mumps(mat, mat_mumps)
!
! Fill mumps structure (based on routine to_mumps_mat)
!
INCLUDE 'mpif.h'
TYPE(cds_mat) :: mat
TYPE(mumps_mat) :: mat_mumps
!
INTEGER :: i, ii, i1, i2, j, k, rank, d, bw0, s, e
INTEGER :: comm, ierr, nnz_loc
!
CALL init(mat%rank, mat%nterms, mat_mumps)
!
comm = mat_mumps%mumps_par%COMM
mat_mumps%mumps_par%ICNTL(18) = 3 ! Distributed assembled matrix
!
! Compute nnz_loc
!
rank = mat_mumps%rank
s = mat_mumps%istart
e = mat_mumps%iend
!
nnz_loc=0
DO k=-mat%kl,mat%ku
d = mat%dists(k)
i1 = MAX(s,1-d)
i2 = MIN(e,rank-d)
nnz_loc = nnz_loc + (i2-i1+1)
END DO
!
! Extra col and row from unicity conditions
!
bw0 = SIZE(mat%rowv)
IF(bw0.GT.0) THEN
IF(mat%ny.GE.s .AND. mat%ny.LE.e) THEN
nnz_loc = nnz_loc + bw0-mat%ny ! rowh(ny+1:bw0)
END IF
nnz_loc = nnz_loc + (MIN(bw0,e)-MAX(mat%ny,s)) ! colh(ny+1:bw0)
END IF
!
mat_mumps%nnz_start = 0
CALL mpi_exscan(nnz_loc, mat_mumps%nnz_start, 1, MPI_INTEGER, MPI_SUM, comm, ierr)
mat_mumps%nnz_start = mat_mumps%nnz_start + 1
mat_mumps%nnz_end = mat_mumps%nnz_start + nnz_loc - 1
mat_mumps%nnz_loc = nnz_loc
CALL mpi_allreduce(nnz_loc, mat_mumps%nnz, 1, MPI_INTEGER, MPI_SUM, comm, ierr)
!
mat_mumps%mumps_par%N = rank
mat_mumps%mumps_par%NZ_loc = nnz_loc
!
! Construct MUMPS (IRN, JCN, A)
!
ALLOCATE(mat_mumps%mumps_par%IRN_loc(nnz_loc))
ALLOCATE(mat_mumps%mumps_par%JCN_loc(nnz_loc))
ALLOCATE(mat_mumps%mumps_par%A_loc(nnz_loc))
!
ii=0
DO k=-mat%kl,mat%ku
d = mat%dists(k)
i1 = MAX(s,1-d)
i2 = MIN(e,rank-d)
DO i=i1,i2
ii = ii+1
mat_mumps%mumps_par%IRN_loc(ii) = i
mat_mumps%mumps_par%JCN_loc(ii) = i+d
mat_mumps%mumps_par%A_loc(ii) = mat%val(i,k)
END DO
END DO
!
IF(bw0.GT.0) THEN
IF(mat%ny.GE.s .AND. mat%ny.LE.e) THEN
DO j=mat%ny+1,bw0 ! rowh(ny+1:bw0)
ii = ii+1
mat_mumps%mumps_par%IRN_loc(ii) = mat%ny
mat_mumps%mumps_par%JCN_loc(ii) = j
mat_mumps%mumps_par%A_loc(ii) = mat%rowv(j)
END DO
END IF
DO i=MAX(mat%ny,s)+1,MIN(bw0,e) ! colh(ny+1:bw0)
ii = ii+1
mat_mumps%mumps_par%IRN_loc(ii) = i
mat_mumps%mumps_par%JCN_loc(ii) = mat%ny
mat_mumps%mumps_par%A_loc(ii) = mat%colh(i)
END DO
END IF
!
CALL destroy(mat_mumps%mat)
NULLIFY(mat_mumps%mat)
END SUBROUTINE cds2mumps
!===========================================================================
DOUBLE PRECISION FUNCTION matnorm_cds(mat, p)
!
! Compute matrix norm
!
TYPE(cds_mat), INTENT(in) :: mat
CHARACTER(len=*), OPTIONAL, INTENT(in) :: p
!
CHARACTER(len=4) :: norm_type
INTEGER :: i, j, k, d
DOUBLE PRECISION :: temp(mat%rank)
!
norm_type = 'fro'
IF(PRESENT(p)) norm_type = p
!
SELECT CASE (norm_type)
CASE ('inf')
DO i=1,mat%rank
temp(i) = SUM(ABS(mat%val(i,:)))
END DO
matnorm_cds = MAXVAL(temp)
CASE ('1')
temp = 0.0d0
DO k=-mat%kl,mat%ku
d = mat%dists(k)
DO i=MAX(1,1-d),MIN(mat%rank,mat%rank-d)
temp(i+d) = temp(i+d) + ABS(mat%val(i,k))
END DO
END DO
matnorm_cds = MAXVAL(temp)
CASE('fro')
matnorm_cds = SQRT(SUM(mat%val**2))
END SELECT
END FUNCTION matnorm_cds
!===========================================================================
!
END MODULE cds
c4science · Help