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