C IMSL ROUTINE NAME - DCADRE
C-----------------------------------------------------------------------
C
C   COMPUTER            - VAX/DOUBLE
C
C   LATEST REVISION     - JUNE 1, 1982
C
C   PURPOSE             - NUMERICAL INTEGRATION OF A FUNCTION USING
C                           CAUTIOUS ADAPTIVE ROMBERG EXTRAPOLATION
C
C   USAGE               - FUNCTION DCADRE (F,A,B,AERR,RERR,ERROR,IER)
C
C   ARGUMENTS    DCADRE - ESTIMATE OF THE INTEGRAL OF F(X) FROM A TO B.
C                           (OUTPUT).
C                F      - A SINGLE-ARGUMENT REAL FUNCTION SUBPROGRAM
C                           SUPPLIED BY THE USER. (INPUT)
C                           F MUST BE DECLARED EXTERNAL IN THE
C                           CALLING PROGRAM.
C                A,B    - THE TWO ENDPOINTS OF THE INTERVAL OF
C                           INTEGRATION. (INPUT)
C                AERR   - DESIRED ABSOLUTE ERROR IN THE ANSWER. (INPUT)
C                RERR   - DESIRED RELATIVE ERROR IN THE ANSWER. (INPUT)
C                ERROR  - ESTIMATED BOUND ON THE ABSOLUTE ERROR OF
C                           THE OUTPUT NUMBER, DCADRE. (OUTPUT)
C                IER    - ERROR PARAMETER. (OUTPUT)
C                         WARNING ERROR(WITH FIX)
C                           IER = 65 IMPLIES THAT ONE OR MORE
C                             SINGULARITIES WERE SUCCESSFULLY HANDLED.
C                           IER = 66 IMPLIES THAT, IN SOME
C                             SUBINTERVAL(S), THE ESTIMATE OF THE
C                             INTEGRAL WAS ACCEPTED MERELY BECAUSE THE
C                             ESTIMATED ERROR WAS SMALL, EVEN THOUGH NO
C                             REGULAR BEHAVIOR WAS RECOGNIZED.
C                         TERMINAL ERROR
C                           IER = 131 INDICATES FAILURE DUE TO
C                             INSUFFICIENT INTERNAL WORKING STORAGE.
C                           IER = 132 INDICATES FAILURE DUE TO
C                             TOO MUCH NOISE IN THE FUNCTION (RELATIVE
C                             TO THE GIVEN ERROR REQUIREMENTS) OR
C                             DUE TO AN ILL-BEHAVED INTEGRAND.
C                           IER = 133 INDICATES THAT RERR IS GREATER
C                             THAN 0.1, OR RERR IS LESS THAN 0.0, OR
C                             RERR IS TOO SMALL FOR THE PRECISION OF
C                             THE MACHINE.
C-----------------------------------------------------------------------
C
     DOUBLE PRECISION FUNCTION DCADRE (F,A,B,AERR,RERR,ERROR,IER)
C                                  SPECIFICATIONS FOR ARGUMENTS
     INTEGER            IER
     DOUBLE PRECISION   F,A,B,AERR,RERR,ERROR

C IMSL ROUTINE NAME - ICSCCU
C-----------------------------------------------------------------------
C
C   PURPOSE             - CUBIC SPLINE INTERPOLATION
C                           (EASY-TO-USE VERSION)
C
C   USAGE               - CALL ICSCCU (X,Y,NX,C,IC,IER)
C
C   ARGUMENTS    X      - VECTOR OF LENGTH NX CONTAINING THE ABSCISSAE
C                           OF THE NX DATA POINTS (X(I),Y(I)) I=1,...,
C                           NX. (INPUT) X MUST BE ORDERED SO THAT
C                           X(I) .LT. X(I+1).
C                Y      - VECTOR OF LENGTH NX CONTAINING THE ORDINATES
C                           (OR FUNCTION VALUES) OF THE NX DATA POINTS.
C                           (INPUT)
C                NX     - NUMBER OF ELEMENTS IN X AND Y. (INPUT) NX
C                           MUST BE .GE. 2.
C                C      - SPLINE COEFFICIENTS. (OUTPUT) C IS AN NX-1 BY
C                           3 MATRIX. THE VALUE OF THE SPLINE
C                           APPROXIMATION AT T IS
C                           S(T) = ((C(I,3)*D+C(I,2))*D+C(I,1))*D+Y(I)
C                           WHERE X(I) .LE. T .LT. X(I+1) AND
C                           D = T-X(I).
C                IC     - ROW DIMENSION OF MATRIX C EXACTLY AS
C                           SPECIFIED IN THE DIMENSION STATEMENT IN
C                           THE CALLING PROGRAM. (INPUT)
C                IER    - ERROR PARAMETER. (OUTPUT)
C                         TERMINAL ERROR
C                           IER = 129, IC IS LESS THAN NX-1.
C                           IER = 130, NX IS LESS THAN 2.
C                           IER = 131, INPUT ABSCISSA ARE NOT ORDERED
C                             SO THAT X(1) .LT. X(2) ... .LT. X(NX).
C-----------------------------------------------------------------------
C
     SUBROUTINE ICSCCU (X,Y,NX,C,IC,IER)
