solve_pbicgstab Subroutine

   subroutine solve_pbicgstab(mesh, diag, coeff, nb, rhs, x, max_iter, tol)
      type(mesh_t), intent(in) :: mesh
      real(rk), intent(in) :: diag(:)
      real(rk), intent(in) :: coeff(:,:)
      integer, intent(in) :: nb(:,:)
      real(rk), intent(in) :: rhs(:)
      real(rk), intent(inout) :: x(:)
      integer, intent(in) :: max_iter
      real(rk), intent(in) :: tol

      real(rk), allocatable :: r(:), r0_star(:), p(:), v(:), s(:), t(:)
      real(rk), allocatable :: y(:), z(:), ax(:), inv_diag(:)
      real(rk) :: rho_prev, rho_curr, alpha, omega, beta, r0_norm, r_norm, denom, tt
      integer :: iter, c
      logical :: converged

      allocate(r(mesh%ncells))
      allocate(r0_star(mesh%ncells))
      allocate(p(mesh%ncells))
      allocate(v(mesh%ncells))
      allocate(s(mesh%ncells))
      allocate(t(mesh%ncells))
      allocate(y(mesh%ncells))
      allocate(z(mesh%ncells))
      allocate(ax(mesh%ncells))
      allocate(inv_diag(mesh%ncells))

      do c = 1, mesh%ncells
         if (diag(c) <= tiny_safe) then
            call fatal_error('dom', 'Non-positive diagonal encountered in DOM BiCGStab solver')
         end if
         inv_diag(c) = 1.0_rk / diag(c)
      end do

      call apply_matrix(mesh, diag, coeff, nb, x, ax)
      r = rhs - ax
      r0_star = r
      r0_norm = sqrt(dot_product(r, r))
      r_norm = r0_norm
      converged = .false.

      if (r0_norm < tol) then
         converged = .true.
         deallocate(r, r0_star, p, v, s, t, y, z, ax, inv_diag)
         return
      end if

      p = zero
      v = zero
      rho_prev = 1.0_rk
      alpha = 1.0_rk
      omega = 1.0_rk

      do iter = 1, max_iter
         rho_curr = dot_product(r0_star, r)
         if (abs(rho_curr) < 1.0e-30_rk) then
            write(*,'(A)') ' [WARNING] mod_radiation_dom: BiCGStab breakdown (rho = 0).'
            exit
         end if

         if (iter == 1) then
            p = r
         else
            if (abs(omega) < 1.0e-30_rk) then
               write(*,'(A)') ' [WARNING] mod_radiation_dom: BiCGStab breakdown (omega = 0).'
               exit
            end if
            beta = (rho_curr / rho_prev) * (alpha / omega)
            p = r + beta * (p - omega * v)
         end if

         ! Left Jacobi preconditioning: y = M^{-1} p.
         y = inv_diag * p
         call apply_matrix(mesh, diag, coeff, nb, y, v)

         denom = dot_product(r0_star, v)
         if (abs(denom) < 1.0e-30_rk) then
            write(*,'(A)') ' [WARNING] mod_radiation_dom: BiCGStab breakdown (denom = 0).'
            exit
         end if

         alpha = rho_curr / denom
         s = r - alpha * v

         r_norm = sqrt(dot_product(s, s))
         if (r_norm < tol * r0_norm .or. r_norm < 1.0e-12_rk) then
            x = x + alpha * y
            converged = .true.
            exit
         end if

         ! z = M^{-1} s.
         z = inv_diag * s
         call apply_matrix(mesh, diag, coeff, nb, z, t)

         tt = dot_product(t, t)
         if (abs(tt) < 1.0e-30_rk) then
            write(*,'(A)') ' [WARNING] mod_radiation_dom: BiCGStab breakdown (t.t = 0).'
            exit
         end if

         omega = dot_product(t, s) / tt
         x = x + alpha * y + omega * z
         r = s - omega * t

         r_norm = sqrt(dot_product(r, r))
         if (r_norm < tol * r0_norm .or. r_norm < 1.0e-12_rk) then
            converged = .true.
            exit
         end if

         rho_prev = rho_curr
      end do

      if (.not. converged) then
         write(*,'(A,I0,A,ES12.5)') ' [WARNING] mod_radiation_dom: BiCGStab failed to converge in ', &
                                      max_iter, ' iterations. Final residual norm: ', r_norm
      end if

      deallocate(r, r0_star, p, v, s, t, y, z, ax, inv_diag)
   end subroutine solve_pbicgstab