scale_sparsemat.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> Module common / sparsemat
!!
!! @par Description
!!          A module to treat sparse matrix and the associated operations
!!
!! @author Yuta Kawai, Xuanzhengbo Ren, and Team SCALE
!!
!<
#include "scaleFElib.h"
module scale_sparsemat
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_const, only:    &
    undef8 => const_undef8, &
    const_eps
  
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !  
 
  !> Derived type to manage a sparse matrix
  type, public :: sparsemat
    integer :: m                      !< Number of row of original matrix
    integer :: n                      !< Number of column of original matrix
    integer :: nnz                    !< Number of nonzero of the matrix
    real(rp), allocatable :: val(:)   !< Array storing the values of the nonzeros
    integer, allocatable :: colidx(:) !< Array saving the column indices of the nonzeros
 
    ! for CSR format
    integer, allocatable :: rowptr(:) !< Array saving the start and end pointers of the nonzeros of the rows.
    integer :: rowptrsize             !< Size of rowPtr
 
    ! for ELL format
    integer :: col_size               !< Size of column of compressed matrix for ELL format (maximum number of nonzeros in a row)
 
    integer, private :: storage_format_id      !< Number of row of original matrix
  contains
    procedure, public :: init => sparsemat_init
    procedure, public :: final => sparsemat_final
    procedure, public :: print => sparsemat_print
    procedure, public :: replaceval => sparsemat_replceval
    procedure, public :: getval => sparsemat_getval
    procedure, public :: getstorageformatid => sparsemat_get_storage_format_id
  end type sparsemat
 
  interface sparsemat_matmul
    module procedure sparsemat_matmul1
    module procedure sparsemat_matmul1_2
    module procedure sparsemat_matmul2
  end interface
  public :: sparsemat_matmul
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  integer, public, parameter :: sparsemat_storage_typeid_csr = 1  !< Storage format ID for CSR format
  integer, public, parameter :: sparsemat_storage_typeid_ell = 2  !< Storage format ID for ELL format
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: sparsemat_matmul_csr_1
  private :: sparsemat_matmul_csr_2
  private :: sparsemat_matmul_ell_1
  private :: sparsemat_matmul_ell_2
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  
  !-----------------------------------------------------------------------------
 
contains
  !> Initialize an object to manage a sparse matrix
