← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/qk15w.f

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-01-16 12:29:42 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130116122942-x86e42akjq31repw
Tags: 0.0.0+20130107.gitd9656f41-1
* New upstream snashot
* No longer try to rebuild helpdb.jl.
   + debian/rules: remove helpdb.jl from build-arch rule
   + debian/control: move back python-sphinx to Build-Depends-Indep
* debian/copyright: reflect upstream changes
* Add Build-Conflicts on libatlas3-base (makes linalg tests fail)
* debian/rules: replace obsolete USE_DEBIAN makeflag by a list of
  USE_SYSTEM_* flags
* debian/rules: on non-x86 systems, use libm instead of openlibm
* dpkg-buildflags.patch: remove patch, applied upstream
* Refreshed other patches

Show diffs side-by-side

added added

removed removed

Lines of Context:
deps/openlibm/slatec/qk15w.f
 
1
*DECK QK15W
 
2
      SUBROUTINE QK15W (F, W, P1, P2, P3, P4, KP, A, B, RESULT, ABSERR,
 
3
     +   RESABS, RESASC)
 
4
C***BEGIN PROLOGUE  QK15W
 
5
C***PURPOSE  To compute I = Integral of F*W over (A,B), with error
 
6
C                           estimate
 
7
C                       J = Integral of ABS(F*W) over (A,B)
 
8
C***LIBRARY   SLATEC (QUADPACK)
 
9
C***CATEGORY  H2A2A2
 
10
C***TYPE      SINGLE PRECISION (QK15W-S, DQK15W-D)
 
11
C***KEYWORDS  15-POINT GAUSS-KRONROD RULES, QUADPACK, QUADRATURE
 
12
C***AUTHOR  Piessens, Robert
 
13
C             Applied Mathematics and Programming Division
 
14
C             K. U. Leuven
 
15
C           de Doncker, Elise
 
16
C             Applied Mathematics and Programming Division
 
17
C             K. U. Leuven
 
18
C***DESCRIPTION
 
19
C
 
20
C           Integration rules
 
21
C           Standard fortran subroutine
 
22
C           Real version
 
23
C
 
24
C           PARAMETERS
 
25
C             ON ENTRY
 
26
C              F      - Real
 
27
C                       Function subprogram defining the integrand
 
28
C                       function F(X). The actual name for F needs to be
 
29
C                       declared E X T E R N A L in the driver program.
 
30
C
 
31
C              W      - Real
 
32
C                       Function subprogram defining the integrand
 
33
C                       WEIGHT function W(X). The actual name for W
 
34
C                       needs to be declared E X T E R N A L in the
 
35
C                       calling program.
 
36
C
 
37
C              P1, P2, P3, P4 - Real
 
38
C                       Parameters in the WEIGHT function
 
39
C
 
40
C              KP     - Integer
 
41
C                       Key for indicating the type of WEIGHT function
 
42
C
 
43
C              A      - Real
 
44
C                       Lower limit of integration
 
45
C
 
46
C              B      - Real
 
47
C                       Upper limit of integration
 
48
C
 
49
C            ON RETURN
 
50
C              RESULT - Real
 
51
C                       Approximation to the integral I
 
52
C                       RESULT is computed by applying the 15-point
 
53
C                       Kronrod rule (RESK) obtained by optimal addition
 
54
C                       of abscissae to the 7-point Gauss rule (RESG).
 
55
C
 
56
C              ABSERR - Real
 
57
C                       Estimate of the modulus of the absolute error,
 
58
C                       which should equal or exceed ABS(I-RESULT)
 
59
C
 
60
C              RESABS - Real
 
61
C                       Approximation to the integral of ABS(F)
 
62
C
 
63
C              RESASC - Real
 
64
C                       Approximation to the integral of ABS(F-I/(B-A))
 
65
C
 
66
C***REFERENCES  (NONE)
 
