hecmw_precond_SAAMG_adapt.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2026 FrontISTR Commons
! This software is released under the MIT License, see License.txt
!-------------------------------------------------------------------------------
module hecmw_precond_saamg_adapt
  use hecmw_util
  use hecmw_precond_saamg_matrix, only: hecmwst_saamg_bcsr, hecmw_saamg_bcsr_free
  use hecmw_precond_saamg_comm,   only: hecmwst_saamg_comm, hecmw_saamg_comm_free, &
       hecmw_saamg_check_alloc
  implicit none
 
  private
  public :: hecmw_saamg_from_hecmat
  public :: hecmw_saamg_comm_from_mesh
 
contains
 
  subroutine hecmw_saamg_from_hecmat(hecMAT, A, include_halo)
    implicit none
    type(hecmwst_matrix),  intent(in)  :: hecmat
    type(hecmwst_saamg_bcsr), intent(out) :: a
    logical, optional,     intent(in)  :: include_halo
    integer(kind=kint) :: ndof, ndof2, nnode, ncolnode, i, ii, kk, j, jn, pos, cnt, b0, astat
    logical :: halo
 
    halo = .false.
    if (present(include_halo)) halo = include_halo
 
    ndof  = hecmat%NDOF
    ndof2 = ndof*ndof
    nnode = hecmat%N                     ! internal rows
    if (halo) then
      ncolnode = hecmat%NP               ! columns reference internal + halo
    else
      ncolnode = hecmat%N
    end if
 
    ! Build the block-CSR directly from the hecMAT BCRS (diagonal D + lower AL +
    ! upper AU, one nb x nb block per neighbor) WITHOUT an intermediate triplet
    ! buffer -- avoids holding the operator twice (critical for large meshes).
    ! hecMAT stores each block row-major ((ii-1)*ndof+kk); transpose to col-major.
    call hecmw_saamg_bcsr_free(a)
    a%n = nnode*ndof; a%ncol = ncolnode*ndof; a%nb = ndof; a%mb = ndof
    a%nbrow = nnode; a%nbcol = ncolnode
    allocate(a%browptr(nnode+1))
 
    ! pass 1: count blocks per row (diagonal + kept L + kept U)
    a%browptr(1) = 1
    do i = 1, nnode
      cnt = 1                                       ! diagonal block
      do j = hecmat%indexL(i-1)+1, hecmat%indexL(i)
        if (halo .or. hecmat%itemL(j) <= nnode) cnt = cnt + 1
      end do
      do j = hecmat%indexU(i-1)+1, hecmat%indexU(i)
        if (halo .or. hecmat%itemU(j) <= nnode) cnt = cnt + 1
      end do
      a%browptr(i+1) = a%browptr(i) + cnt
    end do
    a%nnzb = a%browptr(nnode+1) - 1
    a%nnz  = ndof2 * a%nnzb
    allocate(a%bcol(a%nnzb), a%bval(ndof2*a%nnzb), stat=astat)
    call hecmw_saamg_check_alloc(astat, operator'from_hecmat finest  (bcol/bval)')
 
    ! pass 2: fill blocks in place (diagonal first, then L, then U)
    do i = 1, nnode
      pos = a%browptr(i)
      a%bcol(pos) = i; b0 = (pos-1)*ndof2          ! diagonal block
      do ii = 1, ndof
        do kk = 1, ndof
          a%bval(b0 + (kk-1)*ndof+ii) = hecmat%D(ndof2*(i-1) + (ii-1)*ndof + kk)
        end do
      end do
      pos = pos + 1
      do j = hecmat%indexL(i-1)+1, hecmat%indexL(i)
        jn = hecmat%itemL(j)
        if (.not. halo .and. jn > nnode) cycle
        a%bcol(pos) = jn; b0 = (pos-1)*ndof2
        do ii = 1, ndof
          do kk = 1, ndof
            a%bval(b0 + (kk-1)*ndof+ii) = hecmat%AL(ndof2*(j-1) + (ii-1)*ndof + kk)
          end do
        end do
        pos = pos + 1
      end do
      do j = hecmat%indexU(i-1)+1, hecmat%indexU(i)
        jn = hecmat%itemU(j)
        if (.not. halo .and. jn > nnode) cycle
        a%bcol(pos) = jn; b0 = (pos-1)*ndof2
        do ii = 1, ndof
          do kk = 1, ndof
            a%bval(b0 + (kk-1)*ndof+ii) = hecmat%AU(ndof2*(j-1) + (ii-1)*ndof + kk)
          end do
        end do
        pos = pos + 1
      end do
    end do
  end subroutine hecmw_saamg_from_hecmat
 
  subroutine hecmw_saamg_comm_from_mesh(hecMESH, ndof, cmt)
    implicit none
    type(hecmwst_local_mesh), intent(in)  :: hecmesh
    integer(kind=kint),       intent(in)  :: ndof
    type(hecmwst_saamg_comm),   intent(out) :: cmt
    integer(kind=kint) :: nnb, nimp, nexp
 
    call hecmw_saamg_comm_free(cmt)
    cmt%comm       = hecmesh%MPI_COMM
    cmt%my_rank    = hecmesh%my_rank
    cmt%nint       = hecmesh%nn_internal
    cmt%nnode      = hecmesh%n_node
    cmt%nb         = ndof
    cmt%n_neighbor = hecmesh%n_neighbor_pe
    ! global node id per local node (for partition-invariant Lanczos seeds)
    allocate(cmt%gnode(hecmesh%n_node))
    cmt%gnode(1:hecmesh%n_node) = hecmesh%global_node_ID(1:hecmesh%n_node)
    nnb = cmt%n_neighbor
    if (nnb == 0) return
 
    nimp = hecmesh%import_index(nnb)
    nexp = hecmesh%export_index(nnb)
    allocate(cmt%neighbor(nnb))
    allocate(cmt%import_index(0:nnb), cmt%import_item(nimp))
    allocate(cmt%export_index(0:nnb), cmt%export_item(nexp))
    cmt%neighbor(1:nnb)        = hecmesh%neighbor_pe(1:nnb)
    cmt%import_index(0:nnb)    = hecmesh%import_index(0:nnb)
    cmt%export_index(0:nnb)    = hecmesh%export_index(0:nnb)
    if (nimp > 0) cmt%import_item(1:nimp) = hecmesh%import_item(1:nimp)
    if (nexp > 0) cmt%export_item(1:nexp) = hecmesh%export_item(1:nexp)
  end subroutine hecmw_saamg_comm_from_mesh
 
end module hecmw_precond_saamg_adapt