hecmw_solver_direct_mf.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_solver_direct_mf
  use hecmw_util
  use m_hecmw_comm_f
  use hecmw_matrix_misc
  use hecmw_ordering
  use hecmw_mf_graph
  use hecmw_mf_symbolic
  use hecmw_mf_dist
  use hecmw_mf_numeric
  use hecmw_mf_kernel
  use hecmw_solver_las
  use hecmw_solver_misc
  use hecmw_solver_direct
  !$ use omp_lib
  implicit none
 
  private
  public :: hecmw_solve_direct_mf
  public :: hecmw_mf_compare_builtin
 
  integer(kind=kint), parameter :: MF_RELAX_COLS = 16
  integer(kind=kint), parameter :: MF_TILE = 256
  integer(kind=kint), parameter :: MF_IR_STEPS = 5
  real(kind=kreal), parameter :: mf_ir_tol = 1.0d-8
 
  type(hecmwST_mf_symbolic), save :: SYM
  type(hecmwST_mf_factor), save :: FCT
  type(hecmwST_mf_map), save :: MAP
  type(hecmwST_mf_gmat), save :: GMAT
 
contains
 
  subroutine hecmw_solve_direct_mf(hecMESH, hecMAT, imsg)
    implicit none
    type(hecmwst_local_mesh), intent(inout) :: hecmesh
    type(hecmwst_matrix), intent(inout) :: hecmat
    integer(kind=kint), intent(in) :: imsg
    type(hecmwst_mf_graph) :: graph
    integer(kind=kint), allocatable :: perm(:), invp(:)
    integer(kind=kint) :: loglevel, iterlog, timelog, ordering, n, nerr, relax, tile, ierr, idof, node, nthreads, irmax, it, eta
    real(kind=kreal) :: t1, t2, irtol, bnrm, rnrm, rprev, tcomm
    real(kind=kreal), allocatable :: rr(:), dd(:), xp(:)
    logical :: clustered
 
    ! TIMELOG holds the stage times, ITERLOG the refinement residuals, LOGLEVEL the
    ! diagnostics; an unset LOGLEVEL (-1) means no diagnostics here, with no fallback
    ! to the other two flags
    loglevel = hecmw_mat_get_loglevel(hecmat)
    if (loglevel < 0) loglevel = 0
    iterlog = hecmw_mat_get_iterlog(hecmat)
    timelog = hecmw_mat_get_timelog(hecmat)
    if (hecmw_comm_get_rank() /= 0) then
      loglevel = 0
      iterlog = 0
      timelog = 0
    endif
 
    if (hecmesh%PETOT > 1) then
      call hecmw_solve_direct_mf_dist(hecmesh, hecmat, imsg, loglevel, iterlog, timelog)
      return
    endif
 
    if (hecmat%Iarray(98) == 1) then
      t1 = hecmw_wtime()
      call hecmw_mf_graph_from_hecmat(hecmat, graph)
      n = graph%nnode
      allocate(perm(n), invp(n))
      ordering = hecmat%Iarray(41)
      call hecmw_ordering_gen(n, n + (graph%xadj(n+1)-1)/2, graph%xadj, graph%adjncy, perm, invp, ordering, loglevel)
      relax = hecmat%Iarray(46)
      if (relax <= 0) relax = mf_relax_cols
      relax = max(1, relax / hecmat%NDOF)
      call hecmw_mf_symbolic_finalize(sym)
      call hecmw_mf_symbolic_build(graph, perm, relax, sym)
      clustered = hecmat%Iarray(43) /= 0 .and. hecmw_mf_kernel_blr_available()
      if (clustered) then
        call hecmw_mf_symbolic_cluster(graph, sym, perm)
        call hecmw_mf_symbolic_finalize(sym)
        call hecmw_mf_symbolic_build(graph, perm, relax, sym)
      endif
      tile = hecmat%Iarray(45)
      if (tile <= 0) tile = mf_tile
      eta = hecmat%Iarray(47)
      if (eta == 0) eta = mf_blr_eta
      if (eta < 0) eta = 0
      call hecmw_mf_numeric_finalize(fct)
      call hecmw_mf_numeric_init(graph, sym, tile, clustered, eta, fct)
      t2 = hecmw_wtime()
      if (timelog > 0) write(*,'(a,f10.3,a)') '[DIRECTmf]: symbolic fct done (', t2 - t1, ' sec)'
      if (loglevel > 0) then
        if (clustered) write(*,'(a)') '[DIRECTmf]: separator nodes clustered for BLR'
        call hecmw_mf_symbolic_print(sym)
        call hecmw_mf_numeric_print(fct)
      endif
      if (loglevel > 1) then
        call hecmw_mf_symbolic_check(graph, sym, nerr)
        write(*,'(a,i0)') '[DIRECTmf]: self check violations = ', nerr
        if (.not. clustered) then
          call hecmw_mf_compare_builtin(hecmesh, hecmat, graph, ordering, nerr)
          write(*,'(a,i0)') '[DIRECTmf]: mismatches against built-in DIRECT = ', nerr
        endif
      endif
      deallocate(perm, invp)
      call hecmw_mf_graph_finalize(graph)
    endif
 
    hecmat%Iarray(98) = 0
 
    if (hecmat%Iarray(97) == 1) then
      fct%mode = hecmat%Iarray(42)
      if (fct%mode < 0 .or. fct%mode > 2) then
        if (hecmw_comm_get_rank() == 0) write(*,'(a,i0,a)') &
          '[DIRECTmf]: WARNING: MF_Mode ', fct%mode, ' is outside 0-2; using the automatic mode'
        fct%mode = 0
      endif
      fct%scan = loglevel > 1
      fct%blr = hecmat%Iarray(43) /= 0
      fct%eps = hecmat%Rarray(41)
      fct%pivot_u = hecmat%Rarray(43)
      ! threshold partial pivoting is defined for u in (0, 0.5]
      if (fct%pivot_u > 0.5d0) then
        if (hecmw_comm_get_rank() == 0) write(*,'(a,1pe9.2,a)') &
          '[DIRECTmf]: WARNING: MF_PivotThreshold ', fct%pivot_u, ' is above 0.5; using 0.5'
        fct%pivot_u = 0.5d0
      endif
      fct%pivot_zero = hecmat%Rarray(44)
      fct%blr_beta = hecmat%Rarray(45)
      if (fct%blr_beta <= 0.0d0) fct%blr_beta = 1.0d0
      fct%blr_reuse = hecmat%Iarray(48) /= 0
      if (fct%blr .and. .not. hecmw_mf_kernel_blr_available()) then
        fct%blr = .false.
        ! the user asked for BLR, so the fallback is reported regardless of LOGLEVEL
        if (hecmw_comm_get_rank() == 0) write(*,'(a)') '[DIRECTmf]: BLR disabled (built without LAPACK)'
      endif
      t1 = hecmw_wtime()
      call hecmw_mf_numeric_factor(hecmat, sym, fct, ierr)
      t2 = hecmw_wtime()
      if (ierr /= 0) then
        if (ierr > 0) then
          idof = fct%pdof(ierr)
          ! the node prints as the user's id where the matrix rows are the mesh nodes
          ! (not the expanded saddle-point or coarse-grid systems)
          node = (idof-1)/hecmat%NDOF + 1
          if (hecmat%NP == hecmesh%n_node) node = hecmesh%global_node_ID(node)
          write(imsg,'(a,i0,a,i0,a)') 'ERROR: DIRECTmf: no acceptable pivot at node ', node, &
            ' dof ', mod(idof-1, hecmat%NDOF) + 1, ' (the matrix may be singular)'
          write(*,'(a,i0,a,i0,a)') 'ERROR: DIRECTmf: no acceptable pivot at node ', node, &
            ' dof ', mod(idof-1, hecmat%NDOF) + 1, ' (the matrix may be singular)'
        else if (ierr == -1) then
          write(imsg,*) 'ERROR: DIRECTmf: block size of the matrix does not match the symbolic structure'
          write(*,*) 'ERROR: DIRECTmf: block size of the matrix does not match the symbolic structure'
        else
          write(imsg,*) 'ERROR: DIRECTmf: the nonzero structure of the matrix is not symmetric'
          write(*,*) 'ERROR: DIRECTmf: the nonzero structure of the matrix is not symmetric'
        endif
        call hecmw_abort(hecmw_comm_get_comm())
      endif
      hecmat%Iarray(97) = 0
      if (timelog > 0) write(*,'(a,f10.3,a)') '[DIRECTmf]: numeric fct done (', t2 - t1, ' sec)'
      if (loglevel > 0) then
        nthreads = 1
        !$ nthreads = omp_get_max_threads()
        write(*,'(a,i0)') '[DIRECTmf]: threads = ', nthreads
        write(*,'(a,1pe9.2,a,1pe9.2)') '[DIRECTmf]: options: pivot u = ', fct%pivot_u, &
          ', zero pivot factor = ', fct%pivot_zero
        if (fct%lu) then
          if (fct%mode == 2) then
            write(*,'(a)') '[DIRECTmf]: mode = LU (forced by option)'
          else
            write(*,'(a)') '[DIRECTmf]: mode = LU (unsymmetric matrix)'
          endif
          if (loglevel > 1) write(*,'(a,1pe9.2)') '[DIRECTmf]: asymmetry = ', fct%asym
          write(*,'(a,i0,a,i0,a,i0,a)') '[DIRECTmf]: row swaps = ', fct%n_swap, ', delayed = ', fct%n_delay, &
            ' (max front growth = ', fct%max_growth, ' dofs)'
        else
          if (fct%mode == 1) then
            write(*,'(a)') '[DIRECTmf]: mode = LDLt (forced by option)'
          else
            write(*,'(a)') '[DIRECTmf]: mode = LDLt'
          endif
          if (loglevel > 1) write(*,'(a,1pe9.2)') '[DIRECTmf]: asymmetry = ', fct%asym
          write(*,'(a,i0,a,i0,a,i0,a,i0,a,i0,a,i0,a)') '[DIRECTmf]: inertia (+,-) = (', fct%n_pos, ',', fct%n_neg, &
            '), 2x2 pivots = ', fct%n_2x2, ', swaps = ', fct%n_swap, ', delayed = ', fct%n_delay, &
            ' (max front growth = ', fct%max_growth, ' dofs)'
        endif
        write(*,'(a,i0,a,f10.3,a,i0,a,i0)') '[DIRECTmf]: factor words = ', fct%factor_words_act, ' (', &
          real(fct%factor_words_act, kind=kreal)*8.0d0/1024.0d0**3, ' GB), stack peak words = ', fct%stack_peak_act, &
          ', front words = ', fct%front_words_act
        if (fct%blr) then
          write(*,'(a,1pe9.2,a,i0,a,0pf6.1,a)') '[DIRECTmf]: BLR: eps = ', fct%eps, &
            ', factor words full rank = ', fct%blr_words_fr, ' (', &
            100.0d0*(1.0d0 - real(fct%factor_words_act, kind=kreal)/max(real(fct%blr_words_fr, kind=kreal), 1.0d0)), &
            '% saved)'
          write(*,'(a,i0,a,f5.2,a,l1)') '[DIRECTmf]: BLR: skip eta = ', fct%blr_eta, &
            ', gain cap beta = ', fct%blr_beta, ', rank reuse = ', fct%blr_reuse
          write(*,'(a,i0,a,i0,a,i0,a,i0,a,f8.1)') '[DIRECTmf]: BLR: tiles compressed = ', fct%blr_tiles_lr, ' of ', &
            fct%blr_tiles, ', skipped = ', fct%blr_tiles_skip, ', rank max = ', fct%blr_rank_max, ', rank avg = ', &
            real(fct%blr_rank_sum, kind=kreal)/max(real(fct%blr_tiles_lr, kind=kreal), 1.0d0)
        endif
        write(*,'(a,i0,a,f10.3,a)') '[DIRECTmf]: front words held concurrently (peak) = ', fct%front_peak, ' (', &
          real(fct%front_peak, kind=kreal)*8.0d0/1024.0d0**3, ' GB)'
      endif
    endif
 
    if (.not. fct%factored) then
      write(imsg,*) 'ERROR: DIRECTmf: numeric factorization not performed'
      write(*,*) 'ERROR: DIRECTmf: numeric factorization not performed'
      call hecmw_abort(hecmw_comm_get_comm())
    endif
    t1 = hecmw_wtime()
    call hecmw_mf_numeric_solve(sym, fct, hecmat%B, hecmat%X)
    t2 = hecmw_wtime()
    if (timelog > 0) write(*,'(a,f10.3,a)') '[DIRECTmf]: solve done (', t2 - t1, ' sec)'
 
    irmax = hecmat%Iarray(44)
    if (fct%blr .and. irmax == 0) irmax = mf_ir_steps
    if (irmax > 0) then
      irtol = hecmat%Rarray(42)
      if (.not. (irtol > 0.0d0)) irtol = mf_ir_tol
      n = hecmat%NP * hecmat%NDOF
      allocate(rr(n), dd(n), xp(n))
      tcomm = 0.0d0
      call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, hecmat%B, hecmat%B, bnrm, tcomm)
      t1 = hecmw_wtime()
      rprev = huge(rprev)
      do it = 1, irmax
        call hecmw_matresid(hecmesh, hecmat, hecmat%X, hecmat%B, rr, tcomm)
        call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, rr, rr, rnrm, tcomm)
        rnrm = sqrt(rnrm / max(bnrm, tiny(bnrm)))
        if (iterlog > 0) write(*,'(a,i0,a,1pe11.4)') '[DIRECTmf]: refinement ', it - 1, ': residual = ', rnrm
        if (rnrm <= irtol) exit
        ! a step that grew the residual would hand back a worse solution than it got,
        ! so the step is undone and the sweep stops with a warning
        if (rnrm > rprev) then
          hecmat%X(1:n) = xp(1:n)
          if (hecmw_comm_get_rank() == 0) write(*,'(a,1pe11.4)') &
            '[DIRECTmf]: WARNING: refinement diverged, reverted to residual = ', rprev
          exit
        endif
        rprev = rnrm
        xp(1:n) = hecmat%X(1:n)
        call hecmw_mf_numeric_solve(sym, fct, rr, dd)
        hecmat%X(1:n) = hecmat%X(1:n) + dd(1:n)
      enddo
      t2 = hecmw_wtime()
      if (timelog > 0) write(*,'(a,f10.3,a)') '[DIRECTmf]: refinement done (', t2 - t1, ' sec)'
      deallocate(rr, dd, xp)
    else if (irmax == 0 .and. iterlog > 0) then
      ! a negative irmax asks for no residual work at all (the contact wrapper measures
      ! its saddle-point system itself and disables this path with -1)
      n = hecmat%NP * hecmat%NDOF
      allocate(rr(n))
      tcomm = 0.0d0
      call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, hecmat%B, hecmat%B, bnrm, tcomm)
      call hecmw_matresid(hecmesh, hecmat, hecmat%X, hecmat%B, rr, tcomm)
      call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, rr, rr, rnrm, tcomm)
      rnrm = sqrt(rnrm / max(bnrm, tiny(bnrm)))
      write(*,'(a,1pe11.4)') '[DIRECTmf]: residual = ', rnrm
      deallocate(rr)
    endif
  end subroutine hecmw_solve_direct_mf
 
  subroutine hecmw_solve_direct_mf_dist(hecMESH, hecMAT, imsg, loglevel, iterlog, timelog)
    implicit none
    type(hecmwst_local_mesh), intent(inout) :: hecmesh
    type(hecmwst_matrix), intent(inout) :: hecmat
    integer(kind=kint), intent(in) :: imsg, loglevel, iterlog, timelog
    type(hecmwst_mf_graph) :: graph
    integer(kind=kint), allocatable :: perm(:), invp(:)
    real(kind=kreal), allocatable :: gb(:), gx(:), rr(:), wpr(:), wmr(:), xp(:)
    integer(kind=kint) :: ordering, n, nerr, relax, tile, ierr, idof, nthreads, irmax, it, i, ofs, eta, gnode, d
    integer(kind=8) :: wrepl
    real(kind=kreal) :: t1, t2, irtol, bnrm, rnrm, rprev, tcomm
    logical :: clustered
 
    clustered = .false.
    if (hecmat%Iarray(98) == 1) then
      t1 = hecmw_wtime()
      call hecmw_mf_dist_gmat_build(hecmesh, hecmat, gmat)
      call hecmw_mf_symbolic_finalize(sym)
      relax = hecmat%Iarray(46)
      if (relax <= 0) relax = mf_relax_cols
      relax = max(1, relax / hecmat%NDOF)
      clustered = hecmat%Iarray(43) /= 0 .and. hecmw_mf_kernel_blr_available()
      if (hecmw_comm_get_rank() == 0) then
        call hecmw_mf_graph_from_hecmat(gmat%mat, graph)
        n = graph%nnode
        allocate(perm(n), invp(n))
        ordering = hecmat%Iarray(41)
        call hecmw_ordering_gen(n, n + (graph%xadj(n+1)-1)/2, graph%xadj, graph%adjncy, perm, invp, ordering, loglevel)
        call hecmw_mf_symbolic_build(graph, perm, relax, sym)
        if (clustered) then
          call hecmw_mf_symbolic_cluster(graph, sym, perm)
          call hecmw_mf_symbolic_finalize(sym)
          call hecmw_mf_symbolic_build(graph, perm, relax, sym)
        endif
        if (loglevel > 1) then
          call hecmw_mf_symbolic_check(graph, sym, nerr)
          write(*,'(a,i0)') '[DIRECTmf]: self check violations = ', nerr
        endif
        deallocate(perm, invp)
      endif
      call hecmw_mf_dist_symbolic_bcast(sym, 0)
      call hecmw_mf_dist_map_finalize(map)
      call hecmw_mf_dist_map_build(sym, map)
      call hecmw_mf_dist_gmat_part(sym, map, gmat)
      allocate(wmr(map%nprocs))
      call mf_gather_matwords(gmat, wmr, wrepl)
      tile = hecmat%Iarray(45)
      if (tile <= 0) tile = mf_tile
      eta = hecmat%Iarray(47)
      if (eta == 0) eta = mf_blr_eta
      if (eta < 0) eta = 0
      call hecmw_mf_numeric_finalize(fct)
      call hecmw_mf_numeric_init(graph, sym, tile, clustered, eta, fct)
      call hecmw_mf_graph_finalize(graph)
      if (clustered .and. eta > 0) call hecmw_mf_numeric_admis_bcast(fct, 0)
      t2 = hecmw_wtime()
      if (timelog > 0) write(*,'(a,f10.3,a)') '[DIRECTmf]: symbolic fct done (', t2 - t1, ' sec)'
      if (loglevel > 0) then
        if (clustered) write(*,'(a)') '[DIRECTmf]: separator nodes clustered for BLR'
        write(*,'(a,i0,a,i0)') '[DIRECTmf]: MPI ranks = ', map%nprocs, ', upper fronts = ', map%nupper
        write(*,'(a,i0,a,*(i0,1x))') '[DIRECTmf]: matrix words replicated = ', wrepl, ', held per rank = ', &
          (nint(wmr(i), kind=8), i = 1, map%nprocs)
        call hecmw_mf_symbolic_print(sym)
        call hecmw_mf_numeric_print(fct)
      endif
      deallocate(wmr)
      hecmat%Iarray(98) = 0
    endif
 
    if (hecmat%Iarray(97) == 1) then
      fct%mode = hecmat%Iarray(42)
      if (fct%mode < 0 .or. fct%mode > 2) then
        if (hecmw_comm_get_rank() == 0) write(*,'(a,i0,a)') &
          '[DIRECTmf]: WARNING: MF_Mode ', fct%mode, ' is outside 0-2; using the automatic mode'
        fct%mode = 0
      endif
      fct%scan = loglevel > 1
      fct%blr = hecmat%Iarray(43) /= 0
      fct%eps = hecmat%Rarray(41)
      fct%pivot_u = hecmat%Rarray(43)
      ! threshold partial pivoting is defined for u in (0, 0.5]
      if (fct%pivot_u > 0.5d0) then
        if (hecmw_comm_get_rank() == 0) write(*,'(a,1pe9.2,a)') &
          '[DIRECTmf]: WARNING: MF_PivotThreshold ', fct%pivot_u, ' is above 0.5; using 0.5'
        fct%pivot_u = 0.5d0
      endif
      fct%pivot_zero = hecmat%Rarray(44)
      fct%blr_beta = hecmat%Rarray(45)
      if (fct%blr_beta <= 0.0d0) fct%blr_beta = 1.0d0
      fct%blr_reuse = hecmat%Iarray(48) /= 0
      if (fct%blr .and. .not. hecmw_mf_kernel_blr_available()) then
        fct%blr = .false.
        ! the user asked for BLR, so the fallback is reported regardless of LOGLEVEL
        if (hecmw_comm_get_rank() == 0) write(*,'(a)') '[DIRECTmf]: BLR disabled (built without LAPACK)'
      endif
      call hecmw_mf_dist_gmat_vals(hecmat, gmat)
      t1 = hecmw_wtime()
      call hecmw_mf_numeric_factor_mpi(gmat%mat, sym, map, fct, ierr)
      t2 = hecmw_wtime()
      if (ierr /= 0) then
        if (hecmw_comm_get_rank() == 0) then
          if (ierr > 0) then
            idof = fct%pdof(ierr)
            ! the row is in the gathered numbering; gid maps it back to the user's node id
            write(imsg,'(a,i0,a,i0,a)') 'ERROR: DIRECTmf: no acceptable pivot at node ', &
              gmat%gid((idof-1)/hecmat%NDOF + 1), &
              ' dof ', mod(idof-1, hecmat%NDOF) + 1, ' (the matrix may be singular)'
            write(*,'(a,i0,a,i0,a)') 'ERROR: DIRECTmf: no acceptable pivot at node ', &
              gmat%gid((idof-1)/hecmat%NDOF + 1), &
              ' dof ', mod(idof-1, hecmat%NDOF) + 1, ' (the matrix may be singular)'
          else if (ierr == -1) then
            write(imsg,*) 'ERROR: DIRECTmf: block size of the matrix does not match the symbolic structure'
            write(*,*) 'ERROR: DIRECTmf: block size of the matrix does not match the symbolic structure'
          else
            write(imsg,*) 'ERROR: DIRECTmf: the nonzero structure of the matrix is not symmetric'
            write(*,*) 'ERROR: DIRECTmf: the nonzero structure of the matrix is not symmetric'
          endif
        endif
        call hecmw_abort(hecmw_comm_get_comm())
      endif
      hecmat%Iarray(97) = 0
      allocate(wpr(map%nprocs))
      call mf_reduce_stats(fct, wpr)
      if (timelog > 0) write(*,'(a,f10.3,a)') '[DIRECTmf]: numeric fct done (', t2 - t1, ' sec)'
      if (loglevel > 0) then
        nthreads = 1
        !$ nthreads = omp_get_max_threads()
        write(*,'(a,i0,a,i0)') '[DIRECTmf]: ranks = ', map%nprocs, ', threads per rank = ', nthreads
        write(*,'(a,1pe9.2,a,1pe9.2)') '[DIRECTmf]: options: pivot u = ', fct%pivot_u, &
          ', zero pivot factor = ', fct%pivot_zero
        if (fct%lu) then
          if (fct%mode == 2) then
            write(*,'(a)') '[DIRECTmf]: mode = LU (forced by option)'
          else
            write(*,'(a)') '[DIRECTmf]: mode = LU (unsymmetric matrix)'
          endif
          if (loglevel > 1) write(*,'(a,1pe9.2)') '[DIRECTmf]: asymmetry = ', fct%asym
          write(*,'(a,i0,a,i0,a,i0,a)') '[DIRECTmf]: row swaps = ', fct%n_swap, ', delayed = ', fct%n_delay, &
            ' (max front growth = ', fct%max_growth, ' dofs)'
        else
          if (fct%mode == 1) then
            write(*,'(a)') '[DIRECTmf]: mode = LDLt (forced by option)'
          else
            write(*,'(a)') '[DIRECTmf]: mode = LDLt'
          endif
          if (loglevel > 1) write(*,'(a,1pe9.2)') '[DIRECTmf]: asymmetry = ', fct%asym
          write(*,'(a,i0,a,i0,a,i0,a,i0,a,i0,a,i0,a)') '[DIRECTmf]: inertia (+,-) = (', fct%n_pos, ',', fct%n_neg, &
            '), 2x2 pivots = ', fct%n_2x2, ', swaps = ', fct%n_swap, ', delayed = ', fct%n_delay, &
            ' (max front growth = ', fct%max_growth, ' dofs)'
        endif
        write(*,'(a,i0,a,f10.3,a,i0,a,i0)') '[DIRECTmf]: factor words = ', fct%factor_words_act, ' (', &
          real(fct%factor_words_act, kind=kreal)*8.0d0/1024.0d0**3, ' GB), stack peak words = ', fct%stack_peak_act, &
          ', front words = ', fct%front_words_act
        write(*,'(a,*(i0,1x))') '[DIRECTmf]: factor words per rank = ', (nint(wpr(i), kind=8), i = 1, map%nprocs)
        if (fct%blr) then
          write(*,'(a,1pe9.2,a,i0,a,0pf6.1,a)') '[DIRECTmf]: BLR: eps = ', fct%eps, &
            ', factor words full rank = ', fct%blr_words_fr, ' (', &
            100.0d0*(1.0d0 - real(fct%factor_words_act, kind=kreal)/max(real(fct%blr_words_fr, kind=kreal), 1.0d0)), &
            '% saved)'
          write(*,'(a,i0,a,f5.2,a,l1)') '[DIRECTmf]: BLR: skip eta = ', fct%blr_eta, &
            ', gain cap beta = ', fct%blr_beta, ', rank reuse = ', fct%blr_reuse
          write(*,'(a,i0,a,i0,a,i0,a,i0,a,f8.1)') '[DIRECTmf]: BLR: tiles compressed = ', fct%blr_tiles_lr, ' of ', &
            fct%blr_tiles, ', skipped = ', fct%blr_tiles_skip, ', rank max = ', fct%blr_rank_max, ', rank avg = ', &
            real(fct%blr_rank_sum, kind=kreal)/max(real(fct%blr_tiles_lr, kind=kreal), 1.0d0)
        endif
        write(*,'(a,i0,a,f10.3,a)') '[DIRECTmf]: front words held concurrently (peak) = ', fct%front_peak, ' (', &
          real(fct%front_peak, kind=kreal)*8.0d0/1024.0d0**3, ' GB)'
      endif
      deallocate(wpr)
    endif
 
    if (.not. fct%factored) then
      if (hecmw_comm_get_rank() == 0) then
        write(imsg,*) 'ERROR: DIRECTmf: numeric factorization not performed'
        write(*,*) 'ERROR: DIRECTmf: numeric factorization not performed'
      endif
      call hecmw_abort(hecmw_comm_get_comm())
    endif
    allocate(gb(fct%ndof_tot), gx(fct%ndof_tot))
    ! the global rows follow the gid numbering; ofs locates my nodes in the gathered order
    ofs = gmat%ndisp(hecmw_comm_get_rank()+1)
    call hecmw_mf_dist_gather_vec(gmat, hecmat%B, gb)
    t1 = hecmw_wtime()
    call hecmw_mf_numeric_solve_mpi(sym, map, fct, gb, gx)
    t2 = hecmw_wtime()
    if (timelog > 0) write(*,'(a,f10.3,a)') '[DIRECTmf]: solve done (', t2 - t1, ' sec)'
    do i = 1, hecmat%N
      gnode = gmat%g2s(ofs + i)
      do d = 1, hecmat%NDOF
        hecmat%X((i-1)*hecmat%NDOF + d) = gx((gnode-1)*hecmat%NDOF + d)
      enddo
    enddo
 
    irmax = hecmat%Iarray(44)
    if (fct%blr .and. irmax == 0) irmax = mf_ir_steps
    if (irmax > 0) then
      irtol = hecmat%Rarray(42)
      if (.not. (irtol > 0.0d0)) irtol = mf_ir_tol
      n = hecmat%NP * hecmat%NDOF
      allocate(rr(n), xp(hecmat%N * hecmat%NDOF))
      tcomm = 0.0d0
      call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, hecmat%B, hecmat%B, bnrm, tcomm)
      t1 = hecmw_wtime()
      rprev = huge(rprev)
      do it = 1, irmax
        call hecmw_matresid(hecmesh, hecmat, hecmat%X, hecmat%B, rr, tcomm)
        call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, rr, rr, rnrm, tcomm)
        rnrm = sqrt(rnrm / max(bnrm, tiny(bnrm)))
        if (iterlog > 0) write(*,'(a,i0,a,1pe11.4)') '[DIRECTmf]: refinement ', it - 1, ': residual = ', rnrm
        if (rnrm <= irtol) exit
        ! a step that grew the residual would hand back a worse solution than it got, so
        ! the step is undone and the sweep stops (rnrm is collective, every rank agrees)
        if (rnrm > rprev) then
          hecmat%X(1:hecmat%N * hecmat%NDOF) = xp(1:hecmat%N * hecmat%NDOF)
          if (hecmw_comm_get_rank() == 0) write(*,'(a,1pe11.4)') &
            '[DIRECTmf]: WARNING: refinement diverged, reverted to residual = ', rprev
          exit
        endif
        rprev = rnrm
        xp(1:hecmat%N * hecmat%NDOF) = hecmat%X(1:hecmat%N * hecmat%NDOF)
        call hecmw_mf_dist_gather_vec(gmat, rr, gb)
        call hecmw_mf_numeric_solve_mpi(sym, map, fct, gb, gx)
        do i = 1, hecmat%N
          gnode = gmat%g2s(ofs + i)
          do d = 1, hecmat%NDOF
            hecmat%X((i-1)*hecmat%NDOF + d) = hecmat%X((i-1)*hecmat%NDOF + d) + gx((gnode-1)*hecmat%NDOF + d)
          enddo
        enddo
      enddo
      t2 = hecmw_wtime()
      if (timelog > 0) write(*,'(a,f10.3,a)') '[DIRECTmf]: refinement done (', t2 - t1, ' sec)'
      deallocate(rr, xp)
    else if (irmax == 0 .and. hecmw_mat_get_iterlog(hecmat) > 0) then
      ! every rank enters the residual evaluation (it is collective); the passed iterlog
      ! is zeroed off rank 0 and only gates the print
      n = hecmat%NP * hecmat%NDOF
      allocate(rr(n))
      tcomm = 0.0d0
      call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, hecmat%B, hecmat%B, bnrm, tcomm)
      call hecmw_matresid(hecmesh, hecmat, hecmat%X, hecmat%B, rr, tcomm)
      call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, rr, rr, rnrm, tcomm)
      rnrm = sqrt(rnrm / max(bnrm, tiny(bnrm)))
      if (iterlog > 0) write(*,'(a,1pe11.4)') '[DIRECTmf]: residual = ', rnrm
      deallocate(rr)
    endif
    call hecmw_update_r(hecmesh, hecmat%X, hecmat%NP, hecmat%NDOF)
    deallocate(gb, gx)
  end subroutine hecmw_solve_direct_mf_dist
 
  subroutine mf_gather_matwords(gmat, wmr, wrepl)
    implicit none
    type(hecmwst_mf_gmat), intent(in) :: gmat
    real(kind=kreal), intent(out) :: wmr(:)
    integer(kind=8), intent(out) :: wrepl
    real(kind=kreal) :: w1(1)
    integer(kind=kint), allocatable :: rcs(:), disp(:)
    integer(kind=kint) :: comm, np, r
    integer(kind=8) :: wpart
 
    comm = hecmw_comm_get_comm()
    np = hecmw_comm_get_size()
    call hecmw_mf_dist_gmat_words(gmat, wpart, wrepl)
    allocate(rcs(np), disp(np))
    do r = 1, np
      rcs(r) = 1
      disp(r) = r - 1
    enddo
    w1(1) = real(wpart, kind=kreal)
    call hecmw_allgatherv_real(w1, 1, wmr, rcs, disp, comm)
    deallocate(rcs, disp)
  end subroutine mf_gather_matwords
 
  subroutine mf_reduce_stats(fct, wpr)
    implicit none
    type(hecmwst_mf_factor), intent(inout) :: fct
    real(kind=kreal), intent(out) :: wpr(:)
    real(kind=kreal) :: rs(13), rm(3), w1(1)
    integer(kind=kint), allocatable :: rcs(:), disp(:)
    integer(kind=kint) :: comm, np, r
 
    comm = hecmw_comm_get_comm()
    np = hecmw_comm_get_size()
    allocate(rcs(np), disp(np))
    do r = 1, np
      rcs(r) = 1
      disp(r) = r - 1
    enddo
    w1(1) = real(fct%factor_words_act, kind=kreal)
    call hecmw_allgatherv_real(w1, 1, wpr, rcs, disp, comm)
    deallocate(rcs, disp)
    rs(1) = real(fct%n_pos, kind=kreal)
    rs(2) = real(fct%n_neg, kind=kreal)
    rs(3) = real(fct%n_2x2, kind=kreal)
    rs(4) = real(fct%n_swap, kind=kreal)
    rs(5) = real(fct%n_delay, kind=kreal)
    rs(6) = real(fct%factor_words_act, kind=kreal)
    rs(7) = real(fct%stack_peak_act, kind=kreal)
    rs(8) = real(fct%front_peak, kind=kreal)
    rs(9) = real(fct%blr_words_fr, kind=kreal)
    rs(10) = real(fct%blr_tiles, kind=kreal)
    rs(11) = real(fct%blr_tiles_lr, kind=kreal)
    rs(12) = real(fct%blr_rank_sum, kind=kreal)
    rs(13) = real(fct%blr_tiles_skip, kind=kreal)
    call hecmw_allreduce_r_comm(rs, 13, hecmw_sum, comm)
    rm(1) = real(fct%max_growth, kind=kreal)
    rm(2) = real(fct%front_words_act, kind=kreal)
    rm(3) = real(fct%blr_rank_max, kind=kreal)
    call hecmw_allreduce_r_comm(rm, 3, hecmw_max, comm)
    fct%n_pos = nint(rs(1), kind=kint)
    fct%n_neg = nint(rs(2), kind=kint)
    fct%n_2x2 = nint(rs(3), kind=kint)
    fct%n_swap = nint(rs(4), kind=kint)
    fct%n_delay = nint(rs(5), kind=kint)
    fct%factor_words_act = nint(rs(6), kind=8)
    fct%stack_peak_act = nint(rs(7), kind=8)
    fct%front_peak = nint(rs(8), kind=8)
    fct%blr_words_fr = nint(rs(9), kind=8)
    fct%blr_tiles = nint(rs(10), kind=8)
    fct%blr_tiles_lr = nint(rs(11), kind=8)
    fct%blr_rank_sum = nint(rs(12), kind=8)
    fct%blr_tiles_skip = nint(rs(13), kind=8)
    fct%max_growth = nint(rm(1), kind=kint)
    fct%front_words_act = nint(rm(2), kind=8)
    fct%blr_rank_max = nint(rm(3), kind=kint)
  end subroutine mf_reduce_stats
 
  subroutine hecmw_mf_compare_builtin(hecMESH, hecMAT, graph, ordering, nerr)
    implicit none
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    type(hecmwst_matrix), intent(in) :: hecmat
    type(hecmwst_mf_graph), intent(in) :: graph
    integer(kind=kint), intent(in) :: ordering
    integer(kind=kint), intent(out) :: nerr
    type(cholesky_factor) :: fct
    type(hecmwst_mf_symbolic) :: ref
    integer(kind=kint), allocatable :: cnt(:)
    integer(kind=kint) :: n, k, l, ir, pb, myk, mp, nbad_tree, nbad_cnt
 
    nerr = 0
    call setij(hecmesh, hecmat, fct)
    call matini(fct, ordering, 0, ir)
    if (ir /= 0) then
      write(*,'(a,i0)') '[DIRECTmf] compare: built-in MATINI failed, ir = ', ir
      nerr = 1
      return
    endif
    n = fct%NEQns
    call hecmw_mf_symbolic_build(graph, fct%IPErm(1:n), 0, ref)
 
    nbad_tree = 0
    do k = 1, n
      pb = fct%PARent(k)
      myk = ref%invp(fct%IPErm(k))
      mp = ref%parent(myk)
      if (pb == 0 .or. pb == n + 1) then
        if (mp /= 0) nbad_tree = nbad_tree + 1
      else
        if (mp == 0) then
          nbad_tree = nbad_tree + 1
        else if (ref%perm(mp) /= fct%IPErm(pb)) then
          nbad_tree = nbad_tree + 1
        endif
      endif
    enddo
 
    allocate(cnt(n))
    cnt(1:n) = 1
    do k = 1, n
      do l = fct%XLNzr(k), fct%XLNzr(k+1) - 1
        cnt(fct%COLno(l)) = cnt(fct%COLno(l)) + 1
      enddo
    enddo
    nbad_cnt = 0
    do k = 1, n
      if (cnt(k) /= ref%colcnt(ref%invp(fct%IPErm(k)))) nbad_cnt = nbad_cnt + 1
    enddo
    if (nbad_tree > 0) write(*,'(a,i0)') '[DIRECTmf] compare: elimination tree mismatches = ', nbad_tree
    if (nbad_cnt > 0) write(*,'(a,i0)') '[DIRECTmf] compare: column count mismatches = ', nbad_cnt
    nerr = nbad_tree + nbad_cnt
 
    deallocate(cnt)
    call hecmw_mf_symbolic_finalize(ref)
    deallocate(fct%IROw, fct%JCOl, fct%IPErm, fct%INVp, fct%PARent, fct%NCH, fct%XLNzr, fct%COLno)
  end subroutine hecmw_mf_compare_builtin
 
end module hecmw_solver_direct_mf