subroutine psolrs(neq, t, cc, ccprime, savr, wk, cj, wt, wp, iwp, b, eplin, ier, rpar, ipar)
!! This routine applies the inverse of a product preconditioner matrix to the vector in the
!! array `b`. Depending on the flag `jpre`, this involves a call to [[gs]], for the inverse
!! of the spatial factor, and/or a call to [[DRBDPS]] or [[DRBGPS]] for the inverse of the
!! reaction-based factor (`cj*I_d - dR/dy`). The latter factor uses block-grouping (with a
!! call to [[DRBGPS]]) if `jbg == 1`, and does not (with a call to [[DRBDPS]]) if `jbg == 0`.
!! the flag `jbg` is passed as `ipar(2)`. The array `b` is overwritten with the solution.
integer, intent(in) :: neq
real(rk), intent(in) :: t
real(rk), intent(in) :: cc(*)
real(rk), intent(in) :: ccprime(*)
real(rk), intent(in) :: savr(*)
real(rk), intent(inout) :: wk(*)
real(rk), intent(in) :: cj
real(rk), intent(in) :: wt(*)
real(rk), intent(in) :: wp(*)
integer, intent(in) :: iwp(*)
real(rk), intent(inout) :: b(*)
real(rk), intent(in) :: eplin
integer, intent(inout) :: ier
real(rk), intent(inout) :: rpar(*)
integer, intent(in) :: ipar(*)
integer :: jbg, jpre
real(rk) :: hl0
external :: drbdps, drbgps
jpre = ipar(1)
ier = 0
hl0 = one/cj
jbg = ipar(2)
if (jpre == 2 .or. jpre == 3) call gauss_seidel(neq, hl0, b, wk)
if (jpre /= 2) then
if (jbg == 0) call drbdps(b, wp, iwp)
if (jbg == 1) call drbgps(b, wp, iwp)
end if
if (jpre == 4) call gauss_seidel(neq, hl0, b, wk)
end subroutine psolrs