sparse_matrix_hec.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_hec
  use hecmw_util
  use m_hecmw_comm_f
  use m_sparse_matrix
  implicit none
 
  private
 
  public :: sparse_matrix_hec_init_prof
  public :: sparse_matrix_hec_set_prof
  public :: sparse_matrix_hec_set_conv_ext
  public :: sparse_matrix_hec_set_vals
  public :: sparse_matrix_hec_set_rhs
  public :: sparse_matrix_hec_get_rhs
 
contains
 
  !!! public subroutines
 
  subroutine sparse_matrix_hec_init_prof(spMAT, hecMAT, hecMESH)
    type (sparse_matrix), intent(inout) :: spmat
    type (hecmwst_matrix), intent(in) :: hecmat
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    integer(kind=kint) :: ndof, ndof2, n_loc, nl, nu, nz
    ndof=hecmat%NDOF; ndof2=ndof*ndof
    n_loc=hecmat%N*ndof
    if (sparse_matrix_is_sym(spmat)) then
      nu=hecmat%indexU(hecmat%N)
      nz=hecmat%N*(ndof2+ndof)/2+nu*ndof2
    else
      nl=hecmat%indexL(hecmat%N)
      nu=hecmat%indexU(hecmat%N)
      nz=(hecmat%N+nu+nl)*ndof2
    endif
    call sparse_matrix_init(spmat, n_loc, nz)
    call sparse_matrix_hec_set_conv_ext(spmat, hecmesh, hecmat%NDOF)
    spmat%iterlog = hecmat%Iarray(21)
    spmat%timelog = hecmat%Iarray(22)
    call sparse_matrix_hec_set_prof(spmat, hecmat)
  end subroutine sparse_matrix_hec_init_prof
 
  subroutine sparse_matrix_hec_set_conv_ext(spMAT, hecMESH, ndof)
    type(sparse_matrix), intent(inout) :: spmat
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    integer(kind=kint), intent(in) :: ndof
    integer(kind=kint) :: i,j,nn_external,id,ierr
    integer(kint)   :: pid,lid,j0
    if (hecmesh%n_neighbor_pe==0) return
    ! create conversion list
    nn_external = hecmesh%n_node - hecmesh%nn_internal
    allocate(spmat%conv_ext(nn_external*ndof), stat=ierr)
    if (ierr /= 0) then
      write(*,*) " Allocation error, spMAT%conv_ext"
      call hecmw_abort(hecmw_comm_get_comm())
    endif
    spmat%conv_ext(:) = -1
    do i=1,nn_external
      id = i + hecmesh%nn_internal
      pid = hecmesh%node_ID(id*2)
      lid = hecmesh%node_ID(id*2-1)
      j0 = spmat%DISPLS(pid+1) + (lid-1)*ndof
      do j=1,ndof
        spmat%conv_ext((i-1)*ndof+j) = j0+j
      enddo
    enddo
  end subroutine sparse_matrix_hec_set_conv_ext
 
  subroutine sparse_matrix_hec_set_prof(spMAT, hecMAT)
    type(sparse_matrix), intent(inout) :: spmat
    type(hecmwst_matrix), intent(in) :: hecmat
    integer(kind=kint) :: ndof, ndof2
    integer(kind=kint) :: m, i, idof, i0, ii, ls, le, l, j, j0, jdof, jdofs
    !integer(kind=kint) :: offset_l, offset_d, offset_u
    ! CONVERT TO CSR or COO STYLE
    ndof=hecmat%NDOF; ndof2=ndof*ndof
    m=1
    ii = 0
    do i=1,hecmat%N
      do idof=1,ndof
        i0=spmat%offset+ndof*(i-1)
        ii=i0+idof
        if (spmat%type==sparse_matrix_type_csr) spmat%IRN(ii-spmat%offset)=m
        ! Lower
        if (.not. sparse_matrix_is_sym(spmat)) then
          ls=hecmat%indexL(i-1)+1
          le=hecmat%indexL(i)
          do l=ls,le
            j=hecmat%itemL(l)
            !if (j <= hecMAT%N) then
            j0=spmat%offset+ndof*(j-1)
            !else
            !   j0=spMAT%conv_ext(ndof*(j-hecMAT%N))-ndof
            !endif
            !offset_l=ndof2*(l-1)+ndof*(idof-1)
            do jdof=1,ndof
              if (spmat%type==sparse_matrix_type_coo) spmat%IRN(m)=ii
              spmat%JCN(m)=j0+jdof
              !spMAT%A(m)=hecMAT%AL(offset_l+jdof)
              m=m+1
            enddo
          enddo
        endif
        ! Diag
        !offset_d=ndof2*(i-1)+ndof*(idof-1)
        if (sparse_matrix_is_sym(spmat)) then; jdofs=idof; else; jdofs=1; endif
        do jdof=jdofs,ndof
          if (spmat%type==sparse_matrix_type_coo) spmat%IRN(m)=ii
          spmat%JCN(m)=i0+jdof
          !spMAT%A(m)=hecMAT%D(offset_d+jdof)
          m=m+1
        enddo
        ! Upper
        ls=hecmat%indexU(i-1)+1
        le=hecmat%indexU(i)
        do l=ls,le
          j=hecmat%itemU(l)
          if (j <= hecmat%N) then
            j0=spmat%offset+ndof*(j-1)
          else
            j0=spmat%conv_ext(ndof*(j-hecmat%N))-ndof
            if (sparse_matrix_is_sym(spmat) .and. j0 < i0) cycle
          endif
          !offset_u=ndof2*(l-1)+ndof*(idof-1)
          do jdof=1,ndof
            if (spmat%type==sparse_matrix_type_coo) spmat%IRN(m)=ii
            spmat%JCN(m)=j0+jdof
            !spMAT%A(m)=hecMAT%AU(offset_u+jdof)
            m=m+1
          enddo
        enddo
      enddo
    enddo
    if (spmat%type == sparse_matrix_type_csr) spmat%IRN(ii+1-spmat%offset)=m
    if (sparse_matrix_is_sym(spmat) .and. m-1 < spmat%NZ) spmat%NZ=m-1
    if (m-1 /= spmat%NZ) then
      write(*,*) 'ERROR: sparse_matrix_set_ij on rank ',hecmw_comm_get_rank()
      write(*,*) 'm-1 = ',m-1,', NZ=',spmat%NZ
      stop
    endif
  end subroutine sparse_matrix_hec_set_prof
 
  subroutine sparse_matrix_hec_set_vals(spMAT, hecMAT)
    type(sparse_matrix), intent(inout) :: spmat
    type(hecmwst_matrix), intent(in) :: hecmat
    integer(kind=kint) :: ndof, ndof2
    integer(kind=kint) :: m, i, idof, i0, ii, ls, le, l, j, j0, jdof, jdofs
    integer(kind=kint) :: offset_l, offset_d, offset_u
    ndof=hecmat%NDOF; ndof2=ndof*ndof
    ii = 0
    m=1
    do i=1,hecmat%N
      do idof=1,ndof
        i0=spmat%offset+ndof*(i-1)
        ii=i0+idof
        if (spmat%type == sparse_matrix_type_csr) then
          if (spmat%IRN(ii-spmat%offset)/=m) stop "ERROR: sparse_matrix_set_a1"
        endif
        ! Lower
        if (.not. sparse_matrix_is_sym(spmat)) then
          ls=hecmat%indexL(i-1)+1
          le=hecmat%indexL(i)
          do l=ls,le
            j=hecmat%itemL(l)
            !if (j <= hecMAT%N) then
            j0=spmat%offset+ndof*(j-1)
            !else
            !   j0=spMAT%conv_ext(ndof*(j-hecMAT%N))-ndof
            !endif
            offset_l=ndof2*(l-1)+ndof*(idof-1)
            do jdof=1,ndof
              if (spmat%type==sparse_matrix_type_coo) then
                if (spmat%IRN(m)/=ii) stop "ERROR: sparse_matrix_set_a2"
              endif
              if (spmat%type == sparse_matrix_type_csr) spmat%JCN(m)=j0+jdof
              if (spmat%JCN(m)/=j0+jdof) stop "ERROR: sparse_matrix_set_a3"
              spmat%A(m)=hecmat%AL(offset_l+jdof)
              m=m+1
            enddo
          enddo
        endif
        ! Diag
        offset_d=ndof2*(i-1)+ndof*(idof-1)
        if (sparse_matrix_is_sym(spmat)) then; jdofs=idof; else; jdofs=1; endif
        do jdof=jdofs,ndof
          if (spmat%type==sparse_matrix_type_coo) then
            if (spmat%IRN(m)/=ii) stop "ERROR: sparse_matrix_set_a4"
          endif
          if (spmat%type == sparse_matrix_type_csr) spmat%JCN(m)=i0+jdof
          if (spmat%JCN(m)/=i0+jdof) stop "ERROR: sparse_matrix_set_a5"
          spmat%A(m)=hecmat%D(offset_d+jdof)
          m=m+1
        enddo
        ! Upper
        ls=hecmat%indexU(i-1)+1
        le=hecmat%indexU(i)
        do l=ls,le
          j=hecmat%itemU(l)
          if (j <= hecmat%N) then
            j0=spmat%offset+ndof*(j-1)
          else
            j0=spmat%conv_ext(ndof*(j-hecmat%N))-ndof
            if (sparse_matrix_is_sym(spmat) .and. j0 < i0) cycle
          endif
          offset_u=ndof2*(l-1)+ndof*(idof-1)
          do jdof=1,ndof
            if (spmat%type==sparse_matrix_type_coo) then
              if (spmat%IRN(m)/=ii) stop "ERROR: sparse_matrix_set_a6"
            endif
            if (spmat%type == sparse_matrix_type_csr) spmat%JCN(m)=j0+jdof
            if (spmat%JCN(m)/=j0+jdof) stop "ERROR: sparse_matrix_set_a7"
            spmat%A(m)=hecmat%AU(offset_u+jdof)
            m=m+1
          enddo
        enddo
      enddo
    enddo
    if (spmat%type == sparse_matrix_type_csr) then
      if (spmat%IRN(ii+1-spmat%offset)/=m) stop "ERROR: sparse_matrix_set_a8"
    endif
    if (m-1 /= spmat%NZ) stop "ERROR: sparse_matrix_set_a9"
  end subroutine sparse_matrix_hec_set_vals
 
  subroutine sparse_matrix_hec_set_rhs(spMAT, hecMAT)
    implicit none
    type (sparse_matrix), intent(inout) :: spmat
    type (hecmwst_matrix), intent(in) :: hecmat
    integer(kind=kint) :: ierr,i
    allocate(spmat%rhs(spmat%N_loc), stat=ierr)
    if (ierr /= 0) then
      write(*,*) " Allocation error, spMAT%rhs"
      call hecmw_abort(hecmw_comm_get_comm())
    endif
    do i=1,spmat%N_loc
      spmat%rhs(i)=hecmat%b(i)
    enddo
  end subroutine sparse_matrix_hec_set_rhs
 
  subroutine sparse_matrix_hec_get_rhs(spMAT, hecMAT)
    implicit none
    type (sparse_matrix), intent(inout) :: spmat
    type (hecmwst_matrix), intent(inout) :: hecmat
    integer(kind=kint) :: i
    do i=1,spmat%N_loc
      hecmat%x(i)=spmat%rhs(i)
    enddo
    deallocate(spmat%rhs)
  end subroutine sparse_matrix_hec_get_rhs
 
end module m_sparse_matrix_hec