gmres.f90

Go to the documentation of this file.

! Copyright (c) 2020-2026, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * Redistributions in binary form must reproduce the above
!     copyright notice, this list of conditions and the following
!     disclaimer in the documentation and/or other materials provided
!     with the distribution.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module gmres
  !$ use omp_lib
  use krylov, only : ksp_t, ksp_monitor_t
  use precon, only : pc_t
  use ax_product, only : ax_t
  use num_types, only : rp, xp
  use field, only : field_t
  use coefs, only : coef_t
  use gather_scatter, only : gs_t, gs_op_add
  use scalar_bc_projector, only : scalar_bc_projector_t
  use vector_bc_projector, only : vector_bc_projector_t, &
       vector_bc_projector_components
  use matrix, only : matrix_t
  use scratch_registry, only : neko_scratch_registry
  use math, only : glsc3, rzero, copy, sub2, cmult2, abscmp
  use neko_config, only : neko_blk_size
  use comm, only : neko_comm, mpi_extra_precision
  use mpi_f08, only : mpi_allreduce, mpi_in_place, mpi_sum
  implicit none
  private
 
  type, public, extends(ksp_t) :: gmres_t
     integer :: lgmres = 30
     real(kind=rp), pointer :: w(:) => null()
     real(kind=rp), pointer :: r(:) => null()
     type(matrix_t), pointer :: z => null()
     type(matrix_t), pointer :: v => null()
     real(kind=xp), allocatable :: h(:,:)
     real(kind=xp), allocatable :: s(:)
     real(kind=xp), allocatable :: gam(:)
     real(kind=xp), allocatable :: c(:)
     real(kind=xp), allocatable :: hp(:,:)
   contains
     procedure, pass(this) :: init => gmres_init
     procedure, pass(this) :: free => gmres_free
     procedure, pass(this) :: solve => gmres_solve
     procedure, pass(this) :: solve_coupled => gmres_solve_coupled
  end type gmres_t
 
