program main

c*********************************************************************72
c
cc MAIN is the main program for NMS_PRB.
c
c  Discussion:
c
c    NMS_PRB tests the NMS library.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    14 February 2012
c
c  Author:
c
c    John Burkardt
c
      implicit none

      call timestamp ( )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'NMS_PRB'
      write ( *, '(a)' ) '  FORTRAN77 version'
      write ( *, '(a)' ) '  Test the NMS library.'

      call test01 ( )
      call test02 ( )
      call test03 ( )
      call test04 ( )
      call test05 ( )
      call test06 ( )
      call dpchez_test ( )
      call test08 ( )
      call test09 ( )

      call test10 ( )
      call test11 ( )
      call test12 ( )
      call test13 ( )
      call test14 ( )
      call test15 ( )
      call test16 ( )
      call test17 ( )
      call test18 ( )
      call test19 ( )

      call test20 ( )
      call test21 ( )
c
c  Terminate.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'NMS_PRB'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ' '
      call timestamp ( )

      stop
      end
      subroutine test01 ( )

c*********************************************************************72
c
cc TEST01 tests DCFT2D.
c
C EXAMPLE 11.8:  Plot image and transform of 8 by 8 unit source
C in 64 by 64 (otherwise zero) array.
C
      parameter (n=64) 
      double precision  a(n,n),a2(n,n),err,ermax,data
      complex*16  image(n,n),image2(n,n),w(3*n+8)
      logical   forwd

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST01'
      write ( *, '(a)' ) '  DCFT2D computes a 2D Fourier transform.'
      write ( *, '(a)' ) ' '

      WRITE (*,*) 'COMPUTING...'
      LDA = N
c
C Set up data. IMAGE is original. IMAGE2 is scaled by (-1)**(I+J)
C to place Fourier coefficients in the correct place for viewing.
      DO I = 1,N
        DO J = 1,N
          A(I,J) = 0.0D0
          IF (I.GE.(N/2-4) .AND. I.LE.(N/2+4) .AND. J.GE.(N/2-4)
     *         .AND. J.LE.(N/2+4)) A(I,J) = 1.0D0
          IMAGE(I,J) = DCMPLX(A(I,J),0.0D0)
          IMAGE2(I,J) = IMAGE(I,J)*(-1)**(I-1+J-1)
        end do
      end do
C
C Plot image  with axonometric plotting routine.
C
c      CALL AXNPLT (A,LDA,N)
C
C Compute Fourier transform of IMAGE and IMAGE2
C
      FORWD = .TRUE.
      CALL DCFT2D (N,IMAGE,LDA,W,FORWD)
      CALL DCFT2D (N,IMAGE2,LDA,W,FORWD)
C
C Compute magnitude of components of transforms.
C Actual transforms are unscaled and need to be divided by N*N
C to be correct.
C
      DO 2 I = 1,N
        DO 2 J = 1,N
           A(I,J) = ABS(IMAGE(I,J))
           A2(I,J) = ABS(IMAGE2(I,J))
    2 CONTINUE
C
C PLOT MODULUS OF TRANSFORM
C
c    CALL AXNPLT (A,LDA,N)
c      CALL AXNPLT (A2,LDA,N)
C
C Compute inverse transform of image and see if it agrees with original.
C
      FORWD = .FALSE.
      CALL DCFT2D (N,IMAGE,LDA,W,FORWD)
      ERMAX = 0.0D0
      DO 3 I = 1,N
        DO 3 J = 1,N
           DATA = 0.0D0
           IF (I.GE.(N/2-4) .AND. I.LE.(N/2+4) .AND. J.GE.(N/2-4)
     *         .AND. J.LE.(N/2+4)) DATA = 1.0D0
           ERR = ABS( DATA-ABS(IMAGE(I,J))/(N*N) )
           IF (ERR .GT. ERMAX) ERMAX = ERR
    3 CONTINUE
      WRITE (*,*) 'DCFT2D RESULTS (EX 11.8: PLOTS HAVE BEEN SKIPPED)'
      WRITE (*,'(2X,A,1X,D18.10)') 'MAXIMUM ERROR IS ',ERMAX

      WRITE (*,*)
      WRITE (*,*) 'REFERENCE RESULTS FROM IBM PC/AT'
      WRITE (*,*) ' MAXIMUM ERROR IS    0.1885421720E-15'
      RETURN
      END
      SUBROUTINE AXNPLT(A,LDA,N)

c*********************************************************************72
c
cc AXNPLT creates an axonometric plot of a 2d array.
C
C       N: SIZE OF A, N.LE.350
C       LDA: DIMENSION OF FIRST COMPONENT OF A IN CALLING PROGRAM
C
C       ASSUMPTIONS:    0.0 .LE. A(I,J) .LE 1.0
C
      INTEGER LDA,J,K,IDAMAX,N,IRET,I
      double precision A(LDA,*)
      double precision WX(350)
      double precision WY(350),X(350),Y(350),YMX(350),X0,DEL,YMAX
      double precision YM,YX
C
C  FIND THE LARGEST Y SO CAN SCALE
c
      do j=1,n
        k=idamax(n,a(1,j),1)
        ymx(j)=a(k,j)
      end do

      k=idamax(n,ymx,1)
      ymax=ymx(k)
      x0=0.1
      del=.4
c      call agraf0(iret,0)
c      call minmax(iret,0.0, 1.0, 0.0, 1.0, ym, yx)
c      call howplt(iret,0,1,15)
      do 2 i=1,n
        do j=1,n
           x(j)=(j-1.)/(n-1.)*del+x0+(i-1.)/(n-1.)*del
           y(j)=(a(i,j)/ymax)*del+x0+(i-1.)/(n-1.)*del
        end do
        if(i.eq.1)then
c           call nowait
c           call plot1(iret,n,x,y,wx,wy)
        else
c           if(i.eq.n)call wait
c           call howplt(iret,0,1,15)
c           call plot(iret,n,x,y,wx,wy)
        endif
    2 continue
      return
      END
      subroutine test02 ( )

c*********************************************************************72
c
cc TEST02 tests DFZERO.
c
c  Modified:
c
c    04 February 2012
c
      double precision b, c, ae, re, test02_f
      external test02_f

      b   = 2.0d0
      c   = 3.0d0
      ae  = 1.0d-06
      re  = 1.0d-06

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST02:'
      write ( *, '(a)' ) '  DFZERO finds the root of a function.'
      write ( *, '(a)' ) ' '

      write (*,'(a)') '  Initial interval: '
      write (*,'(8x,2d18.10)') b, c
      write (*,'(a)') '  Tolerances:       '
      write (*,'(8x,2d18.10)') ae, re

      call dfzero ( test02_f, b, c, c, re, ae, iflag)

      write (*,'(a)') ' '
      write (*,'(a)') '  DFZERO results'
      if (iflag .ne. 1) write (*,'(a)') '  Error code =', iflag
      write (*,'(a19,2x,d18.10)') '  Estimate of zero =', b
      write (*,'(a19,2x,d18.10)') '  Function value =  ', test02_f(b)

      write (*,'(a)') ' '
      write (*,'(a)') '  Reference results:'
      write (*,'(a)') '  Estimate of zero =   0.2094551435e+01'    
      write (*,'(a)') '  Function value =    -0.5157851953e-06'

      return
      end
      double precision function test02_f ( x )

c*********************************************************************72
c
cc TEST02_F is a function needed for TEST02_DFZERO.
c
c  Modified:
c
c    04 February 2012
c
      double precision x

      test02_f = x * (x*x-2.0d0) - 5.0d0

      return
      end
      subroutine test03 ( )

c*********************************************************************72
c
cc TEST03 tests DGEFS
c
c  Modified:
c
c    04 February 2012
c
      integer   lda
      parameter (lda=10)
      double precision  a(lda,lda), b(lda), work(lda), rcond
      integer   iwork(lda), i, j, n, itask, ind

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST03:'
      write ( *, '(a)' ) '  DGEFS factors and solves a linear system.'
      write ( *, '(a)' ) ' '
C
C Set up problem
c
      n      =  3
      itask  =  1
      a(1,1) = 10.0d0
      a(2,1) = -3.0d0
      a(3,1) =  5.0d0
      a(1,2) = -7.0d0
      a(2,2) =  2.0d0
      a(3,2) = -1.0d0
      a(1,3) =  0.0d0
      a(2,3) =  6.0d0
      a(3,3) =  5.0d0
      b(1)   =  7.0d0
      b(2)   =  4.0d0
      b(3)   =  6.0d0
C
C Print problem information
C
      write (*,*) '  Coefficient matrix ='
      do i = 1,n
         write (*,800) (a(i,j), j = 1,n)
      end do
      write (*,*) '  Right-hand side ='
      write (*,800) (b(j), j = 1,n)