!OCL SERIAL
  subroutine sparsemat_init( this, mat, &
      EPS, storage_format )
    implicit none
 
    class(sparsemat), intent(inout) :: this
    real(rp), intent(in) :: mat(:,:)
    real(rp), optional, intent(in) :: eps
    character(len=*), optional, intent(in) :: storage_format
 
    integer :: i
    integer :: j
    integer :: l
 
    integer :: val_counter
    integer :: rowptr_counter
    
    real(dp) :: tmp_val(size(mat)+1)
    integer :: tmp_colind(size(mat)+1)
    integer :: tmp_rowptr(0:size(mat,1)+1)
 
    real(rp) :: eps_
    character(len=H_MID) :: storage_format_ = 'CSR' 
 
    integer :: row_nonzero_counter(size(mat,1))
    integer :: col_size_l
    !--------------------------------------------------------------------------- 
 
    !$acc enter data create(this)
 
    this%M = size(mat,1)
    this%N = size(mat,2)
 
    val_counter     = 1
    rowptr_counter  = 1
    tmp_rowptr(1)   = 1
    tmp_val(:)      = 0.0_rp
 
    if (present(eps)) then
      eps_ = eps
    else
      eps_ = const_eps * 500.0_rp ! ~ 1x10^-13
    end if
 
    if (present(storage_format)) then
      storage_format_ = storage_format
    end if
 
    select case (storage_format_)
    case ('CSR')
      this%storage_format_id = sparsemat_storage_typeid_csr
 
      do i=1, this%M
        do j=1, this%N
          if ( abs(mat(i,j)) > eps_ ) then
            tmp_val(val_counter) = mat(i,j)
            tmp_colind(val_counter) = j
            val_counter = val_counter + 1
          end if
        end do
        rowptr_counter = rowptr_counter + 1
        tmp_rowptr(rowptr_counter) = val_counter
      end do
    case ('ELL')
      this%storage_format_id = sparsemat_storage_typeid_ell
 
      do i=1, this%M
        row_nonzero_counter(i) = 0
        do j=1, this%N
          if ( abs(mat(i,j)) >= eps_ ) &
            row_nonzero_counter(i) = row_nonzero_counter(i) + 1
        end do
      end do
      this%col_size = maxval(row_nonzero_counter(:))
      val_counter   = this%M * this%col_size
 
      tmp_val(:)    = 0.0_rp
      tmp_colind(:) = -1
      do i=1, this%M
        col_size_l = 0
        do j=1,this%N
          if ( abs(mat(i,j)) >= eps_ ) then
              col_size_l = col_size_l + 1
              l = i+(col_size_l-1)*this%M
              tmp_val(l) = mat(i,j)   
              tmp_colind(l) = j
          end if
        end do
      end do
    case default
      log_error("SparseMat_Init",*) 'Not appropriate names of storage format. Check!', trim(storage_format_)
      call prc_abort
    end select
 
    this%nnz = val_counter
    allocate( this%val(val_counter) )
    allocate( this%colIdx(val_counter) )
    this%val(:)    = tmp_val(1:val_counter)
    this%colIdx(:) = tmp_colind(1:val_counter)
 
    select case (storage_format_)
    case ('CSR')
      allocate( this%rowPtr(rowptr_counter) )
      this%rowPtr(:) = tmp_rowptr(1:rowptr_counter)
      !$acc enter data copyin( this%rowPtr )
      !$acc enter data attach( this%rowPtr )
      this%rowPtrSize = rowptr_counter
    case('ELL')
      do l=2, val_counter
        if (this%colIdx(l) == -1) then
          this%colIdx(l) = this%colIdx(l-1)
        end if
      end do
    end select
    !$acc enter data copyin( this%val, this%colIdx )
    !$acc enter data attach( this%val, this%colIdx )
 
    !$acc update device( this%M, this%N, this%nnz, this%rowPtrSize, this%storage_format_id, this%col_size )
 
    !write(*,*) "--- Mat ------"
    !write(*,*) "shape:", shape(mat)
    !write(*,*) "size:", size(mat)
    ! do j=1, size(mat,2)
    !    write(*,*) mat(:,j)
    ! end do
    !write(*,*) "--- Compressed Mat ------"
    !write(*,*) "shape (Mat, colInd, rowPtr):", &
    !  & shape(this%val), shape(this%colInd), shape(this%rowPtr)
    ! write(*,*) "val:", this%val(:)
    ! write(*,*) "colInd:", this%colInd(:)
    ! write(*,*) "rowPtr:", this%rowPtr(:)
 
    return
  end subroutine sparsemat_init
 
  !> Finalize an object to manage a sparse matrix