67
C***ROUTINES CALLED  R1MACH
 
68
C***REVISION HISTORY  (YYMMDD)
 
69
C   810101  DATE WRITTEN
 
70
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
71
C   890531  REVISION DATE from Version 3.2
 
72
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
73
C***END PROLOGUE  QK15W
 
74
C
 
75
      REAL A,ABSC,ABSC1,ABSC2,ABSERR,B,CENTR,DHLGTH,
 
76
     1  R1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,
 
77
     2  HLGTH,P1,P2,P3,P4,RESABS,RESASC,RESG,RESK,RESKH,RESULT,UFLOW,
 
78
     3  W,WG,WGK,XGK
 
79
      INTEGER J,JTW,JTWM1,KP
 
80
      EXTERNAL F, W
 
81
C
 
82
      DIMENSION FV1(7),FV2(7),XGK(8),WGK(8),WG(4)
 
83
C
 
84
C           THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1).
 
85
C           BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR
 
86
C           CORRESPONDING WEIGHTS ARE GIVEN.
 
87
C
 
88
C           XGK    - ABSCISSAE OF THE 15-POINT GAUSS-KRONROD RULE
 
89
C                    XGK(2), XGK(4), ... ABSCISSAE OF THE 7-POINT
 
90
C                    GAUSS RULE
 
91
C                    XGK(1), XGK(3), ... ABSCISSAE WHICH ARE OPTIMALLY
 
92
C                    ADDED TO THE 7-POINT GAUSS RULE
 
93
C
 
94
C           WGK    - WEIGHTS OF THE 15-POINT GAUSS-KRONROD RULE
 
95
C
 
96
C           WG     - WEIGHTS OF THE 7-POINT GAUSS RULE
 
97
C
 
98
      SAVE XGK, WGK, WG
 
99
      DATA XGK(1),XGK(2),XGK(3),XGK(4),XGK(5),XGK(6),XGK(7),
 
100
     1  XGK(8)/
 
101
     2     0.9914553711208126E+00,     0.9491079123427585E+00,
 
102
     3     0.8648644233597691E+00,     0.7415311855993944E+00,
 
103
     4     0.5860872354676911E+00,     0.4058451513773972E+00,
 
104
     5     0.2077849550078985E+00,     0.0000000000000000E+00/
 
105
C
 
106
      DATA WGK(1),WGK(2),WGK(3),WGK(4),WGK(5),WGK(6),WGK(7),
 
107
     1  WGK(8)/
 
108
     2     0.2293532201052922E-01,     0.6309209262997855E-01,
 
109
     3     0.1047900103222502E+00,     0.1406532597155259E+00,
 
110
     4     0.1690047266392679E+00,     0.1903505780647854E+00,
 
111
     5     0.2044329400752989E+00,     0.2094821410847278E+00/
 
112
C
 
113
      DATA WG(1),WG(2),WG(3),WG(4)/
 
114
     1     0.1294849661688697E+00,    0.2797053914892767E+00,
 
115
     2     0.3818300505051889E+00,    0.4179591836734694E+00/
 
116
C
 
117
C
 
118
C           LIST OF MAJOR VARIABLES
 
119
C           -----------------------
 
120
C
 
121
C           CENTR  - MID POINT OF THE INTERVAL
 
122
C           HLGTH  - HALF-LENGTH OF THE INTERVAL
 
123
C           ABSC*  - ABSCISSA
 
124
C           FVAL*  - FUNCTION VALUE
 
125
C           RESG   - RESULT OF THE 7-POINT GAUSS FORMULA
 
126
C           RESK   - RESULT OF THE 15-POINT KRONROD FORMULA
 
127
C           RESKH  - APPROXIMATION TO THE MEAN VALUE OF F*W OVER (A,B),
 
128
C                    I.E. TO I/(B-A)
 
129
C
 
130
C           MACHINE DEPENDENT CONSTANTS
 
