Stationary.F
#include "parmdefines.h"
#include "iadefines.h"
c***********************************************************************
#include "author.inc"
c* $Id: Stationary.F,v 1.28 1996/04/29 02:41:35 turner Exp $
c*
c* Computes the solution to a linear system of equations of the
c* form Ax=b by a stationary iterative method (Jacobi, Jacobi with
c* relaxation, Gauss-Seidel, SOR, Symmetric Gauss-Seidel, or SSOR).
c*
c* An error message is written out and control is returned to the
c* calling routine if the maximum number of iterations is exceeded.
c* At that time, err(iparm(_JT_iter_)) can be examined to determine
c* if execution can continue.
c*
c* <PARAMETER LIST>
c*
c* Input:
c* b - source vector
c* a - coefficient matrix
c* ia - integer vector containing info about how "a" is stored
c* NOTE: see description of ia below
c* ja - integer map for coefficient
c*
c* In/Out:
c* iparm - array of integer parameters
c* rparm - array of floating point parameters
c* NOTE: see description of iparm and rparm below
c* x - solution vector (whatever is in x on entry is used as an
c* initial guess)
c*
c* Output:
c* cpu - cumulative cpu time after each iteration
c* rnormt - norm of true residual
c* errt - error estimate based on true residual
c* err - error estimate at each iteration
c* status - return status
c* -4 ==> iteration limit exceeded
c* -3 ==> internal error
c* -2 ==> memory allocation failure
c* -1 ==> invalid argument(s)
c* 0 ==> success
c*
#include "parmdesc.inc"
c*
#include "iadesc.inc"
c*
c* <SUBROUTINES REQUIRED>
c*
c* JT_CheckConvergence
c* JT_Clock
c* JT_Jacobi
c* JT_SOR
c* JT_WriteSystem
c* JT_WriteVectorFloat
c* JT_y_eq_x
c* JT_y_eq_y_minus_Ax
c*
c* <FUNCTIONS REQUIRED>
c*
c* JT_MatrixNorm
c* JT_VectorNorm
c*
#include "arrays-Stationary.inc"
c*
#include "copyright.inc"
c***********************************************************************
subroutine JT_Stationary (b, a, ia, ja, iparm, rparm, x,
& cpu, rnormt, errt, err, status)
implicit none
c
c ... Input:
integer ia(_JT_no_of_storage_parameters_), ja(*)
real a(*)
#ifdef strict_f77
real b(*)
#else
real b(ia(_JT_nrows_))
#endif
c
c ... In/Out:
integer iparm(_JT_no_of_iparms_)
real rparm(_JT_no_of_rparms_)
#ifdef strict_f77
real x(*)
#else
real x(ia(_JT_nrows_))
#endif
c
c ... Output:
integer status
#ifdef strict_f77
real cpu(0:*), err(0:*), rnormt(0:*), errt(0:*)
#else
real cpu(0:iparm(_JT_itmax_))
real err(0:iparm(_JT_itmax_)),
& rnormt(0:iparm(_JT_itmax_)),
& errt(0:iparm(_JT_itmax_))
#endif
c
c ... Local:
integer i
real zero, bnorm, one, rnorm, cpu_tmp
real JT_MatrixNorm, JT_VectorNorm
#include "declare-Stationary.inc"
c
c ... Parameters.
parameter (zero=0.0d0, one=1.0d0)
c
c ... Obtain CPU time at start of solution.
call JT_Clock (1, cpu(0), status)
c
c ... Check to make sure a valid value of iparm(_JT_stop_) has been chosen.
if (ABS(iparm(_JT_stop_)) .gt. 5) then
status = -1
return
endif
c
#include "allocate-Stationary.inc"
c
c ... Initialize iteration counter.
iparm(_JT_iter_) = 0
c
c ... Initialize cpu time and residual and error norm arrays.
call JT_FillVectorFloat (iparm(_JT_itmax_), zero, cpu(1), status)
call JT_FillVectorFloat (iparm(_JT_itmax_)+1, zero, err, status)
call JT_FillVectorFloat (iparm(_JT_itmax_)+1, zero, errt, status)
call JT_FillVectorFloat (iparm(_JT_itmax_)+1, zero, rnormt, status)
c
c ... rnorm is a dummy variable since it is not meaningful in this
c routine, but go ahead and set it to prevent compilers from
c complaining about it being referenced before being set.
rnorm = one
c
c ... Write initial system if desired.
if (iparm(_JT_out_) .ge. 4) then
call JT_WriteSystem (iparm(_JT_luout_), a, ia, ja, x, b,
& 'JT_Stationary: ON ENTRY (a,x,b):', status)
endif
c
c ... Check to make sure that if one of the residual-based stopping
c tests is chosen it is one based on the true residual. If not,
c correct and continue.
if (iparm(_JT_stop_) .gt. 0) iparm(_JT_stop_) = -iparm(_JT_stop_)
c
c ... Initialization dependent on the stopping test.
if (iparm(_JT_stop_).eq.-1 .or. iparm(_JT_stop_).eq.-2) then
c
c Compute || b || if needed.
bnorm = JT_VectorNorm(iparm(_JT_norm_), ia(_JT_nrows_), b, status)
else
bnorm = zero
endif
c
rparm(_JT_anorm_) = zero
if (iparm(_JT_stop_) .eq. _JT_stop_relchg_) then
c
c Initialize xold if using || x - xold || / || x || stopping test.
call JT_y_eq_x (ia(_JT_nrows_), x, xold, status)
elseif (iparm(_JT_stop_) .eq. -_JT_stop_axb_) then
c
c Compute || A || if needed.
rparm(_JT_anorm_) = JT_MatrixNorm(iparm(_JT_norm_), a, ia, ja, status)
if (status .eq. -2) goto 9999
elseif (iparm(_JT_stop_) .eq. -_JT_stop_r0_) then
c
c Compute initial residual if using || r || / || r0 || stopping
c test and write out if desired.
call JT_y_eq_x (ia(_JT_nrows_), b, w, status)
call JT_y_eq_y_minus_Ax (a, ia, ja, x, w, status)
if (status .lt. 0) then
if (status .eq. -1) status = -3
goto 9999
endif
if (iparm(_JT_out_) .ge. 6) then
call JT_WriteVectorFloat (iparm(_JT_luout_), ia(_JT_nrows_), w,
& 'JT_Stationary: INITIAL RESIDUAL (w):', status)
endif
rnormt(0) = JT_VectorNorm(iparm(_JT_norm_), ia(_JT_nrows_), w, status)
if (rnormt(0) .le. rparm(_JT_tiny_)) then
status = 0
goto 9999
endif
endif
c
c ... Check to see if initial guess meets the convergence
c criterion (note that rnorm is a dummy variable since
c it is not meaningful in this routine).
c
c Note that on entry iparm(_JT_iter_) = 0 and JT_CheckConvergence
c increments iparm(_JT_iter_).
call JT_CheckConvergence (bnorm, rnormt(0), rnorm,
& x, xold, b, a, ia, ja, iparm, rparm, cpu_tmp,
& err(iparm(_JT_iter_)), rnormt(iparm(_JT_iter_)),
& errt(iparm(_JT_iter_)), status)
if (status .le. 0) goto 9999
c
c ... Main loop.
10 continue
c
if (iparm(_JT_method_) .eq. _JT_method_Jacobi_) then
c
c .... Jacobi, Jacobi with relaxation.
call JT_Jacobi (rparm(_JT_omega_), b, a, ia, ja, x, status)
else
c
c .... Gauss-Seidel, SOR, Symmetric Gauss-Seidel, SSOR.
call JT_SOR (iparm(_JT_method_), rparm(_JT_omega_), b, a, ia, ja,
& x, status)
endif
c
c ... Check return status of iteration routines.
if (status .lt. 0) then
if (status .eq. -1) status = -3
goto 9999
endif
c
c ... Check convergence (note that rnorm is a dummy variable since
c it is not meaningful in this routine).
call JT_CheckConvergence (bnorm, rnormt(0), rnorm,
& x, xold, b, a, ia, ja, iparm, rparm, cpu(iparm(_JT_iter_)),
& err(iparm(_JT_iter_)), rnormt(iparm(_JT_iter_)), errt(iparm(_JT_iter_)),
& status)
if (status .le. 0) goto 9999
c
goto 10
c
9999 continue
c
c ... Subtract time at start of solution from each element of cpu array
c to obtain true CPU time for each iteration.
do i=1,iparm(_JT_iter_)
cpu(i) = cpu(i) - cpu(0)
enddo
#include "deallocate-Stationary.inc"
c
return
end