← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/dpchsw.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/dpchsw.f
 
1
*DECK DPCHSW
 
2
      SUBROUTINE DPCHSW (DFMAX, IEXTRM, D1, D2, H, SLOPE, IERR)
 
3
C***BEGIN PROLOGUE  DPCHSW
 
4
C***SUBSIDIARY
 
5
C***PURPOSE  Limits excursion from data for DPCHCS
 
6
C***LIBRARY   SLATEC (PCHIP)
 
7
C***TYPE      DOUBLE PRECISION (PCHSW-S, DPCHSW-D)
 
8
C***AUTHOR  Fritsch, F. N., (LLNL)
 
9
C***DESCRIPTION
 
10
C
 
11
C         DPCHSW:  DPCHCS Switch Excursion Limiter.
 
12
C
 
13
C     Called by  DPCHCS  to adjust D1 and D2 if necessary to insure that
 
14
C     the extremum on this interval is not further than DFMAX from the
 
15
C     extreme data value.
 
16
C
 
17
C ----------------------------------------------------------------------
 
18
C
 
19
C  Calling sequence:
 
20
C
 
21
C        INTEGER  IEXTRM, IERR
 
22
C        DOUBLE PRECISION  DFMAX, D1, D2, H, SLOPE
 
23
C
 
24
C        CALL  DPCHSW (DFMAX, IEXTRM, D1, D2, H, SLOPE, IERR)
 
25
C
 
26
C   Parameters:
 
27
C
 
28
C     DFMAX -- (input) maximum allowed difference between F(IEXTRM) and
 
29
C           the cubic determined by derivative values D1,D2.  (assumes
 
30
C           DFMAX.GT.0.)
 
31
C
 
32
C     IEXTRM -- (input) index of the extreme data value.  (assumes
 
33
C           IEXTRM = 1 or 2 .  Any value .NE.1 is treated as 2.)
 
34
C
 
35
C     D1,D2 -- (input) derivative values at the ends of the interval.
 
36
C           (Assumes D1*D2 .LE. 0.)
 
37
C          (output) may be modified if necessary to meet the restriction
 
38
C           imposed by DFMAX.
 
39
C
 
40
C     H -- (input) interval length.  (Assumes  H.GT.0.)
 
41
C
 
42
C     SLOPE -- (input) data slope on the interval.
 
43
C
 
44
C     IERR -- (output) error flag.  should be zero.
 
45
C           If IERR=-1, assumption on D1 and D2 is not satisfied.
 
46
C           If IERR=-2, quadratic equation locating extremum has
 
47
C                       negative discriminant (should never occur).
 
48
C
 
49
C    -------
 
50
C    WARNING:  This routine does no validity-checking of arguments.
 
51
C    -------
 
52
C
 
53
C  Fortran intrinsics used:  ABS, SIGN, SQRT.
 
54
C
 
55
C***SEE ALSO  DPCHCS
 
56
C***ROUTINES CALLED  D1MACH, XERMSG
 
57
C***REVISION HISTORY  (YYMMDD)
 
58
C   820218  DATE WRITTEN
 
59
C   820805  Converted to SLATEC library version.
 
60
C   870707  Corrected XERROR calls for d.p. name(s).
 
61
C   870707  Replaced DATA statement for SMALL with a use of D1MACH.
 
62
C   870813  Minor cosmetic changes.
 
63
C   890206  Corrected XERROR calls.
 
64
C   890411  1. Added SAVE statements (Vers. 3.2).
 
65
C           2. Added DOUBLE PRECISION declaration for D1MACH.
 
66
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
67
C   890531  REVISION DATE from Version 3.2
 
68
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
69
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
 
70
C   900328  Added TYPE section.  (WRB)
 
71
C   910408  Updated AUTHOR and DATE WRITTEN sections in prologue.  (WRB)
 
72
C   920526  Eliminated possible divide by zero problem.  (FNF)
 
73
C   930503  Improved purpose.  (FNF)
 
74
C***END PROLOGUE  DPCHSW
 
75
C
 
76
C**End
 
77
C
 
78
C  DECLARE ARGUMENTS.
 
79
C
 
80
      INTEGER  IEXTRM, IERR
 
81
      DOUBLE PRECISION  DFMAX, D1, D2, H, SLOPE
 
82
C
 
83
C  DECLARE LOCAL VARIABLES.
 
84
C
 
85
      DOUBLE PRECISION  CP, FACT, HPHI, LAMBDA, NU, ONE, PHI, RADCAL,
 
86
     *                  RHO, SIGMA, SMALL, THAT, THIRD, THREE, TWO, ZERO
 
87
      SAVE ZERO, ONE, TWO, THREE, FACT
 
88
      SAVE THIRD
 
89
      DOUBLE PRECISION  D1MACH
 
90
C
 
91
      DATA  ZERO /0.D0/,  ONE /1.D0/,  TWO /2.D0/, THREE /3.D0/,
 
92
     *      FACT /100.D0/
 
93
C        THIRD SHOULD BE SLIGHTLY LESS THAN 1/3.
 
94
      DATA  THIRD /0.33333D0/
 
95
C
 
96
C  NOTATION AND GENERAL REMARKS.
 
97
C
 
98
C     RHO IS THE RATIO OF THE DATA SLOPE TO THE DERIVATIVE BEING TESTED.
 
99
C     LAMBDA IS THE RATIO OF D2 TO D1.
 
100
C     THAT = T-HAT(RHO) IS THE NORMALIZED LOCATION OF THE EXTREMUM.
 
101
C     PHI IS THE NORMALIZED VALUE OF P(X)-F1 AT X = XHAT = X-HAT(RHO),
 
