c===========================================================
c     newtsqrt:  Uses Newton's method to find (positive)
c     square root of number supplied on command line, i.e.
c     solves 
c 
c     f(x) = x^2 - a = 0
c
c     for given 'a'.  Optional second argument specifies 
c     convergence criteria (relative dx). Tracing output
c     includes iteration number, estimated root (xn), 
c     change in estimate (dxn), log(dxn), residual and 
c     log(residual).  Optional second argument
c     specifies convergence criteria (relative dx).
c===========================================================
      program         newtsqrt

      implicit        none

      integer         iargc
      real*8          r8arg,        drelabs

      real*8          r8_never
      parameter     ( r8_never = -1.0d-60 )

      integer         mxiter
      parameter     ( mxiter = 50 )

      real*8          a,            xtol

      integer         iter
      real*8          xn,           resn,           dxn

c-----------------------------------------------------------
c     Argument parsing.
c-----------------------------------------------------------
      if( iargc() .lt. 1 ) go to 900
      a     = r8arg(1,r8_never)
      if( a .eq. r8_never .or. a .lt. 0.0d0 ) go to 900
      xtol  = r8arg(2,1.0d-8)
      if( xtol .le. 0.0d0 ) xtol = 1.0d-8

c-----------------------------------------------------------
c     Un-inspired initial guess: x^(0) = a / 2.
c-----------------------------------------------------------
      xn = 0.5d0 * a
 
c-----------------------------------------------------------
c     Newton loop.
c-----------------------------------------------------------
      write(0,*) 'Iter          xn                 '//
     &           'dxn       log(dxn)    rn         log(rn)'
      write(0,*)
      do iter = 1 , mxiter
         resn = xn**2 - a
         dxn  = resn / (2.0d0 * xn)
         xn   = xn - dxn
         write(0,1000) iter, xn, dxn, log(abs(dxn)),
     &                 resn, log(abs(resn))
1000     format(i2,1p,e26.16,e12.3,0p,f10.2,1p,e12.3,0p,f10.2)
c-----------------------------------------------------------
c        Check for convergence.
c-----------------------------------------------------------
         if( drelabs(dxn,xn,1.0d-6) .le. xtol ) go to 100
      end do
c-----------------------------------------------------------
c     No-convergence exit.
c-----------------------------------------------------------
      write(0,*) 'No convergence after ', mxiter,  
     &           ' iterations'
      stop

c-----------------------------------------------------------
c     Normal exit, write input and estimated square root 
c     to standard output.
c-----------------------------------------------------------
 100  continue
      write(0,*)
      write(*,*) a, xn
      stop

 900  continue
         write(0,*) 'usage: newtsqrt <a> [<xtol>]'
      stop

      end

c-----------------------------------------------------------
c     drelabs:  Function useful for 'relativizing' quantity
c     being monitored for detection of convergence.  
c-----------------------------------------------------------
      real*8 function drelabs(dx,x,xfloor)

         implicit      none

         real*8        dx,     x,    xfloor

         if( abs(x) .lt. abs(xfloor) ) then
            drelabs = abs(dx)
         else
            drelabs = abs(dx/x)
         end if

         return

      end