Trailing-Edge
-
PDP-10 Archives
-
decuslib10-09
-
43,50466/itgtn.for
There are no other files named itgtn.for in the archive.
C WESTERN MICHIGAN UNIVERSITY
C ITGTN.FOR (FILENAME ON LIBRARY DECTAPE)
C ITGTN.FOR IS CALLED (BY RUNUUO FROM INTGR.FOR
C EXTERNAL FUNCTIONS USED (GENERATED BY INTGR.FOR): SEE
C COMMENTS AT BEGINNING OF SOURCE LISTING OF INTGR.FOR
C INTERNAL FUNCTIONS USED: QMULT2, QMULT3,SQUANK
C INTERNAL SUBR. USED: INT2, INT3, GLO16
C ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
C THIS IS THE SECOND PART OF THE INTEGRATION PROGRAM IN WHICH
C THE FOLLOWING VARIABLES HAVE TO BE READ IN FROM THE TEMPORARY
C FILE:
C
C IDIM---DIMENSION
C FF-----THE FUNCTION TO BE INTEGRATED
C FFL1---LOWER LIMIT OF THE INNER INTEGRAL
C FFU1---UPPER LIMIT OF THE INNER INTEGRAL
C FFL2---LOWER LIMIT OF THE MIDDLE INTEGRAL IN THE 3-
C DIMENSION INTEGRATION
C FFU2---UPPER LIMIT OF THE MIDDLE INTEGRAL
C
C
C THE FOLLOWING ARE FUNCTIONS AND/OR LIMITS TO BE USED IN THE
C INTEGRATION:
C
C F----FUNCTION TO BE INTEGRATED
C FL1--LOWER LIMIT OF THE INNER INTEGRAL
C FU1--UPPER LIMIT OF THE INNER INTEGRAL
C FL2--LOWER LIMIT OF THE MIDDLE INTEGRAL IN THE 3-DIMENSION
C INTEGRATION
C FU2--UPPER LIMIT OF THE MIDDLE INTEGRAL
C FL3--LOWER LIMIT OF THE OUTER INTEGRAL
C FU3--UPPER LIMIT OF THEOUTER INTEGRAL
C
C***********************************************************************
C
DIMENSION FF(12),FFU1(12),FFL1(12),FFU2(12),FFL2(12)
EXTERNAL F,FU1,FL1,FU2,FL2
C
C***********************************************************************
C DEVICES USED: SEE PART 1
C***********************************************************************
C
IDLG=-1
ICC=-4
IOUT=30
IDSK=1
ITEMP=20
C
C***********************************************************************
C READ IN INFORMATION THAT WERE TEMPORARILY STORED ON THE
C DISK BY THE FIRST PART OF THE PROGRAM.
C***********************************************************************
C
READ(ITEMP) IDIM,FF
READ(ITEMP) FFL1,FFU1
READ(ITEMP) FFL2,FFU2
CALL RELEAS (ITEMP)
C
C***********************************************************************
C ASK AND OBTAIN OUTER LIMITS--MUST BE REAL NUMBERS
C**********************************************************************
C
IF (IDIM.EQ.1) WRITE(IDLG,61)
IF (IDIM.GT.1) WRITE(IDLG,580)
580 FORMAT('-ENTER OUTER LIMITS:')
WRITE(IDLG,101)
101 FORMAT(' LOWER: ',$)
READ(ICC,53) FL3
WRITE(IDLG,102)
102 FORMAT('+UPPER: ',$)
READ(ICC,53) FU3
IF (IDIM.GT.1) GO TO 540
C
C***********************************************************************
C ONE-DIMENSION INTEGRATION
C***********************************************************************
C
WRITE(IDLG,541)
541 FORMAT(' ENTER VALUE FOR THE REQUIRED TOLERANCE: ',$)
READ(ICC,53) ERROR
53 FORMAT(F)
C---------------SQUANK MAKES USE OF FUNCTION F GENERATED BY INTGR.FOR
CALL SQUANK(FL3,FU3,ERROR,FIFTH,RUM,NO)
S=SQUANK(FL3,FU3,ERROR,FIFTH,RUM,NO)
C
C***********************************************************************
C START OF THE REPORT
C***********************************************************************
C
540 DO 5410 J=12,1,-1
IF (FF(J).NE.' ') GO TO 542
5410 CONTINUE
544 WRITE(IDLG,543) FF
543 FORMAT('- ERROR IN FUNCTION, JOB TERMINATES'/1X,12A5)
CALL EXIT
542 WRITE(IOUT,54) IDIM,(FF(JJ),JJ=1,J)
54 FORMAT('1',I1,'-DIMENSION INTEGRATION'//' FUNCTION: ',12A5)
IF (IDIM.LT.2) GO TO 60
C
C***********************************************************************
C TWO AND THREE-DIMENSION INTEGRATION
C***********************************************************************
C
WRITE(IOUT,30)
30 FORMAT('-LIMITS FOR INNER INTEGRAL:')
DO 550 J=12,1,-1
IF (FFL1(J).NE.' ') GO TO 551
550 CONTINUE
GO TO 544
551 WRITE(IOUT,55) (FFL1(JJ),JJ=1,J)
55 FORMAT(' LOWER: ',12A5)
DO 560 J=12,1,-1
IF (FFU1(J).NE.' ') GO TO 561
560 CONTINUE
GO TO 544
561 WRITE(IOUT,56) (FFU1(JJ),JJ=1,J)
56 FORMAT(' UPPER: ',12A5)
IF (IDIM.LT.3) GO TO 57
WRITE(IOUT,43)
43 FORMAT('-LIMITS FOR MIDDLE INTEGRAL:')
DO 5500 J=12,1,-1
IF (FFL2(J).NE.' ') GO TO 5501
5500 CONTINUE
GO TO 544
5501 WRITE(IOUT,55) (FFL2(JJ),JJ=1,J)
DO 5503 J=12,1,-1
IF (FFU2(J).NE.' ') GO TO 5504
5503 CONTINUE
GO TO 544
5504 WRITE(IOUT,56) (FFU2(JJ),JJ=1,J)
57 WRITE(IOUT,58)
58 FORMAT('-OUTER LIMITS:')
WRITE(IOUT,5510) FL3
WRITE(IOUT,5610) FU3
IF (IDIM.EQ.2) CALL INT2(FU3,FL3,IOUT)
IF (IDIM.EQ.3) CALL INT3(FU3,FL3,IOUT)
GO TO 70
C
C***********************************************************************
C CONTINUATION OF REPORT FOR ONE-DIMENSION INTEGRATION
C***********************************************************************
C
60 WRITE(IOUT,61)
61 FORMAT('-LIMITS OF INTEGRATION')
WRITE(IOUT,5510) FL3
5510 FORMAT(' LOWER: ',G15.7)
WRITE(IOUT,5610) FU3
5610 FORMAT(' UPPER: ', G15.7)
WRITE(IOUT,65) S,FIFTH,NO,RUM
65 FORMAT('-ANSWER: ',G15.7/'-FIFTH-ORDER ADJUSTMENT TERM ',G15.7/
1 ' NUMBER OF FUNCTION EVALUATION REQUIRED: ',G15.7/' THE CLAIMED
1 TOLERANCE (ADJUSTED FOR ROUNDOFF ERROR): ',G15.7)
C
C***********************************************************************
C END OF ONE INTEGRATION, EITHER CALL EXIT OR BRANCH TO FIRST
C PART OF THE PROGRAM
C***********************************************************************
C
70 WRITE(IDLG,71)
71 FORMAT('-END OF INTEGRATION, TYPE 1 TO CONTINUE OR 0 TO
1 TERMINATE'/)
READ(ICC,14) LAST
14 FORMAT(I)
C---------------INTGR.EXE IS IN SYS:
IF (LAST.GT.0) CALL RUNUUO('R INTGR')
CALL EXIT
END
C---------------ALL ARGS. ARE INPUT.
SUBROUTINE INT2(FU3,FL3,IOUT)
C THIS PROGRAM APPROXIMATES A 2-DIMENSIONAL ITERATED
C INTEGRAL OF F(X,Y) AS DISCUSSED IN SECTION 1.6.
C SUBROUTINE GLO16 PROVIDES A TABLE OF THE 16-POINT
C GAUSS-LEGENDRE FORMULA. THE LIMITS ON THE VARIABLES ARE
C AA .LE. X .LE. BB
C FL1(X) .LE. Y. LE. FU1(X)
C
DOUBLE PRECISION DX(16),DA(16)
DIMENSION X(40,2),A(40,2),MM(2)
EXTERNAL F,FL1,FU1
CALL GLO16(DX,DA,-1.D0,1.D0)
DO 2 I = 1,16
XX=DX(I)
AB=DA(I)
X(I,1)=XX
X(I,2)=XX
A(I,1)=AB
2 A(I,2)=AB
MM(1)=16
MM(2)=16
Q = QMULT2(F,FL3,FU3,FU1,FL1,X,A,MM)
WRITE(IOUT,4) Q
4 FORMAT('-ANSWER: ',G15.7)
RETURN
END
C---------------ALL ARGS. ARE INPUT
FUNCTION QMULT2(FCN,AA,BB,FU,FL,X,A,MM)
DIMENSION X(40,2),A(40,2),MM(2)
H1=(BB-AA)/2.
G1=(BB+AA)/2.
Q1=0.
M1=MM(1)
M2=MM(2)
DO 4 I=1,M1
UI=H1*X(I,1)+G1
AI=H1*A(I,1)
D=FU(UI)
C=FL(UI)
H=(D-C)/2.
G=(D+C)/2.
Q=0.
DO 2 J=1,M2
VJ=H*X(J,2)+G
C---------------THERE IS NO FUNCTION FCN; HOWEVER THE EXTERNAL F
C--------------- IN SUBR. INT2 AND THE F IN ARG. LIST OF QMULT2
C--------------- IN SUBR. INT2 ALLOWS FOR SUBST. OF F FOR FCN
2 Q=Q+A(J,2)*FCN(UI,VJ)
4 Q1=Q1+AI*H*Q
QMULT2 = Q1
RETURN
END
C---------------ALL ARGS. ARE INPUT.
SUBROUTINE INT3(FU3,FL3,IOUT)
C
C THIS PROGRAM APPROXIMATES A 3-DIMENSIONAL ITERATED
C INTEGRAL OF F(X,Y,Z) AS DISCUSSED IN SECTION 1.6.
C SUBROUTINE GLO16 PROVIDES A TABLE OF THE 16-POINT
C GAUSS-LEGENDRE FORMULA. THE LIMITS ON THE VARIABLES ARE
C AA .LE. X .LE. BB
C FL2(X) .LE. Y .LE. FU2(X)
C FL1(X,Y) .LE. Z .LE. FU1(X,Y)
C
DOUBLE PRECISION DX(16), DA(16)
DIMENSION X(40,3),A(40,3), MM(3)
EXTERNAL F,FL1,FU1,FL2,FU2
CALL GLO16(DX,DA,-1.D0,1.D0)
DO 2 I=1,16
XX=DX(I)
X(I,1)=XX
X(I,2)=XX
X(I,3)=XX
AB=DA(I)
A(I,1)=AB
A(I,2)=AB
2 A(I,3)=AB
MM(1)=16
MM(2)=16
MM(3)=16
Q = QMULT3(F, FL3, FU3, FU2, FL2, FU1, FL1, X, A, MM)
WRITE(IOUT,4) Q
4 FORMAT('-ANSWER: ',G15.7)
RETURN
END
C---------------ALL ARGS. ARE INPUT.
FUNCTION QMULT3(FCN,AA,BB,FU1,FL1,FU2,FL2,X,A,MM)
DIMENSION X(40,3), A(40,3), MM(3)
H1= (BB-AA)/2.
G1= (BB+AA)/2.
Q1= 0.
M1= MM(1)
M2= MM(2)
M3= MM(3)
DO 6 I = 1,M1
UI = H1*X(I,1) + G1
AI = H1*A(I,1)
D1 = FU1(UI)
C1 = FL1(UI)
H2 = (D1-C1)/2.
G2 = (D1+C1)/2.
Q2 = 0.
DO 4 J=1,M2
VJ = H2*X(J,2) + G2
AJ = H2*A(J,2)
D = FU2(UI,VJ)
C = FL2(UI,VJ)
H = (D-C)/2.
G = (D+C)/2.
Q = 0.
DO 2 K=1,M3
WK = H*X(K,3)+G
2 Q = Q + A(K,3)*FCN(UI,VJ,WK)
4 Q2= Q2+ AJ*H*Q
6 Q1= Q1 + AI*Q2
QMULT3=Q1
RETURN
END
C---------------A, BIG, ERROR ARE INPUT. FIFTH, RUM, NO ARE RETURNED.
FUNCTION SQUANK(A,BIG,ERROR,FIFTH,RUM,NO)
C
C S*Q*U*A*N*K STANDS FOR * SIMPSON QUADRATURE USED
C ADAPTIVELY. NOISE KILLED.*
C
C
C
C
C
C CALLING PROGRAM REQUIRES
C EXTERNAL F
C THIS IS FUNCTION TO BE INTEGRATED
C A THE LOWER LIMIT OF INTEGRATION
C BIG THE UPPER LIMIT
C ERROR THE REQUIRED TOLERANCE (ABSOLUTE ERROR)
C
C
C
C
C OUTPUT
C
C SQUANK THE FIFTH ORDER RESULT = THIRD + FIFTH
C FIFTH THE FIFTH ORDER ADJUSTMENT TERM
C RUM THE CLAIMED TOLERANCE (ADJUSTED FOR ROUNDOFF ERROR)
C NO THE NUMBER OF FUNCTION EVALUATIONS REQUIRED
C
C
C
C
C NOTES ON USE. (1) DISCONTINUOUS FUNCTIONS
C
C THIS ROUTINE IS BASED ON DEGREE 3 AND DEGREE 5 LOCAL
C POLYNOMIAL APPROXIMATION. CONSEQUENTLY IT SHOULD NOT BE USED
C WITH FUNCTIONS WHICH HAVE DISCONTINUITIES IN THE FOURTH OR LOWER
C DERIVATIVES WITHIN THE INTERVAL OF INTEGRATION. IF THERE ARE
C SUCH DISCONTINUITIES, THIS ROUTINE WILL TAKE THIS TO BE EVIDENCE
C OF ROUND OFF ERROR IN FUNCTION VALUES AND WILL ADJUST THE
C TOLERANCE.
C
C IF THE LOCATIONS OF SUCH DISCONTINUITIES ARE KNOWN, THE
C ROUTINE MAY BE USED SEPARATELY FOR EACH INTERVAL BETWEEN
C CONSECUTIVE DISCONTINUITIES. WHILE, LIKE ALL SUCH ROUTINES, IT
C DISLIKES DISCONTINUITIES, IT CAN HANDLE THEM IF THEY ARE LOCATED
C AT THE END POINTS OF THE INTEGRATION INTERVAL.
C
C
C
C
C NOTES ON USE. (2) FUNCTIONS WITH HIGH-FREQUENCY OSCILLATIONS.
C
C THE ROUTINE WILL RETURN UNRELIABLE RESULTS FOR FUNCTIONS LIKE
C G(X) TIMES COS(100*X). IF THE HIGHEST PERIOD LIKELY TO BE
C ENCOUNTERED IS KNOWN, THE INTERVAL SHOULD BE SUB-DIVIDED IN
C SUCH A WAY THAT,
C
C (ONE) THERE ARE NOT MORE THAN THREE PERIODS PER INTERVAL, AND
C (TWO) THE PERIOD *P* DIVIDED BY SUB-INTERVAL, *(B-A)* IS NOT A
C SIMPLE FRACTION *N/M* WITH N OR M LESS THAN 9.
C
C
C
C
C NOTES ON USE. (3) INTERVAL SUB-DIVISION
C
C THE FAILURES DESCRIBED ABOVE ARE GENERALLY WORSE FOR *SQUANK*
C THAN FOR OTHER ROUTINES BECAUSE *SQUANK* TAKES THE INCONVENIENT
C BEHAVIOR AS AN INDICATION OF ROUND OFF ERROR. IN GENERAL SUB-
C DIVISION OF THE INTERVAL IS ADVOCATED. ESSENTIALLY THE USER
C CARRIES OUT, UNDER DRIVING PROGRAM CONTROL, A SEQUENCE OF
C CALCULATIONS WHICH SHOULD HAVE BEEN CARRIED OUT IN THE
C SUBROUTINE IN ANY CASE. IN THIS WAY HE PREVENTS CHANCE LOW
C ORDER FALSE CONVERGENCE AT VIRTUALLY NO ADDITIONAL COST. NOTE
C THAT THE SUM OF THE PARAMETERS *ERROR* FOR THE SUB-INTERVALS
C SHOULD CORRESPOND TO THE VALUE REQUIRED FOR THE WHOLE INTERVAL.
C
C
C
C
C NIM NUMBERING SYSTEM AND LOGIC
C
C THE INTERVAL (A,B) IS DEFINED NIM = 1, LEVEL = 0.
C
C THE INTERVAL NIM = N, LEVEL = L IS BISECTED, IF NECESSARY,
C INTO TWO INTERVALS, NIM = 2*N AND NIM = 2*N+1,
C BOTH AT LEVEL = L+1.
C
C IF INTERVAL NIM = N, LEVEL = L DOES NOT CONVERGE, THE NEXT
C INTERVAL CONSIDERED IS NIM = 2*N, LEVEL = L+1.
C
C IF INTERVAL NIM=N, LEVEL = L DOES CONVERGE, THE NEXT INTERVAL
C CONSIDERED IS NIM = M(R) + 1, LEVEL = L-R, WHERE M(R)
C IS THE FIRST EVEN MEMBER OF THE SEQUENCE M(0) = N,
C M(S+1) =(M(S)-1)/2. IF THIS GIVES LEVEL = 0,
C THE CALCULATION IS COMPLETE.
C
C
C
C SCALING TO AVOIDING EXCESSIVE DIVISION BY TWO.
C
CTHE INTERVAL (X1,X5) IS OF LENGTH H = X5-X1. THE POINTS
C X1,X2,X3,X4,X5 ARE THE POINTS OF QUARTERSECTION OF THIS INTERVAL
C AND FX1,FX2,FX3,FX4,FX5 ARE THE CORRESPONDING FUNCTION VALUES.
C
C EST IS APPROXIMATION TO (6.0/H)* INTEGRAL(X1,X5).
C
C (EST1+EST2) IS APPROXIMATION TO (12.0/H)*INTEGRAL(X1,X5).
C
C SUM IS APPROXIMATION TO (12.0) * INTEGRAL(A,X1).
C
C
C
C
C STORAGE
C
C X3ST(L) =0.5*(X5ST(L)+X1). THUS X3ST(L) COULD BE RECALCULATED
C AT EACH STAGE TO AVOID STORAGE. ESTST(L) IS SAME
C IN THIS RESPECT. THE RESULTS OF ABOVE RECALCULATION
C ARE IDENTICAL MACHINE NUMBERS.
C
C X5ST(L) = X1 + (B-A)*(2**(-L)). THIS COULD ALSO BE RECALCU-
C LATED. BUT IN THIS CASE CALCULATION IS EXCESSIVE
C AND THERE IS A POSSIBILITY OF ROUND OFF ERROR ARISING
C BECAUSE THE SAME POINT IS BEING CALCULATED IN TWO OR
C MORE DIFFERENT WAYS.
C
C
C
C
C AVOIDANCE OF ROUND OFF ERROR TROUBLE
C
C IF INTERVAL DOES NOT CONVERGE, FOLLOWING INTERVAL SHOULD HAVE
C ADIFF VALUE APPROXIMATELY EQUAL TO (1/16) TIMES PREVIOUS ADIFF
C VALUES, CALLED ADIFF1 IN THE CODE. THERE IS A THOREM WHICH
C STATES THAT, UNLESS THE FOURTH DERIVATION OF F(X) VANISHES
C IN THE PREVIOUS INTERVAL, ADIFF IS LESS THAN OR EQUAL TO ADIFF1.
C IF THIS DOES NOT HAPPEN, IT IS TAKEN TO BE AN INDICATION OF
C POSSIBLE ROUND OFF LEVEL. IN THIS CASE, UNLESS LEV IS LESS
C THAN FIVE, THE CURRENT TOLERANCE LEVEL, CEPS, IS APPROPRIATELY
C ADJUSTED. HOWEVER CEPS IS RESET AS AND WHEN APPEARS THAT IT
C SHOULD BE ADJUSTED EITHER UP OR DOWN. IT IS REDUCED IF
C CONVERGENCE OCCURS WITH A NON-ZERO ADIFF STRICTLY LESS THAN
C 0.25*CEPS. AN INVOLVED SECTION OF CODING GUARDS TO SOME
C EXTENT AGAINST AN UNREALISTIC VALUE ARISING AS A RESULT OF
C A ZERO IN THE FOURTH DERIVATIVE. A FACTOR EFACT IS CALCULATED
C WHICH ADJUSTS THE CLAIMED TOLERANCE TO TAKE INTO ACCOUNT
C THESE ALTERATIONS IN THE TOLERANCE LEVEL.
C
C THE ROUTINE ENTERS THESE INVOLVED SECTIONS OF CODING ONLY IF
C ROUND OFF ERROR APPEARS TO BE PRESENT. IN A NORMAL (ROUND OFF
C ERROR FREE) RUN, THESE SECTIONS ARE SKIPPED AT A COST OF A
C SINGLE COMPARISION PER ITERATION (TWO FUNCTION EVALUATIONS).
C
C
C ARBITRARY CONSTANTS
C
C THE FOLLOWING CONSTANTS HAVE BEEN ASSIGNED IN THE LIGHT OF
C EXPERIENCE WITH NO THEORETICAL JUSTIFICATON.
C
C (1) NO CONVERGENCE IS ALLOWED AT LEVEL = **ZERO**. THIS MEANS
C THAT THE ROUTINE IS CONSTRAINED TO BASE THE RESULT ON AT
C LEAST 9 FUNCTION VALUES.
C
C (2) NO UPWARD ADJUSTMENT OF THE TOLERANCE LEVEL IS CONSIDERED
C AT LEVELS LOWER THAN LEVEL = **FIVE**. THE POINT SPACING
C IS THEN (BIG-A)/128.0.
C
C (3) PHYSICAL LIMIT. HIGHEST LEVEL ALLOWED IS LEVEL =
C **THIRTY**. HERE CONVERGENCE IS ASSIGNED WHETHER OR NOT
C THE INTERVAL HAS CONVERGED. THE POINT SPACING IS THEN
C ABOUT (BIG-A)*2.0*10**-10.
C
C (4) UPWARD ADJUSTMENT OF TOLERANCE LEVEL IS LIMITED IN GENERAL
C TO A FACTOR **2.0** OR LESS.
C
C (5) DOWNWARD ADJUSTMENT OF TOLERANCE LEVEL IS INHIBITED IN
C GENERAL UNLESS BY A FACTOR GREATER THAN **4.0**.
C
C
C
C
C SOME NOTATION
C
C SUM AND SIM ARE RUNNING SUMS, INCREASED AT STAGE EIGHT. THEY
C ARE RESPECTIVELY 12.0 * (THIRD ORDER APPROXIMATION TO THE
C INTEGRAL) AND -180.0* (FIFTH ORDER ADJUSTMENT TO THE INTEGRAL).
C
C CEPSF IS THE REQUIRED (SCALED) TOLERANCE.
C
C CEPS IS THE RUNNING VALUE OF THE ADJUSTED TOLERANCE.
C
C QCEPS = 0.25 * CEPS
C
C LEVTAG = -1 OR 0,2,3 INDICATES WHETHER TOLERANCE IS NOT OR IS
C CURRENTLY ADJUSTED. (SEE COMMENT IN STAGE SEVEN).
C
C EFACT IS RUNNING SUM CORRESPONDING TO 180.0 * RUM
C
C FACERR = 15.0 OR 1.0 DEPENDING ON WHETHER TOLERANCE IS OR
C IS NOT CURRENTLY ADJUSTED. IF IT IS, THERE IS NO JUSTIFICATION
C FOR THE DIFFERENCE OF APPROXIMATIONS. FACERR = 15.0 REMOVES
C THE BUILT IN 15.0 FACTOR FOR CALCULATION OF EFACT.
C
C EPMACH THE MACHINE ACCURACY PARAMETER. THE ROUND OFF ERROR
C GUARD DOES NOT REQUIRE THIS NUMBER. IT IS MACHINE INDEPENDENT.
C THIS IS ONLY USED TO HELP IN AN INITIAL GUESS IN STAGE TWO
C IF THE VALUE OF ERROR HAPPENS TO BE ZERO. ANY NON-ZERO NUMBER
C MAY BE USED INSTEAD, WITH A VERY SMALL PENALTY IN NUMBER OF
C FUNCTION EVALUATIONS IF A COMPLETELY UNREASONABLE NUMBER IS
C USED.
C
C
C
DIMENSION FX3ST(30),X3ST(30),ESTST(30),FX5ST(30),X5ST(30)
DIMENSION PREDIF(30)
C---------------EXTERNAL F IS APPARENTLY NOT NECESSARY
EXTERNAL F
DOUBLE PRECISION SUM, SIM
EPMACH = 0.0000000000075
C
C **** STAGE ONE ****
C INITIALISE ALL QUANTITIES REQUIRED FOR CENTRAL CALCULATION
C (STAGE 3).
C
SUM=0.
SIM=0.
CEPSF=180.0*ERROR/(BIG-A)
CEPS=CEPSF
ADIFF = 0.0
LEVTAG = -1
FACERR = 1.0
XZERO = A
EFACT= 0.0
NIM=1
LEV=0
C
C FIRST INTERVAL
C
X1=A
X5=BIG
X3=0.5*(A+BIG)
FX1=F(X1)
FX3=F(X3)
FX5=F(X5)
NO=3
EST= FX1 + FX5 + 4.0*FX3
C
C **** STAGE TWO ****
C SET A STARTING VALUE FOR TOLERANCE IN CASE THAT CEPSF = 0.0
C
IF (CEPSF) 295,205,295
205 LEVTAG = 0
FACERR = 15.0
CEPS = EPMACH*ABS(FX1)
IF (FX1) 295,210,295
210 CEPS = EPMACH*ABS(FX3)
LEVTAG = 3
IF (FX3) 295,215,295
215 CEPS =EPMACH*ABS(FX5)
IF (FX5) 295,220,295
220 CEPS = EPMACH
295 QCEPS = 0.25*CEPS
C
C INITIALISING COMPLETE
C
C **** STAGE THREE ****
C CENTRAL CALCULATION.
C REQUIRES X1,X3,X5,FX1,FX3,FX5,EST,ADIFF.
C
300 CONTINUE
X2=0.5*(X1+X3)
X4=0.5*(X3+X5)
FX2=F(X2)
FX4=F(X4)
NO=NO+2
EST1 = FX1 + 4.0*FX2 + FX3
EST2 = FX3 + 4.0*FX4 + FX5
ADIFF1 = ADIFF
DIFF = EST + EST - EST1 - EST2
IF (LEV-30) 305,800,800
305 ADIFF = ABS(DIFF)
CRIT = ADIFF - CEPS
IF (CRIT) 700,700,400
C
C END OF CENTRAL LOOP
C NEXT STAGE IS STAGE FOUR IN CASE OF NO NATURAL CONVERGENCE
C NEXT STAGE IS STAGE SEVEN IN CASE OF NATURAL CONVERGENCE
C
C **** STAGE FOUR ****
C NO NATURAL CONVERGENCE. A COMPLEX SEQUENCE OF INSTRUCTIONS
C FOLLOWS WHICH ASSIGNS CONVERGENCE AND/OR ALTERS TOLERANCE
C LEVEL IN UPWARD DIRECTION IF THERE ARE INDICATIONS
C OF ROUND OFF ERROR.
C
400 CONTINUE
IF (ADIFF1-ADIFF) 410,410,500
C
C IN A NORMAL RUN WITH NO ROUND OFF ERROR PROBLEM, ADIFF1 IS
C GREATER THAN ADIFF AND THE REST OF STAGE FOUR IS OMITTED.
C
410 IF (LEV-5) 500,415,415
415 EFACT =EFACT + CEPS *(X1-XZERO)*FACERR
XZERO = X1
FACERR = 15.0
C
C THE REST OF STAGE FOUR DEALS WITH UPWARD ADJUSTMENT OF TOLERANCE
C (CEPS) BECAUSE OF SUSPECTED ROUND OFF ERROR TROUBLE.
C
IF (ADIFF-2.0*CEPS) 420,420,425
C
C SMALL JUMP IN CEPS. ASSIGN CONVERGENCE.
C
420 CEPS = ADIFF
LEVTAG = 0
GO TO 780
425 IF (ADIFF1-ADIFF) 435,430,435
C
C LARGE JUMP IN CEPS
C
430 CEPS = ADIFF
GO TO 445
C
C FACTOR TWO JUMP IN CEPS
C
435 CEPS = 2.0 * CEPS
IF (LEVTAG-3) 440,445,445
440 LEVTAG = 2
445 QCEPS = 0.25*CEPS
C
C **** STAGE FIVE ****
C NO ACTUAL CONVERGENCE.
C STORE RIGHT HAND ELEMENTS
C
500 CONTINUE
NIM= 2* NIM
LEV = LEV + 1
ESTST(LEV) = EST2
X3ST(LEV) = X4
X5ST(LEV) = X5
FX3ST(LEV) = FX4
FX5ST(LEV) = FX5
PREDIF(LEV) = ADIFF
C
C **** STAGE SIX ****
C SET UP QUANTITIES FOR CENTRAL CALCULATION.
C
C READY TO GO AHEAD AT LEVEL LOWER WITH LEFT HAND ELEMENTS
C X1 AND FX1 ARE THE SAME AS BEFORE
C
X5=X3
X3=X2
FX5=FX3
FX3=FX2
EST = EST1
GO TO 300
C
C **** STAGE SEVEN ****
C NATURAL CONVERGENCE IN PREVIOUS INTERVAL. THE FOLLOWING
C COMPLEX SEQUENCE CHECKS PRIMARILY THAT TOLERANCE LEVEL IS
C NOT TOO HIGH. UNDER CERTAIN CIRCUMSTANCES NON CONVERGENCE
C IS ASSIGNED AND/OR TOLERANCE LEVEL IS RE-SET.
C
700 CONTINUE
C
C CHECK THAT IT WAS NOT LEVEL ZERO INTERVAL. IF SO ASSIGN
C NON CONVERGENCE
C
IF (LEV) 400,400,705
C
C LEVTAG = -1 CEPS = CEPSF, ITS ORIGINAL VALUE.
C LEVTAG = 0 CEPS IS GREATER THAN CEPSF. REGULAR SITUATION.
C LEVTAG = 2 CEPS IS GREATER THAN CEPSF. CEPS PREVIOUSLY ASKED
C FOR A BIG JUMP, BUT DID NOT GET ONE.
C LEVTAG = 3 CEPS IS GREATER THAN CEPSF. CEPS PREVIOUSLY HAD
C A BIG JUMP.
C
705 IF (LEVTAG) 800,710,710
C
C IN A NORMAL RUN WITH NO ROUND OFF ERROR PROBLEM, LEVTAG = -1
C AND THE REST OF STAGE SEVEN IS OMITTED.
C
710 CEPST = 15.0*CEPS
C
C CEPST HERE IS FACERR*CURRENT VALUE OF CEPS
C
IF (CRIT) 715,800,800
715 IF (LEVTAG-2) 720,740,750
C
C LEVTAG = 0
C
720 IF (ADIFF) 800,800,725
725 IF (ADIFF-QCEPS) 730,800,800
730 IF (ADIFF-CEPSF) 770,770,735
735 LEVTAG = 0
CEPS = ADIFF
EFACT = EFACT + CEPST*(X1-XZERO)
XZERO=X1
GO TO 445
C
C LEVTAG = 2
C
740 LEVTAG = 0
IF (ADIFF) 765,765,725
C
C LEVTAG = 3
C
750 LEVTAG = 0
IF (ADIFF) 775,775,730
765 CEPS = ADIFF1
GO TO 775
770 LEVTAG = -1
FACERR = 1.0
CEPS = CEPSF
775 EFACT = EFACT + CEPST*(X1-XZERO)
XZERO = X1
780 CONTINUE
QCEPS = 0.25*CEPS
C
C **** STAGE EIGHT ****
C ACTUAL CONVERGENCE IN PREVIOUS INTERVAL. INCREMENTS ADDED INTO
C RUNNING SUMS
C
C ADD INTO SUM AND SIM
C
800 CONTINUE
SUM= SUM+ (EST1+EST2)*(X5-X1)
IF (LEVTAG) 805,810,810
C
C WE ADD INTO SIM ONLY IF WE ARE CLEAR OF ROUND OFF LEVEL
C
805 SIM=SIM + DIFF*(X5-X1)
810 CONTINUE
C
C **** STAGE NINE ****
C SORT OUT WHICH LEVEL TO GO TO. THIS INVOLVES NIM NUMBERING
C SYSTEM DESCRIBED BEFORE STAGE ONE.
C
905 NUM = NIM/2
NOM= NIM - 2*NUM
IF (NOM) 910,915,910
910 NIM = NUM
LEV = LEV - 1
GO TO 905
915 NIM = NIM + 1
C
C NEW LEVEL IS SET. IF LEV = 0 WE HAVE FINISHED
C
IF (LEV) 1100,1100,1000
C
C **** STAGE TEN ****
C SET UP QUANTITIES FOR CENTRAL CALCULATION.
C
1000 CONTINUE
X1=X5
FX1=FX5
X3=X3ST(LEV)
X5=X5ST(LEV)
FX3=FX3ST(LEV)
FX5=FX5ST(LEV)
EST=ESTST(LEV)
ADIFF = PREDIF(LEV)
GO TO 300
C
C **** STAGE ELEVENT ****
C
1100 CONTINUE
EFACT = EFACT + CEPS * (BIG-XZERO)*FACERR
RUM = EFACT/180.0
THIRD = SUM/12.0
FIFTH = -SIM/180.0
SQUANK = THIRD + FIFTH
RETURN
C
C END OF SQUANK
C
END
C---------------X, A ARE RETURNED. C, D ARE INPUT.
SUBROUTINE GLO16(X,A,C,D)
DOUBLE PRECISION DMC,DPC,C,D,X(16),A(16),XX(8),AA(8)
DATA (XX(I),AA(I),I=1,8)/
1.9894009349916499325961541734D0 ,.2715245941175409485178057245D-1,
1.9445750230732325760779884155D0 ,.6225352393864789286284383699D-1,
1.8656312023878317438804678977D0 ,.9515851168249278480992510760D-1,
1.7554044083550030338951011948D0 ,.1246289712555338720524762821D0 ,
1.6178762444026437484466717640D0 ,.1495959888165767320815017305D0 ,
1.4580167776572273863424194429D0 ,.1691565193950025381893120790D0 ,
1.2816035507792589132304605014D0 ,.1826034150449235888667636679D0 ,
1.9501250983763744018531933542D-1,.1894506104550684962853967232D0/
DMC = .5D0*(D-C)
DPC = .5D0*(D+C)
DO 2 I=1,8
NI = 17 - I
X(I) = -DMC*XX(I) + DPC
X(NI)= DMC*XX(I) + DPC
A(I) = DMC*AA(I)
2 A(NI)= DMC*AA(I)
RETURN
END