C                                  SPECIFICATIONS FOR ARGUMENTS
     INTEGER            NX,IC,IER
     DOUBLE PRECISION   X(NX),Y(NX),C(IC,3)

C IMSL ROUTINE NAME - ICSEVU
C-----------------------------------------------------------------------
C
C   PURPOSE             - EVALUATION OF A CUBIC SPLINE
C
C   USAGE               - CALL ICSEVU(X,Y,NX,C,IC,U,S,M,IER)
C
C   ARGUMENTS    X      - VECTOR OF LENGTH NX CONTAINING THE ABSCISSAE
C                           OF THE NX DATA POINTS (X(I),Y(I)) I=1,...,
C                           NX (INPUT). X MUST BE ORDERED SO THAT
C                           X(I) .LT. X(I+1).
C                Y      - VECTOR OF LENGTH NX CONTAINING THE ORDINATES
C                           (OR FUNCTION VALUES) OF THE NX DATA POINTS
C                           (INPUT).
C                NX     - NUMBER OF ELEMENTS IN X AND Y (INPUT).
C                           NX MUST BE .GE. 2.
C                C      - SPLINE COEFFICIENTS (INPUT). C IS AN NX-1 BY
C                           3 MATRIX.
C                IC     - ROW DIMENSION OF MATRIX C EXACTLY AS
C                           SPECIFIED IN THE DIMENSION STATEMENT
C                           IN THE CALLING PROGRAM (INPUT).
C                           IC MUST BE .GE. NX-1.
C                U      - VECTOR OF LENGTH M CONTAINING THE ABSCISSAE
C                           OF THE M POINTS AT WHICH THE CUBIC SPLINE
C                           IS TO BE EVALUATED (INPUT).
C                S      - VECTOR OF LENGTH M (OUTPUT).
C                           THE VALUE OF THE SPLINE APPROXIMATION AT
C                           U(I) IS
C                           S(I) = ((C(J,3)*D+C(J,2))*D+C(J,1))*D+Y(J)
C                           WHERE X(J) .LE. U(I) .LT. X(J+1) AND
C                           D = U(I)-X(J).
C                M      - NUMBER OF ELEMENTS IN U AND S (INPUT).
C                IER    - ERROR PARAMETER (OUTPUT).
C                         WARNING ERROR
C                           IER = 33, U(I) IS LESS THAN X(1).
C                           IER = 34, U(I) IS GREATER THAN X(NX).

C IMSL ROUTINE NAME - MDCH
C-----------------------------------------------------------------------
C
C   PURPOSE             - CHI-SQUARED PROBABILITY DISTRIBUTION FUNCTION
C
C   USAGE               - CALL MDCH (CS,DF,P,IER)
C
C   ARGUMENTS    CS     - INPUT VALUE FOR WHICH THE PROBABILITY IS
C                           COMPUTED. CS MUST BE GREATER THAN OR EQUAL
C                           TO ZERO.
C                DF     - INPUT NUMBER OF DEGREES OF FREEDOM OF THE
C                           CHI-SQUARED DISTRIBUTION. DF MUST BE GREATER
C                           THAN OR EQUAL TO .5 AND LESS THAN OR EQUAL
C                           TO 200,000.
C                P      - OUTPUT PROBABILITY THAT A RANDOM VARIABLE
C                           WHICH FOLLOWS THE CHI-SQUARED DISTRIBUTION
C                           WITH DF DEGREES OF FREEDOM IS LESS THAN OR
C                           EQUAL TO CS.
C                IER    - ERROR PARAMETER. (OUTPUT)
C                         TERMINAL ERROR
C                           IER = 129 INDICATES THAT CS OR DF WAS
C                             SPECIFIED INCORRECTLY.
C                         WARNING ERROR
C                           IER = 34 INDICATES THAT THE NORMAL PDF
C                             WOULD HAVE PRODUCED AN UNDERFLOW.
C-----------------------------------------------------------------------
C
     SUBROUTINE MDCH (CS,DF,P,IER)
C                                  SPECIFICATIONS FOR ARGUMENTS
     REAL               CS,DF,P

