| 1 | *DECK PCHSP |
|---|
| 2 | SUBROUTINE PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR) |
|---|
| 3 | C***BEGIN PROLOGUE PCHSP |
|---|
| 4 | C***PURPOSE Set derivatives needed to determine the Hermite represen- |
|---|
| 5 | C tation of the cubic spline interpolant to given data, with |
|---|
| 6 | C specified boundary conditions. |
|---|
| 7 | C***LIBRARY SLATEC (PCHIP) |
|---|
| 8 | C***CATEGORY E1A |
|---|
| 9 | C***TYPE SINGLE PRECISION (PCHSP-S, DPCHSP-D) |
|---|
| 10 | C***KEYWORDS CUBIC HERMITE INTERPOLATION, PCHIP, |
|---|
| 11 | C PIECEWISE CUBIC INTERPOLATION, SPLINE INTERPOLATION |
|---|
| 12 | C***AUTHOR Fritsch, F. N., (LLNL) |
|---|
| 13 | C Lawrence Livermore National Laboratory |
|---|
| 14 | C P.O. Box 808 (L-316) |
|---|
| 15 | C Livermore, CA 94550 |
|---|
| 16 | C FTS 532-4275, (510) 422-4275 |
|---|
| 17 | C***DESCRIPTION |
|---|
| 18 | C |
|---|
| 19 | C PCHSP: Piecewise Cubic Hermite Spline |
|---|
| 20 | C |
|---|
| 21 | C Computes the Hermite representation of the cubic spline inter- |
|---|
| 22 | C polant to the data given in X and F satisfying the boundary |
|---|
| 23 | C conditions specified by IC and VC. |
|---|
| 24 | C |
|---|
| 25 | C To facilitate two-dimensional applications, includes an increment |
|---|
| 26 | C between successive values of the F- and D-arrays. |
|---|
| 27 | C |
|---|
| 28 | C The resulting piecewise cubic Hermite function may be evaluated |
|---|
| 29 | C by PCHFE or PCHFD. |
|---|
| 30 | C |
|---|
| 31 | C NOTE: This is a modified version of C. de Boor's cubic spline |
|---|
| 32 | C routine CUBSPL. |
|---|
| 33 | C |
|---|
| 34 | C ---------------------------------------------------------------------- |
|---|
| 35 | C |
|---|
| 36 | C Calling sequence: |
|---|
| 37 | C |
|---|
| 38 | C PARAMETER (INCFD = ...) |
|---|
| 39 | C INTEGER IC(2), N, NWK, IERR |
|---|
| 40 | C REAL VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(NWK) |
|---|
| 41 | C |
|---|
| 42 | C CALL PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR) |
|---|
| 43 | C |
|---|
| 44 | C Parameters: |
|---|
| 45 | C |
|---|
| 46 | C IC -- (input) integer array of length 2 specifying desired |
|---|
| 47 | C boundary conditions: |
|---|
| 48 | C IC(1) = IBEG, desired condition at beginning of data. |
|---|
| 49 | C IC(2) = IEND, desired condition at end of data. |
|---|
| 50 | C |
|---|
| 51 | C IBEG = 0 to set D(1) so that the third derivative is con- |
|---|
| 52 | C tinuous at X(2). This is the "not a knot" condition |
|---|
| 53 | C provided by de Boor's cubic spline routine CUBSPL. |
|---|
| 54 | C < This is the default boundary condition. > |
|---|
| 55 | C IBEG = 1 if first derivative at X(1) is given in VC(1). |
|---|
| 56 | C IBEG = 2 if second derivative at X(1) is given in VC(1). |
|---|
| 57 | C IBEG = 3 to use the 3-point difference formula for D(1). |
|---|
| 58 | C (Reverts to the default b.c. if N.LT.3 .) |
|---|
| 59 | C IBEG = 4 to use the 4-point difference formula for D(1). |
|---|
| 60 | C (Reverts to the default b.c. if N.LT.4 .) |
|---|
| 61 | C NOTES: |
|---|
| 62 | C 1. An error return is taken if IBEG is out of range. |
|---|
| 63 | C 2. For the "natural" boundary condition, use IBEG=2 and |
|---|
| 64 | C VC(1)=0. |
|---|
| 65 | C |
|---|
| 66 | C IEND may take on the same values as IBEG, but applied to |
|---|
| 67 | C derivative at X(N). In case IEND = 1 or 2, the value is |
|---|
| 68 | C given in VC(2). |
|---|
| 69 | C |
|---|
| 70 | C NOTES: |
|---|
| 71 | C 1. An error return is taken if IEND is out of range. |
|---|
| 72 | C 2. For the "natural" boundary condition, use IEND=2 and |
|---|
| 73 | C VC(2)=0. |
|---|
| 74 | C |
|---|
| 75 | C VC -- (input) real array of length 2 specifying desired boundary |
|---|
| 76 | C values, as indicated above. |
|---|
| 77 | C VC(1) need be set only if IC(1) = 1 or 2 . |
|---|
| 78 | C VC(2) need be set only if IC(2) = 1 or 2 . |
|---|
| 79 | C |
|---|
| 80 | C N -- (input) number of data points. (Error return if N.LT.2 .) |
|---|
| 81 | C |
|---|
| 82 | C X -- (input) real array of independent variable values. The |
|---|
| 83 | C elements of X must be strictly increasing: |
|---|
| 84 | C X(I-1) .LT. X(I), I = 2(1)N. |
|---|
| 85 | C (Error return if not.) |
|---|
| 86 | C |
|---|
| 87 | C F -- (input) real array of dependent variable values to be inter- |
|---|
| 88 | C polated. F(1+(I-1)*INCFD) is value corresponding to X(I). |
|---|
| 89 | C |
|---|
| 90 | C D -- (output) real array of derivative values at the data points. |
|---|
| 91 | C These values will determine the cubic spline interpolant |
|---|
| 92 | C with the requested boundary conditions. |
|---|
| 93 | C The value corresponding to X(I) is stored in |
|---|
| 94 | C D(1+(I-1)*INCFD), I=1(1)N. |
|---|
| 95 | C No other entries in D are changed. |
|---|
| 96 | C |
|---|
| 97 | C INCFD -- (input) increment between successive values in F and D. |
|---|
| 98 | C This argument is provided primarily for 2-D applications. |
|---|
| 99 | C (Error return if INCFD.LT.1 .) |
|---|
| 100 | C |
|---|
| 101 | C WK -- (scratch) real array of working storage. |
|---|
| 102 | C |
|---|
| 103 | C NWK -- (input) length of work array. |
|---|
| 104 | C (Error return if NWK.LT.2*N .) |
|---|
| 105 | C |
|---|
| 106 | C IERR -- (output) error flag. |
|---|
| 107 | C Normal return: |
|---|
| 108 | C IERR = 0 (no errors). |
|---|
| 109 | C "Recoverable" errors: |
|---|
| 110 | C IERR = -1 if N.LT.2 . |
|---|
| 111 | C IERR = -2 if INCFD.LT.1 . |
|---|
| 112 | C IERR = -3 if the X-array is not strictly increasing. |
|---|
| 113 | C IERR = -4 if IBEG.LT.0 or IBEG.GT.4 . |
|---|
| 114 | C IERR = -5 if IEND.LT.0 of IEND.GT.4 . |
|---|
| 115 | C IERR = -6 if both of the above are true. |
|---|
| 116 | C IERR = -7 if NWK is too small. |
|---|
| 117 | C NOTE: The above errors are checked in the order listed, |
|---|
| 118 | C and following arguments have **NOT** been validated. |
|---|
| 119 | C (The D-array has not been changed in any of these cases.) |
|---|
| 120 | C IERR = -8 in case of trouble solving the linear system |
|---|
| 121 | C for the interior derivative values. |
|---|
| 122 | C (The D-array may have been changed in this case.) |
|---|
| 123 | C ( Do **NOT** use it! ) |
|---|
| 124 | C |
|---|
| 125 | C***REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- |
|---|
| 126 | C Verlag, New York, 1978, pp. 53-59. |
|---|
| 127 | C***ROUTINES CALLED PCHDF, XERMSG |
|---|
| 128 | C***REVISION HISTORY (YYMMDD) |
|---|
| 129 | C 820503 DATE WRITTEN |
|---|
| 130 | C 820804 Converted to SLATEC library version. |
|---|
| 131 | C 870707 Minor cosmetic changes to prologue. |
|---|
| 132 | C 890411 Added SAVE statements (Vers. 3.2). |
|---|
| 133 | C 890703 Corrected category record. (WRB) |
|---|
| 134 | C 890831 Modified array declarations. (WRB) |
|---|
| 135 | C 890831 REVISION DATE from Version 3.2 |
|---|
| 136 | C 891214 Prologue converted to Version 4.0 format. (BAB) |
|---|
| 137 | C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) |
|---|
| 138 | C 920429 Revised format and order of references. (WRB,FNF) |
|---|
| 139 | C***END PROLOGUE PCHSP |
|---|
| 140 | C Programming notes: |
|---|
| 141 | C |
|---|
| 142 | C To produce a double precision version, simply: |
|---|
| 143 | C a. Change PCHSP to DPCHSP wherever it occurs, |
|---|
| 144 | C b. Change the real declarations to double precision, and |
|---|
| 145 | C c. Change the constants ZERO, HALF, ... to double precision. |
|---|
| 146 | C |
|---|
| 147 | C DECLARE ARGUMENTS. |
|---|
| 148 | C |
|---|
| 149 | INTEGER IC(2), N, INCFD, NWK, IERR |
|---|
| 150 | REAL VC(2), X(*), F(INCFD,*), D(INCFD,*), WK(2,*) |
|---|
| 151 | C |
|---|
| 152 | C DECLARE LOCAL VARIABLES. |
|---|
| 153 | C |
|---|
| 154 | INTEGER IBEG, IEND, INDEX, J, NM1 |
|---|
| 155 | REAL G, HALF, ONE, STEMP(3), THREE, TWO, XTEMP(4), ZERO |
|---|
| 156 | SAVE ZERO, HALF, ONE, TWO, THREE |
|---|
| 157 | REAL PCHDF |
|---|
| 158 | C |
|---|
| 159 | DATA ZERO /0./, HALF /0.5/, ONE /1./, TWO /2./, THREE /3./ |
|---|
| 160 | C |
|---|
| 161 | C VALIDITY-CHECK ARGUMENTS. |
|---|
| 162 | C |
|---|
| 163 | C***FIRST EXECUTABLE STATEMENT PCHSP |
|---|
| 164 | IF ( N.LT.2 ) GO TO 5001 |
|---|
| 165 | IF ( INCFD.LT.1 ) GO TO 5002 |
|---|
| 166 | DO 1 J = 2, N |
|---|
| 167 | IF ( X(J).LE.X(J-1) ) GO TO 5003 |
|---|
| 168 | 1 CONTINUE |
|---|
| 169 | C |
|---|
| 170 | IBEG = IC(1) |
|---|
| 171 | IEND = IC(2) |
|---|
| 172 | IERR = 0 |
|---|
| 173 | IF ( (IBEG.LT.0).OR.(IBEG.GT.4) ) IERR = IERR - 1 |
|---|
| 174 | IF ( (IEND.LT.0).OR.(IEND.GT.4) ) IERR = IERR - 2 |
|---|
| 175 | IF ( IERR.LT.0 ) GO TO 5004 |
|---|
| 176 | C |
|---|
| 177 | C FUNCTION DEFINITION IS OK -- GO ON. |
|---|
| 178 | C |
|---|
| 179 | IF ( NWK .LT. 2*N ) GO TO 5007 |
|---|
| 180 | C |
|---|
| 181 | C COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO, |
|---|
| 182 | C COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.). |
|---|
| 183 | DO 5 J=2,N |
|---|
| 184 | WK(1,J) = X(J) - X(J-1) |
|---|
| 185 | WK(2,J) = (F(1,J) - F(1,J-1))/WK(1,J) |
|---|
| 186 | 5 CONTINUE |
|---|
| 187 | C |
|---|
| 188 | C SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL. |
|---|
| 189 | C |
|---|
| 190 | IF ( IBEG.GT.N ) IBEG = 0 |
|---|
| 191 | IF ( IEND.GT.N ) IEND = 0 |
|---|
| 192 | C |
|---|
| 193 | C SET UP FOR BOUNDARY CONDITIONS. |
|---|
| 194 | C |
|---|
| 195 | IF ( (IBEG.EQ.1).OR.(IBEG.EQ.2) ) THEN |
|---|
| 196 | D(1,1) = VC(1) |
|---|
| 197 | ELSE IF (IBEG .GT. 2) THEN |
|---|
| 198 | C PICK UP FIRST IBEG POINTS, IN REVERSE ORDER. |
|---|
| 199 | DO 10 J = 1, IBEG |
|---|
| 200 | INDEX = IBEG-J+1 |
|---|
| 201 | C INDEX RUNS FROM IBEG DOWN TO 1. |
|---|
| 202 | XTEMP(J) = X(INDEX) |
|---|
| 203 | IF (J .LT. IBEG) STEMP(J) = WK(2,INDEX) |
|---|
| 204 | 10 CONTINUE |
|---|
| 205 | C -------------------------------- |
|---|
| 206 | D(1,1) = PCHDF (IBEG, XTEMP, STEMP, IERR) |
|---|
| 207 | C -------------------------------- |
|---|
| 208 | IF (IERR .NE. 0) GO TO 5009 |
|---|
| 209 | IBEG = 1 |
|---|
| 210 | ENDIF |
|---|
| 211 | C |
|---|
| 212 | IF ( (IEND.EQ.1).OR.(IEND.EQ.2) ) THEN |
|---|
| 213 | D(1,N) = VC(2) |
|---|
| 214 | ELSE IF (IEND .GT. 2) THEN |
|---|
| 215 | C PICK UP LAST IEND POINTS. |
|---|
| 216 | DO 15 J = 1, IEND |
|---|
| 217 | INDEX = N-IEND+J |
|---|
| 218 | C INDEX RUNS FROM N+1-IEND UP TO N. |
|---|
| 219 | XTEMP(J) = X(INDEX) |
|---|
| 220 | IF (J .LT. IEND) STEMP(J) = WK(2,INDEX+1) |
|---|
| 221 | 15 CONTINUE |
|---|
| 222 | C -------------------------------- |
|---|
| 223 | D(1,N) = PCHDF (IEND, XTEMP, STEMP, IERR) |
|---|
| 224 | C -------------------------------- |
|---|
| 225 | IF (IERR .NE. 0) GO TO 5009 |
|---|
| 226 | IEND = 1 |
|---|
| 227 | ENDIF |
|---|
| 228 | C |
|---|
| 229 | C --------------------( BEGIN CODING FROM CUBSPL )-------------------- |
|---|
| 230 | C |
|---|
| 231 | C **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF |
|---|
| 232 | C F AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM- |
|---|
| 233 | C INATION, WITH S(J) ENDING UP IN D(1,J), ALL J. |
|---|
| 234 | C WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE. |
|---|
| 235 | C |
|---|
| 236 | C CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM |
|---|
| 237 | C WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1) |
|---|
| 238 | C |
|---|
| 239 | IF (IBEG .EQ. 0) THEN |
|---|
| 240 | IF (N .EQ. 2) THEN |
|---|
| 241 | C NO CONDITION AT LEFT END AND N = 2. |
|---|
| 242 | WK(2,1) = ONE |
|---|
| 243 | WK(1,1) = ONE |
|---|
| 244 | D(1,1) = TWO*WK(2,2) |
|---|
| 245 | ELSE |
|---|
| 246 | C NOT-A-KNOT CONDITION AT LEFT END AND N .GT. 2. |
|---|
| 247 | WK(2,1) = WK(1,3) |
|---|
| 248 | WK(1,1) = WK(1,2) + WK(1,3) |
|---|
| 249 | D(1,1) =((WK(1,2) + TWO*WK(1,1))*WK(2,2)*WK(1,3) |
|---|
| 250 | * + WK(1,2)**2*WK(2,3)) / WK(1,1) |
|---|
| 251 | ENDIF |
|---|
| 252 | ELSE IF (IBEG .EQ. 1) THEN |
|---|
| 253 | C SLOPE PRESCRIBED AT LEFT END. |
|---|
| 254 | WK(2,1) = ONE |
|---|
| 255 | WK(1,1) = ZERO |
|---|
| 256 | ELSE |
|---|
| 257 | C SECOND DERIVATIVE PRESCRIBED AT LEFT END. |
|---|
| 258 | WK(2,1) = TWO |
|---|
| 259 | WK(1,1) = ONE |
|---|
| 260 | D(1,1) = THREE*WK(2,2) - HALF*WK(1,2)*D(1,1) |
|---|
| 261 | ENDIF |
|---|
| 262 | C |
|---|
| 263 | C IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND |
|---|
| 264 | C CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH |
|---|
| 265 | C EQUATION READS WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J). |
|---|
| 266 | C |
|---|
| 267 | NM1 = N-1 |
|---|
| 268 | IF (NM1 .GT. 1) THEN |
|---|
| 269 | DO 20 J=2,NM1 |
|---|
| 270 | IF (WK(2,J-1) .EQ. ZERO) GO TO 5008 |
|---|
| 271 | G = -WK(1,J+1)/WK(2,J-1) |
|---|
| 272 | D(1,J) = G*D(1,J-1) |
|---|
| 273 | * + THREE*(WK(1,J)*WK(2,J+1) + WK(1,J+1)*WK(2,J)) |
|---|
| 274 | WK(2,J) = G*WK(1,J-1) + TWO*(WK(1,J) + WK(1,J+1)) |
|---|
| 275 | 20 CONTINUE |
|---|
| 276 | ENDIF |
|---|
| 277 | C |
|---|
| 278 | C CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM |
|---|
| 279 | C (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N) |
|---|
| 280 | C |
|---|
| 281 | C IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK- |
|---|
| 282 | C SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT |
|---|
| 283 | C AT THIS POINT. |
|---|
| 284 | IF (IEND .EQ. 1) GO TO 30 |
|---|
| 285 | C |
|---|
| 286 | IF (IEND .EQ. 0) THEN |
|---|
| 287 | IF (N.EQ.2 .AND. IBEG.EQ.0) THEN |
|---|
| 288 | C NOT-A-KNOT AT RIGHT ENDPOINT AND AT LEFT ENDPOINT AND N = 2. |
|---|
| 289 | D(1,2) = WK(2,2) |
|---|
| 290 | GO TO 30 |
|---|
| 291 | ELSE IF ((N.EQ.2) .OR. (N.EQ.3 .AND. IBEG.EQ.0)) THEN |
|---|
| 292 | C EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT* |
|---|
| 293 | C NOT-A-KNOT AT LEFT END POINT). |
|---|
| 294 | D(1,N) = TWO*WK(2,N) |
|---|
| 295 | WK(2,N) = ONE |
|---|
| 296 | IF (WK(2,N-1) .EQ. ZERO) GO TO 5008 |
|---|
| 297 | G = -ONE/WK(2,N-1) |
|---|
| 298 | ELSE |
|---|
| 299 | C NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR ALSO NOT-A- |
|---|
| 300 | C KNOT AT LEFT END POINT. |
|---|
| 301 | G = WK(1,N-1) + WK(1,N) |
|---|
| 302 | C DO NOT NEED TO CHECK FOLLOWING DENOMINATORS (X-DIFFERENCES). |
|---|
| 303 | D(1,N) = ((WK(1,N)+TWO*G)*WK(2,N)*WK(1,N-1) |
|---|
| 304 | * + WK(1,N)**2*(F(1,N-1)-F(1,N-2))/WK(1,N-1))/G |
|---|
| 305 | IF (WK(2,N-1) .EQ. ZERO) GO TO 5008 |
|---|
| 306 | G = -G/WK(2,N-1) |
|---|
| 307 | WK(2,N) = WK(1,N-1) |
|---|
| 308 | ENDIF |
|---|
| 309 | ELSE |
|---|
| 310 | C SECOND DERIVATIVE PRESCRIBED AT RIGHT ENDPOINT. |
|---|
| 311 | D(1,N) = THREE*WK(2,N) + HALF*WK(1,N)*D(1,N) |
|---|
| 312 | WK(2,N) = TWO |
|---|
| 313 | IF (WK(2,N-1) .EQ. ZERO) GO TO 5008 |
|---|
| 314 | G = -ONE/WK(2,N-1) |
|---|
| 315 | ENDIF |
|---|
| 316 | C |
|---|
| 317 | C COMPLETE FORWARD PASS OF GAUSS ELIMINATION. |
|---|
| 318 | C |
|---|
| 319 | WK(2,N) = G*WK(1,N-1) + WK(2,N) |
|---|
| 320 | IF (WK(2,N) .EQ. ZERO) GO TO 5008 |
|---|
| 321 | D(1,N) = (G*D(1,N-1) + D(1,N))/WK(2,N) |
|---|
| 322 | C |
|---|
| 323 | C CARRY OUT BACK SUBSTITUTION |
|---|
| 324 | C |
|---|
| 325 | 30 CONTINUE |
|---|
| 326 | DO 40 J=NM1,1,-1 |
|---|
| 327 | IF (WK(2,J) .EQ. ZERO) GO TO 5008 |
|---|
| 328 | D(1,J) = (D(1,J) - WK(1,J)*D(1,J+1))/WK(2,J) |
|---|
| 329 | 40 CONTINUE |
|---|
| 330 | C --------------------( END CODING FROM CUBSPL )-------------------- |
|---|
| 331 | C |
|---|
| 332 | C NORMAL RETURN. |
|---|
| 333 | C |
|---|
| 334 | RETURN |
|---|
| 335 | C |
|---|
| 336 | C ERROR RETURNS. |
|---|
| 337 | C |
|---|
| 338 | 5001 CONTINUE |
|---|
| 339 | C N.LT.2 RETURN. |
|---|
| 340 | IERR = -1 |
|---|
| 341 | CALL XERMSG ('SLATEC', 'PCHSP', |
|---|
| 342 | + 'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1) |
|---|
| 343 | RETURN |
|---|
| 344 | C |
|---|
| 345 | 5002 CONTINUE |
|---|
| 346 | C INCFD.LT.1 RETURN. |
|---|
| 347 | IERR = -2 |
|---|
| 348 | CALL XERMSG ('SLATEC', 'PCHSP', 'INCREMENT LESS THAN ONE', IERR, |
|---|
| 349 | + 1) |
|---|
| 350 | RETURN |
|---|
| 351 | C |
|---|
| 352 | 5003 CONTINUE |
|---|
| 353 | C X-ARRAY NOT STRICTLY INCREASING. |
|---|
| 354 | IERR = -3 |
|---|
| 355 | CALL XERMSG ('SLATEC', 'PCHSP', 'X-ARRAY NOT STRICTLY INCREASING' |
|---|
| 356 | + , IERR, 1) |
|---|
| 357 | RETURN |
|---|
| 358 | C |
|---|
| 359 | 5004 CONTINUE |
|---|
| 360 | C IC OUT OF RANGE RETURN. |
|---|
| 361 | IERR = IERR - 3 |
|---|
| 362 | CALL XERMSG ('SLATEC', 'PCHSP', 'IC OUT OF RANGE', IERR, 1) |
|---|
| 363 | RETURN |
|---|
| 364 | C |
|---|
| 365 | 5007 CONTINUE |
|---|
| 366 | C NWK TOO SMALL RETURN. |
|---|
| 367 | IERR = -7 |
|---|
| 368 | CALL XERMSG ('SLATEC', 'PCHSP', 'WORK ARRAY TOO SMALL', IERR, 1) |
|---|
| 369 | RETURN |
|---|
| 370 | C |
|---|
| 371 | 5008 CONTINUE |
|---|
| 372 | C SINGULAR SYSTEM. |
|---|
| 373 | C *** THEORETICALLY, THIS CAN ONLY OCCUR IF SUCCESSIVE X-VALUES *** |
|---|
| 374 | C *** ARE EQUAL, WHICH SHOULD ALREADY HAVE BEEN CAUGHT (IERR=-3). *** |
|---|
| 375 | IERR = -8 |
|---|
| 376 | CALL XERMSG ('SLATEC', 'PCHSP', 'SINGULAR LINEAR SYSTEM', IERR, |
|---|
| 377 | + 1) |
|---|
| 378 | RETURN |
|---|
| 379 | C |
|---|
| 380 | 5009 CONTINUE |
|---|
| 381 | C ERROR RETURN FROM PCHDF. |
|---|
| 382 | C *** THIS CASE SHOULD NEVER OCCUR *** |
|---|
| 383 | IERR = -9 |
|---|
| 384 | CALL XERMSG ('SLATEC', 'PCHSP', 'ERROR RETURN FROM PCHDF', IERR, |
|---|
| 385 | + 1) |
|---|
| 386 | RETURN |
|---|
| 387 | C------------- LAST LINE OF PCHSP FOLLOWS ------------------------------ |
|---|
| 388 | END |
|---|
