← Back to branch summary

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

« back to all changes in this revision

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