C
C Solve linear system
c
      call dgefs (a, lda, n, b, itask, ind, work, iwork, rcond)
C
C Print results
c
      write (*,*)
      write (*,*) '  DGEFS results'
      IF (IND .EQ. -10) THEN
         WRITE (*,*) ' ERROR CODE =', IND
      ELSE IF (IND .LT.0) THEN
         WRITE (*,*) '  ERROR CODE =', IND
         STOP
      ELSE
         WRITE (*,*) '  Number of accurate digits =', IND
      END IF

      WRITE (*,*) '  SOLUTION ='
      WRITE (*,800) (B(J), J = 1,N)

      WRITE (*,*)
      WRITE (*,*) '  REFERENCE RESULTS FROM IBM PC/AT '
      WRITE (*,*) '   NUMBER OF ACCURATE DIGITS =          14'
      WRITE (*,*) '   SOLUTION ='
      WRITE (*,*) '     0.000000000000E+00     -0.100000000000E+01',    
     *'      0.100000000000E+01'

      return
800   FORMAT (4X, D20.12, 4X, D20.12, 4X, D20.12)
      END
      subroutine test04 ( )

c*********************************************************************72
c
cc TEST04 tests DUNI
c
c  Modified:
c
c    04 February 2012
c
      DOUBLE PRECISION  U,DUNI,DUSTAR,USEED
      INTEGER ISEED,I

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST04'
      write ( *, '(a)' ) '  DUNI computes uniform random numbers.'
c
C SET INITIAL SEED
c
      ISEED = 305
      USEED = DUSTAR(ISEED)
c
C USTART RETURNS FLOATING ECHO OF ISEED 
c
      WRITE (*,*) '  DUNI RESULTS '
      WRITE (*,'(5X,I5,5X,D15.8)')ISEED, USEED
      DO I = 1,1000
        U = DUNI()
      end do
      WRITE (*,'(12X,D25.15)') U
      WRITE (*,*) '  REFERENCE RESULTS FROM AN IBM P/C'
      WRITE (*,*) '        305      0.30500000E+03' 
      WRITE (*,*) '                 0.241576136599321E+00'

      return
      END
      subroutine test05 ( )

c*********************************************************************72
c
cc TEST05 tests DCFFTI and DCFFTF.
c
c  Modified:
c
c    04 February 2012
c
C Using complex discrete Fourier transform, find the approximate Fourier
C coefficients to Runge's function on [-1,1] with N=16 and N=17. 
C 
       DOUBLE PRECISION  WSAVE(150),RUNGE,X,X0,PI,DEL,F,XJ
       COMPLEX*16 COEFF(0:16), SQTM1
C
C Note: Complex Double Precision is not Fortran Standard and many
C       compiler vendors may fail to implement it or implement it
C       in some other manner.  If the COMPLEX*16 declaration does
C       not work on your compiler we recommend that you check
C       with the vendor's instructions and modify the above accordingly.
C
C Note 0 subscript, which makes indexing easier, allowed in Fortran 77.
C
C Arithmetic statement function for Runge's function. 
       RUNGE(X) = 1.0D0/(1.0D0+25.0D0*X*X) 

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST05'
      write ( *, '(a)' ) '  DCFFTI initializes a complex FFT.'
      write ( *, '(a)' ) '  DCFFTF computes it.'

       X0 = -1.0D0 
       PI = ASIN(1.0D0)*2.0D0 
       SQTM1 = CMPLX(0.0D0,-1.0D0)

       DO 10 N = 16,17 
          CALL DCFFTI (N,WSAVE) 
C Function assumed to be periodic on [-1,1], of length 2.
          DEL = 2.0D0/N 
          F = 2.0D0*PI/(N*DEL) 
          DO 1 J = 0,N-1
C First sample point at -1, last at 1-DEL
              XJ = (-1.0D0) + J*DEL 
              COEFF(J) = CMPLX(RUNGE(XJ),0.0D0) 
    1     CONTINUE 
          CALL DCFFTF (N,COEFF,WSAVE) 
C Returned coefficients must be divided by N for correct normaliziation.
C 
C Note repetition after N/2 in original coefficients. Scaling because
C X0 not at origin destroys this to some extent. 
C
          WRITE (*,*) ' DCFFTF RESULTS FOR N = ' ,N 
          WRITE (*,'(A,1X,2D17.8)') ' CZERO = ',COEFF(0)/N*2.0D0
          WRITE (*,*)  
     *' J          OUTPUT FROM DCFFTF,         SCALED COEFFICIENTS'   
          DO 11 J = 1,N-1
              WRITE (*,'(I5,2D17.8,5X,2D17.8)') 
     *              J, COEFF(J), EXP(-SQTM1*J*F*X0) * COEFF(J)/N *2
   11     CONTINUE 
          WRITE (*,*)

   10  CONTINUE 
       WRITE (*,*) 'REFERENCE RESULTS (PARTIAL) FROM IBM PC/AT'
       WRITE (*,*) ' DCFFTF RESULTS FOR N =           17'
       WRITE (*,*) ' CZERO =  (0.54916124E+00, 0.00000000E+00)'
       WRITE (*,*) '    .  '
       WRITE (*,*) '    .  '
       WRITE (*,*) '  10  -0.60345452E-01  -0.37837322E-16 '  
       WRITE (*,*) '  11   0.11248162E+00  -0.47424091E-16 '  
       WRITE (*,*) '  12  -0.23457294E+00   0.73973161E-16 '  
       WRITE (*,*) '  13   0.42197537E+00  -0.25268310E-17 '  
       WRITE (*,*) '  14  -0.82441735E+00  -0.14952867E-16 '
       WRITE (*,*) '  15   0.14919178E+01   0.31242621E-16 '
       WRITE (*,*) '  16  -0.29260639E+01  -0.69144627E-16 '
      return
      END  
      subroutine test06 ( )

c*********************************************************************72
c
cc TEST06 tests DRNOR.
c
c  Modified:
c
c    06 February 2012
c
C HISTOGRAM FOR DRNOR
C
      DOUBLE PRECISION A,B
      PARAMETER(NBINS=32,A=-3.0D0,B=3.0D0)
      INTEGER ISEED,I,J,H(NBINS),NR,test06_INBIN
      DOUBLE PRECISION  R,DRNOR,DSTART,RSEED,WIDTH

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST06'
      write ( *, '(a)' ) '  DRNOR computes normal random numbers.'

      ISEED = 305
      RSEED = DSTART(ISEED)
      WIDTH=(B-A)/(NBINS-2)
CCCCCC  100 WRITE (*,*) 'EX 10.1: ENTER NUMBER OF NORMALS : '
      DO 200 I = 1,NBINS
          H(I)=0
  200 CONTINUE
CCCCC      READ (*,*) NR
      WRITE (*,*) 'RUNNING 10,000 NORMALS INTO 32 BINS...'
      NR = 10000
      IF (NR .LE. 0) return
      DO I = 1,NR 
          R = DRNOR()
          J = test06_INBIN (A,B,NBINS,WIDTH,R)
          H(J) = H(J)+1
      end do

      WRITE (*,*) 'HISTOGRAM FOR DRNOR: NUMBER IN BIN 1,...,32'
      WRITE (*,*) '   (-INFINITY,-3],(-3,-2.8],...,(2.8,3],(3,INFINITY)'
      WRITE (*,*) ' (VALUES ARE SLIGHTLY COMPUTER DEPENDENT)'
      WRITE (*,'(9I8)') (H(I),I=1,NBINS)

      WRITE (*,*)
      WRITE (*,*)
     *   'REFERENCE RESULTS FROM IBM PC/AT '
      WRITE(*,*)'     16      14      21      32      46      96     135
     *      198     292'
      WRITE(*,*)'    334     454     549     611     665     743     801
     *      785     751'
      WRITE(*,*)'    747     634     503     421     351     274     184
     *      124      90'
      WRITE(*,*)'     53      35      19      13       9'
      return
      END 
      INTEGER FUNCTION test06_INBIN (XMIN,XMAX,NBINS,WIDTH,DATA)

c*********************************************************************72
c
cc TEST06_INBIN bins data for TEST06.
C
C THIS FUNCTION TAKES A DOUBLE PRECISION VALUE IN DATA, AND FINDS
C THE CORRECT BIN FOR IT.  VALES BELOW XMIN COME BACK
C IN 1.  VALUES ABOVE XMAX COME BACK IN NBINS.
C
      integer nbins 
      double precision xmin,xmax,width,data

      if (data .lt. xmin) then
         test06_inbin = 1
      else if (data .ge. xmax) then
         test06_inbin = nbins
      else
         test06_inbin = 2 + (data-xmin)/width
      endif

      return
      end
      subroutine dpchez_test ( )

