Homec4science
Diffusion SPClibs (master) Browse

spclibs/bsplines/csr_mod.f90
1a965ce17d3bmaster

Last Change
Raw File

csr_mod.f90
View Options

!>
!> @file csr_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 csr
!
! CSR: Implement CSR (Compressed Sparse Row) matrice
!
! T.M. Tran, CRPP-EPFL
! October 2012
!
USE sparse
USE mumps_bsplines
IMPLICIT NONE
!
TYPE, EXTENDS(spmat) :: csr_mat
INTEGER :: mrows, ncols
INTEGER :: nnz = 0 ! Number of non-zeros
INTEGER :: nterms ! Number of terms in weak form
LOGICAL :: nlforce_zero ! Keep exixting nodes with zero value if .true.
INTEGER, POINTER :: irow(:) => NULL() ! points to start of row
INTEGER, POINTER :: idiag(:) => NULL() ! points to diagonal element
INTEGER, POINTER :: cols(:) => NULL() ! Column indices
DOUBLE PRECISION, POINTER :: val(:) => NULL() ! Elelement values
TYPE(mumps_mat), ALLOCATABLE :: mumps
END TYPE csr_mat
!
TYPE, EXTENDS(zspmat) :: zcsr_mat
INTEGER :: mrows, ncols
INTEGER :: nnz = 0 ! Number of non-zeros
INTEGER :: nterms ! Number of terms in weak form
LOGICAL :: nlforce_zero ! Keep exixting nodes with zero value if .true.
INTEGER, POINTER :: irow(:) => NULL() ! points to start of row
INTEGER, POINTER :: idiag(:) => NULL() ! points to diagonal element
INTEGER, POINTER :: cols(:) => NULL() ! Column indices
DOUBLE COMPLEX, POINTER :: val(:) => NULL() ! Elelement values
!
TYPE(zmumps_mat), ALLOCATABLE :: mumps
END TYPE zcsr_mat
!
INTERFACE init
MODULE PROCEDURE init_csr_mat, init_zcsr_mat
END INTERFACE init
INTERFACE clear_mat
MODULE PROCEDURE clear_csr_mat, clear_zcsr_mat
END INTERFACE clear_mat
INTERFACE updtmat
MODULE PROCEDURE updt_csr_mat, updt_zcsr_mat
END INTERFACE updtmat
INTERFACE putele
MODULE PROCEDURE putele_csr_mat, putele_zcsr_mat
END INTERFACE putele
INTERFACE getele
MODULE PROCEDURE getele_csr_mat, getele_zcsr_mat
END INTERFACE getele
INTERFACE putrow
MODULE PROCEDURE putrow_csr_mat, putrow_zcsr_mat
END INTERFACE putrow
INTERFACE getrow
MODULE PROCEDURE getrow_csr_mat, getrow_zcsr_mat
END INTERFACE getrow
INTERFACE getdiag
MODULE PROCEDURE getdiag_csr_mat, getdiag_zcsr_mat
END INTERFACE getdiag
INTERFACE putcol
MODULE PROCEDURE putcol_csr_mat, putcol_zcsr_mat
END INTERFACE putcol
INTERFACE getcol
MODULE PROCEDURE getcol_csr_mat, getcol_zcsr_mat
END INTERFACE getcol
INTERFACE to_mat
MODULE PROCEDURE to_csr_mat, to_zcsr_mat
END INTERFACE to_mat
INTERFACE vmx
MODULE PROCEDURE vmx_csr_mat, vmx_csr_matn, vmx_zcsr_mat, vmx_zcsr_matn
END INTERFACE vmx
INTERFACE destroy
MODULE PROCEDURE destroy_csr_mat, destroy_zcsr_mat
END INTERFACE destroy
INTERFACE putmat
MODULE PROCEDURE put_csr_mat, put_zcsr_mat
END INTERFACE putmat
!>>>>>
!>>>>> CONMAT
!>>>>
INTERFACE conmat
MODULE PROCEDURE conmat_1d_csr, conmat_2d_csr, conmat_1d_zcsr, conmat_2d_zcsr
END INTERFACE conmat
!>>>>
!>>>> MULTIGRID_MOD
!>>>>
INTERFACE femat
MODULE PROCEDURE femat_csr
END INTERFACE femat
INTERFACE matnorm
MODULE PROCEDURE matnorm_csr
END INTERFACE matnorm
INTERFACE kron
MODULE PROCEDURE kron_csr
END INTERFACE kron
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
CONTAINS
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE init_csr_mat(n, nterms, mat, nlforce_zero, ncols)
!
! Initialize an empty CSR matrix
!
INTEGER, INTENT(in) :: n, nterms
TYPE(csr_mat) :: mat
LOGICAL, OPTIONAL, INTENT(in) :: nlforce_zero
INTEGER, OPTIONAL, intent(in) :: ncols
!
CALL init(n, mat%spmat)
mat%mrows = n
mat%ncols = n
IF(PRESENT(ncols)) mat%ncols = ncols
mat%nterms = nterms
mat%nlforce_zero = .TRUE. ! Do not remove existing nodes with zero val
IF(PRESENT(nlforce_zero)) mat%nlforce_zero = nlforce_zero
!
END SUBROUTINE init_csr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE init_zcsr_mat(n, nterms, mat, nlforce_zero, ncols)
!
! Initialize an empty CSR matrix
!
INTEGER, INTENT(in) :: n, nterms
TYPE(zcsr_mat) :: mat
LOGICAL, OPTIONAL, INTENT(in) :: nlforce_zero
INTEGER, OPTIONAL, intent(in) :: ncols
!
CALL init(n, mat%zspmat)
mat%mrows = n
mat%ncols = n
IF(PRESENT(ncols)) mat%ncols = ncols
mat%nterms = nterms
mat%nlforce_zero = .TRUE. ! Do not remove existing nodes with zero val
IF(PRESENT(nlforce_zero)) mat%nlforce_zero = nlforce_zero
!
END SUBROUTINE init_zcsr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE clear_zcsr_mat(mat)
!
! Clear matrix, keeping its sparse structure unchanged
!
TYPE(zcsr_mat) :: mat
!
mat%val = (0.0d0,0.0d0)
END SUBROUTINE clear_zcsr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE clear_csr_mat(mat)
!
! Clear matrix, keeping its sparse structure unchanged
!
TYPE(csr_mat) :: mat
!
mat%val = 0.0d0
END SUBROUTINE clear_csr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE updt_csr_mat(mat, i, j, val)
!
! Update element Aij of csr matrix
!
TYPE(csr_mat) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE PRECISION, INTENT(in) :: val
INTEGER :: k, s, e
!
IF(mat%nnz.EQ.0) THEN ! Still using linked lists
CALL updtmat(mat%spmat, 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 csr_mod ***'
END IF
END IF
END SUBROUTINE updt_csr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE updt_zcsr_mat(mat, i, j, val)
!
! Update element Aij of csr matrix
!
TYPE(zcsr_mat) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE COMPLEX, INTENT(in) :: val
INTEGER :: k, s, e
!
IF(mat%nnz.EQ.0) THEN ! Still using linked lists
CALL updtmat(mat%zspmat, 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 csr_mod ***'
END IF
END IF
END SUBROUTINE updt_zcsr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putele_csr_mat(mat, i, j, val)
!
! Put element (i,j) of matrix
!
TYPE(csr_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%nnz.EQ.0) THEN ! Still using linked lists
CALL putele(mat%spmat, 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
PRINT*, 'val', val
PRINT*, 'matrix m, n', mat%mrows, mat%ncols
STOP '*** Abnormal EXIT in MODULE csr_mod ***'
END IF
END IF
END IF
END SUBROUTINE putele_csr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putele_zcsr_mat(mat, i, j, val)
!
! Put element (i,j) of matrix
!
TYPE(zcsr_mat) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE COMPLEX, INTENT(in) :: val
INTEGER :: iput, jput
INTEGER :: k, s, e
!
iput = i
jput = j
!
IF(mat%nnz.EQ.0) THEN ! Still using linked lists
CALL putele(mat%zspmat, 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
PRINT*, 'val', val
PRINT*, 'matrix m, n', mat%mrows, mat%ncols
STOP '*** Abnormal EXIT in MODULE csr_mod ***'
END IF
END IF
END IF
END SUBROUTINE putele_zcsr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getele_csr_mat(mat, i, j, val)
!
! Get element (i,j) of sparse matrix
!
TYPE(csr_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%nnz.EQ.0) THEN ! Still using linked lists
CALL getele(mat%spmat, 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_csr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getele_zcsr_mat(mat, i, j, val)
!
! Get element (i,j) of sparse matrix
!
TYPE(zcsr_mat) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE COMPLEX, INTENT(out) :: val
INTEGER :: iget, jget
INTEGER :: k, s, e
!
iget = i
jget = j
IF(mat%nnz.EQ.0) THEN ! Still using linked lists
CALL getele(mat%zspmat, 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_zcsr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putrow_csr_mat(amat, i, arr)
!
! Put a row into sparse matrix
!
TYPE(csr_mat), INTENT(inout) :: amat
INTEGER, INTENT(in) :: i
DOUBLE PRECISION, INTENT(in) :: arr(:)
INTEGER :: s, e, j
!
IF(amat%nnz.EQ.0) THEN ! Still using linked lists
DO j=1,amat%ncols
CALL putele(amat, i, j, arr(j))
END DO
ELSE
s = amat%irow(i)
e = amat%irow(i+1)-1
DO j=s,e
amat%val(j) = arr(amat%cols(j))
END DO
END IF
END SUBROUTINE putrow_csr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putrow_zcsr_mat(amat, i, arr)
!
! Put a row into sparse matrix
!
TYPE(zcsr_mat), INTENT(inout) :: amat
INTEGER, INTENT(in) :: i
DOUBLE COMPLEX, INTENT(in) :: arr(:)
INTEGER :: s, e, j
!
IF(amat%nnz.EQ.0) THEN ! Still using linked lists
DO j=1,amat%ncols
CALL putele(amat, i, j, arr(j))
END DO
ELSE
s = amat%irow(i)
e = amat%irow(i+1)-1
DO j=s,e
amat%val(j) = arr(amat%cols(j))
END DO
END IF
END SUBROUTINE putrow_zcsr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getrow_csr_mat(amat, i, arr)
!
! Get a row from sparse matrix
!
TYPE(csr_mat), INTENT(in) :: amat
INTEGER, INTENT(in) :: i
DOUBLE PRECISION, INTENT(out) :: arr(:)
INTEGER :: s, e, j
!
arr = 0.0d0
IF(amat%nnz.EQ.0) THEN ! Still using linked lists
DO j=1,amat%ncols
CALL getele(amat, i, j, arr(j))
END DO
ELSE
s = amat%irow(i)
e = amat%irow(i+1)-1
DO j=s,e
arr(amat%cols(j)) = amat%val(j)
END DO
END IF
END SUBROUTINE getrow_csr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getrow_zcsr_mat(amat, i, arr)
!
! Get a row from sparse matrix
!
TYPE(zcsr_mat), INTENT(in) :: amat
INTEGER, INTENT(in) :: i
DOUBLE COMPLEX, INTENT(out) :: arr(:)
INTEGER :: s, e, j
!
arr = 0.0d0
IF(amat%nnz.EQ.0) THEN ! Still using linked lists
DO j=1,amat%ncols
CALL getele(amat, i, j, arr(j))
END DO
ELSE
s = amat%irow(i)
e = amat%irow(i+1)-1
DO j=s,e
arr(amat%cols(j)) = amat%val(j)
END DO
END IF
END SUBROUTINE getrow_zcsr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getdiag_csr_mat(amat, arr)
!
! Get the diagonal from matrix
!
TYPE(csr_mat), INTENT(in) :: amat
DOUBLE PRECISION, INTENT(out) :: arr(:)
!
! WARNING: assume that CSR matrix has been converted from linked lists
!
arr(:) = amat%val(amat%idiag(:))
END SUBROUTINE getdiag_csr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getdiag_zcsr_mat(amat, arr)
!
! Get the diagonal from matrix
!
TYPE(zcsr_mat), INTENT(in) :: amat
DOUBLE COMPLEX, INTENT(out) :: arr(:)
!
! WARNING: assume that CSR matrix has been converted from linked lists
!
arr(:) = amat%val(amat%idiag(:))
END SUBROUTINE getdiag_zcsr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putcol_csr_mat(amat, j, arr)
!
! Put a column into sparse matrix
!
TYPE(csr_mat), INTENT(inout) :: amat
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, INTENT(in) :: arr(:)
INTEGER :: i
!
DO i=1,amat%mrows
CALL putele(amat, i, j, arr(i))
END DO
END SUBROUTINE putcol_csr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putcol_zcsr_mat(amat, j, arr)
!
! Put a column into sparse matrix
!
TYPE(zcsr_mat), INTENT(inout) :: amat
INTEGER, INTENT(in) :: j
DOUBLE COMPLEX, INTENT(in) :: arr(:)
INTEGER :: i
!
DO i=1,amat%mrows
CALL putele(amat, i, j, arr(i))
END DO
END SUBROUTINE putcol_zcsr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getcol_csr_mat(amat, j, arr)
!
! Get a column from sparse matrix
!
TYPE(csr_mat), INTENT(in) :: amat
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, INTENT(out) :: arr(:)
INTEGER :: i
!
DO i=1,amat%mrows
CALL getele(amat, i, j, arr(i))
END DO
END SUBROUTINE getcol_csr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getcol_zcsr_mat(amat, j, arr)
!
! Get a column from sparse matrix
!
TYPE(zcsr_mat), INTENT(in) :: amat
INTEGER, INTENT(in) :: j
DOUBLE COMPLEX, INTENT(out) :: arr(:)
INTEGER :: i
!
DO i=1,amat%mrows
CALL getele(amat, i, j, arr(i))
END DO
END SUBROUTINE getcol_zcsr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE to_csr_mat(mat, nlkeep)
!
! Convert linked list spmat to csr matrice structure
!
TYPE(csr_mat) :: mat
LOGICAL, INTENT(in), OPTIONAL :: nlkeep
INTEGER :: nnz_arr(mat%rank)
INTEGER :: i, nnz, rank, s, e
LOGICAL :: nlclean
!
nlclean = .TRUE.
IF(PRESENT(nlkeep)) THEN
nlclean = .NOT. nlkeep
END IF
!
! Allocate the csr matrix structure
!
nnz = get_count(mat%spmat, nnz_arr)
rank = mat%rank
mat%nnz = nnz
IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
IF(ASSOCIATED(mat%idiag)) DEALLOCATE(mat%idiag)
IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
ALLOCATE(mat%irow(rank+1))
ALLOCATE(mat%idiag(rank))
ALLOCATE(mat%cols(nnz))
ALLOCATE(mat%val(nnz))
!
! Fill csr structure and optionally deallocate the sparse rows
!
mat%irow = 1
DO i=1,rank
mat%irow(i+1) = mat%irow(i) + nnz_arr(i)
s = mat%irow(i)
e = mat%irow(i+1)-1
CALL getrow(mat%spmat%row(i), mat%val(s:e), mat%cols(s:e))
mat%idiag(i) = isearch(mat%cols(s:e), i) + s
IF(nlclean) CALL destroy(mat%spmat%row(i))
END DO
!!$!
!!$! MUMPS mat for direct solver
!!$!
!!$ ALLOCATE(mat%mumps)
!!$ CALL csr2mumps(mat, mat%mumps)
END SUBROUTINE to_csr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE to_zcsr_mat(mat, nlkeep)
!
! Convert linked list spmat to csr matrice structure
!
TYPE(zcsr_mat) :: mat
LOGICAL, INTENT(in), OPTIONAL :: nlkeep
INTEGER :: nnz_arr(mat%rank)
INTEGER :: i, nnz, rank, s, e
LOGICAL :: nlclean
!
nlclean = .TRUE.
IF(PRESENT(nlkeep)) THEN
nlclean = .NOT. nlkeep
END IF
!
! Allocate the csr matrix structure
!
nnz = get_count(mat%zspmat, nnz_arr)
rank = mat%rank
mat%nnz = nnz
IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
IF(ASSOCIATED(mat%idiag)) DEALLOCATE(mat%idiag)
IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
ALLOCATE(mat%irow(rank+1))
ALLOCATE(mat%idiag(rank))
ALLOCATE(mat%cols(nnz))
ALLOCATE(mat%val(nnz))
!
! Fill csr structure and optionally deallocate the sparse rows
!
mat%irow = 1
DO i=1,rank
mat%irow(i+1) = mat%irow(i) + nnz_arr(i)
s = mat%irow(i)
e = mat%irow(i+1)-1
CALL getrow(mat%zspmat%row(i), mat%val(s:e), mat%cols(s:e))
mat%idiag(i) = isearch(mat%cols(s:e), i) + s
IF(nlclean) CALL destroy(mat%zspmat%row(i))
END DO
!!$!
!!$! MUMPS mat for direct solver
!!$!
!!$ ALLOCATE(mat%mumps)
!!$ CALL csr2mumps(mat, mat%mumps)
END SUBROUTINE to_zcsr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE csr2mumps(mat, mat_mumps)
!
! Fill mumps structure (based on routine to_mumps_mat)
!
INCLUDE 'mpif.h'
TYPE(csr_mat) :: mat
TYPE(mumps_mat) :: mat_mumps
!
INTEGER :: i, rank, 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
!
! 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_mumps%rank
nnz_loc = mat%nnz
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
!
mat_mumps%cols => mat%cols
mat_mumps%irow => mat%irow
mat_mumps%val => mat%val
!
! (A,JCN) picked from CSR mat
mat_mumps%mumps_par%A_loc => mat_mumps%val
mat_mumps%mumps_par%JCN_loc => mat_mumps%cols
!
! Determine IRN array
IF(ASSOCIATED(mat_mumps%mumps_par%IRN_loc)) DEALLOCATE(mat_mumps%mumps_par%IRN_loc)
ALLOCATE(mat_mumps%mumps_par%IRN_loc(nnz_loc))
DO i=mat_mumps%istart,mat_mumps%iend
s = mat_mumps%irow(i) - mat_mumps%nnz_start + 1
e = mat_mumps%irow(i+1) - mat_mumps%nnz_start
mat_mumps%mumps_par%IRN_loc(s:e) = i
END DO
CALL destroy(mat_mumps%mat)
NULLIFY(mat_mumps%mat)
!
END SUBROUTINE csr2mumps
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE zcsr2mumps(mat, mat_mumps)
!
! Fill mumps structure (based on routine to_mumps_mat)
!
INCLUDE 'mpif.h'
TYPE(zcsr_mat) :: mat
TYPE(zmumps_mat) :: mat_mumps
!
INTEGER :: i, rank, 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
!
! 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_mumps%rank
nnz_loc = mat%nnz
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
!
mat_mumps%cols => mat%cols
mat_mumps%irow => mat%irow
mat_mumps%val => mat%val
!
! (A,JCN) picked from CSR mat
mat_mumps%mumps_par%A_loc => mat_mumps%val
mat_mumps%mumps_par%JCN_loc => mat_mumps%cols
!
! Determine IRN array
IF(ASSOCIATED(mat_mumps%mumps_par%IRN_loc)) DEALLOCATE(mat_mumps%mumps_par%IRN_loc)
ALLOCATE(mat_mumps%mumps_par%IRN_loc(nnz_loc))
DO i=mat_mumps%istart,mat_mumps%iend
s = mat_mumps%irow(i) - mat_mumps%nnz_start + 1
e = mat_mumps%irow(i+1) - mat_mumps%nnz_start
mat_mumps%mumps_par%IRN_loc(s:e) = i
END DO
CALL destroy(mat_mumps%mat)
NULLIFY(mat_mumps%mat)
!
END SUBROUTINE zcsr2mumps
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE destroy_csr_mat(mat)
!
! Deallocate csr mat
!
TYPE(csr_mat) :: mat
!
CALL destroy(mat%spmat)
IF(mat%nnz.GT.0) THEN
DEALLOCATE(mat%irow)
DEALLOCATE(mat%idiag)
DEALLOCATE(mat%cols)
DEALLOCATE(mat%val)
END IF
mat%nnz = 0
END SUBROUTINE destroy_csr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE destroy_zcsr_mat(mat)
!
! Deallocate csr mat
!
TYPE(zcsr_mat) :: mat
!
CALL destroy(mat%zspmat)
IF(mat%nnz.GT.0) THEN
DEALLOCATE(mat%irow)
DEALLOCATE(mat%idiag)
DEALLOCATE(mat%cols)
DEALLOCATE(mat%val)
END IF
mat%nnz = 0
END SUBROUTINE destroy_zcsr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION vmx_csr_mat(mat, xarr, transa_in) RESULT(yarr)
!
! Return product mat*x
!
TYPE(csr_mat) :: mat
DOUBLE PRECISION, INTENT(in) :: xarr(:)
CHARACTER(len=1), OPTIONAL, INTENT(in) :: transa_in
DOUBLE PRECISION :: yarr(SIZE(xarr))
!
DOUBLE PRECISION :: alpha=1.0d0, beta=0.0d0
CHARACTER(len=1) :: transa
CHARACTER(len=6) :: matdescra
INTEGER :: n, i, j
!
n = mat%rank
transa = 'N'
IF(PRESENT(transa_in)) transa = transa_in
!
#ifdef MKL
matdescra = 'g'
CALL mkl_dcsrmv(transa, n, n, alpha, matdescra, mat%val, &
& mat%cols, mat%irow(1), mat%irow(2), xarr, &
& beta, yarr)
#else
yarr = 0.0d0
DO i=1,mat%rank
DO j=mat%irow(i), mat%irow(i+1)-1
IF(transa .EQ. 'N') THEN
yarr(i) = yarr(i) + mat%val(j)*xarr(mat%cols(j))
ELSE
yarr(mat%cols(j)) = yarr(mat%cols(j)) + mat%val(j)*xarr(i)
END IF
END DO
END DO
#endif
!
END FUNCTION vmx_csr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION vmx_zcsr_mat(mat, xarr, transa_in) RESULT(yarr)
!
! Return product mat*x
!
TYPE(zcsr_mat) :: mat
DOUBLE COMPLEX, INTENT(in) :: xarr(:)
CHARACTER(len=1), OPTIONAL, INTENT(in) :: transa_in
DOUBLE COMPLEX :: yarr(SIZE(xarr))
!
DOUBLE COMPLEX :: alpha=(1.0d0,0.0d0), beta=(0.0d0,0.0d0)
CHARACTER(len=1) :: transa
CHARACTER(len=6) :: matdescra
INTEGER :: n, i, j
!
n = mat%rank
transa = 'N'
IF(PRESENT(transa_in)) transa = transa_in
!
#ifdef MKL
matdescra = 'g'
CALL mkl_zcsrmv(transa, n, n, alpha, matdescra, mat%val, &
& mat%cols, mat%irow(1), mat%irow(2), xarr, &
& beta, yarr)
#else
yarr = 0.0d0
DO i=1,mat%rank
DO j=mat%irow(i), mat%irow(i+1)-1
IF(transa .EQ. 'N') THEN
yarr(i) = yarr(i) + mat%val(j)*xarr(mat%cols(j))
ELSE
yarr(mat%cols(j)) = yarr(mat%cols(j)) + mat%val(j)*xarr(i)
END IF
END DO
END DO
#endif
!
END FUNCTION vmx_zcsr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION vmx_csr_matn(mat, xarr, transa_in) RESULT(yarr)
!
! Return product mat*x
!
TYPE(csr_mat) :: mat
DOUBLE PRECISION, INTENT(in) :: xarr(:,:)
CHARACTER(len=1), OPTIONAL, INTENT(in) :: transa_in
DOUBLE PRECISION :: yarr(SIZE(xarr,1),SIZE(xarr,2))
!
DOUBLE PRECISION :: alpha=1.0d0, beta=0.0d0
CHARACTER(len=1) :: transa
INTEGER :: n, nrhs, i, j
CHARACTER(len=6) :: matdescra
!
n = mat%rank
nrhs = SIZE(xarr,2)
transa = 'N'
IF(PRESENT(transa_in)) transa = transa_in
!
#ifdef MKL
matdescra = 'g'
CALL mkl_dcsrmm(transa, n, nrhs, n, alpha, matdescra, mat%val,&
& mat%cols, mat%irow(1), mat%irow(2), xarr, &
& n, beta, yarr, n)
#else
yarr = 0.0d0
DO i=1,n
DO j=mat%irow(i), mat%irow(i+1)-1
IF(transa .EQ. 'N') THEN
yarr(i,:) = yarr(i,:) + mat%val(j)*xarr(mat%cols(j),:)
ELSE
yarr(mat%cols(j),:) = yarr(mat%cols(j),:) + mat%val(j)*xarr(i,:)
END IF
END DO
END DO
#endif
!
END FUNCTION vmx_csr_matn
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION vmx_zcsr_matn(mat, xarr, transa_in) RESULT(yarr)
!
! Return product mat*x
!
TYPE(zcsr_mat) :: mat
DOUBLE COMPLEX, INTENT(in) :: xarr(:,:)
CHARACTER(len=1), OPTIONAL, INTENT(in) :: transa_in
DOUBLE COMPLEX :: yarr(SIZE(xarr,1),SIZE(xarr,2))
!
DOUBLE COMPLEX :: alpha=(1.0d0,0.0d0), beta=(0.0d0,0.0d0)
CHARACTER(len=1) :: transa
INTEGER :: n, nrhs, i, j
CHARACTER(len=6) :: matdescra
!
n = mat%rank
nrhs = SIZE(xarr,2)
transa = 'N'
IF(PRESENT(transa_in)) transa = transa_in
!
#ifdef MKL
matdescra = 'g'
CALL mkl_zcsrmm(transa, n, nrhs, n, alpha, matdescra, mat%val,&
& mat%cols, mat%irow(1), mat%irow(2), xarr, &
& n, beta, yarr, n)
#else
yarr = (0.0d0,0.0d0)
DO i=1,n
DO j=mat%irow(i), mat%irow(i+1)-1
IF(transa .EQ. 'N') THEN
yarr(i,:) = yarr(i,:) + mat%val(j)*xarr(mat%cols(j),:)
ELSE
yarr(mat%cols(j),:) = yarr(mat%cols(j),:) + mat%val(j)*xarr(i,:)
END IF
END DO
END DO
#endif
!
END FUNCTION vmx_zcsr_matn
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE conmat_1d_csr(spl, mat, coefeq, maxder)
!
! Construction of FE matrix mat for 1D differential operator
! using spline spl
!
USE bsplines
TYPE(csr_mat) :: mat
TYPE(spline1d), INTENT(in) :: spl
!
INCLUDE 'conmat_1d.tpl'
END SUBROUTINE conmat_1d_csr
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE conmat_1d_zcsr(spl, mat, coefeq, maxder)
!
! Construction of FE matrix mat for 1D differential operator
! using spline spl
!
USE bsplines
TYPE(zcsr_mat) :: mat
TYPE(spline1d), INTENT(in) :: spl
!
INCLUDE 'zconmat_1d.tpl'
END SUBROUTINE conmat_1d_zcsr
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE conmat_2d_csr(spl, mat, coefeq, maxder, nat_order)
!
! Construction of FE matrix mat for 2D differential operator
! using spline spl
!
USE bsplines
TYPE(spline2d), INTENT(in) :: spl
TYPE(csr_mat) :: mat
!
INCLUDE 'conmat.tpl'
END SUBROUTINE conmat_2d_csr
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE conmat_2d_zcsr(spl, mat, coefeq, maxder, nat_order)
!
! Construction of FE matrix mat for 2D differential operator
! using spline spl
!
USE bsplines
TYPE(spline2d), INTENT(in) :: spl
TYPE(zcsr_mat) :: mat
!
INCLUDE 'zconmat.tpl'
END SUBROUTINE conmat_2d_zcsr
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE femat_csr(spl, mat, coefeq, nterms)
!
! Compute fe matrix
!
USE bsplines
TYPE(spline1d), INTENT(in) :: spl
TYPE(csr_mat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: nterms
INTERFACE
SUBROUTINE coefeq(x, idt, idw, c)
DOUBLE PRECISION, INTENT(in) :: x
INTEGER, INTENT(out) :: idt(:), idw(:)
DOUBLE PRECISION, INTENT(out) :: c(:)
END SUBROUTINE coefeq
END INTERFACE
!
INTEGER :: nrank, nx, nidbas
!
CALL get_dim(spl, nrank, nx, nidbas)
IF(spl%period) nrank = nx
IF(mat%nnz.EQ.0) THEN
WRITE(*,'(a,i0,a)') 'FEMAT: Initialize mat with ', &
& nterms, ' terms ...'
CALL init(nrank, nterms, mat)
END IF
CALL conmat(spl, mat, coefeq)
END SUBROUTINE femat_csr
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE put_csr_mat(fid, label, mat, str)
!
! Write matrix to hdf5 file
!
USE futils
!
INTEGER, INTENT(in) :: fid
CHARACTER(len=*), INTENT(in) :: label
TYPE(csr_mat) :: mat
CHARACTER(len=*), OPTIONAL, INTENT(in) :: str
!
IF(PRESENT(str)) THEN
CALL creatg(fid, label, str)
ELSE
CALL creatg(fid, label)
END IF
CALL attach(fid, label, 'RANK', mat%rank)
CALL attach(fid, label, 'NNZ', mat%nnz)
CALL putarr(fid, TRIM(label)//'/irow', mat%irow)
CALL putarr(fid, TRIM(label)//'/cols', mat%cols)
CALL putarr(fid, TRIM(label)//'/idiag', mat%idiag)
CALL putarr(fid, TRIM(label)//'/val', mat%val)
END SUBROUTINE put_csr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE put_zcsr_mat(fid, label, mat, str)
!
! Write matrix to hdf5 file
!
USE futils
!
INTEGER, INTENT(in) :: fid
CHARACTER(len=*), INTENT(in) :: label
TYPE(zcsr_mat) :: mat
CHARACTER(len=*), OPTIONAL, INTENT(in) :: str
!
IF(PRESENT(str)) THEN
CALL creatg(fid, label, str)
ELSE
CALL creatg(fid, label)
END IF
CALL attach(fid, label, 'RANK', mat%rank)
CALL attach(fid, label, 'NNZ', mat%nnz)
CALL putarr(fid, TRIM(label)//'/irow', mat%irow)
CALL putarr(fid, TRIM(label)//'/cols', mat%cols)
CALL putarr(fid, TRIM(label)//'/idiag', mat%idiag)
CALL putarr(fid, TRIM(label)//'/val', mat%val)
END SUBROUTINE put_zcsr_mat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
DOUBLE PRECISION FUNCTION matnorm_csr(mat, p)
!
! Compute matrix norm
!
TYPE(csr_mat), INTENT(in) :: mat
CHARACTER(len=*), OPTIONAL, INTENT(in) :: p
!
CHARACTER(len=4) :: norm_type
INTEGER :: i, j
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(mat%irow(i):mat%irow(i+1)-1)))
END DO
matnorm_csr = MAXVAL(temp)
CASE ('1')
temp = 0.0d0
DO i=1,mat%rank
DO j=mat%irow(i), mat%irow(i+1)-1
temp(mat%cols(j)) = temp(mat%cols(j)) + ABS(mat%val(j))
END DO
END DO
matnorm_csr = MAXVAL(temp)
CASE('fro')
matnorm_csr = SQRT(DOT_PRODUCT(mat%val, mat%val))
END SELECT
END FUNCTION matnorm_csr
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE full_to_csr(fullmat, mat)
!
! Convert full rectangular matrix to csr mat
!
DOUBLE PRECISION, INTENT(inout) :: fullmat(:,:)
TYPE(csr_mat), INTENT(out) :: mat
!
INTEGER :: m, n, nnz
INTEGER :: i, j, k
!
m = SIZE(fullmat,1)
n = SIZE(fullmat,2)
CALL init(m, 0, mat, ncols=n)
!
! Determine nnz of matrix
!
IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
IF(ASSOCIATED(mat%idiag)) DEALLOCATE(mat%idiag)
IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
!
ALLOCATE(mat%irow(m+1))
ALLOCATE(mat%idiag(m))
!
! Clear matrix small elements of fullmat
WHERE( ABS(fullmat) < 1.d-8) fullmat=0.0d0
!
mat%irow(1) = 1
nnz = 0
DO i=1,m
DO j=1,n
IF(ABS(fullmat(i,j)).GT.EPSILON(0.0d0)) THEN
nnz = nnz+1
IF(m.EQ.n .AND. i.EQ.j) THEN ! Only for square matrix
mat%idiag(i) = nnz
END IF
END IF
END DO
mat%irow(i+1) = nnz+1
END DO
!
! Allocate and fill the csr matrix structure
!
mat%nnz = nnz
ALLOCATE(mat%cols(nnz))
ALLOCATE(mat%val(nnz))
k=0
DO i=1,m
DO j=1,n
IF(ABS(fullmat(i,j)).GT.EPSILON(0.0d0)) THEN
k=k+1
mat%cols(k) = j
mat%val(k) = fullmat(i,j)
END IF
END DO
END DO
END SUBROUTINE full_to_csr
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE full_to_zcsr(fullmat, mat)
!
! Convert full rectangular matrix to csr mat
!
DOUBLE COMPLEX, INTENT(inout) :: fullmat(:,:)
TYPE(zcsr_mat), INTENT(out) :: mat
!
INTEGER :: m, n, nnz
INTEGER :: i, j, k
!
m = SIZE(fullmat,1)
n = SIZE(fullmat,2)
CALL init(m, 0, mat, ncols=n)
!
! Determine nnz of matrix
!
IF(ASSOCIATED(mat%irow)) DEALLOCATE(mat%irow)
IF(ASSOCIATED(mat%idiag)) DEALLOCATE(mat%idiag)
IF(ASSOCIATED(mat%cols)) DEALLOCATE(mat%cols)
IF(ASSOCIATED(mat%val)) DEALLOCATE(mat%val)
!
ALLOCATE(mat%irow(m+1))
ALLOCATE(mat%idiag(m))
!
! Clear matrix small elements of fullmat
WHERE( ABS(fullmat) < 1.d-8) fullmat=(0.0d0,0.0d0)
!
mat%irow(1) = 1
nnz = 0
DO i=1,m
DO j=1,n
IF(ABS(fullmat(i,j)).GT.EPSILON(0.0d0)) THEN
nnz = nnz+1
IF(m.EQ.n .AND. i.EQ.j) THEN ! Only for square matrix
mat%idiag(i) = nnz
END IF
END IF
END DO
mat%irow(i+1) = nnz+1
END DO
!
! Allocate and fill the csr matrix structure
!
mat%nnz = nnz
ALLOCATE(mat%cols(nnz))
ALLOCATE(mat%val(nnz))
k=0
DO i=1,m
DO j=1,n
IF(ABS(fullmat(i,j)).GT.EPSILON(0.0d0)) THEN
k=k+1
mat%cols(k) = j
mat%val(k) = fullmat(i,j)
END IF
END DO
END DO
END SUBROUTINE full_to_zcsr
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE check_dom(mat)
!
! Check whether mat is strict diagonal dominabce.
!
TYPE(csr_mat), INTENT(in) :: mat
DOUBLE PRECISION :: arow(mat%rank), asum(mat%rank)
INTEGER :: n, i, j1, j2, jdiag
!
n = mat%rank
DO i=1,n
j1 = mat%irow(i)
jdiag = mat%idiag(i)
j2 = mat%irow(i+1)-1
asum(i) = SUM(ABS(mat%val(j1:j2))) / ABS(mat%val(jdiag)) - 1.0d0
END DO
WRITE(*,'(/a,1pe12.3)') 'Max of sum of off-diag', MAXVAL(asum)
!
END SUBROUTINE check_dom
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE kron_csr(mata, matb, matc)
!
! Kronecker product of 2 CSR matrices
!
USE sparse, ONLY : isearch
TYPE(csr_mat), INTENT(in) :: mata, matb
TYPE(csr_mat), INTENT(out) :: matc
!
INTEGER :: m1, n1, nnz1, m2, n2, nnz2, m, n, nnz
INTEGER :: i,i1,i2,j1,s,s1,s2,e,e1,e2,k,nc2
!
m1 = mata%mrows
n1 = mata%ncols
nnz1 = mata%nnz
m2 = matb%mrows
n2 = matb%ncols
nnz2 = matb%nnz
m = m1*m2
n = n1*n2
nnz = nnz1*nnz2
!
CALL init(m, 0, matc, ncols=n)
matc%nnz = nnz
IF(ASSOCIATED(matc%irow)) DEALLOCATE(matc%irow)
IF(ASSOCIATED(matc%idiag)) DEALLOCATE(matc%idiag)
IF(ASSOCIATED(matc%cols)) DEALLOCATE(matc%cols)
IF(ASSOCIATED(matc%val)) DEALLOCATE(matc%val)
ALLOCATE(matc%irow(m+1))
IF(m.EQ.n) THEN
ALLOCATE(matc%idiag(m)) ! Only for square matrices
END IF
ALLOCATE(matc%cols(nnz))
ALLOCATE(matc%val(nnz))
!
k = 0
matc%irow(1) = 1
DO i1=1,m1
s1=mata%irow(i1)
e1=mata%irow(i1+1)-1
DO i2=1,m2
s2=matb%irow(i2)
e2=matb%irow(i2+1)-1
nc2=e2-s2+1
DO j1=s1,e1
matc%val(k+1:k+nc2) = mata%val(j1)*matb%val(s2:e2)
matc%cols(k+1:k+nc2) = (mata%cols(j1)-1)*n2 + matb%cols(s2:e2)
k = k+nc2
matc%irow((i1-1)*m2+i2+1) = k+1 ! Points to next row
END DO
END DO
END DO
!
! Search the diagonals
DO i=1,matc%mrows
s = matc%irow(i)
e = matc%irow(i+1)-1
matc%idiag(i) = isearch(matc%cols(s:e),i) + s
END DO
END SUBROUTINE kron_csr
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
END MODULE csr
c4science · Help