sparse_matrix.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
module m_sparse_matrix
  use hecmw_util
  use m_hecmw_comm_f
  implicit none
 
  private
 
  public :: sparse_matrix_type_csr
  public :: sparse_matrix_type_coo
 
  public :: sparse_matrix_symtype_asym
  public :: sparse_matrix_symtype_spd
  public :: sparse_matrix_symtype_sym
 
  public :: sparse_matrix
  public :: sparse_matrix_set_type
  public :: sparse_matrix_init
  public :: sparse_matrix_finalize
  public :: sparse_matrix_dump
  public :: sparse_matrix_gather_rhs
  public :: sparse_matrix_scatter_rhs
  public :: sparse_matrix_is_sym
 
  integer(kind=kint), parameter :: sparse_matrix_type_csr=1
  integer(kind=kint), parameter :: sparse_matrix_type_coo=2
 
  integer(kind=kint), parameter :: sparse_matrix_symtype_asym=0
  integer(kind=kint), parameter :: sparse_matrix_symtype_spd=1
  integer(kind=kint), parameter :: sparse_matrix_symtype_sym=2
 
  type sparse_matrix
    integer(kind=kint) :: type
    integer(kind=kint) :: symtype
    integer(kind=kint) :: n, n_loc
    integer(kind=kint) :: nz
    integer(kind=kint), pointer :: irn(:) => null()
    integer(kind=kint), pointer :: jcn(:) => null()
    real(kind=kreal), pointer :: a(:) => null()
    real(kind=kreal), pointer :: rhs(:)
    integer(kind=kint) :: offset
    integer(kind=kint), pointer :: n_counts(:) => null()
    integer(kind=kint), pointer :: displs(:) => null()
    integer(kind=kint), pointer :: conv_ext(:) => null()
    integer(kind=kint) :: iterlog
    integer(kind=kint) :: timelog
    logical :: is_initialized = .false.
  end type sparse_matrix
 