c*********************************************************************72
c
cc DPCHEZ_TEST tests DPCHEZ.
c
c  Modified:
c
c    06 February 2012
c
c  Reference:
c
c    Fred Fritsch, Ralph Carlson,
c    Monotone Piecewise Cubic Interpolation,
c    SIAM Journal on Numerical Analysis,
c    Volume 17, Number 2, April 1980, pages 238-246.
c
      implicit none

      integer n
      parameter ( n = 21 )
      integer ne
      parameter ( ne = 101 )
      integer nwk
      parameter ( nwk = 2 * n )

      double precision a
      double precision b
      double precision d(n)
      double precision dpchqa
      double precision edans(4)
      double precision erans(4)
      double precision error
      double precision errord
      double precision f(n)
      double precision fans(4)
      double precision fd(ne)
      double precision fe(ne) 
      integer i
      integer ians
      integer ierr
      double precision q
      double precision qans
      double precision r
      double precision rp
      logical spline
      double precision u
      double precision wk(nwk)
      double precision x(n)
      double precision xans(4)
      double precision xe(ne)

      save edans
      save erans
      save fans
      save ians
      save qans
      save xans

      data edans /
     &  0.2932883472d0, 0.2677039506d0, 0.1744906562d0, 0.0d0 /
      data erans/
     &  -0.6075110024d-02, -0.3219009901d-02,
     &  -0.946234414d-03, 0.0d0 /
      data fans / 0.971920d0, 0.986880d0, 0.996560d0, 1.0d0 /
      data ians / 0 /
      data qans / 0.274679262701527d0 /
      data xans / -0.3d-01, -0.2d-01, -0.1d-01, 0.0d0 /
c
c  Arithmetic statement functions for Runge's function and derivative.
c
      r ( u ) = 1.0D0 / ( 1.0D0 + 25.0D0 * u * u )
      rp ( u ) = -50.0D0 * u * r ( u ) ** 2

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'DPCHEZ_TEST'
      write ( *, '(a)' ) 
     &  '  DPCHEZ sets up piecewise cubic Hermite splines.'
      write ( *, '(a)' ) ' '
c
c  Compute Runge's function at 21 points in [-1,1].
c
      do i = 1, n
        x(i) = -1.0D+00 + dble ( i - 1 ) / 10.0D+00
        f(i) = r ( x(i) )
      end do
 
      spline = .false.
c
c  Compute cubic Hermite interpolant because SPLINE is .FALSE.
c
      call dpchez ( n, x, f, d, spline, wk, nwk, ierr )

      if ( ierr .lt. 0 ) then
        write ( *, '(a)' ) ''
        write ( *, '(a)' ) 'DPCHEZ_TEST - Fatal error!'
        write ( *, '(a,i6)' ) 
     &    '  DPCHEZ returned error flag IERR= ', ierr
        stop 1
      end if
c
c  Evaluate interpolant and derivative at 101 points from -1 to 0.
c
      do i = 1, ne
        xe(i) = -1.0D0 + dble ( i - 1 ) / dble ( ne - 1 )
      end do

      call dpchev ( n, x, f, d, ne, xe, fe, fd, ierr )

      if ( ierr .ne. 0 ) then
        write ( *, '(a)' ) ''
        write ( *, '(a)' ) 'DPCHEZ_TEST - Fatal error!'
        write ( *, '(a,i6)' ) 
     &    '  DPCHEV returned error flag IERR= ', ierr
        stop 1
      end if

      do i = 1, ne
        error = fe(i) - r(xe(i))
        errord = fd(i) - rp(xe(i))
        write ( *, 3 ) xe(i), fe(i), error, errord
    3   FORMAT(1X,D17.10,3X,D17.10,3X,D17.10,3X,D17.10)
      end do
c
c  Compute integral over the interval [0,1] 
c
      a = 0.0D+00
      b = 1.0D+00
      q = dpchqa ( n, x, f, d, a, b, ierr )

      write ( *, '(2X,A,D20.12,3X,A,I5)' ) 
     &  'Integral from 0 to 1 = ', q, ' IERR = ', ierr

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  Reference results:'

      do i = 1, 4
        write ( *, 5 ) xans(i), fans(i), erans(i), edans(i)
    5   FORMAT(1X,D17.10,3X,D17.10,3X,D17.10,3X,D17.10)
      end do

      write ( *, '(2X,A,D20.12,3X,A,I5)' ) 
     &  'Integral from 0 to 1 = ', qans, ' IERR = ', ians           

      return
      end
      subroutine test08 ( )

c*********************************************************************72
c
cc TEST08 tests DFMIN.
c
      DOUBLE PRECISION A
      double precision B, XSTAR, TOL, DFMIN, test08_f
      EXTERNAL test08_f

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST08'
      write ( *, '(a)' ) '  DFMIN minimizes a function.'
      write ( *, '(a)' ) ' '

      A     = 0.1D0
      B     = 0.9D0
      TOL   = 1.0D-5
      XSTAR = DFMIN ( A, B, test08_f, TOL )
      WRITE (*,*) 'DFMIN RESULTS'
      WRITE (*,'(2X,A,D18.10)') ' XSTAR =', XSTAR
      WRITE (*,*)
      WRITE (*,*) 'REFERENCE RESULTS FROM IBM PC/AT'
      WRITE (*,*) '  XSTAR =  0.8164961769E+00'    
      return
      END
      DOUBLE PRECISION FUNCTION test08_f(X)

c*********************************************************************72
c
cc TEST08_F is the function to be minimized.
c
      DOUBLE PRECISION X
      test08_f = X*(X*X - 2.0D0) - 5.0D0
      RETURN
      END
      subroutine test09 ( )

c*********************************************************************72
c
cc TEST09 tests DSNQE.
c
      PARAMETER   (N = 2, LW = 19)
      DOUBLE PRECISION   X(N), FVEC(N), W(LW), TOL
      EXTERNAL    test09_f

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST09'
      write ( *, '(a)' ) 
     &  '  DNSQE solves a system of nonlinear equations.'
      write ( *, '(a)' ) ' '

      TOL  = 1.0D-6
      X(1) = 2.0D0
      X(2) = 3.0D0
      WRITE (*,800) X(1), X(2)
      IOPT = 2
      NPRINT = 0
C
C SOLVE NONLINEAR EQUATIONS
C
      CALL DNSQE ( test09_f, test09_f, IOPT, N, X, FVEC, TOL, NPRINT, 
     &  INFO, W, LW)
C
C PRINT RESULTS
C
      WRITE (*,*)
      WRITE (*,*) '  DNSQE RESULTS'
      IF (INFO .NE. 1) WRITE (*,810) INFO
      WRITE (*,820) X(1), X(2)
      WRITE (*,830) FVEC(1), FVEC(2)

      WRITE (*,*)
      WRITE (*,*) 'REFERENCE RESULTS FROM IBM PC/AT'
      WRITE (*,*) ' ESTIMATE OF SOLUTION '
      WRITE (*,*) '     0.199999999997E+01  0.100000000001E+01'
      WRITE (*,*) ' VALUES OF NONLINEAR FUNCTIONS'
      WRITE (*,*) '    -0.415920631269E-10 -0.940900690694E-10'

      return
800   FORMAT (' INITIAL GUESS', /, 4X,2D20.12)
810   FORMAT (' INFO =', I3)
820   FORMAT (' ESTIMATE OF SOLUTION ', /, 4X,2D20.12)
830   FORMAT (' VALUES OF NONLINEAR FUNCTIONS', /, 4X,2D20.12)
      END
      SUBROUTINE test09_f (N, X, FVEC, IFLAG)

c*********************************************************************72
c
cc TEST09 evaluates the nonlinear equations.
c
      DOUBLE PRECISION   X(N), FVEC(N)
C
C COMPUTE NONLINEAR FUNCTIONS
C
      FVEC(1) = X(1)*X(2) - X(2)**3 - 1.0D0
      FVEC(2) = X(1)**2*X(2) + X(2) - 5.0D0

      RETURN
      END
      subroutine test10 ( )

c*********************************************************************72
c
cc TEST10 tests DQRLS.
c
      PARAMETER (MM = 5, NN = 3)
      DOUBLE PRECISION A(MM,NN), B(MM), X(NN), QRAUX(NN), WORK(NN), TOL
      INTEGER   JPVT(NN)
      DATA      B / 1.0D0, 2.3D0, 4.6D0, 3.1D0, 1.2D0/

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST10'
      write ( *, '(a)' ) '  DQRLS solves a least squares problem.'
      write ( *, '(a)' ) ' '