131
C           ---------------------------
 
132
C
 
133
C           EPMACH IS THE LARGEST RELATIVE SPACING.
 
134
C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
 
135
C
 
136
C***FIRST EXECUTABLE STATEMENT  QK15W
 
137
      EPMACH = R1MACH(4)
 
138
      UFLOW = R1MACH(1)
 
139
C
 
140
      CENTR = 0.5E+00*(A+B)
 
141
      HLGTH = 0.5E+00*(B-A)
 
142
      DHLGTH = ABS(HLGTH)
 
143
C
 
144
C           COMPUTE THE 15-POINT KRONROD APPROXIMATION TO THE
 
145
C           INTEGRAL, AND ESTIMATE THE ERROR.
 
146
C
 
147
      FC = F(CENTR)*W(CENTR,P1,P2,P3,P4,KP)
 
148
      RESG = WG(4)*FC
 
149
      RESK = WGK(8)*FC
 
150
      RESABS = ABS(RESK)
 
151
      DO 10 J=1,3
 
152
        JTW = J*2
 
153
        ABSC = HLGTH*XGK(JTW)
 
154
        ABSC1 = CENTR-ABSC
 
155
        ABSC2 = CENTR+ABSC
 
156
        FVAL1 = F(ABSC1)*W(ABSC1,P1,P2,P3,P4,KP)
 
157
        FVAL2 = F(ABSC2)*W(ABSC2,P1,P2,P3,P4,KP)
 
158
        FV1(JTW) = FVAL1
 
159
        FV2(JTW) = FVAL2
 
160
        FSUM = FVAL1+FVAL2
 
161
        RESG = RESG+WG(J)*FSUM
 
162
        RESK = RESK+WGK(JTW)*FSUM
 
163
        RESABS = RESABS+WGK(JTW)*(ABS(FVAL1)+ABS(FVAL2))
 
164
   10 CONTINUE
 
165
      DO 15 J=1,4
 
166
        JTWM1 = J*2-1
 
167
        ABSC = HLGTH*XGK(JTWM1)
 
168
        ABSC1 = CENTR-ABSC
 
169
        ABSC2 = CENTR+ABSC
 
170
        FVAL1 = F(ABSC1)*W(ABSC1,P1,P2,P3,P4,KP)
 
171
        FVAL2 = F(ABSC2)*W(ABSC2,P1,P2,P3,P4,KP)
 
172
        FV1(JTWM1) = FVAL1
 
173
        FV2(JTWM1) = FVAL2
 
174
        FSUM = FVAL1+FVAL2
 
175
        RESK = RESK+WGK(JTWM1)*FSUM
 
176
        RESABS = RESABS+WGK(JTWM1)*(ABS(FVAL1)+ABS(FVAL2))
 
177
   15 CONTINUE
 
178
      RESKH = RESK*0.5E+00
 
179
      RESASC = WGK(8)*ABS(FC-RESKH)
 
180
      DO 20 J=1,7
 
181
        RESASC = RESASC+WGK(J)*(ABS(FV1(J)-RESKH)+ABS(FV2(J)-RESKH))
 
182
   20 CONTINUE
 
183
      RESULT = RESK*HLGTH
 
184
      RESABS = RESABS*DHLGTH
 
185
      RESASC = RESASC*DHLGTH
 
186
      ABSERR = ABS((RESK-RESG)*HLGTH)
 
187
      IF(RESASC.NE.0.0E+00.AND.ABSERR.NE.0.0E+00)
 
188
     1  ABSERR = RESASC*MIN(0.1E+01,
 
189
     2  (0.2E+03*ABSERR/RESASC)**1.5E+00)
 
190
      IF(RESABS.GT.UFLOW/(0.5E+02*EPMACH)) ABSERR = MAX((EPMACH*
 
191
     1  0.5E+02)*RESABS,ABSERR)
 
192
      RETURN
 
193
      END