C IMSL ROUTINE NAME - DREBS
C-----------------------------------------------------------------------
C
C   PURPOSE             - DIFFERENTIAL EQUATION SOLVER -
C                           EXTRAPOLATION METHOD
C
C   USAGE               - CALL DREBS (FCN,Y,X,N,JM,IND,JSTART,H,HMIN,
C                           TOL,R,S,WK,IER)
C
C   ARGUMENTS    FCN    - NAME OF SUBROUTINE FOR EVALUATING FUNCTIONS.
C                           (INPUT)
C                           THE SUBROUTINE ITSELF MUST ALSO BE PROVIDED
C                             BY THE USER AND IT SHOULD BE OF THE
C                             FOLLOWING FORM
C                               SUBROUTINE FCN(N,X,Y,YPRIME)
C                               REAL Y(N),YPRIME(N)
C                                    .
C                                    .
C                                    .
C                           FCN SHOULD EVALUATE YPRIME(1),...,YPRIME(N)
C                             GIVEN N,X, AND Y(1),...,Y(N).  YPRIME(I)
C                             IS THE FIRST DERIVATIVE OF Y(I) WITH
C                             RESPECT TO X.
C                           FCN MUST APPEAR IN AN EXTERNAL STATEMENT IN
C                             THE CALLING PROGRAM AND N,X,Y(1),...,Y(N)
C                             MUST NOT BE ALTERED BY FCN.
C                Y      - DEPENDENT VARIABLES, VECTOR OF LENGTH N.
C                           (INPUT AND OUTPUT)
C                         ON INPUT, Y(1),...,Y(N) SUPPLY INITIAL
C                           VALUES.
C                         ON OUTPUT, Y(1),...,Y(N) ARE REPLACED WITH
C                           AN APPROXIMATE SOLUTION AT X (AS SET ON
C                           OUTPUT).
C                X      - INDEPENDENT VARIABLE. (INPUT AND OUTPUT)
C                           ON INPUT, X SUPPLIES THE INITIAL VALUE.
C                           ON OUTPUT, X IS REPLACED WITH THE UPDATED
C                             VALUE OF THE INDEPENDENT VARIABLE.
C                N      - THE NUMBER OF EQUATIONS. (INPUT)
C                JM     - THE MAXIMUM ORDER OF THE RATIONAL APPROX-
C                           IMATION. (INPUT) JM MUST BE LESS THAN 7.
C                           A SUGGESTED VALUE IS JM=6. SEE REMARKS.
C                IND    - CONVERGENCE TYPE INDICATOR. (INPUT)
C                           IND = 1 SPECIFIES THE STANDARD ERROR TEST
C                           IND = 2 SPECIFIES THE RELATIVE ERROR TEST
C                           IND = 3 SPECIFIES THE ABSOLUTE ERROR TEST
C                           SEE REMARKS FOR FURTHER DETAILS.
C                JSTART - INDICATOR. (INPUT)
C                           THE USER MUST SET JSTART TO 0 OR -1
C                           JSTART = 0 IMPLIES PERFORM A STEP.
C                             THE FIRST STEP MUST BE DONE WITH THIS
C                             VALUE OF JSTART SO THAT THE SUBROUTINE
C                             CAN INITIALIZE ITSELF.
C                           JSTART = -1 IMPLIES REPEAT THE LAST STEP
C                             WITH A NEW VALUE OF H OR JM.
C                             THE INITIAL VALUES OF Y, S,
C                             AND X ARE SET TO THE INITIAL VALUES OF Y,
C                             S, AND X FROM THE MOST RECENT CALL TO
C                             DREBS WITH JSTART = 0.
C                H      - STEP SIZE. (INPUT AND OUTPUT)
C                         ON INPUT, H IS AN INITIAL GUESS FOR THE STEP
C                           SIZE.
C                         ON OUTPUT, H IS REPLACED BY A SUGGESTED STEP
C                           SIZE FOR THE NEXT STEP. THE SUGGESTED VALUE
C                           MAY BE LARGER OR SMALLER THAN THE ORIGINAL
C                           STEP SIZE.
C                HMIN   - THE SMALLEST PERMISSIBLE STEP SIZE. (INPUT)
C                           DREBS WILL DECREASE THE STEP SIZE
C                           UNTIL CONVERGENCE CAN BE OBTAINED.
C                TOL    - TOLERANCE FOR ERROR CONTROL. (INPUT)
C                R      - VECTOR OF LENGTH N. (OUTPUT)
C                         ON OUTPUT, R CONTAINS THE ABSOLUTE
C                           ERRORS IN EACH COMPONENT FOR
C                           THE CURRENT STEP.
C                S      - VECTOR OF LENGTH N. (INPUT AND OUTPUT)
C                         IF IND = 1,
C                           BEFORE THE FIRST CALL TO THE ROUTINE,
C                           S(I) SHOULD BE SET TO Y(I), I=1,...,N.
C                           ON OUTPUT, S CONTAINS THE LARGEST VALUE
C                           OF EACH Y COMPUTED SINCE THE START OF THE
C                           INTEGRATION.
C                         IF IND = 2,
C                           BEFORE THE FIRST CALL TO THE ROUTINE,
C                           S(I) SHOULD BE SET TO Y(I), I=1,...,N.
C                           ON OUTPUT, S CONTAINS THE LARGEST VALUE
C                           OF EACH Y COMPUTED DURING THE CURRENT STEP.
C                         IF IND = 3,
C                           BEFORE THE FIRST CALL TO THE ROUTINE,
C                           S(I) SHOULD BE SET TO 1, I=1,...,N.
C                           ON OUTPUT, S IS UNCHANGED.
C                WK     - WORK VECTOR OF LENGTH 29*N.
C                           WK MUST REMAIN UNCHANGED BETWEEN SUCCESSIVE
C                           CALLS DURING INTEGRATION.
C                IER    - ERROR PARAMETER. (OUTPUT)
C                         TERMINAL ERROR
C                           IER = 129 INDICATES CONVERGENCE COULD NOT BE
C                             OBTAINED WITH CURRENT VALUES OF H AND
C                             HMIN. Y,X, AND H HAVE BEEN UPDATED.
C                         WARNING ERROR (WITH FIX)
C                           IER = 66 INDICATES JM IS LESS THAN 1 OR
C                             GREATER THAN 6. JM IS RESET TO 6.
C-----------------------------------------------------------------------
C
     SUBROUTINE DREBS  (FCN,Y,X,N,JM,IND,JSTART,H,HMIN,TOL,R,S,WK,IER)
