Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/src/sparse_mod.f90
f1da25c60953feature/refactoring

Last Change
Raw File

sparse_mod.f90
View Options

!>
!> @file sparse_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/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
MODULE sparse
!
! SPARSE: Implement sparse matrix using dynamic linked lists
! as matrix rows.
!
! T.M. Tran, CRPP-EPFL
! October 2010
!
IMPLICIT NONE
!
TYPE elt
INTEGER :: index=0
DOUBLE PRECISION :: val=0.0d0
TYPE(elt), POINTER :: next => NULL()
END TYPE elt
!
TYPE zelt
INTEGER :: index=0
DOUBLE COMPLEX :: val=(0.0d0, 0.0d0)
TYPE(zelt), POINTER :: next => NULL()
END TYPE zelt
!
TYPE sprow
INTEGER :: nnz=0 ! Number of non zeros in a row
TYPE(elt), POINTER :: row0 => NULL() ! Points to head of a (sparse) row
END TYPE sprow
!
TYPE zsprow
INTEGER :: nnz=0 ! Number of non zeros in a row
TYPE(zelt), POINTER :: row0 => NULL() ! Points to head of a (sparse) row
END TYPE zsprow
!
TYPE spmat
INTEGER :: rank
TYPE(sprow), POINTER :: row(:) => NULL()
END TYPE spmat
!
TYPE zspmat
INTEGER :: rank
TYPE(zsprow), POINTER :: row(:) => NULL()
END TYPE zspmat
!
INTERFACE init
MODULE PROCEDURE init_spmat, init_zspmat
END INTERFACE init
!
INTERFACE updtmat
MODULE PROCEDURE updt_sp, updt_zsp, updt_spmat, updt_zspmat
END INTERFACE updtmat
!
INTERFACE putele
MODULE PROCEDURE putele_sp, putele_zsp, putele_spmat, putele_zspmat
END INTERFACE putele
!
INTERFACE getele
MODULE PROCEDURE getele_sp, getele_zsp, getele_spmat, getele_zspmat
END INTERFACE getele
!
INTERFACE putrow
MODULE PROCEDURE putrow_csr, putrow_full, putrow_spmat, &
& putrow_zcsr, putrow_zfull, putrow_zspmat
END INTERFACE putrow
!
INTERFACE getrow
MODULE PROCEDURE getrow_csr, getrow_full, getrow_spmat, &
& getrow_zcsr, getrow_zfull, getrow_zspmat
END INTERFACE getrow
!
INTERFACE putcol
MODULE PROCEDURE putcol_spmat, putcol_zspmat
END INTERFACE putcol
!
INTERFACE getcol
MODULE PROCEDURE getcol_spmat, getcol_zspmat
END INTERFACE getcol
!
INTERFACE get_count
MODULE PROCEDURE get_count_sp, get_count_spmat, &
& get_count_zsp, get_count_zspmat
END INTERFACE get_count
!
INTERFACE destroy
MODULE PROCEDURE destroy_spmat, destroy_row, destroy_node, &
& destroy_zspmat, destroy_zrow, destroy_znode
END INTERFACE destroy
!
CONTAINS
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE init_spmat(n, mat, istart, iend)
!
! Initial an empty sparse matrix
!
INTEGER, INTENT(in) :: n
INTEGER, INTENT(in), OPTIONAL :: istart, iend
TYPE(spmat) :: mat
!
mat%rank = n
IF(ASSOCIATED(mat%row)) DEALLOCATE(mat%row)
IF(PRESENT(istart)) THEN
ALLOCATE(mat%row(istart:iend))
ELSE
ALLOCATE(mat%row(n))
END IF
!
END SUBROUTINE init_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE init_zspmat(n, mat, istart, iend)
!
! Initial an empty sparse matrix
!
INTEGER, INTENT(in) :: n
INTEGER, INTENT(in), OPTIONAL :: istart, iend
TYPE(zspmat) :: mat
!
mat%rank = n
IF(ASSOCIATED(mat%row)) DEALLOCATE(mat%row)
IF(PRESENT(istart)) THEN
ALLOCATE(mat%row(istart:iend))
ELSE
ALLOCATE(mat%row(n))
END IF
!
END SUBROUTINE init_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE updt_sp(arow, j, val)
!
! Update element j of row arow or insert it in an increasing "index"
!
TYPE(sprow), TARGET :: arow
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, INTENT(in) :: val
!
TYPE(elt), TARGET :: pre_root
TYPE(elt), POINTER :: t, p
!
pre_root%next => arow%row0 ! pre_root is linked to the head of the list.
t => pre_root
DO WHILE( ASSOCIATED(t%next) )
p => t%next
IF( p%index .EQ. j ) THEN
p%val = p%val+val
RETURN
END IF
IF( p%index .GT. j ) EXIT
t => t%next
END DO
ALLOCATE(p)
p = elt(j, val, t%next)
t%next => p
!
arow%nnz = arow%nnz+1
arow%row0 => pre_root%next ! In case the head is altered
END SUBROUTINE updt_sp
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE updt_zsp(arow, j, val)
!
! Update element j of row arow or insert it in an increasing "index"
!
TYPE(zsprow), TARGET :: arow
INTEGER, INTENT(in) :: j
DOUBLE COMPLEX, INTENT(in) :: val
!
TYPE(zelt), TARGET :: pre_root
TYPE(zelt), POINTER :: t, p
!
pre_root%next => arow%row0 ! pre_root is linked to the head of the list.
t => pre_root
DO WHILE( ASSOCIATED(t%next) )
p => t%next
IF( p%index .EQ. j ) THEN
p%val = p%val+val
RETURN
END IF
IF( p%index .GT. j ) EXIT
t => t%next
END DO
ALLOCATE(p)
p = zelt(j, val, t%next)
t%next => p
!
arow%nnz = arow%nnz+1
arow%row0 => pre_root%next ! In case the head is altered
END SUBROUTINE updt_zsp
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE updt_spmat(mat, i, j, val)
!
! Update element Aij of sparse matrix
!
TYPE(spmat) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE PRECISION, INTENT(in) :: val
!
CALL updt_sp(mat%row(i), j, val)
END SUBROUTINE updt_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE updt_zspmat(mat, i, j, val)
!
! Update element Aij of sparse matrix
!
TYPE(zspmat) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE COMPLEX, INTENT(in) :: val
!
CALL updt_zsp(mat%row(i), j, val)
END SUBROUTINE updt_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getele_sp(arow, j, val, found)
!
! Get element j from row arow
!
TYPE(sprow), TARGET :: arow
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, INTENT(out) :: val
LOGICAL, INTENT(out), OPTIONAL :: found
!
TYPE(elt), POINTER :: t
INTEGER :: i
!
val = 0.0d0
t => arow%row0 ! Start of a row
DO WHILE( ASSOCIATED(t) )
IF(t%index .EQ. j) THEN
val = t%val
IF(PRESENT(found)) found = .TRUE.
RETURN
END IF
t => t%next
END DO
IF(PRESENT(found)) found = .FALSE.
END SUBROUTINE getele_sp
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getele_zsp(arow, j, val, found)
!
! Get element j from row arow
!
TYPE(zsprow), TARGET :: arow
INTEGER, INTENT(in) :: j
DOUBLE COMPLEX, INTENT(out) :: val
LOGICAL, INTENT(out), OPTIONAL :: found
!
TYPE(zelt), POINTER :: t
INTEGER :: i
!
val = 0.0d0
t => arow%row0 ! Start of a row
DO WHILE( ASSOCIATED(t) )
IF(t%index .EQ. j) THEN
val = t%val
IF(PRESENT(found)) found = .TRUE.
RETURN
END IF
t => t%next
END DO
IF(PRESENT(found)) found = .FALSE.
END SUBROUTINE getele_zsp
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getele_spmat(mat, i, j, val)
!
! Get element (i,j) of sparse matrix
!
TYPE(spmat) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE PRECISION, INTENT(out) :: val
!
CALL getele(mat%row(i), j, val)
END SUBROUTINE getele_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getele_zspmat(mat, i, j, val)
!
! Get element (i,j) of sparse matrix
!
TYPE(zspmat) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE COMPLEX, INTENT(out) :: val
!
CALL getele(mat%row(i), j, val)
END SUBROUTINE getele_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putele_sp(arow, j, val, nlforce_zero)
!
! Put (overwrite) element j of row arow or insert it in an increasing "index"
!
TYPE(sprow), TARGET :: arow
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, INTENT(in) :: val
LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
TYPE(elt), TARGET :: pre_root
TYPE(elt), POINTER :: t, p
LOGICAL :: rmnode
!
pre_root%next => arow%row0 ! pre_root is linked to the head of the list.
t => pre_root
!
! Remove node which has zero val or not?
! But never create new node with zero val
!
rmnode = .TRUE.
IF(PRESENT(nlforce_zero)) rmnode = .NOT.nlforce_zero
!
DO WHILE( ASSOCIATED(t%next) )
p => t%next
IF( p%index .EQ. j ) THEN
IF(ABS(val).LE.EPSILON(0.0d0) .AND. rmnode) THEN ! Remove the node for zero val!
t%next => p%next
arow%nnz = arow%nnz-1
arow%row0 => pre_root%next ! In case the head is altered
DEALLOCATE(p)
ELSE
p%val = val
END IF
RETURN
END IF
IF( p%index .GT. j ) EXIT
t => t%next
END DO
!
! Never create new node with zero val
!
IF(ABS(val).GT.EPSILON(0.0d0)) THEN
ALLOCATE(p)
p = elt(j, val, t%next)
t%next => p
arow%nnz = arow%nnz+1
arow%row0 => pre_root%next
END IF
END SUBROUTINE putele_sp
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putele_zsp(arow, j, val, nlforce_zero)
!
! Put (overwrite) element j of row arow or insert it in an increasing "index"
!
TYPE(zsprow), TARGET :: arow
INTEGER, INTENT(in) :: j
DOUBLE COMPLEX, INTENT(in) :: val
LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
TYPE(zelt), TARGET :: pre_root
TYPE(zelt), POINTER :: t, p
LOGICAL :: rmnode
!
pre_root%next => arow%row0 ! pre_root is linked to the head of the list.
t => pre_root
!
! Remove node which has zero val or not?
! But never create new node with zero val
!
rmnode = .TRUE.
IF(PRESENT(nlforce_zero)) rmnode = .NOT.nlforce_zero
!
DO WHILE( ASSOCIATED(t%next) )
p => t%next
IF( p%index .EQ. j ) THEN
IF(ABS(val).LE.EPSILON(0.0d0) .AND. rmnode) THEN ! Remove the node for zero val!
t%next => p%next
arow%nnz = arow%nnz-1
arow%row0 => pre_root%next ! In case the head is altered
DEALLOCATE(p)
ELSE
p%val = val
END IF
RETURN
END IF
IF( p%index .GT. j ) EXIT
t => t%next
END DO
!
! Never create new node with zero val
!
IF(ABS(val).GT.EPSILON(0.0d0)) THEN
ALLOCATE(p)
p = zelt(j, val, t%next)
t%next => p
arow%nnz = arow%nnz+1
arow%row0 => pre_root%next
END IF
END SUBROUTINE putele_zsp
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putele_spmat(mat, i, j, val, nlforce_zero)
!
! Put element (i,j) of matrix
!
TYPE(spmat) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE PRECISION, INTENT(in) :: val
LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
CALL putele(mat%row(i), j, val, nlforce_zero)
END SUBROUTINE putele_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putele_zspmat(mat, i, j, val, nlforce_zero)
!
! Put element (i,j) of matrix
!
TYPE(zspmat) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE COMPLEX, INTENT(in) :: val
LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
CALL putele(mat%row(i), j, val, nlforce_zero)
END SUBROUTINE putele_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
INTEGER FUNCTION get_count_sp(arow)
!
! Number of elements in arow
!
TYPE(sprow), INTENT(in) :: arow
TYPE(elt), POINTER :: t
INTEGER :: i
!
t => arow%row0 ! Start of a row
i = 0
DO WHILE( ASSOCIATED(t) )
i=i+1
t => t%next
END DO
get_count_sp = i
END FUNCTION get_count_sp
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
INTEGER FUNCTION get_count_zsp(arow)
!
! Number of elements in arow
!
TYPE(zsprow), INTENT(in) :: arow
TYPE(zelt), POINTER :: t
INTEGER :: i
!
t => arow%row0 ! Start of a row
i = 0
DO WHILE( ASSOCIATED(t) )
i=i+1
t => t%next
END DO
get_count_zsp = i
END FUNCTION get_count_zsp
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
INTEGER FUNCTION get_count_spmat(mat, nnz)
!
! Number of non-zeros in sparse matrix
!
TYPE(spmat) :: mat
INTEGER, INTENT(out), OPTIONAL :: nnz(:)
!
INTEGER :: i, c(LBOUND(mat%row,1):UBOUND(mat%row,1))
DO i=LBOUND(mat%row,1),UBOUND(mat%row,1)
c(i) = get_count_sp(mat%row(i))
END DO
IF(PRESENT(nnz)) nnz = c
get_count_spmat = SUM(c)
END FUNCTION get_count_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
INTEGER FUNCTION get_count_zspmat(mat, nnz)
!
! Number of non-zeros in sparse matrix
!
TYPE(zspmat) :: mat
INTEGER, INTENT(out), OPTIONAL :: nnz(:)
!
INTEGER :: i, c(LBOUND(mat%row,1):UBOUND(mat%row,1))
DO i=LBOUND(mat%row,1),UBOUND(mat%row,1)
c(i) = get_count_zsp(mat%row(i))
END DO
IF(PRESENT(nnz)) nnz = c
get_count_zspmat = SUM(c)
END FUNCTION get_count_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getrow_csr(arow, arr, col, count)
!
! Get a row from sparse row arow and put it in a CSR format
!
TYPE(sprow), INTENT(in) :: arow
DOUBLE PRECISION, INTENT(out) :: arr(:)
INTEGER, INTENT(out) :: col(:)
INTEGER, OPTIONAL, INTENT(out) :: count
!
TYPE(elt), POINTER :: t
INTEGER :: i
!
t => arow%row0 ! Start of a row
i = 0
DO WHILE( ASSOCIATED(t) )
i=i+1
col(i) = t%index
arr(i) = t%val
t => t%next
END DO
IF(PRESENT(count)) count = i
END SUBROUTINE getrow_csr
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getrow_zcsr(arow, arr, col, count)
!
! Get a row from sparse row arow and put it in a CSR format
!
TYPE(zsprow), INTENT(in) :: arow
DOUBLE COMPLEX, INTENT(out) :: arr(:)
INTEGER, INTENT(out) :: col(:)
INTEGER, OPTIONAL, INTENT(out) :: count
!
TYPE(zelt), POINTER :: t
INTEGER :: i
!
t => arow%row0 ! Start of a row
i = 0
DO WHILE( ASSOCIATED(t) )
i=i+1
col(i) = t%index
arr(i) = t%val
t => t%next
END DO
IF(PRESENT(count)) count = i
END SUBROUTINE getrow_zcsr
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getrow_full(arow, arr, count)
!
! Get a row from sparse row arow and put it in an full row
!
TYPE(sprow), INTENT(in) :: arow
DOUBLE PRECISION, INTENT(out) :: arr(:)
INTEGER, OPTIONAL, INTENT(out) :: count
!
TYPE(elt), POINTER :: t
INTEGER :: n, i, j
!
n = SIZE(arr)
arr = 0.0d0
t => arow%row0 ! Start of a row
i = 0
DO WHILE( ASSOCIATED(t) )
i=i+1
j = t%index
IF(j.LE.n) THEN
arr(j) = t%val
t => t%next
ELSE
WRITE(*,'(a)') 'GETROW_FULL: size of input ARR too small!'
STOP
END IF
END DO
IF(PRESENT(count)) count = i
END SUBROUTINE getrow_full
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getrow_zfull(arow, arr, count)
!
! Get a row from sparse row arow and put it in an full row
!
TYPE(zsprow), INTENT(in) :: arow
DOUBLE COMPLEX, INTENT(out) :: arr(:)
INTEGER, OPTIONAL, INTENT(out) :: count
!
TYPE(zelt), POINTER :: t
INTEGER :: n, i, j
!
n = SIZE(arr)
arr = 0.0d0
t => arow%row0 ! Start of a row
i = 0
DO WHILE( ASSOCIATED(t) )
i=i+1
j = t%index
IF(j.LE.n) THEN
arr(j) = t%val
t => t%next
ELSE
WRITE(*,'(a)') 'GETROW_FULL: size of input ARR too small!'
STOP
END IF
END DO
IF(PRESENT(count)) count = i
END SUBROUTINE getrow_zfull
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getrow_spmat(mat, i, arr, col)
!
! Get a row from sparse matrix
!
TYPE(spmat), INTENT(in) :: mat
INTEGER, INTENT(in) :: i
DOUBLE PRECISION, INTENT(out) :: arr(:)
INTEGER, INTENT(out), OPTIONAL :: col(:)
!
IF(PRESENT(col)) THEN ! The output row is defined by (col, arr)
CALL getrow_csr(mat%row(i), arr, col)
ELSE
CALL getrow_full(mat%row(i), arr)
END IF
END SUBROUTINE getrow_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getrow_zspmat(mat, i, arr, col)
!
! Get a row from sparse matrix
!
TYPE(zspmat), INTENT(in) :: mat
INTEGER, INTENT(in) :: i
DOUBLE COMPLEX, INTENT(out) :: arr(:)
INTEGER, INTENT(out), OPTIONAL :: col(:)
!
IF(PRESENT(col)) THEN ! The output row is defined by (col, arr)
CALL getrow_zcsr(mat%row(i), arr, col)
ELSE
CALL getrow_zfull(mat%row(i), arr)
END IF
END SUBROUTINE getrow_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putrow_csr(arow, arr, col, nlforce_zero)
!
! Put a row from sparse row arow and put it in a CSR format
!
TYPE(sprow), INTENT(inout) :: arow
DOUBLE PRECISION, INTENT(in) :: arr(:)
INTEGER, INTENT(in) :: col(:)
LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
INTEGER :: n, i
!
n=SIZE(arr)
DO i=1,n
CALL putele(arow, col(i), arr(i), nlforce_zero)
END DO
END SUBROUTINE putrow_csr
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putrow_zcsr(arow, arr, col, nlforce_zero)
!
! Put a row from sparse row arow and put it in a CSR format
!
TYPE(zsprow), INTENT(inout) :: arow
DOUBLE COMPLEX, INTENT(in) :: arr(:)
INTEGER, INTENT(in) :: col(:)
LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
INTEGER :: n, i
!
n=SIZE(arr)
DO i=1,n
CALL putele(arow, col(i), arr(i), nlforce_zero)
END DO
END SUBROUTINE putrow_zcsr
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putrow_full(arow, arr, nlforce_zero)
!
! Put a row from sparse row arow and put it in a full row
!
TYPE(sprow), INTENT(inout) :: arow
DOUBLE PRECISION, INTENT(in) :: arr(:)
LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
INTEGER :: n, i
!
n=SIZE(arr)
DO i=1,n
CALL putele(arow, i, arr(i), nlforce_zero)
END DO
END SUBROUTINE putrow_full
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putrow_zfull(arow, arr, nlforce_zero)
!
! Put a row from sparse row arow and put it in a full row
!
TYPE(zsprow), INTENT(inout) :: arow
DOUBLE COMPLEX, INTENT(in) :: arr(:)
LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
INTEGER :: n, i
!
n=SIZE(arr)
DO i=1,n
CALL putele(arow, i, arr(i), nlforce_zero)
END DO
END SUBROUTINE putrow_zfull
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putrow_spmat(mat, i, arr, col, nlforce_zero)
!
! Put a row to matrix
!
TYPE(spmat) :: mat
INTEGER, intent(in) :: i
DOUBLE PRECISION, INTENT(in) :: arr(:)
INTEGER, INTENT(in), OPTIONAL :: col(:)
LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
IF(PRESENT(col)) THEN ! The input row is defined by (col, arr)
CALL putrow_csr(mat%row(i), arr, col, nlforce_zero)
ELSE
CALL putrow_full(mat%row(i), arr, nlforce_zero)
END IF
END SUBROUTINE putrow_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putrow_zspmat(mat, i, arr, col, nlforce_zero)
!
! Put a row to matrix
!
TYPE(zspmat) :: mat
INTEGER, intent(in) :: i
DOUBLE COMPLEX, INTENT(in) :: arr(:)
INTEGER, INTENT(in), OPTIONAL :: col(:)
LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
IF(PRESENT(col)) THEN ! The input row is defined by (col, arr)
CALL putrow_zcsr(mat%row(i), arr, col, nlforce_zero)
ELSE
CALL putrow_zfull(mat%row(i), arr, nlforce_zero)
END IF
END SUBROUTINE putrow_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putcol_spmat(mat, j, arr, nlforce_zero)
!
! Put a column to mtarix
!
TYPE(spmat) :: mat
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, INTENT(in) :: arr(:)
LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
INTEGER :: i
DO i=1,mat%rank
CALL putele(mat, i, j, arr(i), nlforce_zero)
END DO
END SUBROUTINE putcol_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE putcol_zspmat(mat, j, arr, nlforce_zero)
!
! Put a column to mtarix
!
TYPE(zspmat) :: mat
INTEGER, INTENT(in) :: j
DOUBLE COMPLEX, INTENT(in) :: arr(:)
LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
INTEGER :: i
DO i=1,mat%rank
CALL putele(mat, i, j, arr(i), nlforce_zero)
END DO
END SUBROUTINE putcol_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getcol_spmat(mat, j, arr)
!
! Get column j of matrix
!
TYPE(spmat) :: mat
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, INTENT(out) :: arr(:)
INTEGER :: i
DO i=1,mat%rank
CALL getele(mat, i, j, arr(i))
END DO
END SUBROUTINE getcol_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE getcol_zspmat(mat, j, arr)
!
! Get column j of matrix
!
TYPE(zspmat) :: mat
INTEGER, INTENT(in) :: j
DOUBLE COMPLEX, INTENT(out) :: arr(:)
INTEGER :: i
DO i=1,mat%rank
CALL getele(mat, i, j, arr(i))
END DO
END SUBROUTINE getcol_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE destroy_row(arow)
!
! Deallocate the sparse row
!
TYPE(sprow), INTENT(inout) :: arow
!
IF(ASSOCIATED(arow%row0)) CALL destroy_node(arow%row0)
arow%nnz = get_count(arow)
END SUBROUTINE destroy_row
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE destroy_zrow(arow)
!
! Deallocate the sparse row
!
TYPE(zsprow), INTENT(inout) :: arow
!
IF(ASSOCIATED(arow%row0)) CALL destroy_znode(arow%row0)
arow%nnz = get_count(arow)
END SUBROUTINE destroy_zrow
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
RECURSIVE SUBROUTINE destroy_node(p)
!
! Deallocate recursively the linked list
!
TYPE(elt), POINTER :: p
!
IF(ASSOCIATED(p%next)) CALL destroy_node(p%next)
DEALLOCATE(p)
END SUBROUTINE destroy_node
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
RECURSIVE SUBROUTINE destroy_znode(p)
!
! Deallocate recursively the linked list
!
TYPE(zelt), POINTER :: p
!
IF(ASSOCIATED(p%next)) CALL destroy_znode(p%next)
DEALLOCATE(p)
END SUBROUTINE destroy_znode
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE destroy_spmat(mat)
!
! Deallocate the sparse matrix
!
TYPE(spmat) :: mat
INTEGER :: n, i
!
n = mat%rank
DO i=LBOUND(mat%row,1),UBOUND(mat%row,1)
CALL destroy(mat%row(i))
END DO
DEALLOCATE(mat%row)
END SUBROUTINE destroy_spmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE destroy_zspmat(mat)
!
! Deallocate the sparse matrix
!
TYPE(zspmat) :: mat
INTEGER :: n, i
!
n = mat%rank
DO i=LBOUND(mat%row,1),UBOUND(mat%row,1)
CALL destroy(mat%row(i))
END DO
DEALLOCATE(mat%row)
END SUBROUTINE destroy_zspmat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
INTEGER FUNCTION isearch(karr, k)
!
! Sequential search an ordered table of integers
!
INTEGER, INTENT(in) :: karr(0:)
INTEGER, INTENT(in) :: k
INTEGER :: n
!
n=SIZE(karr)
isearch = -1 ! Failure
IF( k.GT.karr(n-1)) RETURN
!
isearch=0
DO
IF( k.LE.karr(isearch)) EXIT
isearch = isearch+1
END DO
IF( k.NE.karr(isearch)) isearch = -1 ! Failure
END FUNCTION isearch
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
INTEGER FUNCTION isearch_bin(karr, k)
!
! Binary search an ordered table of integers
!
INTEGER, INTENT(in) :: karr(0:)
INTEGER, INTENT(in) :: k
INTEGER :: n
INTEGER :: l, u
!
n=SIZE(karr)
isearch_bin = -1 ! Failure
IF( k.LT.karr(0) .OR. k.GT.karr(n-1)) RETURN
!
l=0; u=n-1
DO WHILE(l.LE.u)
isearch_bin = (l+u)/2
IF(k.EQ.karr(isearch_bin)) THEN
RETURN
ELSE IF(k.LT.karr(isearch_bin)) THEN
u = isearch_bin-1
ELSE
l = isearch_bin+1
END IF
END DO
isearch_bin = -1 ! Failure
END FUNCTION isearch_bin
!
END MODULE sparse
c4science · Help