C
C SET UP LEAST-SQUARES PROBLEM
C     QUADRATIC MODEL, EQUALLY-SPACED POINTS
C
      M = 5
      N = 3
      DO I = 1,M
         A(I,1) = 1.0D0
         DO J = 2,N
            A(I,J) = A(I,J-1)*I
         end do
      end do

      TOL = 1.D-15
      ITASK = 1
      WRITE (*,*)   ' COEFFICIENT MATRIX'
      WRITE (*,800) ((A(I,J), J = 1,N), I = 1,M)
      WRITE (*,*)   ' RIGHT-HAND SIDE'
      WRITE (*,810) (B(I), I = 1,M)
C
C
C SOLVE LEAST-SQUARES PROBLEM
C
      CALL DQRLS (A, MM, M, N, TOL, KR, B, X, B, JPVT, QRAUX, WORK,
     *            ITASK, IND)
      IF (IND .NE. 0) THEN
         WRITE (*,*) ' ERROR CODE =', IND
         STOP
      END IF
      WRITE (*,*)    ' RANK OF MATRIX =', KR
C
C PRINT RESULTS
C
      WRITE (*,*)   ' PARAMETERS'
      WRITE (*,800) (X(J), J = 1,N)
      WRITE (*,*)   ' RESIDUALS'
      WRITE (*,810) (B(I), I = 1,M)
      WRITE (*,*)
      WRITE (*,*) 'REFERENCE RESULTS FROM IBM PC/AT'
      WRITE (*,*) ' RANK OF MATRIX =           3'
      WRITE (*,*) ' PARAMETERS'
      WRITE (*,*) '   -0.30200000E+01   0.44914286E+01  -0.72857143E+00'
      WRITE (*,*) ' RESIDUALS'
      WRITE (*,*) 
     *'   0.2571429E+00 -0.7485714E+00  0.7028571E+00 -0.1885714E+00
     *-0.2285714E-01'

      return
800   FORMAT (2X,3D17.8)
810   FORMAT (2X,5D15.7)
      END
      subroutine test11 ( )

c*********************************************************************72
c
cc TEST11 tests DFZERO.
c
      DOUBLE PRECISION     B, C, AE, RE, test11_f
      EXTERNAL test11_f

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST11'
      write ( *, '(a)' ) '  DFZERO solves a nonlinear equation.'
      write ( *, '(a)' ) ' '

      B   = 2.0D0
      C   = 3.0D0
      AE  = 1.0D-06
      RE  = 1.0D-06
      WRITE (*,*) ' INITIAL INTERVAL: '
      WRITE (*,'(8X,2D18.10)') B, C
      WRITE (*,*) ' TOLERANCES:       '
      WRITE (*,'(8X,2D18.10)') AE, RE

      CALL DFZERO ( test11_f, B, C, C, RE, AE, IFLAG)

      WRITE (*,*)
      WRITE (*,*) '  DFZERO RESULTS'
      IF (IFLAG .NE. 1) WRITE (*,*) '  ERROR CODE =', IFLAG
      WRITE (*,'(A19,2X,D18.10)') '  ESTIMATE OF ZERO =', B
      WRITE (*,'(A19,2X,D18.10)') '  FUNCTION VALUE =', test11_f (B)

      WRITE (*,*)
      WRITE (*,*) '  REFERENCE RESULTS FROM IBM PC/AT'
      WRITE (*,*) '   ESTIMATE OF ZERO =   0.2094551435E+01'    
      WRITE (*,*) '     FUNCTION VALUE =  -0.5157851953E-06'

      return
      END
      DOUBLE PRECISION FUNCTION test11_f ( X )

c*********************************************************************72
c
cc TEST11_F evaluates a nonlinear equation.
c
      DOUBLE PRECISION X

      test11_f = X * (X*X-2.0D0) - 5.0D0

      RETURN
      END
      subroutine test12 ( )

c*********************************************************************72
c
cc TEST12 tests UNCMND.
c
C MAIN PROGRAM FOR NONLINEAR LEAST-SQUARES DATA FITTING
C
      PARAMETER       (N = 2, LWORK = N*(N+10), MD = 4)
      DOUBLE PRECISION  X0(N), X(N), F, WORK(LWORK), T(MD), B(MD)
      COMMON /EXPDAT/ T, B, M
      EXTERNAL test12_dcalcf

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST12'
      write ( *, '(a)' ) 
     &  '  UNCMND performs nonlinear least squares data fitting.'
      write ( *, '(a)' ) ' '
C
C DATA FOR DATA FITTING
C
      T(1) =  0.0D0
      T(2) =  1.0D0
      T(3) =  2.0D0
      T(4) =  3.0D0
      B(1) = 20.0D0
      B(2) =  9.0D0
      B(3) =  3.0D0
      B(4) =  1.0D0
C
C SPECIFY INITIAL ESTIMATE OF THE SOLUTION
C
      M     = 4
      X0(1) = 1.0D0
      X0(2) = 1.0D0
C
C MINIMIZE FUNCTION
C
      CALL UNCMND (N, X0, test12_dcalcf, X, F, IERROR, WORK, LWORK)
C
C PRINT RESULTS
C
      WRITE (*,*) 'UNCMND FOR NONLINEAR LEAST SQUARES RESULTS'
      IF (IERROR .NE. 0) WRITE (*,*) ' ERROR CODE =', IERROR
      WRITE (*,'(A,1X,D15.8)') ' F(X*) =', F
      WRITE (*,'(1X,A)') ' X* =' 
      WRITE (*,'(5X,D20.12)') (X(I), I = 1,N)

      WRITE (*,*) 
      WRITE (*,*) 
     * 'REFERENCE RESULTS (PARTIAL-LAST 8 LINES) FROM IBM PC/AT'
      WRITE (*,*) '-0.52716742E+01  -0.20203044E+02  -0.57189483E+01'
      WRITE (*,*) ' 0.20000000E+02  -0.20605341E+02  -0.52716742E+01'
      WRITE (*,*) ' 0.20000001E+02  -0.20605341E+02  -0.52716742E+01'    
      WRITE (*,*) ' 0.20000000E+02  -0.20605340E+02  -0.52716742E+01'
      WRITE (*,*) 'UNCMND WARNING -- INFO = 1: PROBABLY CONVERGED, GRADI
     *ENT SMALL'
      WRITE (*,*) 'UNCMND FOR NONLINEAR LEAST SQUARES RESULTS'
      WRITE (*,*) ' ERROR CODE =           1'
      WRITE (*,*) ' F(X*) =  0.91000000E+02'
      WRITE (*,*) ' X* = '
      WRITE (*,*) '     0.200000001134E+02'
      WRITE (*,*) '    -0.206053408470E+02'    

      return
      END
      SUBROUTINE test12_dcalcf (N, X, F)

c*********************************************************************72
c
cc TEST12_DCALCF evaluates the objective function.
c
      DOUBLE PRECISION  X(N), F, T(4), B(4)
      COMMON /EXPDAT/ T, B, M

      F = 0.0D0
      WRITE (*,'(1X,D15.8,2X,D15.8,2X,D15.8)') X(1),X(2),X(3)
      DO J = 1,M
         F = F + (B(J) - X(1)*EXP(X(2)*T(J)))**2.0D0
      end do

      RETURN
      END
      subroutine test13 ( )

c*********************************************************************72
c
cc TEST13 tests DQK15.
c
C  Compute erf(1), i.e. integral of 2/sqrt(pi) * exp(-x*x) from 0 to 1.0
c
      EXTERNAL test13_func
      DOUBLE PRECISION A,B,RESULT,ABSERR,RESABS,RESASC,PI

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST13'
      write ( *, '(a)' ) '  DQK15 estimates an integral.'
      write ( *, '(a)' ) ' '

      A = 0.0D0
      B = 1.0D0
      PI = 4.0D0*ATAN(1.0D0)
      CALL DQK15(test13_func,A,B,RESULT,ABSERR,RESABS,RESASC)
      WRITE(*,*) ' DQK15 ESTIMATE OF ERF(1)'
      WRITE(*,*) ' 2.0/SQRT(PI)*RESULT'
      WRITE(*,'(3X,D20.12,3X,D20.12)') 2.0D0/SQRT(PI)*RESULT,ABSERR

      WRITE(*,*)
      WRITE(*,*) 'REFERENCE RESULTS COMPUTED ON IBM PC/AT '
      WRITE(*,*) '    0.842700792950E+00     0.828974787422E-14'
      return
      END
      DOUBLE PRECISION FUNCTION test13_func (X)

c*********************************************************************72
c
cc TEST13_FUNC evaluates the integrand.
c
      DOUBLE PRECISION X
      test13_func = EXP(-X*X)
      RETURN
      END
      subroutine test14 ( )