C                                  SPECIFICATIONS FOR ARGUMENTS
     INTEGER            N,JM,IND,JSTART,IER
     DOUBLE PRECISION   Y(N),X,H,HMIN,TOL,R(N),S(N),WK(1)

C IMSL ROUTINE NAME - FFTCC
C-----------------------------------------------------------------------
C
C   PURPOSE             - COMPUTE THE FAST FOURIER TRANSFORM OF A
C                           COMPLEX VALUED SEQUENCE
C
C   USAGE               - CALL FFTCC (A,N,IWK,WK)
C
C   ARGUMENTS    A      - COMPLEX VECTOR OF LENGTH N. ON INPUT A
C                           CONTAINS THE COMPLEX VALUED SEQUENCE TO BE
C                           TRANSFORMED. ON OUTPUT A IS REPLACED BY THE
C                           FOURIER TRANSFORM.
C                N      - INPUT NUMBER OF DATA POINTS TO BE
C                           TRANSFORMED. N MAY BE ANY POSITIVE
C                           INTEGER.
C                IWK    - INTEGER WORK VECTOR OF LENGTH 6*N+150.
C                           (SEE PROGRAMMING NOTES FOR FURTHER DETAILS)
C                WK     - REAL WORK VECTOR OF LENGTH 6*N+150.
C                           (SEE PROGRAMMING NOTES FOR FURTHER DETAILS)
C
C   PRECISION/HARDWARE  - SINGLE AND DOUBLE/H32
C                       - SINGLE/H36,H48,H60
C
C   REQD. IMSL ROUTINES - NONE REQUIRED
C
C   NOTATION            - INFORMATION ON SPECIAL NOTATION AND
C                           CONVENTIONS IS AVAILABLE IN THE MANUAL
C                           INTRODUCTION OR THROUGH IMSL ROUTINE UHELP
C
C   REMARKS  1.  FFTCC COMPUTES THE FOURIER TRANSFORM, X, ACCORDING
C                TO THE FOLLOWING FORMULA;
C
C                  X(K+1) = SUM FROM J = 0 TO N-1 OF
C                           A(J+1)*CEXP((0.0,(2.0*PI*J*K)/N))
C                  FOR K=0,1,...,N-1 AND PI=3.1415...
C
C                NOTE THAT X OVERWRITES A ON OUTPUT.
C            2.  FFTCC CAN BE USED TO COMPUTE
C
C                  X(K+1) = (1/N)*SUM FROM J = 0 TO N-1 OF
C                           A(J+1)*CEXP((0.0,(-2.0*PI*J*K)/N))
C                  FOR K=0,1,...,N-1 AND PI=3.1415...
C
C                BY PERFORMING THE FOLLOWING STEPS;
C
C                     DO 10 I=1,N
C                        A(I) = CONJG(A(I))
C                  10 CONTINUE
C                     CALL FFTCC (A,N,IWK,WK)
C                     DO 20 I=1,N
C                        A(I) = CONJG(A(I))/N
C                  20 CONTINUE
C-----------------------------------------------------------------------
C
     SUBROUTINE FFTCC (A,N,IWK,WK)
C                                  SPECIFICATIONS FOR ARGUMENTS
     INTEGER            N,IWK(1)
     DOUBLE PRECISION   WK(1)
     DOUBLE COMPLEX     A(N)

C IMSL ROUTINE NAME - FFT2C
C-----------------------------------------------------------------------
C
C   PURPOSE             - COMPUTE THE FAST FOURIER TRANSFORM OF A
C                           COMPLEX VALUED SEQUENCE OF LENGTH EQUAL TO
C                           A POWER TWO
C
C   USAGE               - CALL FFT2C (A,M,IWK)
C
C   ARGUMENTS    A      - COMPLEX VECTOR OF LENGTH N, WHERE N=2**M.
C                           ON INPUT A CONTAINS THE COMPLEX VALUED
C                           SEQUENCE TO BE TRANSFORMED.
C                           ON OUTPUT A IS REPLACED BY THE
C                           FOURIER TRANSFORM.
C                M      - INPUT EXPONENT TO WHICH 2 IS RAISED TO
C                           PRODUCE THE NUMBER OF DATA POINTS, N
C                           (I.E. N = 2**M).
C                IWK    - WORK VECTOR OF LENGTH M+1.
C-----------------------------------------------------------------------
C
     SUBROUTINE  FFT2C (A,M,IWK)