!OCL SERIAL
  subroutine sparsemat_final(this)
    implicit none
    class(sparsemat), intent(inout) :: this
    !--------------------------------------------------------------------------- 
 
    !$acc exit data detach( this%val, this%colIdx )
    !$acc exit data delete( this%val, this%colIdx )
    deallocate( this%val )
    deallocate( this%colIdx )
    
    select case( this%storage_format_id )
    case( sparsemat_storage_typeid_csr )
      !$acc exit data detach( this%rowPtr )
      !$acc exit data delete( this%rowPtr )
      deallocate( this%rowPtr )
    end select
 
    !$acc exit data delete( this )
 
    !write(*,*) "--- Finalize sparsemat ---"
    return
  end subroutine sparsemat_final
 
  function sparsemat_getval(A, i, j) result(v)
    class(sparsemat), intent(in) :: a
    integer, intent(in) :: i, j
 
    real(dp) ::v
    integer :: n
    integer :: jj
    !--------------------------------------------------------------------------- 
 
    v = undef8
    select case( a%storage_format_id )
    case( sparsemat_storage_typeid_csr )
      do n=a%rowPtr(i),a%rowPtr(i+1)-1
        if (a%colIdx(n) == j) then
          v = a%val(n); exit
        end if
      end do
    case ( sparsemat_storage_typeid_ell )
      do jj=1, a%col_size
        n = i + (jj-1)*a%M
        if (a%colIdx(n) == j) then
          v = a%val(n); exit
        end if
      end do
    end select
 
    return
  end function sparsemat_getval
 
  subroutine sparsemat_replceval(A, i, j, v)
    
    class(sparsemat), intent(inout) :: a
    integer, intent(in) :: i, j
    real(rp), intent(in) ::v
 
    integer :: n
    integer :: jj
    !--------------------------------------------------------------------------- 
 
    select case ( a%storage_format_id )
    case ( sparsemat_storage_typeid_csr )
      do n=a%rowPtr(i),a%rowPtr(i+1)-1
        if (a%colIdx(n) == j) then
          a%val(n) = v; exit
        end if
      end do
    case ( sparsemat_storage_typeid_ell )
      do jj=1, a%col_size
        n = i + (jj-1)*a%M
        if (a%colIdx(n) == j) then
          a%val(n) = v; exit
        end if
      end do
    end select
 
    return
  end subroutine sparsemat_replceval
 
  subroutine sparsemat_print(A)
    implicit none
    class(sparsemat), intent(in) :: a
 
    real(rp) :: row_val(a%n)
    integer :: p
    integer :: j1, j2
    !--------------------------------------------------------------------------- 
 
    write(*,*) "-- print matrix:"
    write(*,'(a,i5,a,i5)') "orginal matrix shape:", a%M, 'x', a%N
    write(*,'(a,i5)') "size of compressed matrix:", a%nnz
    select case ( a%storage_format_id )
    case ( sparsemat_storage_typeid_csr )
      write(*,*) "rowPtr:", a%rowPtr(:)
    end select
    write(*,*) "colInd:", a%colIdx(:)
    write(*,*) "val:"
 
    select case ( a%storage_format_id )
    case ( sparsemat_storage_typeid_csr )
      j1 = a%rowPtr(1)
      do p=1, a%rowPtrSize-1
        j2 = a%rowPtr(p+1)
        row_val(:) = 0.0_rp
        row_val(a%colIdx(j1:j2-1)) = a%val(j1:j2-1)
        write(*,*) row_val(:)
        j1 = j2
      end do
    case ( sparsemat_storage_typeid_ell )
      do p=1, a%M
        row_val(:) = 0.0_rp
        do j1=1, a%col_size
          j2 = p + (j1-1)*a%M
         row_val(a%colIdx(j2)) = a%val(j2)
        end do
        write(*,*) row_val(:)
      end do
    end select
 
    return
  end subroutine sparsemat_print
 
  function sparsemat_get_storage_format_id( A ) result(id)
    implicit none
    class(sparsemat), intent(in) :: a
    real(rp) :: id
    !------------------------------------------------------
 
    id = a%storage_format_id
    return 
  end function sparsemat_get_storage_format_id
 
  !> Matrix-vector multiplication for a sparse matrix
!OCL SERIAL  
  subroutine sparsemat_matmul1(A, b, c)
    implicit none
 
    type(sparsemat), intent(in) :: A
    real(RP), intent(in ) :: b(:)
    real(RP), intent(out) :: c(A%M)
 
    !--------------------------------------------------------------------------- 
 
    !$acc routine vector
 
    select case( a%storage_format_id )
    case( sparsemat_storage_typeid_csr )
      call sparsemat_matmul_csr_1( a%val, a%colIdx, a%rowPtr, b, c, &
        a%M, a%N, a%nnz, a%rowPtrSize                          )
    case( sparsemat_storage_typeid_ell )
      call sparsemat_matmul_ell_1( a%val, a%colIdx, b, c,           &
        a%M, a%N, a%nnz, a%col_size                            )
    end select
 
    return
  end subroutine sparsemat_matmul1
 
  !> Matrix-vector multiplication for a sparse matrix with two vectors
  !! This routine computes c = A * (b1 .* b2), where .* is the element-wise multiplication.
!OCL SERIAL  
  subroutine sparsemat_matmul1_2(A, b1, b2, c)
    implicit none
 
    type(sparsemat), intent(in) :: A
    real(RP), intent(in ) :: b1(:)
    real(RP), intent(in ) :: b2(:)
    real(RP), intent(out) :: c(A%M)
 
    !--------------------------------------------------------------------------- 
 
    !$acc routine vector
 
    select case( a%storage_format_id )
    case( sparsemat_storage_typeid_csr )
      call sparsemat_matmul_csr_1_2( a%val, a%colIdx, a%rowPtr, b1, b2, c, &
        a%M, a%N, a%nnz, a%rowPtrSize                          )
    case( sparsemat_storage_typeid_ell )
      call sparsemat_matmul_ell_1_2( a%val, a%colIdx, b1, b2, c, &
        a%M, a%N, a%nnz, a%col_size                              )
    end select
 
    return
  end subroutine sparsemat_matmul1_2
 
  !> Matrix-matrix multiplication for a sparse matrix
