← Back to branch summary

~ubuntu-branches/ubuntu/lucid/pdl/lucid

« back to all changes in this revision

Viewing changes to Lib/Slatec/slatec/pchsw.f

  • Committer: Bazaar Package Importer
  • Author(s): Ben Gertzfield
  • Date: 2002-04-08 18:47:16 UTC
  • Revision ID: james.westby@ubuntu.com-20020408184716-0hf64dc96kin3htp
Tags: upstream-2.3.2
Import upstream version 2.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/Slatec/slatec/pchsw.f
 
1
*DECK PCHSW
 
2
      SUBROUTINE PCHSW (DFMAX, IEXTRM, D1, D2, H, SLOPE, IERR)
 
3
C***BEGIN PROLOGUE  PCHSW
 
4
C***SUBSIDIARY
 
5
C***PURPOSE  Limits excursion from data for PCHCS
 
6
C***LIBRARY   SLATEC (PCHIP)
 
7
C***TYPE      SINGLE PRECISION (PCHSW-S, DPCHSW-D)
 
8
C***AUTHOR  Fritsch, F. N., (LLNL)
 
9
C***DESCRIPTION
 
10
C
 
11
C         PCHSW:  PCHCS Switch Excursion Limiter.
 
12
C
 
13
C     Called by  PCHCS  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        REAL  DFMAX, D1, D2, H, SLOPE
 
23
C
 