C                                  SPECIFICATIONS FOR ARGUMENTS
     INTEGER            M,IWK(1)
     DOUBLE COMPLEX     A(1)

C IMSL ROUTINE NAME - FFTRC
C-----------------------------------------------------------------------
C
C   PURPOSE             - COMPUTE THE FAST FOURIER TRANSFORM OF A
C                           REAL VALUED SEQUENCE
C
C   USAGE               - CALL FFTRC (A,N,X,IWK,WK)
C
C   ARGUMENTS    A      - INPUT REAL VECTOR OF LENGTH N WHICH
C                           CONTAINS THE DATA TO BE TRANSFORMED.
C                N      - INPUT NUMBER OF DATA POINTS TO BE TRANSFORMED.
C                            N MUST BE A POSITIVE EVEN INTEGER.
C                X      - OUTPUT COMPLEX VECTOR OF LENGTH N/2+1
C                           CONTAINING THE FIRST N/2+1 COEFFICIENTS OF
C                           THE FOURIER TRANSFORM. THE REMAINING
C                           COEFFICIENTS MAY BE DETERMINED BY
C                             X(N+2-I) = CONJG(X(I)), FOR I=2,...,N/2.
C                IWK    - INTEGER WORK VECTOR.
C                           IF N IS A POWER OF 2, THEN IWK SHOULD BE OF
C                           LENGTH M WHERE N=2**M.
C                           OTHERWISE, IWK SHOULD BE OF LENGTH
C                           6*(N/2)+150.
C                           (SEE PROGRAMMING NOTES FOR FURTHER DETAILS)
C                WK     - REAL WORK VECTOR OF LENGTH 6*(N/2)+150.
C                           WK IS NOT USED IF N IS A POWER OF 2.
C                           (SEE PROGRAMMING NOTES FOR FURTHER DETAILS)
C   REMARKS  1.  FFTRC COMPUTES THE FOURIER TRANSFORM, X, ACCORDING
C                TO THE FOLLOWING FORMULA;
C
C                  X(K+1) = SUM FROM J = 0 TO N-1 OF
C                           A(J+1)*CEXP((0.0,(2.0*PI*J*K)/N))
C                  FOR K=0,1,...,N/2 AND PI=3.1415...
C
C                THE USER CAN COMPUTE THE REMAINING X VALUES BY
C                PERFORMING THE FOLLOWING STEPS;
C
C                     ND2 = N/2
C                     DO 10 I=2,ND2
C                        X(N+2-I) = CONJG(X(I))
C                  10 CONTINUE
C            2.  FFTRC CAN BE USED TO COMPUTE
C
C                  X(K+1) = (1/N)*SUM FROM J = 0 TO N-1 OF
C                           A(J+1)*CEXP((0.0,(-2.0*PI*J*K)/N))
C                  FOR K=0,1,...,N/2 AND PI=3.1415...
C
C                BY PERFORMING THE FOLLOWING STEPS;
C
C                     CALL FFTRC (A,N,X,IWK,WK)
C                     ND2P1 = N/2+1
C                     DO 10 I=1,ND2P1
C                        X(I) = CONJG(X(I))/N
C                  10 CONTINUE
C-----------------------------------------------------------------------
C
     SUBROUTINE FFTRC  (A,N,X,IWK,WK)
C                                  SPECIFICATIONS FOR ARGUMENTS
     INTEGER            N,IWK(1)
     DOUBLE PRECISION   A(N),WK(1)
     DOUBLE COMPLEX     X(1)

