solve_pcg Subroutine

   subroutine solve_pcg(mesh, diag, coeff, nb, rhs, G, 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) :: G(:)
      integer, intent(in) :: max_iter
      real(rk), intent(in) :: tol

      real(rk), allocatable :: r(:), z(:), p(:), ap(:), inv_diag(:)
      real(rk) :: rz, rz_new, alpha, beta, pap, r_norm, r0_norm
      integer :: iter, c
      logical :: converged

      allocate(r(mesh%ncells))
      allocate(z(mesh%ncells))
      allocate(p(mesh%ncells))
      allocate(ap(mesh%ncells))
      allocate(inv_diag(mesh%ncells))

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

      ! r = rhs - A*G
      call apply_matrix(mesh, diag, coeff, nb, G, ap)
      r = rhs - ap

      r0_norm = sqrt(dot_product(r, r))
      r_norm = r0_norm
      converged = .false.

      if (r0_norm < tol) then
         converged = .true.
         deallocate(r, z, p, ap, inv_diag)
         return
      end if

      z = inv_diag * r
      p = z
      rz = dot_product(r, z)

      do iter = 1, max_iter
         call apply_matrix(mesh, diag, coeff, nb, p, ap)

         pap = dot_product(p, ap)
         if (abs(pap) < 1.0e-30_rk) then
            write(*,'(A)') ' [WARNING] mod_radiation_p1: PCG breakdown (pAp = 0).'
            exit
         end if

         alpha = rz / pap
         G = G + alpha * p
         r = r - alpha * ap

         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

         z = inv_diag * r
         rz_new = dot_product(r, z)
         if (abs(rz) < 1.0e-30_rk) then
            write(*,'(A)') ' [WARNING] mod_radiation_p1: PCG breakdown (rMz = 0).'
            exit
         end if

         beta = rz_new / rz
         p = z + beta * p
         rz = rz_new
      end do

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

      deallocate(r, z, p, ap, inv_diag)
   end subroutine solve_pcg