* * $Id: chebqu.F,v 1.1.1.1 1996/04/01 15:03:13 mclareni Exp $ * * $Log: chebqu.F,v $ * Revision 1.1.1.1 1996/04/01 15:03:13 mclareni * Mathlib gen * * #include "sys/CERNLIB_machine.h" #include "_gen/pilot.h" DOUBLE PRECISION FUNCTION CHEBQU(A,B,EPSIN,EPSOUT,JOP,FUNC) C C BASIC VERSION AUGUST 1971. C PURPOSE = TO OBTAIN THE VALUE OF A REAL AND DEFINITE INTEGRAL AP- C PROXIMATING THE INTEGRAND F(X) BY C F(X)=.5*C(1)+C(2)*T(X,1)+............+C(M)*T(X,M-1) C WHERE T(X,K) IS THE CHEBYSHEV POLYNOMIAL OF ORDER K. C MODIFIED CLENSHAW/CURTIS METHOD REWRITTEN AS AN ORDINARY C QUADRATURE FORMULA. C THE BASIC THEORY FOR THE MODIFIED CLENSHAW-CURTIS ALGORITHM IS GI- C VEN IN T.HAAVIE , ON A MODIFICATION OF THE C355081W 3URT90 QUA- C DRATURE FORMULA , BIT 9(1969),338-350. C THE FORMULAS RELATED TO REWRITING THE ALGORITHM AS AN ORDINARY C QUADRATURE FORMULA IS GIVEN IN , T.HAAVIE , SOME METHODS FOR AUTO- C MATIC INTEGRATION AND THEIR IMPLEMENTATION ON THE CERN 65/6600 C COMPUTERS , (YELLOW REPORT) CERN 71- , GENEVA 1971. C PARAMETERS C A = LOWER BOUNDARY C B = UPPER BOUNDARY C EPSIN = ACCURACY REQUIRED FOR THE APPROXINATION C EPSOUT = IMPROVED ERROR ESTIMATE FOR THE APPROXIMATION C JOP = OPTION PARAMETER , JOP=0 , NO PRINTING OF INTERMEDIATE C CALCULATIONS C JOP=1 , PRINT INTERMEDIATE CALCULA- C TIONS C FUNC = FUNCTION ROUTINE FOR THE FUNCTION FUNC(X).TO BE DE- C CLARED EXTERNAL IN THE CALLING ROUTINE C PARAMETERS IN COMMON BLOCK / CHEINT / C TEND = UPPER BOUND FOR VALUE OF INTEGRAL,I.E. TEND=TN*=TN+DN. C UMID = LOWER BOUND FOR VALUE OF INTEGRAL,I.E. UMID=UN*=UN+DN. C DELN = DN , DEFINITION FOLLOWS FROM DESCRIPTION OF TEND,UMID. C N = THE NUMBER OF INTEGRAND VALUES USED IN THE CALCULATION EXTERNAL FUNC DOUBLE PRECISION A,B,EPSIN,EPSOUT,FUNC C VARIABLES DEPENDING ON STEPSIZE DOUBLE PRECISION ALF,BET,RN,HNSTEP,TEND,UMID,WMEAN,DELN,TNEW,TSUM 1 , ACOFO , ACOFN , BCOFO , VALUE C OTHER VARIABLES USED DOUBLE PRECISION CONST1,CONST2,XPLUS,XMIN,BOUNDS , A0 , A1 , A2 , 1 A3 , A4 COMMON /CHEINT/ TEND , UMID , DELN , N C PARAMETERS IN THE COMMON BLOCK / DPCHEB / C (THESE PARAMETERS ARE PRESET IN COMMON BY MEANS OF A BLOCK DATA C SUBPROGRAM AND CONTAIN THE WEIGHTS AND ABSCISSAE FOR QUADRATURE C FORMULAS OF INCREASING ORDER.) C NUPPER = N CORRESPONDS TO 2**(N+1) SUB-INTERVALS FOR THE UNFOLDED C INTEGRAL.THE MAX NO OF FUNCTION EVALUATIONS THUS BEING C 1+2**(N+1).THE HIGHEST END-POINT APPROXIMATION IS THUS C USING 2**(N+1) INTERVALS WHILE THE HIGHEST MID-POINT AP- C PROXIMATION IS USING 2**N INTERVALS. C BETQJ = AN ARRAY CONTAINING THE WEIGHTS OF THE MID-POINT FORMU- C LAS 1 THROUGH NUPPER. C DIMENSION BETQJ(2**N-1) C KSIQJ = AN ARRAY CONTAINING THE ABSCISSAE OF THE MID-POINT FOR- C MULAS OF ORDER 1 THROUGH NUPPER. C DIMENSION KSIQJ(2**N-1) C ALFQJ = AN ARRAY CONTAINING THE WEIGHTS OF THE END-POINT FORMU- C LAS OF ORDER 2 THROUGH NUPPER+1. C DIMENSION ALFQJ(2**(N+1)+N-2) DOUBLE PRECISION BETQJ , KSIQJ , ALFQJ , FVAL C IN THIS VERSION OF CHEBQU NUPPER=9 , THE DIMENSIONS OF THE ARRAYS C IN THE COMMON BLOCK / DPCHEB / THUS BEING COMMON /DPCHEB/ BETQJ(511),KSIQJ(511),ALFQJ(1031),NUPPER C SOME INTERNAL PARAMETERS. C FVAL = AN INTERNAL ARRAY USED FOR STORING THE INTEGRAND VALU- C ES. C DIMENSION FVAL(2**N+1) C RNDERR = 1.D-28 , THE RELATIVE MACHINE ACCURACY IN DOUBLE PRECI- C SION (CDC-6000 SERIE). C IN THIS VERSION OF CHEBQU NUPPER=9 , THE DIMENSION OF FVAL THUS C BEING. DIMENSION FVAL(513) DOUBLE PRECISION ZERO,HALF,ONE,TWO,TWTHIR DATA ZERO , HALF , ONE , TWO / 0.D0 , .5D0 , 1.D0 , 2.D0 / DATA TWTHIR/.666666666666666666666666666667D0/ C SET INITIAL VALUES OF CONSTANTS DATA IBD/0/ IF(IBD.EQ.0) CALL D115BD(IBD) C THIS IS A CALL TO A SUBROUTINE WHICH ONLY GIVES VALUES TO THE C ARRAYS IN THE COMMON BLOCK DPCHEB BY DATA STATEMENTS C INTEGRATION INTERVAL PARAMETERS C*UL 1000 ALF=HALF*(B-A) ALF=HALF*(B-A) BET=HALF*(B+A) JNCR=2**(NUPPER-1) INCR=2*JNCR C C PARAMETERS FOR INTEGRATION STEPSIZE AND LOOPS RN =TWO N =2 NHALF =1 HNSTEP=ONE C C INITIAL CALCULATION FOR THE END-POINT APPROXIMATION CONST1=HALF*(FUNC(A)+FUNC(B)) CONST2=FUNC(BET) ACOFO=HALF*(CONST1+CONST2) ACOFN=HALF*(CONST1-CONST2) TEND =TWTHIR*(CONST1+TWO*CONST2) FVAL(1)=CONST2+CONST2 FVAL(INCR+1)=CONST1+CONST1 C C START ACTUAL CALCULATION C C--- Transform this DO-loop into a GOTO to avoid illegal jumps into it C C DO 1150 I=1,NUPPER I=0 1021 I=I+1 IF(I.GT.NUPPER) GOTO 1150 C C COMPUTE FUNCTION VALUES,CORRECTION TO MID-POINT APPROXIMATION AND C MID-POINT APPROXIMATION. BCOFO =ZERO UMID =ZERO DO 1030 J=1,NHALF INDEX= NHALF+J-1 XPLUS= ALF*KSIQJ(INDEX)+BET XMIN =-ALF*KSIQJ(INDEX)+BET VALUE= FUNC(XPLUS)+FUNC(XMIN) BCOFO= BCOFO+VALUE UMID = UMID+BETQJ(INDEX)*VALUE INDEX= (2*J-1)*JNCR+1 FVAL(INDEX)=VALUE 1030 CONTINUE C UMID = UMID-TWO*ACOFN/(RN**2-ONE) BCOFO= HALF*HNSTEP*BCOFO C C COMPUTE NEW END-POINT APPROXIMATION WHEN THE INTERVAL OF INTEGRA- C TION IS DIVIDED IN 2N EQUAL SUB INTERVALS WMEAN =HALF*(TEND+UMID) BOUNDS=HALF*(TEND-UMID) IPOINT=I+N-2 TSUM =ALFQJ(IPOINT)*FVAL(1) DO 1110 J=1,N INDEX=1+J*JNCR TSUM =TSUM+ALFQJ(IPOINT+J)*FVAL(INDEX) 1110 CONTINUE DELN =TSUM-WMEAN ACOFN =HALF*(ACOFO-BCOFO) ACOFO =ACOFO-ACOFN C ****************************************************************** C * * C PRINT INTERMEDIATE RESULTS IF WANTED IF (JOP.EQ.0) GO TO 1120 GO TO 2000 C * * C ****************************************************************** 1120 TNEW =TSUM EPSOUT=ABS(BOUNDS/TNEW) IF (EPSOUT.GT.EPSIN) GO TO 1130 C C REQUIRED ACCURACY OBTAINED OR THE MAXIMUM NUMBER OF FUNTION VAL- C UES USED WITHOUT OBTAINING THE REQUIRED ACCURACY. 1124 N=2*N+1 C*UL 1125 TEND=ALF*(TEND+DELN) TEND=ALF*(TEND+DELN) UMID=ALF*(UMID+DELN) DELN=ALF*DELN C*UL 1126 CHEBQU=ALF*TNEW C*UL 1128 RETURN CHEBQU=ALF*TNEW RETURN C C CALCULATION FINISHED. C 1130 IF (I.EQ.NUPPER) GO TO 1124 C C IF I=NUPPER THEN THE REQUIRED ACCURACY IS NOT OBTAINED. C PREPARE PARAMETERS FOR NEXT INTERVAL HALVING. TEND=TNEW NHALF =N N =2*N RN =TWO*RN INCR=JNCR JNCR=JNCR/2 HNSTEP=HALF*HNSTEP GOTO 1021 1150 CONTINUE C ****************************************************************** C * * C PRINT INTERMEDIATE RESULTS 2000 A0 = ALF*TEND A1 = ALF*WMEAN A2 = ALF*UMID CONST1 = ALF*BOUNDS CONST2 = ALF*DELN A3 = A0+CONST2 A4 = A2+CONST2 C IF (I.GT.1) GO TO 2010 WRITE(6,10) 2010 WRITE(6,20) N , A0 WRITE(6,30) A3 WRITE(6,30) A1 , CONST2 , CONST1 WRITE(6,30) A4 WRITE(6,30) A2 C GO TO 1120 C 10 FORMAT(6X,'N',17X,'TEND ', 1 /,24X,'TEND* ', 2 /,20X,'(TEND+UMID)/2',28X,'DELTAN',34X,'(TEND-UMID)/2' 3 , /,24X,'UMID* ', 4 /,24X,'UMID ') 20 FORMAT(/,4X,I3,2X,D36.29) 30 FORMAT(9X,D36.29,3X,D36.29,3X,D36.29) C * * C ****************************************************************** END