RelativeDifference.F

c**********************************************************************
#include "author.inc"
c*    $Id: RelativeDifference.F,v 1.1 1996/02/27 19:29:12 turner Exp $
c*
c*    Computes relative difference in two vectors.
c*
c*    WARNING: Note that JT_RelativeDifference *must* be declared real
c*             in routines that use it.
c*
c*    <PARAMETER LIST>
c*
c*     Input:
c*      norm - determines which vector norm to use
c*         0  ==>  infinity norm
c*         1  ==>  1-norm
c*         2  ==>  2-norm
c*      n - number of elements
c*      x - x-vector
c*      y - y-vector
c*
c*     In/Out:
c*
c*     Output:
c*      JT_RelativeDifference - relative difference
c*      status - return status
c*
c*    <SUBROUTINES REQUIRED>
c*
c*     JT_x_eq_sx
c*     JT_y_eq_y_div_x
c*     JT_z_eq_y_minus_x
c*     JT_z_eq_y_plus_x
c*     
c*    <FUNCTIONS REQUIRED>
c*
c*     JT_VectorNorm
c*
#include "arrays-RelativeDifference.inc"
c*
#include "copyright.inc"
c**********************************************************************
      real function JT_RelativeDifference (norm, n, x, y, status)
      implicit none
c
c ... Input:
      integer n, norm
      real x(n), y(n)
c
c ... Output:
      integer status
c
c ... Local:
      real JT_VectorNorm, two
#include "declare-RelativeDifference.inc"
c
      parameter (two=2.0d0)
c
#include "allocate-RelativeDifference.inc"
c
c ... Compute z = 2(x-y).
      call JT_z_eq_y_minus_x (n, x, y, z, status)
      call JT_x_eq_sx (n, two, z, status)
c
c ... Compute w = x+y.
      call JT_z_eq_y_plus_x (n, x, y, w, status)
c
c ... Compute z = z/w.
      call JT_y_eq_y_div_x (n, w, z, status)
c
c ... Find max. abs. value of z.
      JT_RelativeDifference = JT_VectorNorm(norm, n, z, status)
c
 9999 continue
#include "deallocate-RelativeDifference.inc"
      return
      end