CheckConvergence.F
#include "parmdefines.h"
#include "iadefines.h"
c**********************************************************************
#include "author.inc"
c* $Id: CheckConvergence.F,v 1.30 1996/08/22 20:49:08 turner Exp $
c*
c* Determines whether the iteration has converged and computes
c* true residual and error estimate based on true residual, if
c* required.
c*
c* <PARAMETER LIST>
c*
c* Input:
c* bnorm - norm of right-hand-side
c* rnormt0 - norm of initial residual vector
c* rnorm - norm of residual vector
c* x - x-vector
c* b - source vector
c* a - coefficient
c* ia - integer vector containing info about how "a" is stored
c* NOTE: see description of ia below
c* ja - integer array that maps columns in a to true columns
c*
c* In/Out:
c* rnormt - norm of true residual
c* iparm - array of integer parameters
c* rparm - array of floating point parameters
c* NOTE: see description of iparm and rparm below
c* xold - x-vector at last iteration
c*
c* Output:
c* cpu - cumulative CPU time
c* err - error estimate
c* errt - error estimate based on true residual
c* status - return status:
c* -4 ==> iparm(_JT_itmax_) exceeded
c* -3 ==> internal error
c* -2 ==> memory allocation failure
c* 0 ==> convergence was achieved
c* 1 ==> convergence not achieved
c*
#include "parmdesc.inc"
c*
#include "iadesc.inc"
c*
c* <SUBROUTINES REQUIRED>
c*
c* JT_Clock
c* JT_FlushUnit
c* JT_StoppingTest
c* JT_WriteVectorFloat
c*
#include "copyright.inc"
c**********************************************************************
subroutine JT_CheckConvergence (bnorm, rnormt0, rnorm, x, xold,
& b, a, ia, ja, iparm, rparm, cpu, err, rnormt, errt,
& status)
implicit none
c
c ... Input:
integer ia(_JT_no_of_storage_parameters_), ja(*)
real bnorm, rnormt0, rnorm
#ifdef strict_f77
real b(*), x(*)
#else
real b(ia(_JT_nrows_)), x(ia(_JT_nrows_))
#endif
real a(*)
c
c ... In/Out:
integer iparm(_JT_no_of_iparms_)
real rparm(_JT_no_of_rparms_)
real rnormt
#ifdef strict_f77
real xold(*)
#else
real xold(ia(_JT_nrows_))
#endif
c
c ... Output:
integer status
real cpu, err, errt
c
c ... Local:
integer dummy
real err_test
c
c ... Write out current iterate, if appropriate.
if (iparm(_JT_out_) .ge. 5) then
call JT_WriteVectorFloat (iparm(_JT_luout_), ia(_JT_nrows_), x,
& 'JT_CheckConvergence: CURRENT ITERATE (x):', status)
endif
c
c ... Compute residual norm(s) and error estimate(s) as appropriate.
call JT_StoppingTest (bnorm, rnormt0, rnorm, x, xold, b,
& a, ia, ja, iparm, rparm, err, rnormt, errt, status)
c
c ... Write out current information, as appropriate.
if (iparm(_JT_out_) .ge. 3) then
if (iparm(_JT_resid_).eq.1 .or.
& iparm(_JT_resid_).eq.3 .or.
& iparm(_JT_stop_).lt.0) then
if (iparm(_JT_iter_).eq.1 .or. iparm(_JT_out_).ge.5) write(iparm(_JT_luout_),102)
write(iparm(_JT_luout_),100) iparm(_JT_iter_), rnorm, err, rnormt, errt
else
if (iparm(_JT_iter_).eq.1 .or. iparm(_JT_out_).ge.5) write(iparm(_JT_luout_),101)
write(iparm(_JT_luout_),100) iparm(_JT_iter_), rnorm, err
endif
call JT_FlushUnit (iparm(_JT_luout_), status)
endif
c
c ... Test for convergence using appropriate estimate of
c error. If converged, set status, write final solution
c estimate if desired, and return. If not, increment
c iteration counter and check to see whether iparm(_JT_itmax_) has
c been exceeded.
if (iparm(_JT_stop_) .lt. 0) then
err_test = errt
else
err_test = err
endif
c
if (err_test .le. rparm(_JT_eps_)) then
if (iparm(_JT_out_) .ge. 4) then
call JT_WriteVectorFloat (iparm(_JT_luout_), ia(_JT_nrows_), x,
& 'JT_CheckConvergence: AFTER CONVERGENCE (x):', status)
endif
status = 0
else
status = 1
iparm(_JT_iter_) = iparm(_JT_iter_)+1
if (iparm(_JT_iter_) .gt. iparm(_JT_itmax_)) then
iparm(_JT_iter_) = iparm(_JT_iter_)-1
if (iparm(_JT_out_) .ge. 1) then
write(iparm(_JT_luerr_),*)
write(iparm(_JT_luerr_),901) '************************************'//
& '*************************************'
write(iparm(_JT_luerr_),901) 'JT_CheckConvergence:'
write(iparm(_JT_luerr_),901)
& ' iteration failed to converge in itmax iterations'
write(iparm(_JT_luerr_),905) ' itmax = ', iparm(_JT_itmax_)
write(iparm(_JT_luerr_),940) ' convergence criterion = ', rparm(_JT_eps_)
write(iparm(_JT_luerr_),940) ' parameters at termination:'
if (iparm(_JT_stop_) .ne. 0) then
write(iparm(_JT_luerr_),940) ' norm of recursively-'//
& 'computed residual = ', rnorm
write(iparm(_JT_luerr_),940) ' error estimate using recursively-'//
& 'computed residual = ', err
else
write(iparm(_JT_luerr_),940) ' relative change in successive'//
& ' iterations = ', err
endif
if (iparm(_JT_resid_).eq.1 .or. iparm(_JT_resid_).eq.3 .or. iparm(_JT_stop_).lt.0) then
write(iparm(_JT_luerr_),940) ' norm of b-Ax = ', rnormt
write(iparm(_JT_luerr_),940) ' error estimate using b-Ax = ', errt
endif
if (iparm(_JT_stop_) .eq. 0) then
write(iparm(_JT_luerr_),901) ' convergence based on relative'//
& ' change in successive iterations '
elseif (iparm(_JT_stop_) .lt. 0) then
write(iparm(_JT_luerr_),901) ' convergence based on error'//
& ' estimate using b-Ax'
else
write(iparm(_JT_luerr_),901) ' convergence based on error estimate'//
& ' using recursively-computed residual'
endif
write(iparm(_JT_luerr_),901) '************************************'//
& '*************************************'
endif
status = -4
endif
endif
c
call JT_Clock (1, cpu, dummy)
c
100 format (5x, i6, 2x, 2(3x,e12.5), 2x, 2(3x,e12.5))
101 format (t8, 'iter', t21, 'rnorm', t37, 'err')
102 format (t8, 'iter', t21, 'rnorm', t37, 'err',
& t52, 'rnormt', t68, 'errt')
#include "formats.inc"
return
end