pure subroutine gauss_seidel(n, hl0, z, x)
!! This routine provides the inverse of the spatial factor for a product preconditoner in an
!! ns-species reaction-diffusion problem. It performs `itmax` Gauss-Seidel iterations to
!! compute an approximation to \(A_S^{-1} z\), where \(A_S = I - h_{l0} J_ d \), and \(J_d\)
!! represents the diffusion contributions to the Jacobian. The solution vector is returned
!! in `z`.
use web_par, only: cox, coy, ns, mx, my, mxns
integer, intent(in) :: n
real(rk), intent(in) :: hl0
real(rk), intent(inout) :: z(*)
real(rk), intent(inout) :: x(*)
integer, parameter :: itmax = 5
integer :: i, ic, ici, iter, iyoff, jx, jy
real(rk) :: elamda, beta1(ns), gamma1(ns), beta2(ns), gamma2(ns), dinv(ns)
! Write matrix as A = D - L - U.
! Load local arrays BETA, BETA2, GAMMA1, GAMMA2, and DINV.
do i = 1, ns
elamda = one/(one + 2*hl0*(cox(i) + coy(i)))
beta1(i) = hl0*cox(i)*elamda
beta2(i) = 2*beta1(i)
gamma1(i) = hl0*coy(i)*elamda
gamma2(i) = 2*gamma1(i)
dinv(i) = elamda
end do
! Zero X in all its components, since X is added to Z at the end.
x(1:n) = zero
! Load array X with (D-inverse)*Z for first iteration.
do jy = 1, my
iyoff = mxns*(jy - 1)
do jx = 1, mx
ic = iyoff + ns*(jx - 1)
do i = 1, ns
ici = ic + i
x(ici) = dinv(i)*z(ici)
z(ici) = zero
end do
end do
end do
do iter = 1, itmax
! Calculate (D-inverse)*U*X
if (iter > 1) then
jy = 1
jx = 1
ic = ns*(jx - 1)
do i = 1, ns
ici = ic + i
x(ici) = beta2(i)*x(ici + ns) + gamma2(i)*x(ici + mxns)
end do
do jx = 2, mx - 1
ic = ns*(jx - 1)
do i = 1, ns
ici = ic + i
x(ici) = beta1(i)*x(ici + ns) + gamma2(i)*x(ici + mxns)
end do
end do
jx = mx
ic = ns*(jx - 1)
do i = 1, ns
ici = ic + i
x(ici) = gamma2(i)*x(ici + mxns)
end do
do jy = 2, my - 1
iyoff = mxns*(jy - 1)
jx = 1
ic = iyoff
do i = 1, ns
ici = ic + i
x(ici) = beta2(i)*x(ici + ns) + gamma1(i)*x(ici + mxns)
end do
do jx = 2, mx - 1
ic = iyoff + ns*(jx - 1)
do i = 1, ns
ici = ic + i
x(ici) = beta1(i)*x(ici + ns) + gamma1(i)*x(ici + mxns)
end do
end do
jx = mx
ic = iyoff + ns*(jx - 1)
do i = 1, ns
ici = ic + i
x(ici) = gamma1(i)*x(ici + mxns)
end do
end do
jy = my
iyoff = mxns*(jy - 1)
jx = 1
ic = iyoff
do i = 1, ns
ici = ic + i
x(ici) = beta2(i)*x(ici + ns)
end do
do jx = 2, mx - 1
ic = iyoff + ns*(jx - 1)
do i = 1, ns
ici = ic + i
x(ici) = beta1(i)*x(ici + ns)
end do
end do
jx = mx
ic = iyoff + ns*(jx - 1)
do i = 1, ns
ici = ic + i
x(ici) = zero
end do
end if
! Calculate [(I - (D-inverse)*L)]-inverse * X
jy = 1
do jx = 2, mx - 1
ic = ns*(jx - 1)
do i = 1, ns
ici = ic + i
x(ici) = x(ici) + beta1(i)*x(ici - ns)
end do
end do
jx = mx
ic = ns*(jx - 1)
do i = 1, ns
ici = ic + i
x(ici) = x(ici) + beta2(i)*x(ici - ns)
end do
do jy = 2, my - 1
iyoff = mxns*(jy - 1)
jx = 1
ic = iyoff
do i = 1, ns
ici = ic + i
x(ici) = x(ici) + gamma1(i)*x(ici - mxns)
end do
do jx = 2, mx - 1
ic = iyoff + ns*(jx - 1)
do i = 1, ns
ici = ic + i
x(ici) = (x(ici) + beta1(i)*x(ici - ns)) + gamma1(i)*x(ici - mxns)
end do
end do
jx = mx
ic = iyoff + ns*(jx - 1)
do i = 1, ns
ici = ic + i
x(ici) = (x(ici) + beta2(i)*x(ici - ns)) + gamma1(i)*x(ici - mxns)
end do
end do
jy = my
iyoff = mxns*(jy - 1)
jx = 1
ic = iyoff
do i = 1, ns
ici = ic + i
x(ici) = x(ici) + gamma2(i)*x(ici - mxns)
end do
do jx = 2, mx - 1
ic = iyoff + ns*(jx - 1)
do i = 1, ns
ici = ic + i
x(ici) = (x(ici) + beta1(i)*x(ici - ns)) + gamma2(i)*x(ici - mxns)
end do
end do
jx = mx
ic = iyoff + ns*(jx - 1)
do i = 1, ns
ici = ic + i
x(ici) = (x(ici) + beta2(i)*x(ici - ns)) + gamma2(i)*x(ici - mxns)
end do
! Add increment X to Z
do i = 1, n
z(i) = z(i) + x(i)
end do
end do
end subroutine gauss_seidel