!OCL SERIAL  
  subroutine sparsemat_matmul2(A, b, c)
    implicit none
 
    type(sparsemat), intent(in) :: A
    real(RP), intent(in ) :: b(:,:)
    real(RP), intent(out) :: c(size(b,1),A%M)
 
    !--------------------------------------------------------------------------- 
 
    select case( a%storage_format_id )
    case( sparsemat_storage_typeid_csr )
      call sparsemat_matmul_csr_2( a%val, a%colIdx, a%rowPtr, b, c, &
        a%M, a%N, a%nnz, a%rowPtrSize, size(b,1)               )
    case( sparsemat_storage_typeid_ell )
      call sparsemat_matmul_ell_2( a%val, a%colIdx, b, c,           &
        a%M, a%N, a%nnz, a%col_size, size(b,1)                 )
    end select
 
    return
  end subroutine sparsemat_matmul2
 
!--- private ----------------------------------------------
 
  !> Matrix-vector multiplication for a sparse matrix in CSR format
!OCL SERIAL
  subroutine sparsemat_matmul_csr_1(A, col_Ind, rowPtr, b, c, M, N, buf_size, rowPtr_size)
    implicit none
 
    integer, intent(in) :: M
    integer, intent(in) :: N
    integer, intent(in) :: buf_size
    integer, intent(in) :: rowPtr_size
    real(RP), intent(in) :: A(buf_size)
    integer, intent(in) :: col_Ind(buf_size)
    integer, intent(in) :: rowPtr(rowPtr_size)
    real(RP), intent(in ) :: b(N)
    real(RP), intent(out) :: c(M)
 
    integer :: p, j
    integer :: j1, j2
    real(RP) :: s
 
    !--------------------------------------------------------------------------- 
 
    !$acc routine vector
  
    !call mkl_dcsrgemv( 'N', rowPtr_size-1, A, rowPtr, col_Ind, b, c)
 
    !$acc loop vector
    do p=1, rowptr_size-1
       j1 = rowptr(p)
       j2 = rowptr(p+1) 
       s = 0.0_rp
       do j=j1, j2-1
          s = s + a(j) * b(col_ind(j))
       end do
       c(p) = s
    end do
 
    return
  end subroutine sparsemat_matmul_csr_1
 
  !> Matrix-vector multiplication for a sparse matrix in CSR format with two vectors
  !! This routine computes c = A * (b1 .* b2), where .* is the element-wise multiplication.
!OCL SERIAL
  subroutine sparsemat_matmul_csr_1_2(A, col_Ind, rowPtr, b1, b2, c, M, N, buf_size, rowPtr_size)
    implicit none
 
    integer, intent(in) :: M
    integer, intent(in) :: N
    integer, intent(in) :: buf_size
    integer, intent(in) :: rowPtr_size
    real(RP), intent(in) :: A(buf_size)
    integer, intent(in) :: col_Ind(buf_size)
    integer, intent(in) :: rowPtr(rowPtr_size)
    real(RP), intent(in ) :: b1(N)
    real(RP), intent(in ) :: b2(N)
    real(RP), intent(out) :: c(M)
 
    integer :: p, j
    integer :: j1, j2
    real(RP) :: s
 
    !--------------------------------------------------------------------------- 
 
    !$acc routine vector
  
    !call mkl_dcsrgemv( 'N', rowPtr_size-1, A, rowPtr, col_Ind, b, c)
 
    !$acc loop vector
    do p=1, rowptr_size-1
       j1 = rowptr(p)
       j2 = rowptr(p+1) 
       s = 0.0_rp
       do j=j1, j2-1
          s = s + a(j) * b1(col_ind(j)) * b2(col_ind(j))
       end do
       c(p) = s
    end do
 
    return
  end subroutine sparsemat_matmul_csr_1_2
 
  !> Matrix-matrix multiplication for a sparse matrix in CSR format