C IMSL ROUTINE NAME - ZBRENT
C-----------------------------------------------------------------------
C
C   PURPOSE             - ZERO OF A FUNCTION WHICH CHANGES SIGN IN A
C                           GIVEN INTERVAL (BRENT ALGORITHM)
C
C   USAGE               - CALL ZBRENT (F,EPS,NSIG,A,B,MAXFN,IER)
C
C   ARGUMENTS    F      - AN EXTERNAL FUNCTION SUBPROGRAM F(X)
C                           PROVIDED BY THE USER WHICH COMPUTES F FOR
C                           ANY X IN THE INTERVAL (A,B). (INPUT)
C                           F MUST APPEAR IN AN EXTERNAL STATEMENT IN
C                           THE CALLING PROGRAM
C                EPS    - FIRST CONVERGENCE CRITERION (INPUT).  A ROOT,
C                           B, IS ACCEPTED IF ABS(F(B)) IS LESS THAN OR
C                           EQUAL TO EPS.  EPS MAY BE SET TO ZERO.
C                NSIG   - SECOND CONVERGENCE CRITERION (INPUT).  A ROOT,
C                           B, IS ACCEPTED IF THE CURRENT APPROXIMATION
C                           AGREES WITH THE TRUE SOLUTION TO NSIG
C                           SIGNIFICANT DIGITS.
C                A,B    - ON INPUT, THE USER MUST SUPPLY TWO POINTS, A
C                           AND B, SUCH THAT F(A) AND F(B) ARE OPPOSITE
C                           IN SIGN.
C                           ON OUTPUT, BOTH A AND B ARE ALTERED.  B
C                           WILL CONTAIN THE BEST APPROXIMATION TO THE
C                           ROOT OF F. SEE REMARK 1.
C                MAXFN  - ON INPUT, MAXFN SHOULD CONTAIN AN UPPER BOUND
C                           ON THE NUMBER OF FUNCTION EVALUATIONS
C                           REQUIRED FOR CONVERGENCE.  ON OUTPUT, MAXFN
C                           WILL CONTAIN THE ACTUAL NUMBER OF FUNCTION
C                           EVALUATIONS USED.
C                IER    - ERROR PARAMETER. (OUTPUT)
C                         TERMINAL ERROR
C                           IER = 129 INDICATES THE ALGORITHM FAILED TO
C                             CONVERGE IN MAXFN EVALUATIONS.
C                           IER = 130 INDICATES F(A) AND F(B) HAVE THE
C                             SAME SIGN.
C
C   REMARKS  1.  ON EXIT FROM ZBRENT, WHEN IER=0, A AND B SATISFY THE
C                FOLLOWING,
C                F(A)*F(B) .LE.0,
C                ABS(F(B)) .LE. ABS(F(A)), AND
C                EITHER ABS(F(B)) .LE. EPS OR
C                ABS(A-B) .LE. MAX(ABS(B),0.1)*10.0**(-NSIG).
C                THE PRESENCE OF 0.1 IN THIS ERROR CRITERION CAUSES
C                LEADING ZEROES TO THE RIGHT OF THE DECIMAL POINT TO BE
C                COUNTED AS SIGNIFICANT DIGITS. SCALING MAY BE REQUIRED
C                IN ORDER TO ACCURATELY DETERMINE A ZERO OF SMALL
C                MAGNITUDE.
C            2.  ZBRENT IS GUARANTEED TO REACH CONVERGENCE WITHIN
C                K = (ALOG((B-A)/D)+1.0)**2 FUNCTION EVALUATIONS WHERE
C                  D=MIN(OVER X IN (A,B) OF
C                    MAX(ABS(X),0.1)*10.0**(-NSIG)).
C                THIS IS AN UPPER BOUND ON THE NUMBER OF EVALUATIONS.
C                RARELY DOES THE ACTUAL NUMBER OF EVALUATIONS USED BY
C                ZBRENT EXCEED SQRT(K). D CAN BE COMPUTED AS FOLLOWS,
C                  P = AMIN1(ABS(A),ABS(B))
C                  P = AMAX1(0.1,P)
C                  IF ((A-0.1)*(B-0.1).LT.0.0) P = 0.1
C                  D = P*10.0**(-NSIG)
C
C-----------------------------------------------------------------------
C
     SUBROUTINE ZBRENT (F,EPS,NSIG,A,B,MAXFN,IER)
C                                  SPECIFICATIONS FOR ARGUMENTS
     INTEGER            NSIG,MAXFN,IER
     DOUBLE PRECISION   F,EPS,A,B

