PROGRAM dcitest ! DESCRIPTION: ! This program tests the performance of DOUBLE-PRECISION codes which ! calculate the cosine-integral special function, defined as ! DCOSINT(x) = gamma + log(abs(x)) + ! integral from 0 to x of {[cos(t)-1]/t} dt ! The program compares the value of DCOSINT(x+h) with the Taylor ! series value about x. The derivatives of DCOSINT at x are ! generated by a standard recurrence relation. ! The program assumes the cosine-integral is calculated in a ! function with the name DCOSINT. Thus users should append their ! Ci code at the end of this file and edit the small DCOSINT ! function to include the name of their code. ! The program uses Cody's MACHAR subroutine to determine certain ! machine parameters. This can be avoided if values for the following ! parameters are available: EPS, IT, IBETA, MACHEP, XMIN, XMAX. ! Values for certain machine-compiler combinations are given in ! W.J. CODY Algorithm 665: MACHAR: A subroutine to dynamically ! determine machine parameters, ACM Trans. Math. Soft. 14 (1988) 303-311. ! The MACHAR code and a random generator code written by Cody are included ! in the file SICISUBD.FOR, though with the MACHAR code called DMACHR. ! To ensure that MACHAR and the argument purification section work as ! intended, this code should be compiled with any optimisation switched OFF. ! Users should note that MACHAR tests the machine arithmetic ! at its very limits. Thus it is possible that warning messages ! might appear about underflows and similar events. These DO ! NOT affect the test results. Such warnings can be avoided ! by explicitly declaring the necessary machine parameters. ! MODULES USED: ! This code calls the following modules: ! (a) DCIDER contained in this file ! (b) DCOSINT contained in the file DCOSINT.F ! (c) DMACHR, DREN and DTERMS contained in the file SICISUBD.F ! AUTHOR: Allan MacLeod ! Dept. of Mathematics and Statistics ! University of Paisley ! Scotland ! (e-mail: macl_ms0@paisley.ac.uk) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60) INTEGER :: i, ibeta, ipt, it, itest, ndec, nterms, nt1, nt2, numeq, & numla, numsm, tstart(8), tend(8) REAL (dp) :: cider(0:100), ait, al, albeta, au, cix, ciy, delta, div, erloss, & ERR, hint, maxerr, maxpt, relerr, rmserr, sum, temp, x, y, z REAL :: dren, time_taken ! Define constants REAL (dp), PARAMETER :: alow(5) = (/ 0.25_dp, 2.0_dp, 4.0_dp, 8.0_dp, & 13.5_dp /) REAL (dp), PARAMETER :: aup(5) = (/ 0.375_dp, 2.5_dp, 4.5_dp, 8.5_dp, & 14.0_dp /) REAL (dp), PARAMETER :: pt0625 = 0.0625_dp, sixten = 16.0_dp, half = 0.5_dp REAL (dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp, eps = EPSILON(one) REAL (dp), PARAMETER :: gamma = 0.57721566490153286061_dp REAL (dp), PARAMETER :: xmin = TINY(one), xmax = HUGE(one) ! IOUT = output channel number ! NPTS = number of test points per interval INTEGER, PARAMETER :: iout = 6, npts = 5000, ntest = 5 INTERFACE FUNCTION dcosint(xvalue) RESULT(fn_val) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14,60) REAL (dp), INTENT(IN) :: xvalue REAL (dp) :: fn_val END FUNCTION dcosint END INTERFACE ! DETERMINE MACHINE PARAMETERS ibeta = RADIX( one ) it = DIGITS( one ) ait = DBLE(it) albeta = LOG(DBLE(ibeta)) ndec = ABS( INT( LOG10( DBLE(ibeta) * eps ) - half ) ) CALL DATE_AND_TIME( VALUES=tstart ) ! START TESTS DO itest = 1 , ntest al = alow(itest) au = aup(itest) hint = ( au - al ) / DBLE(npts) ! DETERMINE HOW MANY DERIVATIVES ARE NEEDED CALL dterms(ABS(al), eps, nt1) CALL dterms(ABS(au), eps, nt2) nterms = nt1 IF ( nt2 > nt1 ) nterms = nt2 ! SET UP STATISTICS COLLECTION maxerr = zero rmserr = zero maxpt = zero numsm = 0 numla = 0 delta = pt0625 IF ( itest == 2 .OR. itest == 4 ) delta = -pt0625 DO ipt = 1 , npts CALL RANDOM_NUMBER( dren ) x = al + hint * dren ! PURIFY ARGUMENT z = sixten * x temp = z + x x = temp - z y = x + delta ! CALCULATE CI(Y) AND CI(X+DELTA) ciy = dcosint(y) cix = dcosint(x) CALL dcider(x, nterms, ndec, cider) sum = cider(nterms) div = DBLE(nterms) + one DO i = 0 , nterms-1 sum = sum * delta / div + cider(nterms-1-i) div = div - one END DO sum = delta * sum ! DETERMINE ERROR AND UPDATE STATS ERR = ( ciy - cix ) - sum IF ( ERR < zero ) numsm = numsm + 1 IF ( ERR > zero ) numla = numla + 1 relerr = ABS(ERR/ciy) rmserr = rmserr + relerr * relerr IF ( relerr > maxerr ) THEN maxerr = relerr maxpt = x END IF al = al + hint END DO ! OUTPUT RESULTS FOR EACH TEST INTERVAL rmserr = SQRT(rmserr/DBLE(npts)) WRITE(iout, 1000) WRITE(iout, 1010) alow(itest), aup(itest), npts IF ( rmserr > zero ) THEN erloss = ait + LOG(rmserr) / albeta ELSE erloss = zero END IF IF ( erloss < zero ) erloss = zero WRITE(iout, 1020) rmserr, erloss WRITE(iout, 1030) maxerr, maxpt IF ( maxerr > zero ) THEN erloss = ait + LOG(maxerr) / albeta ELSE erloss = zero END IF IF ( erloss < zero ) erloss = zero WRITE(iout, 1040) erloss numeq = npts - ( numsm + numla ) WRITE(iout, 1050) numsm, numeq, numla END DO ! PRINT VALUES OF CI(XMAX) = ZERO AND CI(XMIN) = LN(XMIN) WRITE(iout, 1100) WRITE(iout, 1200) y = dcosint(xmax) WRITE(iout, 1210) y x = LOG(xmin) + gamma WRITE(iout, 1220) x y = dcosint(xmin) WRITE(iout, 1230) y CALL DATE_AND_TIME( VALUES=tend ) time_taken = 3600.*(tend(5) - tstart(5)) + 60.*(tend(6) - tstart(6)) + & (tend(7) - tstart(7)) + 0.001*(tend(8) - tstart(8)) WRITE(iout, 1240) time_taken STOP ! FORMAT STATEMENTS 1000 FORMAT(/////t11, 'TEST OF CI(X+DELTA) AGAINST TAYLOR SERIES') 1010 FORMAT(/t11, 'TEST CARRIED OUT ON ', f9.4,' TO ', f9.4, ' WITH ', & i5, ' ARGUMENTS') 1020 FORMAT(//t11, 'RMS ERROR = ', e15.6,' IS A LOSS OF', & f8.3, ' SIGNIFICANT DIGITS') 1030 FORMAT(//t11, 'MAX ERROR = ', e15.6, ' OCCURRED AT X = ', f10.4) 1040 FORMAT(/t11, 'THIS IS A LOSS OF ', f8.3, ' SIGNIFICANT DIGITS') 1050 FORMAT(///t11, i6, ' ARGUMENTS GAVE POSITIVE ERROR', / & t11, i6, ' ARGUMENTS GAVE NO ERROR', / & t11, i6, ' ARGUMENTS GAVE NEGATIVE ERROR') 1100 FORMAT(/////t11, 'SPECIAL TESTS') 1200 FORMAT(/////t11, 'CALCULATE CI(XMAX) - SHOULD RETURN ZERO') 1210 FORMAT(/t11, 'VALUE OF CI(XMAX) = ', e15.5) 1220 FORMAT(///t11, 'CALCULATE CI(XMIN) - SHOULD RETURN ', e15.5) 1230 FORMAT(/t11, 'VALUE OF CI(XMIN) = ', e15.5) 1240 FORMAT(/t11, 'Time taken = ', f8.2, 'sec.') CONTAINS SUBROUTINE dterms(x, eps, mval) ! This subroutine determines how many derivatives of Si or Ci ! are needed at the point X, so that the Taylor series with this no. ! of terms has a truncation error less than EPS in the ! relative sense. The no. of terms needed is stored in MVAL. ! INPUT PARAMETERS: ! X - DOUBLE PRECISION - The value at which the derivative is wanted. ! EPS - DOUBLE PRECISION - The relative error wanted in the Taylor series. ! OUTPUT PARAMETERS: ! MVAL - INTEGER - The required no. of terms ! AUTHOR: Allan MacLeod ! Dept. of Mathematics and Statistics ! University of Paisley ! High St. ! Paisley ! SCOTLAND ! PA1 2BE ! (e-mail: macl_ms0@paisley.ac.uk) REAL (dp), INTENT(IN) :: x REAL (dp), INTENT(IN) :: eps INTEGER, INTENT(OUT) :: mval ! Local variables REAL (dp) :: con1, con2, fmid, lower, mid, upper REAL (dp), PARAMETER :: zero = 0.0_dp, tol = 0.1_dp, one = 1.0_dp REAL (dp), PARAMETER :: two = 2.0_dp, four = 4.0_dp, sixten = 16.0_dp con1 = LOG ( four * EXP(x) / eps ) con2 = LOG ( sixten * x ) upper = con1 / con2 IF ( upper <= one ) THEN mval = 1 RETURN END IF lower = one 100 mid = ( upper + lower ) / two IF ( ABS(mid-lower) < tol ) THEN mval = INT(mid) + 1 GO TO 200 END IF fmid = LOG(mid) + mid * con2 - con1 IF ( fmid >= zero ) THEN upper = mid ELSE lower = mid END IF GO TO 100 200 RETURN END SUBROUTINE dterms SUBROUTINE dcider(x, nmax, ndigs, cider) ! This subroutine calculates the necessary number of derivatives of ! Ci by using the recurrence relation code of Gautschi and Klein ! from Comm. ACM., vol. 13 1970 pp53-54. ! INPUT PARAMETERS: ! X - DOUBLE PRECISION - Value of argument ! NMAX - INTEGER - Max. no of derivatives ! NDIGS - INTEGER - No. of decimal digits of accuracy wanted ! OUTPUT PARAMETERS: ! CIDER - DOUBLE PRECISION ARRAY - The values of the derivatives ! MODULES CALLED: ! This subroutine calls the function DTLNTI ! AUTHOR: Allan MacLeod ! Dept. of Mathematics and Statistics ! University of Paisley ! Scotland ! (e-mail: macl_ms0@paisley.ac.uk) REAL (dp), INTENT(IN) :: x INTEGER, INTENT(IN) :: nmax INTEGER, INTENT(IN) :: ndigs REAL (dp), INTENT(OUT) :: cider(0:) ! Local variables INTEGER :: minn, n, nlim, n0, n1 REAL (dp) :: d1, expone, q, sigma(4), s1, x1 REAL (dp), PARAMETER :: one = 1.0_dp, two = 2.0_dp, three = 3.0_dp, & ten = 10.0_dp x1 = ABS(x) sigma(1) = -SIN(x) sigma(2) = -COS(x) sigma(3) = -sigma(1) sigma(4) = -sigma(2) n0 = INT(x1) cider(0) = sigma(4) / x IF ( n0 < nmax ) THEN minn = n0 ELSE minn = nmax END IF IF ( x1 <= three ) THEN nlim = nmax ELSE nlim = minn END IF DO n = 1 , nlim cider(n) = (sigma(n-4*((n-1)/4)) - DBLE(n)*cider(n-1))/x END DO IF ( x1 > three .AND. n0 < nmax ) THEN expone = EXP(one) s1 = expone * x1 d1 = DBLE(ndigs) * LOG(ten) + LOG(two) n1 = INT(s1*dtlnti(d1/s1) - one) IF ( n1 < nmax ) n1 = nmax q = one / x DO n = 1 , n1+1 q = -DBLE(n) * q / x END DO DO n = n1 , n0+1 , -1 q = (sigma(n+1-4*(n/4)) - x*q) / DBLE(n+1) IF ( n <= nmax ) cider(n) = q END DO END IF RETURN END SUBROUTINE dcider FUNCTION dtlnti(y) RESULT(fn_val) ! This function solves the equation ! y = t ln(t) ! for a given value of y. It is not intended to be very accurate. ! N.B. Halley's method has been substituted for Newton's method. ! INPUT PARAMETER: ! Y - DOUBLE PRECISION - Given value of y. ! AUTHOR: Allan MacLeod ! Dept. of Mathematics and Statistics ! University of Paisley ! High St. ! Paisley ! SCOTLAND ! PA1 2BE ! (e-mail: macl_ms0@paisley.ac.uk) REAL (dp), INTENT(IN) :: y REAL (dp) :: fn_val ! Local variables REAL (dp) :: f0, fd, t0, t1 REAL (dp), PARAMETER :: one = 1.0_dp, tol = 0.05_dp, two = 2.0_dp t0 = ( y + y ) / ( one + LOG(y) ) DO f0 = t0 * LOG(t0) - y fd = one + LOG(t0) t1 = t0 - f0 / (fd - f0 / (two*t0*fd) ) IF ( ABS(t0-t1) > tol ) THEN t0 = t1 ELSE EXIT END IF END DO fn_val = t1 RETURN END FUNCTION dtlnti END PROGRAM dcitest