c*********************************************************************72
c
cc TEST14 demonstates DSVDC.
c
C  Compute erf(1), i.e. integral of 2/sqrt(pi) * exp(-x*x) from 0 to 1.0
c
C       CENSUS DATA ANALYSIS WITH THE SVD
C       FROM SECTION 8.1
C     From the book "Numerical Methods and Software"
C          by  D. Kahaner, C. Moler, S. Nash
C               Prentice Hall 1988
C
      INTEGER LDX,N,P,LDU,LDV,JOB,INFO,I,J,JB
      PARAMETER(LDX=8,N=8,P=3,LDU=N,LDV=P,JOB=11)
      DOUBLE PRECISION POP(N),Y(N),X(LDX,P),S(P),E(P),
     *U(LDU,LDU),V(LDV,P),W(N),C(P),YEAR,POP80,
     *TOL,RELERR,SUM,R,RSQ,RI
C
C       C CONTAINS COEFFICIENTS OF POLYNOMIAL C(1)*1+C(2)*T+C(3)*T*T
C         T=YEAR (1900 ETC.)
C
      DATA POP/
     *   75.994575D0,
     *   91.972266D0,
     *  105.710620D0,
     *  122.775046D0,
     *  131.669275D0,
     *  150.697361D0,
     *  179.323175D0,
     *  203.235298D0/

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST14'
      write ( *, '(a)' ) 
     &  '  DSVDC computes the singular value decomposition.'
      write ( *, '(a)' ) ' '

      DO 1 I=1,8
C        Y(I)=(1935.0D0-1900D0*(I-1))/1935.0D0
        Y(I)=1900.0D0+(I-1)*10D0
        X(I,1)=1.0D0
        X(I,2)=Y(I)
        X(I,3)=Y(I)**2
    1 CONTINUE

      CALL DSVDC(X,LDX,N,P,S,E,U,LDU,V,LDV,W,JOB,INFO)
C       NOTE: X IS DESTROYED BY ABOVE CALL!!!
C
C       WRITE SINGULAR VALUES (DESCENDING ORDER)
      WRITE(*,*) 'SINGULAR VALUES ARE: '
      WRITE(*,'(2X,D16.8)') (S(I),I=1,P)
      DO 2 I=1,P
        C(I)=0.0D0
    2 CONTINUE
C
C       RELERR REFLECTS NUMBER OF ACCURATE DIGITS IN DATA
C        E.G. 6 DIGITS ==> RELERR=1.D-6, ETC.
C       MAKING RELERR LARGER INCREASES RESIDUALS
      RELERR=1.D-6
      TOL=RELERR*S(1)
C
C       MULTIPLY U-TRANS * POP, AND SOLVE FOR COEFFICIENTS C(I)
C
      DO 60 J=1,P
         IF(S(J).LE.TOL)GO TO 60
         SUM=0.0D0
         DO I=1,N
            SUM=SUM+U(I,J)*POP(I)
         end do
         SUM=SUM/S(J)
         DO 50 I=1,P
            C(I)=C(I)+SUM*V(I,J)
   50    CONTINUE
   60 CONTINUE
C
      WRITE(*,*) 'COEFFICIENTS (ASSUMING DATA GOOD TO 6 DIGITS) ARE:'
      WRITE(*,'(2X,D16.8)') (C(I),I=1,P)
      WRITE(*,*)