!OCL SERIAL
  subroutine sparsemat_matmul_csr_2(A, col_Ind, rowPtr, b, c, M, N, buf_size, rowPtr_size, NQ)
    implicit none
 
    integer, intent(in) :: M
    integer, intent(in) :: N
    integer, intent(in) :: NQ
    integer, intent(in) :: buf_size
    integer, intent(in) :: rowPtr_size
    real(RP), intent(in) :: A(buf_size)
    integer, intent(in) :: col_Ind(buf_size)
    integer, intent(in) :: rowPtr(rowPtr_size)
    real(RP), intent(in ) :: b(NQ,N)
    real(RP), intent(out) :: c(NQ,M)
 
    integer :: p, j
    integer :: j1, j2
 
    !--------------------------------------------------------------------------- 
 
    !call mkl_dcsrgemv( 'N', rowPtr_size-1, A, rowPtr, col_Ind, b, c)
    j1 = rowptr(1)
    c(:,:) = 0.0_rp
    do p=1, rowptr_size-1
      j2 = rowptr(p+1) 
      do j=j1, j2-1
        c(:,p) = c(:,p) + a(j) * b(:,col_ind(j))
      end do
      j1 = j2
    end do
 
    return
  end subroutine sparsemat_matmul_csr_2
 
  !> Matrix-vector multiplication for a sparse matrix in ELL format
!OCL SERIAL
  subroutine sparsemat_matmul_ell_1(A, col_Ind, b, c, M, N, buf_size, col_size)
    implicit none
 
    integer, intent(in) :: M
    integer, intent(in) :: N
    integer, intent(in) :: buf_size
    integer, intent(in) :: col_size
    real(RP), intent(in) :: A(buf_size)
    integer, intent(in) :: col_Ind(buf_size)
    real(RP), intent(in ) :: b(N)
    real(RP), intent(out) :: c(M)
 
    integer :: k, kk, i
    integer :: j_ptr
    !--------------------------------------------------------------------------- 
 
    !$acc routine vector
 
#ifdef _OPENACC
    !$acc loop vector
    do i=1, m
      c(i) = 0.0_rp
    end do
#else
    c(:) = 0.0_rp
#endif
 
    do k=1, col_size
      kk = m * (k-1)
      !$acc loop vector
      do i=1, m
        j_ptr = kk + i        
        c(i) = c(i) + a(j_ptr) * b(col_ind(j_ptr))
      end do
    end do
 
    return
  end subroutine sparsemat_matmul_ell_1
 
  !> Matrix-vector multiplication for a sparse matrix in ELL format with two vectors
  !! This routine computes c = A * (b1 .* b2), where .* is the element-wise multiplication.
!OCL SERIAL
  subroutine sparsemat_matmul_ell_1_2(A, col_Ind, b1, b2, c, M, N, buf_size, col_size)
    implicit none
 
    integer, intent(in) :: M
    integer, intent(in) :: N
    integer, intent(in) :: buf_size
    integer, intent(in) :: col_size
    real(RP), intent(in) :: A(buf_size)
    integer, intent(in) :: col_Ind(buf_size)
    real(RP), intent(in ) :: b1(N)
    real(RP), intent(in ) :: b2(N)
    real(RP), intent(out) :: c(M)
 
    integer :: k, kk, i
    integer :: j_ptr
    !--------------------------------------------------------------------------- 
 
    !$acc routine vector
 
#ifdef _OPENACC
    !$acc loop vector
    do i=1, m
      c(i) = 0.0_rp
    end do
#else
    c(:) = 0.0_rp
#endif
 
    do k=1, col_size
      kk = m * (k-1)
      !$acc loop vector
      do i=1, m
        j_ptr = kk + i        
        c(i) = c(i) + a(j_ptr) * b1(col_ind(j_ptr)) * b2(col_ind(j_ptr))
      end do
    end do
 
    return
  end subroutine sparsemat_matmul_ell_1_2
 
  !> Matrix-matrix multiplication for a sparse matrix in ELL format
!OCL SERIAL
  subroutine sparsemat_matmul_ell_2(A, col_Ind, b, c, M, N, buf_size, col_size, NQ)
    implicit none
 
    integer, intent(in) :: M
    integer, intent(in) :: N
    integer, intent(in) :: buf_size
    integer, intent(in) :: col_size
    integer, intent(in) :: NQ
    real(RP), intent(in) :: A(buf_size)
    integer, intent(in) :: col_Ind(buf_size)
    real(RP), intent(in ) :: b(NQ,N)
    real(RP), intent(out) :: c(NQ,M)
 
    integer :: k, kk, i
    integer :: j_ptr
    !--------------------------------------------------------------------------- 
 
    c(:,:) = 0.0_rp
    do k=1, col_size
      kk = m * (k-1)
      do i=1, m
        j_ptr = kk + i        
        c(:,i) = c(:,i) + a(j_ptr) * b(:,col_ind(j_ptr))
      end do
    end do
 
    return
  end subroutine sparsemat_matmul_ell_2
 
end module scale_sparsemat