contains
 
  !!! public subroutines
 
  subroutine sparse_matrix_set_type(spMAT, type, symtype)
    type (sparse_matrix), intent(inout) :: spmat
    integer(kind=kint), intent(in) :: type
    integer(kind=kint), intent(in) :: symtype
    if (type == sparse_matrix_type_csr .or. &
        type == sparse_matrix_type_coo) then
      spmat%type = type
    else
      write(*,*) 'ERROR: unknown sparse matrix type'
      stop
    endif
    if (symtype == sparse_matrix_symtype_asym .or. &
        symtype == sparse_matrix_symtype_spd .or. &
        symtype == sparse_matrix_symtype_sym) then
      spmat%symtype = symtype
    else
      write(*,*) 'ERROR: unknown symmetry type for sparse matrix'
      stop
    endif
  end subroutine sparse_matrix_set_type
 
  subroutine sparse_matrix_init(spMAT, N_loc, NZ)
    type (sparse_matrix), intent(inout) :: spmat
    integer(kind=kint), intent(in) :: n_loc
    integer(kind=kint), intent(in) :: nz
    if (spmat%is_initialized) then
      !write(*,*) "WARNING: sparse_matrix_init_prof: already initialized; freeing"
      call sparse_matrix_finalize(spmat)
    endif
    call sparse_matrix_set_dimension(spmat, n_loc)
    call sparse_matrix_set_offsets(spmat)
    call sparse_matrix_allocate_arrays(spmat, nz)
    spmat%is_initialized = .true.
  end subroutine sparse_matrix_init
 
  subroutine sparse_matrix_finalize(spMAT)
    type (sparse_matrix), intent(inout) :: spmat
    call sparse_matrix_free_arrays(spmat)
    call sparse_matrix_clear(spmat)
    spmat%is_initialized = .false.
  end subroutine sparse_matrix_finalize
 
  subroutine sparse_matrix_dump(spMAT)
    type(sparse_matrix), intent(in) :: spmat
    integer(kind=kint),save :: num_call
    character(len=128) :: fname
    integer(kind=kint) :: i,myrank
    myrank=hecmw_comm_get_rank()
    write(fname,'(A,I0,A,I0)') 'sparse-matrix-',num_call,'.coo.',myrank
    num_call = num_call + 1
    write(*,*) 'Dumping sparse matrix to ', fname
    open(91,file=fname,iostat=i)
    if( i /= 0 ) return
    if (spmat%type == sparse_matrix_type_csr) then
      if (spmat%symtype == sparse_matrix_symtype_asym) &
        write(91,*) '%SPARSE MATRIX CSR ASYM'
      if (spmat%symtype == sparse_matrix_symtype_spd) &
        write(91,*) '%SPARSE MATRIX CSR SPD'
      if (spmat%symtype == sparse_matrix_symtype_sym) &
        write(91,*) '%SPARSE MATRIX CSR SYM'
      write(91,*) spmat%N_loc, spmat%N_loc, spmat%NZ
      do i=1,spmat%N_loc+1
        write(91,*) spmat%IRN(i)
      enddo
      do i=1,spmat%NZ
        write(91,*) spmat%JCN(i), spmat%A(i)
      enddo
    elseif (spmat%type == sparse_matrix_type_coo) then
      if (spmat%symtype == sparse_matrix_symtype_asym) &
        write(91,*) '%SPARSE MATRIX COO ASYM'
      if (spmat%symtype == sparse_matrix_symtype_spd) &
        write(91,*) '%SPARSE MATRIX COO SPD'
      if (spmat%symtype == sparse_matrix_symtype_sym) &
        write(91,*) '%SPARSE MATRIX COO SYM'
      write(91,*) spmat%N_loc, spmat%N_loc, spmat%NZ
      do i=1,spmat%NZ
        write(91,*) spmat%IRN(i), spmat%JCN(i), spmat%A(i)
      enddo
    endif
    close(91)
  end subroutine sparse_matrix_dump
 
  subroutine sparse_matrix_gather_rhs(spMAT, rhs_all)
    type (sparse_matrix), intent(in) :: spmat
    real(kind=kreal), intent(out) :: rhs_all(:)
    integer(kind=kint) :: ierr,i
    if (hecmw_comm_get_size() == 1) then
      do i=1,spmat%N_loc
        rhs_all(i) = spmat%rhs(i)
      enddo
    else
      call hecmw_gatherv_real(spmat%rhs, spmat%N_loc, &
        rhs_all, spmat%N_COUNTS, spmat%DISPLS, &
        0, hecmw_comm_get_comm())
    endif
  end subroutine sparse_matrix_gather_rhs
 
  subroutine sparse_matrix_scatter_rhs(spMAT, rhs_all)
    type (sparse_matrix), intent(inout) :: spmat
    real(kind=kreal), intent(in) :: rhs_all(:)
    integer(kind=kint) :: ierr,i
    if (hecmw_comm_get_size() == 1) then
      do i=1,spmat%N_loc
        spmat%rhs(i) = rhs_all(i)
      enddo
    else
      call hecmw_scatterv_real( &
        rhs_all, spmat%N_COUNTS, spmat%DISPLS, &
        spmat%rhs, spmat%N_loc, &
        0, hecmw_comm_get_comm())
    endif
  end subroutine sparse_matrix_scatter_rhs
 
  function sparse_matrix_is_sym(spMAT)
    type(sparse_matrix), intent(inout) :: spmat
    logical :: sparse_matrix_is_sym
    sparse_matrix_is_sym = (spmat%symtype > 0)
  end function sparse_matrix_is_sym
 
  !!! private subroutines
 
  subroutine sparse_matrix_set_dimension(spMAT, N_loc)
    type (sparse_matrix), intent(inout) :: spmat
    integer(kind=kint), intent(in) :: n_loc
    integer(kind=kint) :: ierr
    spmat%N_loc = n_loc
    if (hecmw_comm_get_size() == 1) then
      spmat%N = spmat%N_loc
    else
      call hecmw_allreduce_int_1(spmat%N_loc, spmat%N, hecmw_sum, &
        hecmw_comm_get_comm())
    endif
  end subroutine sparse_matrix_set_dimension
 
  subroutine sparse_matrix_set_offsets(spMAT)
    type (sparse_matrix), intent(inout) :: spmat
    integer(kind=kint) :: i,ierr,nprocs,myrank
    nprocs=hecmw_comm_get_size()
    myrank=hecmw_comm_get_rank()
    ! OFFSET
    !if (myrank == 0) then
    if (associated(spmat%N_COUNTS)) deallocate(spmat%N_COUNTS)
    if (associated(spmat%DISPLS)) deallocate(spmat%DISPLS)
    allocate(spmat%N_COUNTS(nprocs), spmat%DISPLS(nprocs), stat=ierr)
    if (ierr /= 0) then
      write(*,*) " Allocation error, spMAT%N_COUNTS, spMAT%DISPLS"
      call hecmw_abort(hecmw_comm_get_comm())
    endif
    !endif
    if (nprocs > 1) then
      call hecmw_allgather_int_1(spmat%N_loc, &
        spmat%N_COUNTS, hecmw_comm_get_comm())
    endif
    spmat%DISPLS(1)=0
    do i=1,nprocs-1
      spmat%DISPLS(i+1)=spmat%DISPLS(i)+spmat%N_COUNTS(i)
    enddo
    spmat%offset = spmat%DISPLS(myrank+1)
  end subroutine sparse_matrix_set_offsets
 
  subroutine sparse_matrix_allocate_arrays(spMAT, NZ)
    type(sparse_matrix), intent(inout) :: spmat
    integer(kind=kint), intent(in) :: nz
    integer(kind=kint) :: n_loc
    integer(kind=kint) :: ierr
    if (associated(spmat%IRN)) deallocate(spmat%IRN)
    if (associated(spmat%JCN)) deallocate(spmat%JCN)
    if (associated(spmat%A)) deallocate(spmat%A)
    ierr = -1
    n_loc=spmat%N_loc
    if (spmat%type == sparse_matrix_type_csr) then
      allocate(spmat%IRN(n_loc+1), spmat%JCN(nz), spmat%A(nz), stat=ierr)
    elseif (spmat%type == sparse_matrix_type_coo) then
      allocate(spmat%IRN(nz), spmat%JCN(nz), spmat%A(nz), stat=ierr)
    endif
    if (ierr /= 0) then
      write(*,*) " Allocation error, spMAT%IRN, spMAT%JCN, spMAT%A"
      call hecmw_abort(hecmw_comm_get_comm())
    endif
    spmat%NZ = nz
  end subroutine sparse_matrix_allocate_arrays
 
  subroutine sparse_matrix_free_arrays(spMAT)
    type (sparse_matrix), intent(inout) :: spmat
    if (associated(spmat%IRN)) deallocate(spmat%IRN)
    if (associated(spmat%JCN)) deallocate(spmat%JCN)
    if (associated(spmat%A)) deallocate(spmat%A)
    if (associated(spmat%N_COUNTS)) deallocate(spmat%N_COUNTS)
    if (associated(spmat%DISPLS)) deallocate(spmat%DISPLS)
    if (associated(spmat%conv_ext)) deallocate(spmat%conv_ext)
  end subroutine sparse_matrix_free_arrays
 
  subroutine sparse_matrix_clear(spMAT)
    type(sparse_matrix), intent(inout) :: spmat
    spmat%type = 0
    spmat%symtype = 0
    spmat%N = 0
    spmat%N_loc = 0
    spmat%NZ = 0
    spmat%IRN => null()
    spmat%JCN => null()
    spmat%A => null()
    spmat%offset = 0
    spmat%N_COUNTS => null()
    spmat%DISPLS => null()
    spmat%conv_ext => null()
    spmat%iterlog = 0
    spmat%timelog = 0
  end subroutine sparse_matrix_clear
 
end module m_sparse_matrix