StoppingTest.F

#include "parmdefines.h"
#include "iadefines.h"
c**********************************************************************
#include "author.inc"
c*    $Id: StoppingTest.F,v 1.9 1996/04/28 20:15:27 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*      err - error estimate
c*      errt - error estimate based on true residual
c*      status - return status:
c*        -3  ==>  internal error
c*        -2  ==>  memory allocation failure
c*         0  ==>  convergence was achieved
c*
#include "parmdesc.inc"
c*
#include "iadesc.inc"
c*
c*    <SUBROUTINES REQUIRED>
c*
c*     JT_EstimateError
c*     JT_WriteVectorFloat
c*     JT_y_eq_x
c*     JT_y_eq_y_minus_Ax
c*     JT_y_eq_y_minus_x
c*
c*    <FUNCTIONS REQUIRED>
c*
c*     JT_VectorNorm
c*
#include "arrays-StoppingTest.inc"
c*
#include "copyright.inc"
c**********************************************************************
      subroutine JT_StoppingTest (bnorm, rnormt0, rnorm, x, xold, b,
     &           a, ia, ja, iparm, rparm, 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 err, errt
c
c ... Local:
      real JT_VectorNorm, xnorm, zero
#include "declare-StoppingTest.inc"
c
      parameter (zero=0.0d0)
c
#include "allocate-StoppingTest.inc"
c
c ... Calculate || x || if needed.
      if (iparm(_JT_stop_).eq._JT_stop_relchg_ .or.
     &    ABS(iparm(_JT_stop_)).eq._JT_stop_axb_ .or.
     &    ABS(iparm(_JT_stop_)).eq._JT_stop_x_) then
       xnorm = JT_VectorNorm(iparm(_JT_norm_), ia(_JT_nrows_), x, status)
      endif
c
c ... Calculate error estimate based on || b - Ax || as appropriate.
      if (iparm(_JT_resid_).eq.1 .or.
     &    iparm(_JT_resid_).eq.3 .or.
     &    iparm(_JT_stop_).lt.0) then
c
       if (rnormt .eq. zero) then
c
c       || b - Ax || was not input, so compute.
        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 .ne. -2) status = -3
         goto 9999
        endif
        rnormt = JT_VectorNorm(iparm(_JT_norm_), ia(_JT_nrows_), w, status)
c
c       Write out b-Ax, if appropriate.
        if (iparm(_JT_out_) .ge. 5) then
         call JT_WriteVectorFloat (iparm(_JT_luout_), ia(_JT_nrows_), w,
     &       'JT_StoppingTest: b - Ax (w):', status)
        endif
       endif
c
       if (iparm(_JT_stop_) .eq. _JT_stop_relchg_) then
c
c       Stopping test is || x - xold || / || x ||, so do not
c       compute an error estimate based on rnormt.
        errt = zero
       else
c
        call JT_EstimateError (ABS(iparm(_JT_stop_)), rparm(_JT_anorm_),
     &       xnorm, bnorm, rnormt, rnormt0, errt, status)
       endif
      endif
c
c ... Compute standard error estimate, based on relative change in
c     unknowns or recursively-generated residual.
      if (iparm(_JT_stop_) .eq. _JT_stop_relchg_) then
c
c      Use || x - xold || / || x || as stopping test.
       call JT_y_eq_y_minus_x (ia(_JT_nrows_), x, xold, status)
       err = JT_VectorNorm(iparm(_JT_norm_), ia(_JT_nrows_), xold, status)
       if (xnorm .ne. zero) err = err/xnorm
       call JT_y_eq_x (ia(_JT_nrows_), x, xold, status)
      else
c
c      Compute error estimate based on || r ||.
       call JT_EstimateError (ABS(iparm(_JT_stop_)), rparm(_JT_anorm_),
     &      xnorm, bnorm, rnorm, rnormt0, err, status)
      endif
c
 9999 continue
#include "deallocate-StoppingTest.inc"
c
      return
      end