source: LMDZ5/branches/LMDZ5_SPLA/libf/bibio/pchsp.F @ 3064

Last change on this file since 3064 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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