contains
 
  subroutine gmres_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
    class(gmres_t), target, intent(inout) :: this
    integer, intent(in) :: n
    integer, intent(in) :: max_iter
    class(pc_t), optional, intent(in), target :: M
    real(kind=rp), optional, intent(in) :: rel_tol
    real(kind=rp), optional, intent(in) :: abs_tol
    logical, optional, intent(in) :: monitor
    integer :: nthrds
 
    call this%free()
 
    if (present(m)) then
       this%M => m
    end if
 
    allocate(this%c(this%lgmres))
    allocate(this%s(this%lgmres))
    allocate(this%gam(this%lgmres + 1))
 
    allocate(this%h(this%lgmres, this%lgmres))
 
    nthrds = 1
    !$ nthrds = omp_get_max_threads()
    allocate(this%hp(this%lgmres, nthrds))
 
    if (present(rel_tol) .and. present(abs_tol) .and. present(monitor)) then
       call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
    else if (present(rel_tol) .and. present(abs_tol)) then
       call this%ksp_init(max_iter, rel_tol, abs_tol)
    else if (present(monitor) .and. present(abs_tol)) then
       call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
    else if (present(rel_tol) .and. present(monitor)) then
       call this%ksp_init(max_iter, rel_tol, monitor = monitor)
    else if (present(rel_tol)) then
       call this%ksp_init(max_iter, rel_tol = rel_tol)
    else if (present(abs_tol)) then
       call this%ksp_init(max_iter, abs_tol = abs_tol)
    else if (present(monitor)) then
       call this%ksp_init(max_iter, monitor = monitor)
    else
       call this%ksp_init(max_iter)
    end if
 
  end subroutine gmres_init
 
  subroutine gmres_free(this)
    class(gmres_t), intent(inout) :: this
 
    call this%ksp_free()
 
    nullify(this%w)
 
    if (allocated(this%c)) then
       deallocate(this%c)
    end if
 
    nullify(this%r)
    nullify(this%z)
 
    if (allocated(this%h)) then
       deallocate(this%h)
    end if
 
    nullify(this%v)
 
    if (allocated(this%s)) then
       deallocate(this%s)
    end if
 
 
    if (allocated(this%gam)) then
       deallocate(this%gam)
    end if
 
    if (allocated(this%hp)) then
       deallocate(this%hp)
    end if
 
    nullify(this%M)
 
  end subroutine gmres_free
 
  function gmres_solve(this, Ax, x, f, n, coef, bc_projector, gs_h, niter) &
       result(ksp_results)
    class(gmres_t), intent(inout) :: this
    class(ax_t), intent(in) :: ax
    type(field_t), intent(inout) :: x
    integer, intent(in) :: n
    real(kind=rp), dimension(n), intent(in) :: f
    type(coef_t), intent(inout) :: coef
    class(scalar_bc_projector_t), intent(inout) :: bc_projector
    type(gs_t), intent(inout) :: gs_h
    type(ksp_monitor_t) :: ksp_results
    integer, optional, intent(in) :: niter
    integer :: iter, max_iter
    integer :: i, j, k, l, ierr, tid, nthrds
    real(kind=xp) :: w_plus(neko_blk_size), x_plus(neko_blk_size)
    real(kind=xp) :: hl(this%lgmres)
    real(kind=xp) :: alpha, lr, alpha2, norm_fac, tmp, acc
    real(kind=rp) :: temp, rnorm
    logical :: conv
    integer :: temp_indices(4)
 
    conv = .false.
    iter = 0
    rnorm = 0.0_rp
 
    if (present(niter)) then
       max_iter = niter
    else
       max_iter = this%max_iter
    end if
 
    nthrds = 1
    !$ nthrds = omp_get_max_threads()
 
    call neko_scratch_registry%request(this%w, temp_indices(1), n, .false.)
    call neko_scratch_registry%request(this%r, temp_indices(2), n, .false.)
    call neko_scratch_registry%request(this%z, temp_indices(3), &
         n, this%lgmres, .false.)
    call neko_scratch_registry%request(this%v, temp_indices(4), &
         n, this%lgmres, .false.)
 
    associate(w => this%w, c => this%c, r => this%r, z => this%z%x, &
         h => this%h, v => this%v%x, s => this%s, gam => this%gam, &
         hp => this%hp)
 
      norm_fac = 1.0_rp / sqrt(coef%volume)
      call rzero(x%x, n)
      gam = 0.0_xp
      s = 1.0_xp
      c = 1.0_xp
      h = 0.0_xp
      call this%monitor_start('GMRES')
      do while (.not. conv .and. iter .lt. max_iter)
 
         if (iter .eq. 0) then
            call copy(r, f, n)
         else
            call copy(r, f, n)
            call ax%compute(w, x%x, coef, x%msh, x%Xh)
            call gs_h%op(w, n, gs_op_add)
            call bc_projector%apply(w, n)
            call sub2(r, w, n)
         end if
 
         gam(1) = sqrt(glsc3(r, r, coef%mult, n))
         if (iter .eq. 0) then
            ksp_results%res_start = gam(1) * norm_fac
         end if
 
         if (abscmp(gam(1), 0.0_xp)) exit
 
         rnorm = 0.0_rp
         temp = 1.0_rp / gam(1)
         call cmult2(v(1,1), r, temp, n)
         do j = 1, this%lgmres
            iter = iter+1
 
            call this%M%solve(z(1,j), v(1,j), n)
 
            call ax%compute(w, z(1,j), coef, x%msh, x%Xh)
            call gs_h%op(w, n, gs_op_add)
            call bc_projector%apply(w, n)
 
            ! Classical Gram-Schmidt orthogonalization: accumulate
            ! <w, v_l>_mult for l=1..j in a thread-private array hl,
            ! publish it once per thread into a column of hp, merge
            ! across threads, then one MPI_Allreduce of length j.
            !$omp parallel private(i, k, l, tid, acc, hl)
            tid = 1
            !$ tid = omp_get_thread_num() + 1
            do l = 1, j
               hl(l) = 0.0_xp
            end do
            !$omp do
            do i = 0, n-1, neko_blk_size
               if (i + neko_blk_size .le. n) then
                  do l = 1, j
                     acc = hl(l)
                     !$omp simd reduction(+:acc)
                     do k = 1, neko_blk_size
                        acc = acc + &
                             w(i+k) * v(i+k,l) * coef%mult(i+k,1,1,1)
                     end do
                     hl(l) = acc
                  end do
               else
                  do l = 1, j
                     do k = 1, n - i
                        hl(l) = hl(l) + &
                             w(i+k) * v(i+k,l) * coef%mult(i+k,1,1,1)
                     end do
                  end do
               end if
            end do
            !$omp end do nowait
            do l = 1, j
               hp(l,tid) = hl(l)
            end do
            !$omp end parallel
 
            ! Cross-thread merge into hp(:, 1), then one Allreduce of j.
            do k = 2, nthrds
               do l = 1, j
                  hp(l,1) = hp(l,1) + hp(l,k)
               end do
            end do
            call mpi_allreduce(mpi_in_place, hp(1,1), j, &
                 mpi_extra_precision, mpi_sum, neko_comm, ierr)
 
            do l = 1, j
               h(l,j) = hp(l,1)
            end do
 
            ! Projection w = w - sum h(l,j) v_l, with fused <w,w>_mult
            ! reduction for the post-projection norm. alpha2 is computed
            ! directly (not by Pythagoras) so accuracy is preserved when
            ! the Krylov subspace is becoming invariant.
            alpha2 = 0.0_xp
            !$omp parallel do private(k, l, w_plus, tmp) reduction(+:alpha2)
            do i = 0, n-1, neko_blk_size
               if (i + neko_blk_size .le. n) then
                  !$omp simd
                  do k = 1, neko_blk_size
                     w_plus(k) = 0.0_xp
                  end do
                  do l = 1, j
                     !$omp simd
                     do k = 1, neko_blk_size
                        w_plus(k) = w_plus(k) - h(l,j) * v(i+k,l)
                     end do
                  end do
                  !$omp simd reduction(+:alpha2)
                  do k = 1, neko_blk_size
                     w(i+k) = w(i+k) + w_plus(k)
                     alpha2 = alpha2 + w(i+k)**2 * coef%mult(i+k,1,1,1)
                  end do
               else
                  do k = 1, n - i
                     tmp = 0.0_xp
                     do l = 1, j
                        tmp = tmp - h(l,j) * v(i+k,l)
                     end do
                     w(i+k) = w(i+k) + tmp
                     alpha2 = alpha2 + w(i+k)**2 * coef%mult(i+k,1,1,1)
                  end do
               end if
            end do
            !$omp end parallel do
            call mpi_allreduce(mpi_in_place, alpha2, 1, &
                 mpi_extra_precision, mpi_sum, neko_comm, ierr)
            alpha = sqrt(alpha2)
            do i = 1, j-1
               temp = h(i,j)
               h(i,j) = c(i)*temp + s(i) * h(i+1,j)
               h(i+1,j) = -s(i)*temp + c(i) * h(i+1,j)
            end do
 
            rnorm = 0.0_rp
            if (abscmp(alpha, 0.0_xp)) then
               conv = .true.
               exit
            end if
 
            lr = sqrt(h(j,j) * h(j,j) + alpha2)
            temp = 1.0_rp / lr
            c(j) = h(j,j) * temp
            s(j) = alpha * temp
            h(j,j) = lr
            gam(j+1) = -s(j) * gam(j)
            gam(j) = c(j) * gam(j)
            rnorm = abs(gam(j+1)) * norm_fac
            call this%monitor_iter(iter, rnorm)
            if (rnorm .lt. this%abs_tol) then
               conv = .true.
               exit
            end if
 
            if (iter + 1 .gt. max_iter) exit
 
            if (j .lt. this%lgmres) then
               temp = 1.0_rp / alpha
               call cmult2(v(1,j+1), w, temp, n)
            end if
 
         end do
 
         j = min(j, this%lgmres)
         do k = j, 1, -1
            temp = gam(k)
            do i = j, k+1, -1
               temp = temp - h(k,i) * c(i)
            end do
            c(k) = temp / h(k,k)
         end do
 
         !$omp parallel do private(k, l, x_plus, tmp)
         do i = 0, n-1, neko_blk_size
            if (i + neko_blk_size .le. n) then
               !$omp simd
               do k = 1, neko_blk_size
                  x_plus(k) = 0.0_xp
               end do
               do l = 1, j
                  !$omp simd
                  do k = 1, neko_blk_size
                     x_plus(k) = x_plus(k) + c(l) * z(i+k,l)
                  end do
               end do
               !$omp simd
               do k = 1, neko_blk_size
                  x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(k)
               end do
            else
               do k = 1, n - i
                  tmp = 0.0_xp
                  do l = 1, j
                     tmp = tmp + c(l) * z(i+k,l)
                  end do
                  x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + tmp
               end do
            end if
         end do
         !$omp end parallel do
      end do
 
    end associate
    nullify(this%w, this%r, this%z, this%v)
    call neko_scratch_registry%relinquish(temp_indices)
    call this%monitor_stop()
    ksp_results%res_final = rnorm
    ksp_results%iter = iter
    ksp_results%converged = this%is_converged(iter, rnorm)
 
  end function gmres_solve
 
  function gmres_solve_coupled(this, Ax, x, y, z, fx, fy, fz, &
       n, coef, bc_projector, gs_h, niter) result(ksp_results)
    class(gmres_t), intent(inout) :: this
    class(ax_t), intent(in) :: ax
    type(field_t), intent(inout) :: x
    type(field_t), intent(inout) :: y
    type(field_t), intent(inout) :: z
    integer, intent(in) :: n
    real(kind=rp), dimension(n), intent(in) :: fx
    real(kind=rp), dimension(n), intent(in) :: fy
    real(kind=rp), dimension(n), intent(in) :: fz
    type(coef_t), intent(inout) :: coef
    class(vector_bc_projector_t), intent(inout) :: bc_projector
    type(gs_t), intent(inout) :: gs_h
    type(ksp_monitor_t), dimension(3) :: ksp_results
    integer, optional, intent(in) :: niter
    type(scalar_bc_projector_t), pointer :: bc_x, bc_y, bc_z
 
    call vector_bc_projector_components(bc_projector, bc_x, bc_y, bc_z)
    ksp_results(1) = this%solve(ax, x, fx, n, coef, bc_x, gs_h, niter)
    ksp_results(2) = this%solve(ax, y, fy, n, coef, bc_y, gs_h, niter)
    ksp_results(3) = this%solve(ax, z, fz, n, coef, bc_z, gs_h, niter)
 
  end function gmres_solve_coupled
 
end module gmres