C IMSL ROUTINE NAME - ZXMIN
C-----------------------------------------------------------------------
C
C   PURPOSE             - MINIMUM OF A FUNCTION OF N VARIABLES USING
C                           A QUASI-NEWTON METHOD
C
C   USAGE               - CALL ZXMIN (FUNCT,N,NSIG,MAXFN,IOPT,X,H,G,F,
C                           W,IER)
C
C   ARGUMENTS    FUNCT  - A USER SUPPLIED SUBROUTINE WHICH CALCULATES
C                           THE FUNCTION F FOR GIVEN PARAMETER VALUES
C                           X(1),X(2),...,X(N).
C                           THE CALLING SEQUENCE HAS THE FOLLOWING FORM
C                           CALL FUNCT(N,X,F)
C                           WHERE X IS A VECTOR OF LENGTH N.
C                           FUNCT MUST APPEAR IN AN EXTERNAL STATEMENT
C                           IN THE CALLING PROGRAM. FUNCT MUST NOT
C                           ALTER THE VALUES OF X(I),I=1,...,N OR N.
C                N      - THE NUMBER OF PARAMETERS (I.E., THE LENGTH
C                           OF X) (INPUT)
C                NSIG   - CONVERGENCE CRITERION. (INPUT). THE NUMBER
C                           OF DIGITS OF ACCURACY REQUIRED IN THE
C                           PARAMETER ESTIMATES.
C                           THIS CONVERGENCE CONDITION IS SATISIFIED IF
C                           ON TWO SUCCESSIVE ITERATIONS, THE PARAMETER
C                           ESTIMATES (I.E.,X(I), I=1,...,N) AGREE,
C                           COMPONENT BY COMPONENT, TO NSIG DIGITS.
C                MAXFN  - MAXIMUM NUMBER OF FUNCTION EVALUATIONS (I.E.,
C                           CALLS TO SUBROUTINE FUNCT) ALLOWED. (INPUT)
C                IOPT   - OPTIONS SELECTOR. (INPUT)
C                         IOPT = 0 CAUSES ZXMIN TO INITIALIZE THE
C                           HESSIAN MATRIX H TO THE IDENTITY MATRIX.
C                         IOPT = 1 INDICATES THAT H HAS BEEN INITIALIZED
C                           BY THE USER TO A POSITIVE DEFINITE MATRIX.
C                         IOPT = 2 CAUSES ZXMIN TO COMPUTE THE DIAGONAL
C                           VALUES OF THE HESSIAN MATRIX AND SET H TO
C                           A DIAGONAL MATRIX CONTAINING THESE VALUES.
C                         IOPT = 3 CAUSES ZXMIN TO COMPUTE AN ESTIMATE
C                           OF THE HESSIAN IN H.
C                X      - VECTOR OF LENGTH N CONTAINING PARAMETER
C                           VALUES.
C                         ON INPUT, X MUST CONTAIN THE INITIAL
C                           PARAMETER ESTIMATES.
C                         ON OUTPUT, X CONTAINS THE FINAL PARAMETER
C                           ESTIMATES AS DETERMINED BY ZXMIN.
C                H      - VECTOR OF LENGTH N*(N+1)/2 CONTAINING AN
C                           ESTIMATE OF THE HESSIAN MATRIX
C                           D**2F/(DX(I)DX(J)), I,J=1,...,N.
C                           H IS STORED IN SYMMETRIC STORAGE MODE.
C                         ON INPUT, IF IOPT = 0, 2, OR 3 ZXMIN INITIA-
C                           LIZES H. AN INITIAL SETTING OF H BY THE
C                           USER IS INDICATED BY IOPT=1.
C                           H MUST BE POSITIVE DEFINITE. IF IT IS NOT,
C                           A TERMINAL ERROR OCCURS.
C                         ON OUTPUT, H CONTAINS AN ESTIMATE OF THE
C                           HESSIAN AT THE FINAL PARAMETER ESTIMATES
C                           (I.E., AT X(1),X(2),...,X(N))
C                G      - A VECTOR OF LENGTH N CONTAINING AN ESTIMATE
C                           OF THE GRADIENT DF/DX(I),I=1,...,N AT THE
C                           FINAL PARAMETER ESTIMATES. (OUTPUT)
C                F      - A SCALAR CONTAINING THE VALUE OF THE FUNCTION
C                           AT THE FINAL PARAMETER ESTIMATES. (OUTPUT)
C                W      - A VECTOR OF LENGTH 3*N USED AS WORKING SPACE.
C                         ON OUTPUT, WORK(I), CONTAINS FOR
C                           I = 1, THE NORM OF THE GRADIENT (I.E.,
C                             SQRT(G(1)**2+G(2)**2+...+G(N)**2))
C                           I = 2, THE NUMBER OF FUNCTION EVALUATIONS
C                             PERFORMED.
C                           I = 3, AN ESTIMATE OF THE NUMBER OF
C                             SIGNIFICANT DIGITS IN THE FINAL
C                             PARAMETER ESTIMATES.
C                IER    - ERROR PARAMETER (OUTPUT)
C                         TERMINAL ERROR
C                           IER = 129 IMPLIES THAT THE INITIAL HESSIAN
C                             USED BY ZXMIN IS NOT POSITIVE DEFINITE,
C                             EVEN AFTER ADDING A MULTIPLE OF THE
C                             IDENTITY TO MAKE ALL DIAGONAL ELEMENTS
C                             POSITIVE.
C                           IER = 130 IMPLIES THAT THE ITERATION WAS
C                             TERMINATED DUE TO ROUNDING ERRORS
C                             BECOMING DOMINANT. THE PARAMETER
C                             ESTIMATES HAVE NOT BEEN DETERMINED TO
C                             NSIG DIGITS.
C                           IER = 131 IMPLIES THAT THE ITERATION WAS
C                             TERMINATED BECAUSE MAXFN WAS EXCEEDED.
C
C-----------------------------------------------------------------------
C
     SUBROUTINE ZXMIN (FUNCT,N,NSIG,MAXFN,IOPT,X,H,G,F,W,IER)
C                                  SPECIFICATIONS FOR ARGUMENTS
     INTEGER            N,NSIG,MAXFN,IOPT,IER
     DOUBLE PRECISION   X(N),G(N),H(1),F,W(1)

