1*DECK PCHQK3
2      SUBROUTINE PCHQK3 (LUN, KPRINT, IPASS)
3C***BEGIN PROLOGUE  PCHQK3
4C***PURPOSE  Test the PCHIP interpolators PCHIC, PCHIM, PCHSP.
5C***LIBRARY   SLATEC (PCHIP)
6C***TYPE      SINGLE PRECISION (PCHQK3-S, DPCHQ3-D)
7C***KEYWORDS  PCHIP INTERPOLATOR QUICK CHECK
8C***AUTHOR  Fritsch, F. N., (LLNL)
9C***DESCRIPTION
10C
11C              PCHIP QUICK CHECK NUMBER 3
12C
13C     TESTS THE INTERPOLATORS:  PCHIC, PCHIM, PCHSP.
14C *Usage:
15C
16C        INTEGER  LUN, KPRINT, IPASS
17C
18C        CALL PCHQK3 (LUN, KPRINT, IPASS)
19C
20C *Arguments:
21C
22C     LUN   :IN  is the unit number to which output is to be written.
23C
24C     KPRINT:IN  controls the amount of output, as specified in the
25C                SLATEC Guidelines.
26C
27C     IPASS:OUT  will contain a pass/fail flag.  IPASS=1 is good.
28C                IPASS=0 indicates one or more tests failed.
29C
30C *Description:
31C
32C   This routine interpolates a constructed data set with all three
33C    PCHIP interpolators and compares the results with those obtained
34C   on a Cray X/MP.  Two different values of the  PCHIC parameter SWITCH
35C   are used.
36C
37C *Remarks:
38C     1. The Cray results are given only to nine significant figures,
39C        so don't expect them to match to more.
40C     2. The results will depend to some extent on the accuracy of
41C        the EXP function.
42C
43C***ROUTINES CALLED  COMP, PCHIC, PCHIM, PCHSP, R1MACH
44C***REVISION HISTORY  (YYMMDD)
45C   900309  DATE WRITTEN
46C   900314  Converted to a subroutine and added a SLATEC 4.0 prologue.
47C   900315  Revised prologue and improved some output formats.  (FNF)
48C   900316  Made TOLD machine-dependent and added extra output when
49C           KPRINT=3.  (FNF)
50C   900320  Added E0's to DATA statement for X to reduce single/double
51C           differences, and other minor cosmetic changes.
52C   900321  Removed IFAIL from call sequence for SLATEC standards and
53C           made miscellaneous cosmetic changes.  (FNF)
54C   900322  Minor changes to reduce single/double differences.  (FNF)
55C   900530  Tolerance (TOLD) changed.  (WRB)
56C   900802  Modified TOLD formula and constants in PCHIC calls to be
57C           compatible with DPCHQ3.  (FNF)
58C   901130  Several significant changes:  (FNF)
59C           1. Changed comparison between PCHIM and PCHIC to only
60C              require agreement to machine precision.
61C           2. Revised to print more output when KPRINT=3.
62C           3. Added 1P's to formats.
63C   910708  Minor modifications in use of KPRINT.  (WRB)
64C   930317  Improved output formats.  (FNF)
65C***END PROLOGUE  PCHQK3
66C
67C*Internal Notes:
68C
69C     TOLD is used to compare with stored Cray results.  Its value
70C          should be consistent with significance of stored values.
71C     TOLZ is used for cases in which exact equality is expected.
72C     TOL  is used for cases in which agreement to machine precision
73C          is expected.
74C**End
75C
76C  Declare arguments.
77C
78      INTEGER  LUN, KPRINT, IPASS
79      LOGICAL  COMP
80      REAL  R1MACH
81C
82C  Declare variables.
83C
84      INTEGER  I, IC(2), IERR, IFAIL, N, NBAD, NBADZ, NWK
85      PARAMETER  (N = 9,  NWK = 2*N)
86      REAL  D(N), DC(N), DC5, DC6, DM(N), DS(N), ERR, F(N), MONE, TOL,
87     *      TOLD, TOLZ, VC(2), X(N), WK(NWK), ZERO
88      PARAMETER  (ZERO = 0.0E0,  MONE = -1.0E0)
89      CHARACTER*6  RESULT
90C
91C  Initialize.
92C
93C       Data.
94      DATA  IC /0, 0/
95      DATA  X /-2.2E0,-1.2E0,-1.0E0,-0.5E0,-0.01E0, 0.5E0, 1.0E0,
96     *          2.0E0, 2.2E0/
97C
98C       Results generated on Cray X/MP (9 sign. figs.)
99      DATA  DM / 0.            , 3.80027352E-01, 7.17253009E-01,
100     *           5.82014161E-01, 0.            ,-5.68208031E-01,
101     *          -5.13501618E-01,-7.77910977E-02,-2.45611117E-03/
102      DATA  DC5,DC6 / 1.76950158E-02,-5.69579814E-01/
103      DATA  DS /-5.16830792E-02, 5.71455855E-01, 7.40530225E-01,
104     *           7.63864934E-01, 1.92614386E-02,-7.65324380E-01,
105     *          -7.28209035E-01,-7.98445427E-02,-2.85983446E-02/
106C
107C***FIRST EXECUTABLE STATEMENT  PCHQK3
108      IFAIL = 0
109C
110C        Set tolerances.
111      TOL  = 10*R1MACH(4)
112      TOLD = MAX( 1.0E-7, 10*TOL )
113      TOLZ = ZERO
114C
115      IF (KPRINT .GE. 3)  WRITE (LUN, 1000)
116      IF (KPRINT .GE. 2)  WRITE (LUN, 1001)
117C
118C  Set up data.
119C
120      DO 10  I = 1, N
121         F(I) = EXP(-X(I)**2)
122   10 CONTINUE
123C
124      IF (KPRINT .GE. 3)  THEN
125         WRITE (LUN, 1002)
126         DO 12  I = 1, 4
127            WRITE (LUN, 1010)  X(I), F(I), DM(I), DS(I)
128   12    CONTINUE
129            WRITE (LUN, 1011)  X(5), F(5), DM(5), DC5, DS(5)
130            WRITE (LUN, 1011)  X(6), F(6), DM(6), DC6, DS(6)
131         DO 15  I = 7, N
132            WRITE (LUN, 1010)  X(I), F(I), DM(I), DS(I)
133   15    CONTINUE
134      ENDIF
135C
136C  Test PCHIM.
137C
138      IF (KPRINT.GE.3)  WRITE (LUN, 2000) 'IM'
139C     --------------------------------
140      CALL PCHIM (N, X, F, D, 1, IERR)
141C     --------------------------------
142C        Expect IERR=1 (one monotonicity switch).
143      IF ( KPRINT.GE.3 )  WRITE (LUN, 2001) 1
144      IF ( .NOT.COMP (IERR, 1, LUN, KPRINT) )  THEN
145         IFAIL = IFAIL + 1
146      ELSE
147         IF ( KPRINT.GE.3 )  WRITE (LUN, 2002)
148         NBAD = 0
149         NBADZ = 0
150         DO 20  I = 1, N
151            RESULT = '  OK'
152C             D-values should agree with stored values.
153C               (Zero values should agree exactly.)
154            IF ( DM(I).EQ.ZERO )  THEN
155               ERR = ABS( D(I) )
156               IF ( ERR.GT.TOLZ )  THEN
157                  NBADZ = NBADZ + 1
158                  RESULT = '**BADZ'
159               ENDIF
160            ELSE
161               ERR = ABS( (D(I)-DM(I))/DM(I) )
162               IF ( ERR.GT.TOLD )  THEN
163                  NBAD = NBAD + 1
164                  RESULT = '**BAD'
165               ENDIF
166            ENDIF
167            IF (KPRINT.GE.3)
168     *         WRITE (LUN, 2003)  I, X(I), D(I), ERR, RESULT
169   20    CONTINUE
170         IF ( (NBADZ.NE.0).OR.(NBAD.NE.0) )  THEN
171            IFAIL = IFAIL + 1
172            IF ((NBADZ.NE.0).AND.(KPRINT.GE.2))
173     *         WRITE (LUN, 2004)  NBAD
174            IF ((NBAD.NE.0).AND.(KPRINT.GE.2))
175     *         WRITE (LUN, 2005)  NBAD, 'IM', TOLD
176         ELSE
177            IF (KPRINT.GE.2)  WRITE (LUN, 2006)  'IM'
178         ENDIF
179      ENDIF
180C
181C  Test PCHIC -- options set to reproduce PCHIM.
182C
183      IF (KPRINT.GE.3)  WRITE (LUN, 2000) 'IC'
184C     --------------------------------------------------------
185      CALL PCHIC (IC, VC, ZERO, N, X, F, DC, 1, WK, NWK, IERR)
186C     --------------------------------------------------------
187C        Expect IERR=0 .
188      IF ( KPRINT.GE.3 )  WRITE (LUN, 2001) 0
189      IF ( .NOT.COMP (IERR, 0, LUN, KPRINT) )  THEN
190         IFAIL = IFAIL + 1
191      ELSE
192         IF ( KPRINT.GE.3 )  WRITE (LUN, 2002)
193         NBAD = 0
194         DO 30  I = 1, N
195            RESULT = '  OK'
196C           D-values should agree exactly with those computed by PCHIM.
197C            (To be generous, will only test to machine precision.)
198            ERR = ABS( D(I)-DC(I) )
199            IF ( ERR.GT.TOL )  THEN
200               NBAD = NBAD + 1
201               RESULT = '**BAD'
202            ENDIF
203            IF (KPRINT.GE.3)
204     *         WRITE (LUN, 2003)  I, X(I), DC(I), ERR, RESULT
205   30    CONTINUE
206         IF ( NBAD.NE.0 )  THEN
207            IFAIL = IFAIL + 1
208            IF (KPRINT.GE.2)  WRITE (LUN, 2005)  NBAD, 'IC', TOL
209         ELSE
210            IF (KPRINT.GE.2)  WRITE (LUN, 2006)  'IC'
211         ENDIF
212      ENDIF
213C
214C  Test PCHIC -- default nonzero switch derivatives.
215C
216      IF (KPRINT.GE.3)  WRITE (LUN, 2000) 'IC'
217C     -------------------------------------------------------
218      CALL PCHIC (IC, VC, MONE, N, X, F, D, 1, WK, NWK, IERR)
219C     -------------------------------------------------------
220C        Expect IERR=0 .
221      IF ( KPRINT.GE.3 )  WRITE (LUN, 2001) 0
222      IF ( .NOT.COMP (IERR, 0, LUN, KPRINT) )  THEN
223         IFAIL = IFAIL + 1
224      ELSE
225         IF ( KPRINT.GE.3 )  WRITE (LUN, 2002)
226         NBAD = 0
227         NBADZ = 0
228         DO 40  I = 1, N
229            RESULT = '  OK'
230C            D-values should agree exactly with those computed in
231C            previous call, except at points 5 and 6.
232            IF ( (I.LT.5).OR.(I.GT.6) )  THEN
233               ERR = ABS( D(I)-DC(I) )
234               IF ( ERR.GT.TOLZ )  THEN
235                  NBADZ = NBADZ + 1
236                  RESULT = '**BADA'
237               ENDIF
238            ELSE
239               IF ( I.EQ.5 )  THEN
240                  ERR = ABS( (D(I)-DC5)/DC5 )
241               ELSE
242                  ERR = ABS( (D(I)-DC6)/DC6 )
243               ENDIF
244               IF ( ERR.GT.TOLD )  THEN
245                  NBAD = NBAD + 1
246                  RESULT = '**BAD'
247               ENDIF
248            ENDIF
249            IF (KPRINT.GE.3)
250     *         WRITE (LUN, 2003)  I, X(I), D(I), ERR, RESULT
251   40    CONTINUE
252         IF ( (NBADZ.NE.0).OR.(NBAD.NE.0) )  THEN
253            IFAIL = IFAIL + 1
254            IF ((NBADZ.NE.0).AND.(KPRINT.GE.2))
255     *         WRITE (LUN, 2007)  NBAD
256            IF ((NBAD.NE.0).AND.(KPRINT.GE.2))
257     *         WRITE (LUN, 2005)  NBAD, 'IC', TOLD
258         ELSE
259            IF (KPRINT.GE.2)  WRITE (LUN, 2006)  'IC'
260         ENDIF
261      ENDIF
262C
263C  Test PCHSP.
264C
265      IF (KPRINT.GE.3)  WRITE (LUN, 2000) 'SP'
266C     -------------------------------------------------
267      CALL PCHSP (IC, VC, N, X, F, D, 1, WK, NWK, IERR)
268C     -------------------------------------------------
269C        Expect IERR=0 .
270      IF ( KPRINT.GE.3 )  WRITE (LUN, 2001) 0
271      IF ( .NOT.COMP (IERR, 0, LUN, KPRINT) )  THEN
272         IFAIL = IFAIL + 1
273      ELSE
274         IF ( KPRINT.GE.3 )  WRITE (LUN, 2002)
275         NBAD = 0
276         DO 50  I = 1, N
277            RESULT = '  OK'
278C             D-values should agree with stored values.
279            ERR = ABS( (D(I)-DS(I))/DS(I) )
280            IF ( ERR.GT.TOLD )  THEN
281               NBAD = NBAD + 1
282               RESULT = '**BAD'
283            ENDIF
284            IF (KPRINT.GE.3)
285     *         WRITE (LUN, 2003)  I, X(I), D(I), ERR, RESULT
286   50    CONTINUE
287         IF ( NBAD.NE.0 )  THEN
288            IFAIL = IFAIL + 1
289            IF (KPRINT.GE.2)  WRITE (LUN, 2005)  NBAD, 'SP', TOLD
290         ELSE
291            IF (KPRINT.GE.2)  WRITE (LUN, 2006)  'SP'
292         ENDIF
293      ENDIF
294C
295C  PRINT SUMMARY AND TERMINATE.
296C
297      IF ((KPRINT.GE.2).AND.(IFAIL.NE.0))  WRITE (LUN, 3001)  IFAIL
298C
299      IF (IFAIL.EQ.0)  THEN
300         IPASS = 1
301         IF (KPRINT.GE.2) WRITE(LUN,99998)
302      ELSE
303         IPASS = 0
304         IF (KPRINT.GE.1) WRITE(LUN,99999)
305      ENDIF
306C
307      RETURN
308C
309C  FORMATS.
310C
311 1000 FORMAT ('1'//10X,'TEST PCHIP INTERPOLATORS')
312 1001 FORMAT (//10X,'PCHQK3 RESULTS'/10X,'--------------')
313 1002 FORMAT (// 5X,'DATA:'
314     *         /39X,'---------- EXPECTED D-VALUES ----------'
315     *         /12X,'X',9X,'F',18X,'DM',13X,'DC',13X,'DS')
316 1010 FORMAT (5X,F10.2,1P,E15.5,4X,E15.5,15X,E15.5)
317 1011 FORMAT (5X,F10.2,1P,E15.5,4X,3E15.5)
318 2000 FORMAT (/5X,'PCH',A2,' TEST:')
319 2001 FORMAT (15X,'EXPECT  IERR =',I5)
320 2002 FORMAT (/9X,'I',7X,'X',9X,'D',13X,'ERR')
321 2003 FORMAT (5X,I5,F10.2,1P,2E15.5,2X,A)
322 2004 FORMAT (/'    **',I5,'  PCHIM RESULTS FAILED TO BE EXACTLY ZERO.')
323 2005 FORMAT (/'    **',I5,'  PCH',A2,' RESULTS FAILED TOLERANCE TEST.',
324     *                     '  TOL =',1P,E10.3)
325 2006 FORMAT (/5X,'  ALL  PCH',A2,' RESULTS OK.')
326 2007 FORMAT (/'    **',I5,'  PCHIC RESULTS FAILED TO AGREE WITH',
327     *                      ' PREVIOUS CALL.')
328 3001 FORMAT (/' *** TROUBLE ***',I5,' INTERPOLATION TESTS FAILED.')
32999998 FORMAT (/' ------------  PCHIP PASSED  ALL INTERPOLATION TESTS',
330     *        ' ------------')
33199999 FORMAT (/' ************  PCHIP FAILED SOME INTERPOLATION TESTS',
332     *        ' ************')
333C------------- LAST LINE OF PCHQK3 FOLLOWS -----------------------------
334      END
335