C
C       EVALUATE MODEL (HORNER'S RULE) AND RESIDUALS AT YEAR =1900,...,1980
C
      RSQ=0.0D0
      DO 75 I=1,9
         YEAR=1900.0D0+(I-1)*10.0D0
         POP80=0.0D0
         DO JB=1,P
            J=P+1-JB
            POP80=YEAR*POP80+C(J)
         end do
         IF(I.LT.9) THEN
            RI=POP(I)-POP80
            RSQ=RSQ+RI*RI
            WRITE(*,'(A,I6,A)') 
     *       ' FOR YEAR',IDINT(YEAR),
     *' POP ESTM, MEAS, AND RESIDUAL ARE'
            WRITE(*,'(3D16.8)')POP80,POP(I),RI
         ELSE
            WRITE(*,'(A,I6,A,D16.8)') 
     *        ' FOR YEAR',IDINT(YEAR),' POP ESTMATE IS ',POP80
         ENDIF
   75 CONTINUE

      R=SQRT(RSQ)
      WRITE(*,'(1X,A,D16.8)')
     *'SQUARE ROOT OF RESIDUAL SUM OF SQUARES IS: ',R

      WRITE(*,*)
      WRITE(*,*)'REFERENCE RESULTS (PARTIAL) FROM IBM PC/AT'
      WRITE(*,*)'SINGULAR VALUES ARE: '
      WRITE(*,*)'   0.10594723E+08'
      WRITE(*,*)'   0.64774566E+02'
      WRITE(*,*)'   0.34620247E-03'
      WRITE(*,*)'COEFFICIENTS (ASSUMING DATA GOOD TO 6 DIGITS) ARE:'
      WRITE(*,*)'  -0.16714353E-02'
      WRITE(*,*)'  -0.16169731E+01'
      WRITE(*,*)'   0.87095700E-03'
      WRITE(*,*)'      .  '
      WRITE(*,*)'      .  '
      WRITE(*,*)'      .  '    
      WRITE(*,*)'FOR YEAR  1980 POP ESTMATE IS    0.21289144E+03'
      WRITE(*,*)'SQUARE ROOT OF RESIDUAL SUM OF SQUARES IS:  ',
     *'  0.16096596E+02'
      return
      END
      subroutine test15 ( )

c*********************************************************************72
c
cc TEST15 demonstates DQ1DA.
c 
       DOUBLE PRECISION A, B, EPS, R, E
       double precision f_15
       external f_15

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST15'
      write ( *, '(a)' ) 
     &  '  DQ1DA estimates the integral of a function.'
      write ( *, '(a)' ) ' '

         A = 0.0D0
         B = 1.0D0          
C Set interval endpoints to [0,1]
         EPS = 0.1D-02       
C Set accuracy request for three digits
         CALL DQ1DA (f_15,A,B,EPS,R,E,KF,IFLAG)
         WRITE(*,*)'         INTERVAL: '
         WRITE(*,'(3X,D18.10,3X,D18.10)') A,B
         WRITE(*,*)  
         WRITE(*,*)'         DQ1DA RESULTS: '
         WRITE(*,1) EPS,R,E,KF,IFLAG
    1    FORMAT(3X,D18.10,3X,D18.10,3X,D18.10,3X,I5,2X,I5)
         WRITE(*,*)
         WRITE(*,*)'REFERENCE RESULTS ON IBM PC/AT (CALLED DUNI)'
         WRITE(*,*)'    0.0000000000E+00     0.1000000000E+01',
     *'    0.1000000000E-02     0.4140675161E-01     0.4286574121E-11',
     *'      30      0'
      return
      END
      DOUBLE PRECISION FUNCTION F_15(X)

c*********************************************************************72
c
cc F_15 evaluates the integrand.
c
      double precision x

      f_15 = dsin(2.0d0*x)-dsqrt(x)

      return 
      end
      subroutine test16 ( )

c*********************************************************************72
c
cc TEST16 demonstates DEZFTF.
c 
C Using the D.P. discrete Fourier transform, find the approximate 
C Fourier coefficients to Runge's function on [-1,1] with N=16 and 
C N=17. 
C 
       PARAMETER (MCOEF=17)
       DOUBLE PRECISION  A(MCOEF/2),B(MCOEF/2),R(MCOEF),
     *WSAVE(3*MCOEF+15),TN,ER,PI,DEL,XJ,DFTA(MCOEF/2),
     *DFTB(MCOEF/2),F,C(MCOEF/2),S(MCOEF/2),X,RUNGE,X0,AZERO 
C
C Arithmetic statement function for Runge's function. 
       RUNGE(X) = 1.0D0/(1.0D0+25.0D0*X*X) 

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST16'
      write ( *, '(a)' ) 
     &  '  DEZFTF computes the discrete Fourier transform.'
      write ( *, '(a)' ) ' '

       X0 = -1.0D0 
       PI = ASIN(1.0D0)*2.0D0
C
       DO 10 N = MCOEF-1,MCOEF 
          CALL DEZFTI (N,WSAVE) 
C Function assumed to be periodic on [-1,1], of length 2.
          DEL = 2.0D0/N 
          F = 2.0D0*PI/(N*DEL) 
          DO 1 J = 1,N 
C First sample point at -1, last at 1-DEL
              XJ = (-1.0D0) + (J-1)*DEL 
              R(J) = RUNGE(XJ) 
C Compute sines and cosines to adjust output of DEZFTF to give
C approximate Fourier coefficients.
              IF (J .LE. N/2) THEN 
                  C(J) = COS(J*F*X0) 
                  S(J) = SIN(J*F*X0) 
              END IF 
    1     CONTINUE 
          CALL DEZFTF (N,R,AZERO,A,B,WSAVE) 
C 
C As a convenience this loop can go to N/2. If N is even last B is zero.
c
          DO J = 1,N/2 
              DFTA(J) = A(J)*C(J) - B(J)*S(J) 
              DFTB(J) = A(J)*S(J) + B(J)*C(J)  
          end do

          WRITE (*,'(A,I3,A,D18.10)') ' DEZFTF RESULTS FOR N= ' ,
     *N, ' AZERO = ',AZERO
          WRITE (*,*) '    J          DFTA(J)              DFTB(J) '
          DO 12 J = 1,N/2
              WRITE(*,'(1X,I5,3X,D18.10,3X,D18.10)') J,
     *DFTA(J),DFTB(J)
   12     CONTINUE

          M = 101 
C
          WRITE (*,*) 'FOR BREVITY 101 EVALUATION POINTS OMITTED'
          IF (M.GT.0)GOTO 10
C Evaluate interpolant at 101 points on [-1,1]     
          WRITE (*,*) ' RESULTS FOR N= ',N 
          DO 20 K = 1,M 
              X = -1.0D0 + 2.0D0*(K-1.0D0)/(M-1.0D0) 
              TN = AZERO 
              DO 19 J = 1,N/2 
                  TN = TN + DFTA(J)*COS(J*F*X) + DFTB(J)*SIN(J*F*X) 
   19         CONTINUE 
              ER = TN - RUNGE(X) 
              WRITE (*,'(2X,D15.8,2X,D15.8,2X,D15.8)') X,TN,ER 
   20     CONTINUE 

          WRITE (*,*)

   10  CONTINUE 
       WRITE (*,*)
       WRITE (*,*) 'REFERENCE RESULTS FROM IBM PC/AT'
       WRITE (*,*)
     *  ' DEZFTF RESULTS FOR N=  17    AZERO =    0.274581'    
       WRITE (*,*)'    J         DFTA(J)              DFTB(J) '
       WRITE (*,*)'    1     0.3442428079E+00     0.5030951647E-16'
       WRITE (*,*)'    2     0.1755197395E+00     0.4668315002E-16'
       WRITE (*,*)'    3     0.9699027603E-01     0.3740740688E-16'
       WRITE (*,*)'    4     0.4964416088E-01     0.2403131256E-16'
       WRITE (*,*)'    5     0.2759681609E-01     0.8202373493E-17'
       WRITE (*,*)'    6     0.1323313221E-01     0.4148226049E-17'
       WRITE (*,*)'    7     0.7099464912E-02     0.1053997863E-16'
       WRITE (*,*)'    8     0.1413251200E-02     0.5512226224E-17'

      return
      END
      subroutine test17 ( )

c*********************************************************************72
c
cc TEST17 demonstates DDRIV2.
C
C Notice that NROOT is set to 1
C
      DOUBLE PRECISION  H
      PARAMETER (N=2, H=10.0, NROOT=1, MINT=2, 
     *           LW=N*N+10*N+2*NROOT+204, LIW=23)
      DOUBLE PRECISION  Y(N+1),W(LW),GFUN,EPS,T,TOUT,EWT,MASS
      INTEGER   IW(LIW)
      EXTERNAL  test17_FSUB,test17_GFUN
      DATA      MASS /0.125D0/

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST17'
      write ( *, '(a)' ) 
     &  '  DDRIV2 solves a differential equation.'
      write ( *, '(a)' ) ' '

      EPS = 1.D-5
c
C Set initial point
c
      T = 0.0D0
      TOUT = T
C Set for pure relative error
      EWT = 0.0D0
C Set initial conditions
      Y(1) = H
      Y(2) = 0.0D0
C Set parameter value
      Y(3) = MASS
      MS = 1
      WRITE (*,*) 'DDRIV2 RESULTS'
      WRITE (*,5) T,Y(1),Y(2),MS
    5 FORMAT(2X,D18.10,2X,D18.10,2X,D18.10,1X,I3)
C
   10 CALL DDRIV2 (N,T,Y,test17_FSUB,TOUT,MS,NROOT,EPS,EWT,MINT,W,LW,
     &  IW,LIW,test17_GFUN)

      TOUT = TOUT+0.1D0
      IF (MS .EQ. 5) THEN
        WRITE (*,'(2X,D18.10,2X,D18.10,2X,D18.10,I4)') 
     *        T,Y(1),Y(2),MS
        WRITE (*,'(2X,A,D18.10)') ' <-- Y=0 AT T= ',T
        GOTO 100
      ELSE
        WRITE (*,'(2X,D18.10,2X,D18.10,2X,D18.10,I4)') T,Y(1),Y(2),MS
C Stop if any output code but 1 or 2.
        IF (MS .GT. 2) GOTO 100
      END IF
      GOTO 10
  100 WRITE (*,*)
      WRITE (*,*) 'REFERENCE RESULTS (LAST LINE) FROM IBM PC/AT'
      WRITE (*,*)
     *'   0.2624999995E+01   -0.5648151218E-15   -0.3999999828E+01   5', 
     *'  <-- Y=0 AT T=   0.2624999995E+01'
      return
      END 
      SUBROUTINE test17_FSUB (N,T,Y,YDOT)

c*********************************************************************72
c
cc TEST17_FSUB evaluates right hand sides of equations.
c
        DOUBLE PRECISION T, Y(*), YDOT(*), MASS, G
        DATA G/32.0D0/
c
C Retrieve parameter
c
        MASS = Y(3)

        YDOT(1) = Y(2)
        YDOT(2) = -(G+1.0D0/MASS*Y(2))
        RETURN
      END
      DOUBLE PRECISION FUNCTION test17_GFUN(N,T,Y,IROOT)

c*********************************************************************72
c
cc TEST17_GFUN is used for root finding.
c
C Integration will stop when the function changes sign.
c
      DOUBLE PRECISION Y(*),T
      test17_GFUN = Y(1)
      RETURN
      END
      subroutine test18 ( )

c*********************************************************************72
c
cc TEST18 finds autocorrelation to el Nino data using direct and FFT methods.
C
      INTEGER N
      PARAMETER(N=168)
      DOUBLE PRECISION EL(0:2*N-1),WSAVE(4*(2*N)+15),ACOV(0:N-1),
     *           A(2*N),B(2*N),ACOVR(0:2*N-1),SUM,AZERO
      COMPLEX*16 CEL(0:2*N-1),CORR(0:2*N-1)
      LOGICAL EX

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST18'
      write ( *, '(a)' ) 
     &  '  Find autocorrelation in el Nino data.'
      write ( *, '(a)' ) ' '
C
C  NOTE:  DOUBLE PRECISION COMPLEX IS NOT STANDARD FORTRAN.
C         SOME COMPILER VENDORS MAY FAIL TO IMPLEMENT IT
C         OR MAY IMPLEMENT IN A DIFFERENT MANNER.  IF YOUR
C         COMPILER DOES NOT ACCEPT THE COMPLEX*16 STATEMENT
C         WE RECOMMEND THAT YOU CHECK WITH YOUR VENDOR AND
C         MODIFY THE ABOVE ACCORDINGLY.
C
C needs N locations for el, 2N for acov, cel, corr
C and 4*(2N)+15 for wsave in complex case. N has maximum of 168.
c
      INQUIRE (FILE='elnino.dat',EXIST=EX,ERR=1000)
      IF(.NOT.EX)GOTO 1000
      OPEN (UNIT=8,FILE='elnino.dat',ERR=1000)
C
C Read data, find mean.
c
      SUM = 0.0D0
      DO I = 0,N-1
         READ (8,*) EL(I)
         SUM = SUM + EL(I)
      end do

      DO 101 I = 0,N-1
         EL(I) = EL(I) - SUM/N
C Subtract mean, and add N zeros for either complex or d.p. FFT usage
         CEL(I) = CMPLX(EL(I),0.0D0)
         EL(I+N) = 0.0D0
         CEL(I+N) = 0.0D0
  101 CONTINUE
C
C Direct calculation. Only sum as far as there is data. 
C Simple, but slow.
C      
      DO 110 J = 0,N-1
         ACOV(J) = 0.0D0
         DO 110 M = 0,N-1-J
            ACOV(J) = ACOV(J)+ EL(M) * EL(M+J)
  110 CONTINUE
C
C Write, scaled correlation
      WRITE (*,*) 
     * 'EX 11.6: AUTOCORRELATION (DIRECT) OUTPUT SUPRESSED'
CCCCC      WRITE (*,'(5D14.6)') ( (ACOV(I)/ACOV(0)), I = 0,N-1)
Cc
C
C FFT approach (complex).
C    Compute FFT of data of length 2N.
C    Compute square of magnitude of transform components and place  
C       in complex array as d.p. parts.
C    Compute inverse transform, throwing away second half and
C       imaginary parts (which are zero), and multiply by length of 
C       sequence, 2N.  
c
      CALL DCFFTI (2*N,WSAVE)
      CALL DCFFTF (2*N,CEL,WSAVE)
c
C DCFFTF returns unscaled transforms. Actual transforms are output
C divided by (2N).
c
      DO I = 0,2*N-1
         CORR(I) = ABS(CEL(I) / (2*N)) **2 
      end do
c
C Since we compute transform times its conjugate, must divide by
C (2N) for each, i.e., (2N)**2.
      CALL DCFFTB (2*N,CORR,WSAVE)

      DO I = 0,N-1
         ACOV(I) = DBLE(CORR(I))*(2*N)
      end do
c
C Autocovariance is inverse transform times sequence length, 2N.
C Normally, all the scaling would be done  only once
C by dividing by 2N. We've broken it up for exposition.
c
      WRITE (*,*) 
     *  'EX 11.6: AUTOCORRELATION (COMPLEX FFT) OUTPUT SUPRESSED'
CCCCC      WRITE (*,'(5D14.6)') ( (ACOV(I)/ACOV(0) ), I = 0,N-1)
C
C FFT approach (d.p.).
C    Compute FFT of data of length 2N.
C    DEZFTF produces correctly scaled A's and B's so no extra scaling
C       needed to get transform.     
C    Compute array of square of each frequency component and place
C       in cosine array (A's) to be back transformed. Set B's to 0.
C    There are N A's, and N B's.
C    Note that care must be taken to compute magnitude correctly, 
C       0.5*(A(I)**2+B(I)**2) for I < N, twice that for I=N.
C    Compute back transform throwing away its second half.
C
      CALL DEZFTI (2*N,WSAVE)
      CALL DEZFTF (2*N,EL,AZERO,A,B,WSAVE)
      AZERO = AZERO*AZERO
C
      DO 150 I = 1,N
         IF(I.NE.N) THEN
            A(I) = (A(I)**2 + B(I)**2) / 2.0D0
         ELSE
            A(I) = (A(I)**2 + B(I)**2)
         ENDIF
         B(I) = 0.0D0
  150 CONTINUE
      CALL DEZFTB (2*N,ACOVR,AZERO,A,B,WSAVE)
      WRITE (*,*) 
     * 'EX 11.6: AUTOCORRELATION (DBLE FFT) OUTPUT REDUCED'
      WRITE (*,'(1X,3D18.10)') ( (ACOVR(I)/ACOVR(0) ), I = 0,19)
CCCCC      WRITE (*,'(5D14.6)') ( (ACOVR(I)/ACOVR(0) ), I = 0,N-1)
      WRITE (*,*)
      WRITE (*,*) 
     * 'REFERENCE RESULTS (EX 11.6 PARTIAL-LAST 7 LINES) FROM IBM PC/AT'
      WRITE (*,*)    
     *'  0.1000000000E+01  0.6067398062E+00  0.3538571074E+00'
      WRITE (*,*)
     *'  0.1843027185E+00 -0.1529114172E-01 -0.2198717698E+00'
      WRITE (*,*)
     *' -0.2962181307E+00 -0.2914580887E+00 -0.1550562043E+00'
      WRITE (*,*)
     *'  0.3680803130E-01  0.1735588664E+00  0.2665974779E+00'
      WRITE (*,*)
     *'  0.3049896546E+00  0.2011691615E+00  0.1780362880E-01'
      WRITE (*,*)
     *' -0.2103672514E+00 -0.3773868663E+00 -0.4379599788E+00'
      WRITE (*,*)
     *' -0.4603327597E+00 -0.4255892823E+00'

      return
 1000 WRITE (*,*) 'CANNOT FIND THE DATA FILE: elnino.dat '
      END
      subroutine test19 ( )

c*********************************************************************72
c
cc TEST19 minimizes a function using UNCMND.
c
      PARAMETER (N = 10, LWORK = N*(N+10))
      DOUBLE PRECISION X0(N), X(N), F, WORK(LWORK), S, W
      EXTERNAL test19_calcf

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST19'
      write ( *, '(a)' ) 
     &  '  UNCMND finds the minimum of a function.'
      write ( *, '(a)' ) ' '
C
C SPECIFY INITIAL ESTIMATE OF THE SOLUTION
C
      WRITE (*,*) 'COMPUTING...'
      DO I = 1,N
         X0(I) = I / FLOAT(N+1)
      end do
C
C MINIMIZE FUNCTION
C
      CALL UNCMND (N, X0, test19_calcf, X, F, IERROR, WORK, LWORK)
C
C PRINT RESULTS
C
      WRITE (*,*) 'UNCMND RESULTS'
      IF (IERROR .NE. 0) WRITE (*,*) ' ERROR CODE =', IERROR
      WRITE (*,'(1X,A,1X,D20.12)') ' F(X*) = ', F
      WRITE (*,*) ' X* ='
      WRITE (*,800) (X(I), I = 1,N)

      WRITE (*,*)
      WRITE (*,*) 'REFERENCE RESULTS FROM IBM PC/AT'
      WRITE (*,*) 'UNCMND WARNING -- INFO = 1: PROBABLY CONVERGED, GRADI
     *ENT SMALL'
      WRITE (*,*) 'UNCMND RESULTS'
      WRITE (*,*) ' ERROR CODE =           1'
      WRITE (*,*) ' F(X*) =    0.100000000000E+01'
      S=0.1D+1
      W=0.9999999D0
      WRITE (*,*) ' X* ='
      WRITE (*,799) S,S,S,S,S
      WRITE (*,799) S,S,S,S,W

      return
799   FORMAT (1X,D14.7,2X,D14.7,2X,D14.7,2X,D14.7,2X,D14.7)
800   FORMAT (1X,D14.7,2X,D14.7,2X,D14.7,2X,D14.7,2X,D14.7)
      END
      SUBROUTINE test19_calcf (N, X, F)

c*********************************************************************72
c
cc TEST19_CALCF evaluates the function to be minimized.
c
      DOUBLE PRECISION X(N), F, T1, T2

      T1 = 0.0
      T2 = 0.0
      DO I = 2,N
         T1 = T1 + (X(I)-X(I-1)**2)**2
         T2 = T2 + (1.0-X(I-1))**2
      end do
      F = 1.0 + 100.0*T1 + T2

      RETURN
      END
      subroutine test20 ( )

c*********************************************************************72
c
cc TEST20 tests DQAGI.
c
C          COMPUTE INTEGRAL OF EXP(-X)*COS(X*X)**2 ON [0,INFINITY)
C               RESULT IS 0.70260...
C
      INTEGER LIMIT,LENW
      PARAMETER(LIMIT=100,LENW=LIMIT*4)
      DOUBLE PRECISION BOUND,EPSABS,EPSREL,RESULT,ABSERR,WORK(LENW)
      INTEGER INF,NEVAL,IER,LAST,IWORK(LIMIT)
      EXTERNAL test20_f

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST20'
      write ( *, '(a)' ) 
     &  '  DQAGI estimates an integral.'
      write ( *, '(a)' ) ' '

      BOUND = 0.0D0
      INF = 1
      EPSABS = 0.0D0
      EPSREL = 1.D-5
      CALL DQAGI(test20_f,BOUND,INF,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,
     *    IER,LIMIT,LENW,LAST,IWORK,WORK)
      WRITE(*,*)'DQAGI RESULT, ABSERR, NEVAL, IER: '
      WRITE(*,'(2X,2D20.10,I6,I5)') RESULT, ABSERR, NEVAL, IER

      WRITE(*,*)
      WRITE(*,*)'REFERENCE RESULTS FROM IBM PC/AT'
      WRITE(*,*)'   BE SURE THAT UNDERFLOWS ARE SET TO ZERO...'
      WRITE(*,*) 'DQAGI RESULT, ABSERR, NEVAL, IER: '
      WRITE(*,*)'     0.7026028726E+00    0.6401518223E-05  1005    0'

      return
      END
      DOUBLE PRECISION FUNCTION test20_f(X)

c*********************************************************************72
c
cc TEST20_F evaluates the integrand.
c
      double precision x,exp,cos
      test20_f = exp(-x)*cos(x*x)**2
      return
      end
      subroutine test21 ( )

c*********************************************************************72
c
cc TEST21 runs the REACTOR SHIELDING PROBLEM
C
      DOUBLE PRECISION MU,E,AZM,D,DIST2C,X,Y,Z,EA,ER,ET,SA,SR,ST,THICK
      INTEGER particle_exit,I,NA,NT,NR,NTOT,NPART
      LOGICAL ABSORB
      COMMON THICK

      ntot = 1000
      thick = 2.0D+00

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST21'
      write ( *, '(a)' ) '  Reactor shielding problem'
      write ( *, '(a,i8)' ) '  Number of particles = ', ntot
      write ( *, '(a,g14.6)' ) '  Slab thickness = ', thick
c
C  INITIALIZE PARTICLE COUNTS AND ENERGY TALLIES
c
      CALL INIT(NA,NT,NR,NPART,EA,ET,ER,SA,SR,ST)
C
C   MAIN LOOP OVER NTOT PARTICLES.
c
    1 CONTINUE

      NPART=NPART+1
c
C FINISHED, COMPUTE AND PRINT AVERAGES, STANDARD DEVIATIONS
c
      IF(NPART.EQ.NTOT)THEN
        CALL OUTPUT(NA,EA,SA,NR,ER,SR,NT,ET,ST,NTOT)
        GOTO 3
      ENDIF
C
C  SOURCE GENERATES A NEW PARTICLE WITH POSITION, DIRECTION, ENERGY
c
      CALL SOURCE(E,MU,AZM,X,Y,Z)

    2 CONTINUE
C
C  COMPUTE-DISTANCE-TO-COLLISION
c
      D=DIST2C(E)
C
C  UPDATE POSITION TO DECIDE IF PARTICLE EXITS OR COLLIDES
c
      CALL UPDATE(MU,AZM,D,X,Y,Z)

      I=particle_exit(X,Y,Z)
c
C  RETURNS -1, 0 , 1 FOR 'OUT ON LEFT', 'IN', 'OUT ON RIGHT'
c
C  EXIT ON LEFT (REFLECTED), TALLY AND GOTO SOURCE
c
      IF(I.LT.0)THEN

           NR=NR+1
           ER=ER+E  
           SR=SR+E*E
           GOTO 1
c
C  EXIT ON RIGHT (TRANSMITTED), TALLY AND GOTO SOURCE
c
      ELSEIF(I.GT.0)THEN

           NT=NT+1
           ET=ET+E
           ST=ST+E*E
           GOTO 1
c
C  COLLISION & ABSORBED, TALLY AND GOTO SOURCE
c
      ELSEIF(I.EQ.0 .AND. ABSORB())THEN

           NA=NA+1
           EA=EA+E
           SA=SA+E*E
           GOTO 1
c
C  COLLISION & SCATTERED, FIND SCATTERING ANGLE AND ENERGY
c
      ELSE
        CALL SCATT(MU,AZM,E)
        GOTO 2
      ENDIF

3     continue

      RETURN
      END
      SUBROUTINE INIT(NA,NT,NR,NPART,EA,ET,ER,SA,SR,ST)

c*********************************************************************72
c
cc INIT initializes the data.
c
      DOUBLE PRECISION EA,ER,ET,SA,SR,ST
      INTEGER NA,NT,NR,NPART
      NA=0
      NT=0
      NR=0
      NPART=0
      EA=0.0D0
      ET=0.0D0
      ER=0.0D0
      SA=0.0D0
      SR=0.0D0
      ST=0.0D0
      RETURN
      END
      SUBROUTINE SOURCE(E,MU,AZM,X,Y,Z)

c*********************************************************************72
c
cc SOURCE models the emitter.
c
C          EMITTER  
C          RETURNS MU UNIFORM ON [0,1] AS THE COSINE OF AN ANGLE, 
C              AND AZM UNIFORM ON [0,2*PI] AS ASIMUTHAL ANGLE 
c
      DOUBLE PRECISION E,MU,AZM,X,Y,Z,PI,ENERGY,DUNI

      PI=4*ATAN(1.0D0)
      MU=DUNI()
      AZM=DUNI()*(2.0D0*PI)
      X=0.0D0
      Y=0.0D0
      Z=0.0D0
c
C  RETURNS ENERGY FROM SOURCE ENERGY DISTRIBUTION
c
      E=ENERGY()
      RETURN
      END
      SUBROUTINE UPDATE(MU,AZM,D,X,Y,Z)

c*********************************************************************72
c
cc UPDATE
c
      DOUBLE PRECISION MU,AZM,D,X,Y,Z,R     
 
      X=X+D*MU
      R=D*SQRT(1.0D0-MU*MU)
      Y=Y+R*COS(AZM)
      Z=Z+R*SIN(AZM)
      RETURN
      END
      SUBROUTINE SCATT(MU,AZM,E)

c*********************************************************************72
c
cc SCATT returns scattering angle and energy.
c
      DOUBLE PRECISION MU,AZM,E,PI,DUNI
C
C  ISOTROPIC SCATTERING, I.E., UNIFORM ON SPHERE 
C
      PI=4.0D0*ATAN(1.0D0)
      MU=-1.0D0+2.0D0*DUNI()
      AZM=DUNI()*2D0*PI
C
C  FIND SCATTERING ENERGY, UNIFORM IN [.3D0*E,E]
C
      E=(DUNI()*0.7D0+0.3D0)*E
      RETURN
      END
      LOGICAL FUNCTION ABSORB()

c*********************************************************************72
c
cc ABSORB returns true if particle is absorbed. 
c
c  OCCURS WITH PROB PA.
c
      DOUBLE PRECISION PA,DUNI
      PARAMETER(PA=0.1D0)

      IF(DUNI().LE.PA)THEN
        ABSORB=.TRUE.
      ELSE
        ABSORB=.FALSE.
      ENDIF

      RETURN
      END
      INTEGER FUNCTION particle_exit (X,Y,Z)

c*********************************************************************72
c
cc PARTICLE_EXIT reports the location of the particle.
c
c  RETURNS -1, 0, +1 AS PARTICLE IS   OUTSIDE ON LEFT, INSIDE, 
C                            OR OUTSIDE ON RIGHT.
c
      DOUBLE PRECISION X,Y,Z,THICK
      COMMON THICK

      IF(X.GT.THICK)THEN
        particle_exit=1
      ELSEIF(X.LT.0.0D0)THEN
        particle_exit=-1
      ELSE
        particle_exit=0
      ENDIF

      RETURN
      END
      DOUBLE PRECISION FUNCTION DIST2C(E)

c*********************************************************************72
c
cc DIST2C returns distance to collision.
c
c  It uses an exponential distribution with parameter  `cross section'
c
      DOUBLE PRECISION LOG,CROSS,E,DUNI

      DIST2C=-LOG(DUNI())/CROSS(E)

      RETURN
      END
      DOUBLE PRECISION FUNCTION ENERGY()

c*********************************************************************72
c
cc ENERGY returns energy, with distribution const/sqrt(energy) over [EMIN,EMAX].
c
C               USE INVERSE FUNCTION APPROACH TO COMPUTE THIS
C
C      ENERGY MIN, MAX IN MEV
c
      DOUBLE PRECISION C,SQRT,EMIN,EMAX,DUNI
      PARAMETER(EMIN=1.0D-3,EMAX=2.5D0)

      C=1.0D0/(2.0D0*(SQRT(EMAX)-SQRT(EMIN)))
      ENERGY=( DUNI()/(2*C)+SQRT(EMIN) )**2
      RETURN
      END
      DOUBLE PRECISION FUNCTION CROSS(E)

c*********************************************************************72
c
cc CROSS returns cross section (fictional) for energy in range [EMIN,EMAX]
c
      DOUBLE PRECISION E,S,ABS,SIN,Y,EXP

      S=ABS(SIN(100.0D0*(EXP(E)-1.0D0))+SIN(18.81D0*(EXP(E)-1.0D0)))
      Y=MAX(0.02D0,S)
      CROSS=10.0D0*EXP(-0.1D0/Y)
      RETURN
      END
      SUBROUTINE OUTPUT(NA,EA,SA,NR,ER,SR,NT,ET,ST,NTOT)

c*********************************************************************72
c
cc OUTPUT prints information.
c
      DOUBLE PRECISION EA,SA,ER,SR,ET,ST
      INTEGER NA,NR,NT,NTOT

      WRITE(*,*) 'TALLIES '

      IF(NA.GT.0)THEN
           EA=EA/NA
           SA=SQRT(SA/NA-EA*EA)
      ENDIF

      WRITE(*,'(1X,A,D10.4,2D20.10)') '% ABSORBED, ENERGY, SD: ',
     *       FLOAT(NA)/NTOT*100,EA,SA

      IF(NR.GT.0)THEN
           ER=ER/NR
           SR=SQRT(SR/NR-ER*ER)
      ENDIF

      WRITE(*,'(1X,A,D10.4,2D20.10)') '% REFLECTED, ENERGY, SD: ',
     *       FLOAT(NR)/NTOT*100,ER,SR

      IF(NT.GT.0)THEN
        ET=ET/NT
        ST=SQRT(ST/NT-ET*ET)
      ENDIF

      WRITE(*,'(1X,A,D10.4,2D20.10)') '% TRANSMITTED, ENERGY, SD: ',
     *       FLOAT(NT)/NTOT*100,ET,ST
      WRITE(*,*)

      RETURN
      END