1 2 SUBROUTINE DQAWS(F,A,B,ALFA,BETA,INTEGR,EPSABS,EPSREL,RESULT, 3 1 ABSERR,NEVAL,IER,LIMIT,LENW,LAST,IWORK,WORK) 4C***BEGIN PROLOGUE DQAWS 5C***DATE WRITTEN 800101 (YYMMDD) 6C***REVISION DATE 830518 (YYMMDD) 7C***CATEGORY NO. H2A2A1 8C***KEYWORDS ALGEBRAICO-LOGARITHMIC END-POINT SINGULARITIES, 9C AUTOMATIC INTEGRATOR,CLENSHAW-CURTIS,GLOBALLY ADAPTIVE, 10C SPECIAL-PURPOSE 11C***AUTHOR PIESSENS, ROBERT, APPLIED MATH. AND PROGR. DIV. - 12C K. U. LEUVEN 13C DE DONCKER, ELISE, APPLIED MATH. AND PROGR. DIV. - 14C K. U. LEUVEN 15C***PURPOSE The routine calculates an approximation result to a given 16C definite integral I = Integral of F*W over (A,B), 17C (where W shows a singular behaviour at the end points 18C see parameter INTEGR). 19C Hopefully satisfying following claim for accuracy 20C ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)). 21C***DESCRIPTION 22C 23C Integration of functions having algebraico-logarithmic 24C end point singularities 25C Standard fortran subroutine 26C Double precision version 27C 28C PARAMETERS 29C ON ENTRY 30C F - Double precision 31C Function subprogram defining the integrand 32C function F(X). The actual name for F needs to be 33C declared E X T E R N A L in the driver program. 34C 35C A - Double precision 36C Lower limit of integration 37C 38C B - Double precision 39C Upper limit of integration, B.GT.A 40C If B.LE.A, the routine will end with IER = 6. 41C 42C ALFA - Double precision 43C Parameter in the integrand function, ALFA.GT.(-1) 44C If ALFA.LE.(-1), the routine will end with 45C IER = 6. 46C 47C BETA - Double precision 48C Parameter in the integrand function, BETA.GT.(-1) 49C If BETA.LE.(-1), the routine will end with 50C IER = 6. 51C 52C INTEGR - Integer 53C Indicates which WEIGHT function is to be used 54C = 1 (X-A)**ALFA*(B-X)**BETA 55C = 2 (X-A)**ALFA*(B-X)**BETA*LOG(X-A) 56C = 3 (X-A)**ALFA*(B-X)**BETA*LOG(B-X) 57C = 4 (X-A)**ALFA*(B-X)**BETA*LOG(X-A)*LOG(B-X) 58C If INTEGR.LT.1 or INTEGR.GT.4, the routine 59C will end with IER = 6. 60C 61C EPSABS - Double precision 62C Absolute accuracy requested 63C EPSREL - Double precision 64C Relative accuracy requested 65C If EPSABS.LE.0 66C and EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28), 67C the routine will end with IER = 6. 68C 69C ON RETURN 70C RESULT - Double precision 71C Approximation to the integral 72C 73C ABSERR - Double precision 74C Estimate of the modulus of the absolute error, 75C Which should equal or exceed ABS(I-RESULT) 76C 77C NEVAL - Integer 78C Number of integrand evaluations 79C 80C IER - Integer 81C IER = 0 Normal and reliable termination of the 82C routine. It is assumed that the requested 83C accuracy has been achieved. 84C IER.GT.0 Abnormal termination of the routine 85C The estimates for the integral and error 86C are less reliable. It is assumed that the 87C requested accuracy has not been achieved. 88C ERROR MESSAGES 89C IER = 1 Maximum number of subdivisions allowed 90C has been achieved. One can allow more 91C subdivisions by increasing the value of 92C LIMIT (and taking the according dimension 93C adjustments into account). However, if 94C this yields no improvement it is advised 95C to analyze the integrand, in order to 96C determine the integration difficulties 97C which prevent the requested tolerance from 98C being achieved. In case of a jump 99C discontinuity or a local singularity 100C of algebraico-logarithmic type at one or 101C more interior points of the integration 102C range, one should proceed by splitting up 103C the interval at these points and calling 104C the integrator on the subranges. 105C = 2 The occurrence of roundoff error is 106C detected, which prevents the requested 107C tolerance from being achieved. 108C = 3 Extremely bad integrand behaviour occurs 109C at some points of the integration 110C interval. 111C = 6 The input is invalid, because 112C B.LE.A or ALFA.LE.(-1) or BETA.LE.(-1) or 113C or INTEGR.LT.1 or INTEGR.GT.4 or 114C (EPSABS.LE.0 and 115C EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28)) 116C or LIMIT.LT.2 or LENW.LT.LIMIT*4. 117C RESULT, ABSERR, NEVAL, LAST are set to 118C zero. Except when LENW or LIMIT is invalid 119C IWORK(1), WORK(LIMIT*2+1) and 120C WORK(LIMIT*3+1) are set to zero, WORK(1) 121C is set to A and WORK(LIMIT+1) to B. 122C 123C DIMENSIONING PARAMETERS 124C LIMIT - Integer 125C Dimensioning parameter for IWORK 126C LIMIT determines the maximum number of 127C subintervals in the partition of the given 128C integration interval (A,B), LIMIT.GE.2. 129C If LIMIT.LT.2, the routine will end with IER = 6. 130C 131C LENW - Integer 132C Dimensioning parameter for WORK 133C LENW must be at least LIMIT*4. 134C If LENW.LT.LIMIT*4, the routine will end 135C with IER = 6. 136C 137C LAST - Integer 138C On return, LAST equals the number of 139C subintervals produced in the subdivision process, 140C which determines the significant number of 141C elements actually in the WORK ARRAYS. 142C 143C WORK ARRAYS 144C IWORK - Integer 145C Vector of dimension LIMIT, the first K 146C elements of which contain pointers 147C to the error estimates over the subintervals, 148C such that WORK(LIMIT*3+IWORK(1)), ..., 149C WORK(LIMIT*3+IWORK(K)) form a decreasing 150C sequence with K = LAST if LAST.LE.(LIMIT/2+2), 151C and K = LIMIT+1-LAST otherwise 152C 153C WORK - Double precision 154C Vector of dimension LENW 155C On return 156C WORK(1), ..., WORK(LAST) contain the left 157C end points of the subintervals in the 158C partition of (A,B), 159C WORK(LIMIT+1), ..., WORK(LIMIT+LAST) contain 160C the right end points, 161C WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) 162C contain the integral approximations over 163C the subintervals, 164C WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST) 165C contain the error estimates. 166C***REFERENCES (NONE) 167C***ROUTINES CALLED DQAWSE,XERROR 168C***END PROLOGUE DQAWS 169 170 171