C IMSL ROUTINE NAME - ZXCGR
C-----------------------------------------------------------------------
C
C   PURPOSE             - A CONJUGATE GRADIENT ALGORITHM FOR FINDING
C                           THE MINIMUM OF A FUNCTION OF N VARIABLES
C
C   USAGE               - CALL ZXCGR (FUNCT,N,ACC,MAXFN,DFPRED,X,G,F,W,
C                           IER)
C
C   ARGUMENTS    FUNCT  - A USER SUPPLIED SUBROUTINE WHICH CALCULATES
C                           THE OBJECTIVE FUNCTION AND ITS GRADIENT
C                           FOR GIVEN PARAMETER VALUES
C                           X(1),X(2),...,X(N).
C                           THE CALLING SEQUENCE HAS THE FOLLOWING FORM
C                           CALL FUNCT (N,X,F,G)
C                           WHERE X AND G ARE VECTORS OF LENGTH N.
C                           THE SCALAR F IS FOR THE OBJECTIVE FUNCTION.
C                           G(1), G(2), ..., G(N) ARE FOR THE COMPONENTS
C                           OF THE GRADIENT OF F.
C                           FUNCT MUST APPEAR IN AN EXTERNAL STATEMENT
C                           IN THE CALLING PROGRAM. FUNCT MUST NOT
C                           ALTER THE VALUES OF X(I),I=1,...,N OR N.
C                N      - THE NUMBER OF PARAMETERS OF THE OBJECTIVE
C                           FUNCTION. (INPUT) (I.E.,THE LENGTH OF X)
C                ACC    - CONVERGENCE CRITERION. (INPUT)
C                           THE CALCULATION ENDS WHEN THE SUM OF SQUARES
C                           OF THE COMPONENTS OF G IS LESS THAN ACC.
C                MAXFN  - MAXIMUM NUMBER OF FUNCTION EVALUATIONS (I.E.,
C                           CALLS TO SUBROUTINE FUNCT) ALLOWED. (INPUT)
C                           IF MAXFN IS SET TO ZERO, THEN THERE IS
C                           NO RESTRICTION ON THE NUMBER OF FUNCTION
C                           EVALUATIONS.
C                DFPRED - A ROUGH ESTIMATE OF THE EXPECTED REDUCTION
C                           IN F, WHICH IS USED TO DETERMINE THE SIZE
C                           OF THE INITIAL CHANGE TO X. (INPUT)
C                           NOTE THAT DFPRED IS THE EXPECTED REDUCTION
C                           ITSELF, AND DOES NOT DEPEND ON ANY RATIOS.
C                           A BAD VALUE OF DFPRED CAUSES AN ERROR
C                           MESSAGE, WITH IER=129, AND A RETURN ON THE
C                           FIRST ITERATION. (SEE THE DESCRIPTION OF
C                           IER BELOW)
C                X      - VECTOR OF LENGTH N CONTAINING PARAMETER
C                           VALUES.
C                         ON INPUT, X MUST CONTAIN THE INITIAL
C                           PARAMETER ESTIMATES.
C                         ON OUTPUT, X CONTAINS THE FINAL PARAMETER
C                           ESTIMATES AS DETERMINED BY ZXCGR.
C                G      - A VECTOR OF LENGTH N CONTAINING THE
C                           COMPONENTS OF THE GRADIENT OF F AT THE
C                           FINAL PARAMETER ESTIMATES. (OUTPUT)
C                F      - A SCALAR CONTAINING THE VALUE OF THE FUNCTION
C                           AT THE FINAL PARAMETER ESTIMATES. (OUTPUT)
C                W      - WORK VECTOR OF LENGTH 6*N.
C                IER    - ERROR PARAMETER. (OUTPUT)
C                         IER = 0 IMPLIES THAT CONVERGENCE WAS
C                           ACHIEVED AND NO ERRORS OCCURRED.
C                         TERMINAL ERROR
C                           IER = 129 IMPLIES THAT THE LINE SEARCH OF
C                             AN INTEGRATION WAS ABANDONED. THIS
C                             ERROR MAY BE CAUSED BY AN ERROR IN THE
C                             GRADIENT.
C                           IER = 130 IMPLIES THAT THE CALCULATION
C                             CANNOT CONTINUE BECAUSE THE SEARCH
C                             DIRECTION IS UPHILL.
C                           IER = 131 IMPLIES THAT THE ITERATION WAS
C                             TERMINATED BECAUSE MAXFN WAS EXCEEDED.
C                           IER = 132 IMPLIES THAT THE CALCULATION
C                             WAS TERMINATED BECAUSE TWO CONSECUTIVE
C                             ITERATIONS FAILED TO REDUCE F.
C
C-----------------------------------------------------------------------
C
     SUBROUTINE ZXCGR  (FUNCT,N,ACC,MAXFN,DFPRED,X,G,F,W,IER)
C                                  SPECIFICATIONS FOR ARGUMENTS
     INTEGER            N,MAXFN,IER
     DOUBLE PRECISION   ACC,DFPRED,X(N),G(N),F,W(1)