102
C           WHERE  THAT = (XHAT - X1)/H .
 
103
C        THAT IS, P(XHAT)-F1 = D*H*PHI,  WHERE D=D1 OR D2.
 
104
C     SIMILARLY,  P(XHAT)-F2 = D*H*(PHI-RHO) .
 
105
C
 
106
C      SMALL SHOULD BE A FEW ORDERS OF MAGNITUDE GREATER THAN MACHEPS.
 
107
C***FIRST EXECUTABLE STATEMENT  DPCHSW
 
108
      SMALL = FACT*D1MACH(4)
 
109
C
 
110
C  DO MAIN CALCULATION.
 
111
C
 
112
      IF (D1 .EQ. ZERO)  THEN
 
113
C
 
114
C        SPECIAL CASE -- D1.EQ.ZERO .
 
115
C
 
116
C          IF D2 IS ALSO ZERO, THIS ROUTINE SHOULD NOT HAVE BEEN CALLED.
 
117
         IF (D2 .EQ. ZERO)  GO TO 5001
 
118
C
 
119
         RHO = SLOPE/D2
 
120
C          EXTREMUM IS OUTSIDE INTERVAL WHEN RHO .GE. 1/3 .
 
121
         IF (RHO .GE. THIRD)  GO TO 5000
 
122
         THAT = (TWO*(THREE*RHO-ONE)) / (THREE*(TWO*RHO-ONE))
 
123
         PHI = THAT**2 * ((THREE*RHO-ONE)/THREE)
 
124
C
 
125
C          CONVERT TO DISTANCE FROM F2 IF IEXTRM.NE.1 .
 
126
         IF (IEXTRM .NE. 1)  PHI = PHI - RHO
 
127
C
 
128
C          TEST FOR EXCEEDING LIMIT, AND ADJUST ACCORDINGLY.
 
129
         HPHI = H * ABS(PHI)
 
130
         IF (HPHI*ABS(D2) .GT. DFMAX)  THEN
 
131
C           AT THIS POINT, HPHI.GT.0, SO DIVIDE IS OK.
 
132
            D2 = SIGN (DFMAX/HPHI, D2)
 
133
         ENDIF
 
134
      ELSE
 
135
C
 
136
         RHO = SLOPE/D1
 
137
         LAMBDA = -D2/D1
 
138
         IF (D2 .EQ. ZERO)  THEN
 
139
C
 
140
C           SPECIAL CASE -- D2.EQ.ZERO .
 
141
C
 
142
C             EXTREMUM IS OUTSIDE INTERVAL WHEN RHO .GE. 1/3 .
 
143
            IF (RHO .GE. THIRD)  GO TO 5000
 
144
            CP = TWO - THREE*RHO
 
145
            NU = ONE - TWO*RHO
 
146
            THAT = ONE / (THREE*NU)
 
147
         ELSE
 
148
            IF (LAMBDA .LE. ZERO)  GO TO 5001
 
149
C
 
150
C           NORMAL CASE -- D1 AND D2 BOTH NONZERO, OPPOSITE SIGNS.
 
151
C
 
152
            NU = ONE - LAMBDA - TWO*RHO
 
153
            SIGMA = ONE - RHO
 
154
            CP = NU + SIGMA
 
155
            IF (ABS(NU) .GT. SMALL)  THEN
 
156
               RADCAL = (NU - (TWO*RHO+ONE))*NU + SIGMA**2
 
157
               IF (RADCAL .LT. ZERO)  GO TO 5002
 
158
               THAT = (CP - SQRT(RADCAL)) / (THREE*NU)
 
159
            ELSE
 
160
               THAT = ONE/(TWO*SIGMA)
 
161
            ENDIF
 
162
         ENDIF
 
163
         PHI = THAT*((NU*THAT - CP)*THAT + ONE)
 
164
C
 
165
C          CONVERT TO DISTANCE FROM F2 IF IEXTRM.NE.1 .
 
166
         IF (IEXTRM .NE. 1)  PHI = PHI - RHO
 
167
C
 
168
C          TEST FOR EXCEEDING LIMIT, AND ADJUST ACCORDINGLY.
 
169
         HPHI = H * ABS(PHI)
 
170
         IF (HPHI*ABS(D1) .GT. DFMAX)  THEN
 
171
C           AT THIS POINT, HPHI.GT.0, SO DIVIDE IS OK.
 
172
            D1 = SIGN (DFMAX/HPHI, D1)
 
173
            D2 = -LAMBDA*D1
 
174
         ENDIF
 
175
      ENDIF
 
176
C
 
177
C  NORMAL RETURN.
 
178
C
 
179
 5000 CONTINUE
 
180
      IERR = 0
 
181
      RETURN
 
182
C
 
183
C  ERROR RETURNS.
 
184
C
 
185
 5001 CONTINUE
 
186
C     D1 AND D2 BOTH ZERO, OR BOTH NONZERO AND SAME SIGN.
 
187
      IERR = -1
 
188
      CALL XERMSG ('SLATEC', 'DPCHSW', 'D1 AND/OR D2 INVALID', IERR, 1)
 
189
      RETURN
 
190
C
 
191
 5002 CONTINUE
 
192
C     NEGATIVE VALUE OF RADICAL (SHOULD NEVER OCCUR).
 
193
      IERR = -2
 
194
      CALL XERMSG ('SLATEC', 'DPCHSW', 'NEGATIVE RADICAL', IERR, 1)
 
195
      RETURN
 
196
C------------- LAST LINE OF DPCHSW FOLLOWS -----------------------------
 
197
      END