DCADRE : integration
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-----------------------------------------------------------------------
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
ICSCCU : cubic spline interpolation
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-----------------------------------------------------------------------
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)
ICSEVU : cubic spline evaluation
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-----------------------------------------------------------------------
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).
MDCH : chi-suared probability distribution
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-----------------------------------------------------------------------
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
DREBS : differential equation solver
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-----------------------------------------------------------------------
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)
FFTCC : FFT for complex valued sequence
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-----------------------------------------------------------------------
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)
FFT2C : FFT for complex valued sequence of length equal to a power two
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-----------------------------------------------------------------------
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)
FFTRC : FFT for real valued sequence
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-----------------------------------------------------------------------
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)
ZBRENT : zero of a function : Brent algorithm
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-----------------------------------------------------------------------
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
ZXMIN : minimum of a function : quasi-Newton method
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-----------------------------------------------------------------------
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)
ZXCGR : minimum of a function : commjugate gradient algorithm
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)
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)
Comments List
안봐져 ㅠ_ㅠ; 낚인거? @_@;;
음.. FF에서는 문제가 없는데 말이지. 'ㅁ'a 태터툴즈 업데이트하면서 어딘가에서 에러가 났나봐. 예전 글에 있던 펼침 목록도 IE에서는 작동이 안 되네. 조만간 고쳐 놓도록 서버 관리자에게 압박을 가해야겠군. 꼭 지금 보고 싶다면 소스보기를..;;;
안보여.. 가 어법에 맞는 것 아니겠소.
플러그인 충돌에 한 표 입니당.
(왜냐면 나는 잘 되거든)
켜놓은 플러그인 이야기해 주세요 :)