24
C        CALL  PCHSW (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  PCHCS
 
56
C***ROUTINES CALLED  R1MACH, XERMSG
 
57
C***REVISION HISTORY  (YYMMDD)
 
58
C   820218  DATE WRITTEN
 
59
C   820805  Converted to SLATEC library version.
 
60
C   870707  Replaced DATA statement for SMALL with a use of R1MACH.
 
61
C   890411  1. Added SAVE statements (Vers. 3.2).
 
62
C           2. Added REAL R1MACH for consistency with D.P. version.
 
63
C   890411  REVISION DATE from Version 3.2
 
64
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
65
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
 
66
C   900328  Added TYPE section.  (WRB)
 
67
C   910408  Updated AUTHOR and DATE WRITTEN sections in prologue.  (WRB)
 
68
C   920526  Eliminated possible divide by zero problem.  (FNF)
 
69
C   930503  Improved purpose.  (FNF)
 
70
C***END PROLOGUE  PCHSW
 
71
C
 
72
C**End
 
73
C
 
74
C  DECLARE ARGUMENTS.
 
75
C
 
76
      INTEGER  IEXTRM, IERR
 
77
      REAL  DFMAX, D1, D2, H, SLOPE
 
78
C
 
79
C  DECLARE LOCAL VARIABLES.
 
80
C
 
81
      REAL  CP, FACT, HPHI, LAMBDA, NU, ONE, PHI, RADCAL, RHO, SIGMA,
 
82
     *      SMALL, THAT, THIRD, THREE, TWO, ZERO
 
83
      SAVE ZERO, ONE, TWO, THREE, FACT
 
84
      SAVE THIRD
 
85
      REAL R1MACH
 
86
C
 
87
      DATA  ZERO /0./,  ONE /1./,  TWO /2./,  THREE /3./, FACT /100./
 
88
C        THIRD SHOULD BE SLIGHTLY LESS THAN 1/3.
 
89
      DATA  THIRD /0.33333/
 
90
C
 
91
C  NOTATION AND GENERAL REMARKS.
 
92
C
 
93
C     RHO IS THE RATIO OF THE DATA SLOPE TO THE DERIVATIVE BEING TESTED.
 
94
C     LAMBDA IS THE RATIO OF D2 TO D1.
 
95
C     THAT = T-HAT(RHO) IS THE NORMALIZED LOCATION OF THE EXTREMUM.
 
96
C     PHI IS THE NORMALIZED VALUE OF P(X)-F1 AT X = XHAT = X-HAT(RHO),
 
97
C           WHERE  THAT = (XHAT - X1)/H .
 
98
C        THAT IS, P(XHAT)-F1 = D*H*PHI,  WHERE D=D1 OR D2.
 
99
C     SIMILARLY,  P(XHAT)-F2 = D*H*(PHI-RHO) .
 
100
C
 
101
C      SMALL SHOULD BE A FEW ORDERS OF MAGNITUDE GREATER THAN MACHEPS.
 
102
C***FIRST EXECUTABLE STATEMENT  PCHSW
 
103
      SMALL = FACT*R1MACH(4)
 
104
C
 
105
C  DO MAIN CALCULATION.
 
106
C
 
107
      IF (D1 .EQ. ZERO)  THEN
 
108
C
 
109
C        SPECIAL CASE -- D1.EQ.ZERO .
 
110
C
 
111
C          IF D2 IS ALSO ZERO, THIS ROUTINE SHOULD NOT HAVE BEEN CALLED.
 
112
         IF (D2 .EQ. ZERO)  GO TO 5001
 
113
C
 
114
         RHO = SLOPE/D2
 
115
C          EXTREMUM IS OUTSIDE INTERVAL WHEN RHO .GE. 1/3 .
 
116
         IF (RHO .GE. THIRD)  GO TO 5000
 
117
         THAT = (TWO*(THREE*RHO-ONE)) / (THREE*(TWO*RHO-ONE))
 
118
         PHI = THAT**2 * ((THREE*RHO-ONE)/THREE)
 
119
C
 
120
C          CONVERT TO DISTANCE FROM F2 IF IEXTRM.NE.1 .
 
121
         IF (IEXTRM .NE. 1)  PHI = PHI - RHO
 
122
C
 
123
C          TEST FOR EXCEEDING LIMIT, AND ADJUST ACCORDINGLY.
 
124
         HPHI = H * ABS(PHI)
 
125
         IF (HPHI*ABS(D2) .GT. DFMAX)  THEN
 
126
C           AT THIS POINT, HPHI.GT.0, SO DIVIDE IS OK.
 
127
            D2 = SIGN (DFMAX/HPHI, D2)
 
128
         ENDIF
 
129
      ELSE
 
130
C
 
131
         RHO = SLOPE/D1
 
132
         LAMBDA = -D2/D1
 
133
         IF (D2 .EQ. ZERO)  THEN
 
134
C
 
135
C           SPECIAL CASE -- D2.EQ.ZERO .
 
136
C
 
137
C             EXTREMUM IS OUTSIDE INTERVAL WHEN RHO .GE. 1/3 .
 
138
            IF (RHO .GE. THIRD)  GO TO 5000
 
139
            CP = TWO - THREE*RHO
 
140
            NU = ONE - TWO*RHO
 
141
            THAT = ONE / (THREE*NU)
 
142
         ELSE
 
143
            IF (LAMBDA .LE. ZERO)  GO TO 5001
 
144
C
 
145
C           NORMAL CASE -- D1 AND D2 BOTH NONZERO, OPPOSITE SIGNS.
 
146
C
 
147
            NU = ONE - LAMBDA - TWO*RHO
 
148
            SIGMA = ONE - RHO
 
149
            CP = NU + SIGMA
 
150
            IF (ABS(NU) .GT. SMALL)  THEN
 
151
               RADCAL = (NU - (TWO*RHO+ONE))*NU + SIGMA**2
 
152
               IF (RADCAL .LT. ZERO)  GO TO 5002
 
153
               THAT = (CP - SQRT(RADCAL)) / (THREE*NU)
 
154
            ELSE
 
155
               THAT = ONE/(TWO*SIGMA)
 
156
            ENDIF
 
157
         ENDIF
 
158
         PHI = THAT*((NU*THAT - CP)*THAT + ONE)
 
159
C
 
160
C          CONVERT TO DISTANCE FROM F2 IF IEXTRM.NE.1 .
 
161
         IF (IEXTRM .NE. 1)  PHI = PHI - RHO
 
162
C
 
163
C          TEST FOR EXCEEDING LIMIT, AND ADJUST ACCORDINGLY.
 
164
         HPHI = H * ABS(PHI)
 
165
         IF (HPHI*ABS(D1) .GT. DFMAX)  THEN
 
166
C           AT THIS POINT, HPHI.GT.0, SO DIVIDE IS OK.
 
167
            D1 = SIGN (DFMAX/HPHI, D1)
 
168
            D2 = -LAMBDA*D1
 
169
         ENDIF
 
170
      ENDIF
 
171
C
 
172
C  NORMAL RETURN.
 
173
C
 
174
 5000 CONTINUE
 
175
      IERR = 0
 
176
      RETURN
 
177
C
 
178
C  ERROR RETURNS.
 
179
C
 
180
 5001 CONTINUE
 
181
C     D1 AND D2 BOTH ZERO, OR BOTH NONZERO AND SAME SIGN.
 
182
      IERR = -1
 
183
      CALL XERMSG ('SLATEC', 'PCHSW', 'D1 AND/OR D2 INVALID', IERR, 1)
 
184
      RETURN
 
185
C
 
186
 5002 CONTINUE
 
187
C     NEGATIVE VALUE OF RADICAL (SHOULD NEVER OCCUR).
 
188
      IERR = -2
 
189
      CALL XERMSG ('SLATEC', 'PCHSW', 'NEGATIVE RADICAL', IERR, 1)
 
190
      RETURN
 
191
C------------- LAST LINE OF PCHSW FOLLOWS ------------------------------
 
192
      END