| 1 | MODULE MODCUBGCV |
|---|
| 2 | CONTAINS |
|---|
| 3 | C ALGORITHM 642 COLLECTED ALGORITHMS FROM ACM. |
|---|
| 4 | C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.12, NO. 2, |
|---|
| 5 | C JUN., 1986, P. 150. |
|---|
| 6 | C SUBROUTINE NAME - CUBGCV |
|---|
| 7 | C |
|---|
| 8 | C-------------------------------------------------------------------------- |
|---|
| 9 | C |
|---|
| 10 | C COMPUTER - VAX/DOUBLE |
|---|
| 11 | C |
|---|
| 12 | C AUTHOR - M.F.HUTCHINSON |
|---|
| 13 | C CSIRO DIVISION OF MATHEMATICS AND STATISTICS |
|---|
| 14 | C P.O. BOX 1965 |
|---|
| 15 | C CANBERRA, ACT 2601 |
|---|
| 16 | C AUSTRALIA |
|---|
| 17 | C |
|---|
| 18 | C LATEST REVISION - 15 AUGUST 1985 |
|---|
| 19 | C |
|---|
| 20 | C PURPOSE - CUBIC SPLINE DATA SMOOTHER |
|---|
| 21 | C |
|---|
| 22 | C USAGE - CALL CUBGCV (X,F,DF,N,Y,C,IC,VAR,JOB,SE,WK,IER) |
|---|
| 23 | C |
|---|
| 24 | C ARGUMENTS X - VECTOR OF LENGTH N CONTAINING THE |
|---|
| 25 | C ABSCISSAE OF THE N DATA POINTS |
|---|
| 26 | C (X(I),F(I)) I=1..N. (INPUT) X |
|---|
| 27 | C MUST BE ORDERED SO THAT |
|---|
| 28 | C X(I) .LT. X(I+1). |
|---|
| 29 | C F - VECTOR OF LENGTH N CONTAINING THE |
|---|
| 30 | C ORDINATES (OR FUNCTION VALUES) |
|---|
| 31 | C OF THE N DATA POINTS (INPUT). |
|---|
| 32 | C DF - VECTOR OF LENGTH N. (INPUT/OUTPUT) |
|---|
| 33 | C DF(I) IS THE RELATIVE STANDARD DEVIATION |
|---|
| 34 | C OF THE ERROR ASSOCIATED WITH DATA POINT I. |
|---|
| 35 | C EACH DF(I) MUST BE POSITIVE. THE VALUES IN |
|---|
| 36 | C DF ARE SCALED BY THE SUBROUTINE SO THAT |
|---|
| 37 | C THEIR MEAN SQUARE VALUE IS 1, AND UNSCALED |
|---|
| 38 | C AGAIN ON NORMAL EXIT. |
|---|
| 39 | C THE MEAN SQUARE VALUE OF THE DF(I) IS RETURNED |
|---|
| 40 | C IN WK(7) ON NORMAL EXIT. |
|---|
| 41 | C IF THE ABSOLUTE STANDARD DEVIATIONS ARE KNOWN, |
|---|
| 42 | C THESE SHOULD BE PROVIDED IN DF AND THE ERROR |
|---|
| 43 | C VARIANCE PARAMETER VAR (SEE BELOW) SHOULD THEN |
|---|
| 44 | C BE SET TO 1. |
|---|
| 45 | C IF THE RELATIVE STANDARD DEVIATIONS ARE UNKNOWN, |
|---|
| 46 | C SET EACH DF(I)=1. |
|---|
| 47 | C N - NUMBER OF DATA POINTS (INPUT). |
|---|
| 48 | C N MUST BE .GE. 3. |
|---|
| 49 | C Y,C - SPLINE COEFFICIENTS. (OUTPUT) Y |
|---|
| 50 | C IS A VECTOR OF LENGTH N. C IS |
|---|
| 51 | C AN N-1 BY 3 MATRIX. THE VALUE |
|---|
| 52 | C OF THE SPLINE APPROXIMATION AT T IS |
|---|
| 53 | C S(T)=((C(I,3)*D+C(I,2))*D+C(I,1))*D+Y(I) |
|---|
| 54 | C WHERE X(I).LE.T.LT.X(I+1) AND |
|---|
| 55 | C D = T-X(I). |
|---|
| 56 | C IC - ROW DIMENSION OF MATRIX C EXACTLY |
|---|
| 57 | C AS SPECIFIED IN THE DIMENSION |
|---|
| 58 | C STATEMENT IN THE CALLING PROGRAM. (INPUT) |
|---|
| 59 | C VAR - ERROR VARIANCE. (INPUT/OUTPUT) |
|---|
| 60 | C IF VAR IS NEGATIVE (I.E. UNKNOWN) THEN |
|---|
| 61 | C THE SMOOTHING PARAMETER IS DETERMINED |
|---|
| 62 | C BY MINIMIZING THE GENERALIZED CROSS VALIDATION |
|---|
| 63 | C AND AN ESTIMATE OF THE ERROR VARIANCE IS |
|---|
| 64 | C RETURNED IN VAR. |
|---|
| 65 | C IF VAR IS NON-NEGATIVE (I.E. KNOWN) THEN THE |
|---|
| 66 | C SMOOTHING PARAMETER IS DETERMINED TO MINIMIZE |
|---|
| 67 | C AN ESTIMATE, WHICH DEPENDS ON VAR, OF THE TRUE |
|---|
| 68 | C MEAN SQUARE ERROR, AND VAR IS UNCHANGED. |
|---|
| 69 | C IN PARTICULAR, IF VAR IS ZERO, THEN AN |
|---|
| 70 | C INTERPOLATING NATURAL CUBIC SPLINE IS CALCULATED. |
|---|
| 71 | C VAR SHOULD BE SET TO 1 IF ABSOLUTE STANDARD |
|---|
| 72 | C DEVIATIONS HAVE BEEN PROVIDED IN DF (SEE ABOVE). |
|---|
| 73 | C JOB - JOB SELECTION PARAMETER. (INPUT) |
|---|
| 74 | C JOB = 0 SHOULD BE SELECTED IF POINT STANDARD ERROR |
|---|
| 75 | C ESTIMATES ARE NOT REQUIRED IN SE. |
|---|
| 76 | C JOB = 1 SHOULD BE SELECTED IF POINT STANDARD ERROR |
|---|
| 77 | C ESTIMATES ARE REQUIRED IN SE. |
|---|
| 78 | C SE - VECTOR OF LENGTH N CONTAINING BAYESIAN STANDARD |
|---|
| 79 | C ERROR ESTIMATES OF THE FITTED SPLINE VALUES IN Y. |
|---|
| 80 | C SE IS NOT REFERENCED IF JOB=0. (OUTPUT) |
|---|
| 81 | C WK - WORK VECTOR OF LENGTH 7*(N + 2). ON NORMAL EXIT THE |
|---|
| 82 | C FIRST 7 VALUES OF WK ARE ASSIGNED AS FOLLOWS:- |
|---|
| 83 | C |
|---|
| 84 | C WK(1) = SMOOTHING PARAMETER (= RHO/(RHO + 1)) |
|---|
| 85 | C WK(2) = ESTIMATE OF THE NUMBER OF DEGREES OF |
|---|
| 86 | C FREEDOM OF THE RESIDUAL SUM OF SQUARES |
|---|
| 87 | C WK(3) = GENERALIZED CROSS VALIDATION |
|---|
| 88 | C WK(4) = MEAN SQUARE RESIDUAL |
|---|
| 89 | C WK(5) = ESTIMATE OF THE TRUE MEAN SQUARE ERROR |
|---|
| 90 | C AT THE DATA POINTS |
|---|
| 91 | C WK(6) = ESTIMATE OF THE ERROR VARIANCE |
|---|
| 92 | C WK(7) = MEAN SQUARE VALUE OF THE DF(I) |
|---|
| 93 | C |
|---|
| 94 | C IF WK(1)=0 (RHO=0) AN INTERPOLATING NATURAL CUBIC |
|---|
| 95 | C SPLINE HAS BEEN CALCULATED. |
|---|
| 96 | C IF WK(1)=1 (RHO=INFINITE) A LEAST SQUARES |
|---|
| 97 | C REGRESSION LINE HAS BEEN CALCULATED. |
|---|
| 98 | C WK(2) IS AN ESTIMATE OF THE NUMBER OF DEGREES OF |
|---|
| 99 | C FREEDOM OF THE RESIDUAL WHICH REDUCES TO THE |
|---|
| 100 | C USUAL VALUE OF N-2 WHEN A LEAST SQUARES REGRESSION |
|---|
| 101 | C LINE IS CALCULATED. |
|---|
| 102 | C WK(3),WK(4),WK(5) ARE CALCULATED WITH THE DF(I) |
|---|
| 103 | C SCALED TO HAVE MEAN SQUARE VALUE 1. THE |
|---|
| 104 | C UNSCALED VALUES OF WK(3),WK(4),WK(5) MAY BE |
|---|
| 105 | C CALCULATED BY DIVIDING BY WK(7). |
|---|
| 106 | C WK(6) COINCIDES WITH THE OUTPUT VALUE OF VAR IF |
|---|
| 107 | C VAR IS NEGATIVE ON INPUT. IT IS CALCULATED WITH |
|---|
| 108 | C THE UNSCALED VALUES OF THE DF(I) TO FACILITATE |
|---|
| 109 | C COMPARISONS WITH A PRIORI VARIANCE ESTIMATES. |
|---|
| 110 | C |
|---|
| 111 | C IER - ERROR PARAMETER. (OUTPUT) |
|---|
| 112 | C TERMINAL ERROR |
|---|
| 113 | C IER = 129, IC IS LESS THAN N-1. |
|---|
| 114 | C IER = 130, N IS LESS THAN 3. |
|---|
| 115 | C IER = 131, INPUT ABSCISSAE ARE NOT |
|---|
| 116 | C ORDERED SO THAT X(I).LT.X(I+1). |
|---|
| 117 | C IER = 132, DF(I) IS NOT POSITIVE FOR SOME I. |
|---|
| 118 | C IER = 133, JOB IS NOT 0 OR 1. |
|---|
| 119 | C |
|---|
| 120 | C PRECISION/HARDWARE - DOUBLE |
|---|
| 121 | C |
|---|
| 122 | C REQUIRED ROUTINES - SPINT1,SPFIT1,SPCOF1,SPERR1 |
|---|
| 123 | C |
|---|
| 124 | C REMARKS THE NUMBER OF ARITHMETIC OPERATIONS REQUIRED BY THE |
|---|
| 125 | C SUBROUTINE IS PROPORTIONAL TO N. THE SUBROUTINE |
|---|
| 126 | C USES AN ALGORITHM DEVELOPED BY M.F. HUTCHINSON AND |
|---|
| 127 | C F.R. DE HOOG, 'SMOOTHING NOISY DATA WITH SPLINE |
|---|
| 128 | C FUNCTIONS', NUMER. MATH. (IN PRESS) |
|---|
| 129 | C |
|---|
| 130 | C----------------------------------------------------------------------- |
|---|
| 131 | C |
|---|
| 132 | SUBROUTINE CUBGCV(X,F,DF,N,Y,C,IC,VAR,JOB,SE,WK,IER) |
|---|
| 133 | . BIND(C, NAME='cubgcv') |
|---|
| 134 | USE ISO_C_BINDING |
|---|
| 135 | !DEC$ ATTRIBUTES DLLEXPORT::CUBGCV |
|---|
| 136 | C |
|---|
| 137 | C---SPECIFICATIONS FOR ARGUMENTS--- |
|---|
| 138 | INTEGER(KIND=4) N,IC,JOB,IER |
|---|
| 139 | DOUBLE PRECISION X(N),F(N),DF(N),Y(N),C(IC,3),SE(N),VAR, |
|---|
| 140 | . WK(0:N+1,7) |
|---|
| 141 | C |
|---|
| 142 | C---SPECIFICATIONS FOR LOCAL VARIABLES--- |
|---|
| 143 | DOUBLE PRECISION DELTA,ERR,GF1,GF2,GF3,GF4,R1,R2,R3,R4,TAU,RATIO, |
|---|
| 144 | . AVH,AVDF,AVAR,ZERO,ONE,STAT(6),P,Q |
|---|
| 145 | C |
|---|
| 146 | DATA RATIO/2.0D0/ |
|---|
| 147 | DATA TAU/1.618033989D0/ |
|---|
| 148 | DATA ZERO,ONE/0.0D0,1.0D0/ |
|---|
| 149 | C |
|---|
| 150 | C---INITIALIZE--- |
|---|
| 151 | IER = 133 |
|---|
| 152 | IF (JOB.LT.0 .OR. JOB.GT.1) GO TO 140 |
|---|
| 153 | CALL SPINT1(X,AVH,F,DF,AVDF,N,Y,C,IC,WK,WK(0,4),IER) |
|---|
| 154 | IF (IER.NE.0) GO TO 140 |
|---|
| 155 | AVAR = VAR |
|---|
| 156 | IF (VAR.GT.ZERO) AVAR = VAR*AVDF*AVDF |
|---|
| 157 | C |
|---|
| 158 | C---CHECK FOR ZERO VARIANCE--- |
|---|
| 159 | IF (VAR.NE.ZERO) GO TO 10 |
|---|
| 160 | R1 = ZERO |
|---|
| 161 | GO TO 90 |
|---|
| 162 | C |
|---|
| 163 | C---FIND LOCAL MINIMUM OF GCV OR THE EXPECTED MEAN SQUARE ERROR--- |
|---|
| 164 | 10 R1 = ONE |
|---|
| 165 | R2 = RATIO*R1 |
|---|
| 166 | CALL SPFIT1(X,AVH,DF,N,R2,P,Q,GF2,AVAR,STAT,Y,C,IC,WK,WK(0,4), |
|---|
| 167 | . WK(0,6),WK(0,7)) |
|---|
| 168 | 20 CALL SPFIT1(X,AVH,DF,N,R1,P,Q,GF1,AVAR,STAT,Y,C,IC,WK,WK(0,4), |
|---|
| 169 | . WK(0,6),WK(0,7)) |
|---|
| 170 | IF (GF1.GT.GF2) GO TO 30 |
|---|
| 171 | C |
|---|
| 172 | C---EXIT IF P ZERO--- |
|---|
| 173 | IF (P.LE.ZERO) GO TO 100 |
|---|
| 174 | R2 = R1 |
|---|
| 175 | GF2 = GF1 |
|---|
| 176 | R1 = R1/RATIO |
|---|
| 177 | GO TO 20 |
|---|
| 178 | |
|---|
| 179 | 30 R3 = RATIO*R2 |
|---|
| 180 | 40 CALL SPFIT1(X,AVH,DF,N,R3,P,Q,GF3,AVAR,STAT,Y,C,IC,WK,WK(0,4), |
|---|
| 181 | . WK(0,6),WK(0,7)) |
|---|
| 182 | IF (GF3.GT.GF2) GO TO 50 |
|---|
| 183 | C |
|---|
| 184 | C---EXIT IF Q ZERO--- |
|---|
| 185 | IF (Q.LE.ZERO) GO TO 100 |
|---|
| 186 | R2 = R3 |
|---|
| 187 | GF2 = GF3 |
|---|
| 188 | R3 = RATIO*R3 |
|---|
| 189 | GO TO 40 |
|---|
| 190 | |
|---|
| 191 | 50 R2 = R3 |
|---|
| 192 | GF2 = GF3 |
|---|
| 193 | DELTA = (R2-R1)/TAU |
|---|
| 194 | R4 = R1 + DELTA |
|---|
| 195 | R3 = R2 - DELTA |
|---|
| 196 | CALL SPFIT1(X,AVH,DF,N,R3,P,Q,GF3,AVAR,STAT,Y,C,IC,WK,WK(0,4), |
|---|
| 197 | . WK(0,6),WK(0,7)) |
|---|
| 198 | CALL SPFIT1(X,AVH,DF,N,R4,P,Q,GF4,AVAR,STAT,Y,C,IC,WK,WK(0,4), |
|---|
| 199 | . WK(0,6),WK(0,7)) |
|---|
| 200 | C |
|---|
| 201 | C---GOLDEN SECTION SEARCH FOR LOCAL MINIMUM--- |
|---|
| 202 | 60 IF (GF3.GT.GF4) GO TO 70 |
|---|
| 203 | R2 = R4 |
|---|
| 204 | GF2 = GF4 |
|---|
| 205 | R4 = R3 |
|---|
| 206 | GF4 = GF3 |
|---|
| 207 | DELTA = DELTA/TAU |
|---|
| 208 | R3 = R2 - DELTA |
|---|
| 209 | CALL SPFIT1(X,AVH,DF,N,R3,P,Q,GF3,AVAR,STAT,Y,C,IC,WK,WK(0,4), |
|---|
| 210 | . WK(0,6),WK(0,7)) |
|---|
| 211 | GO TO 80 |
|---|
| 212 | |
|---|
| 213 | 70 R1 = R3 |
|---|
| 214 | GF1 = GF3 |
|---|
| 215 | R3 = R4 |
|---|
| 216 | GF3 = GF4 |
|---|
| 217 | DELTA = DELTA/TAU |
|---|
| 218 | R4 = R1 + DELTA |
|---|
| 219 | CALL SPFIT1(X,AVH,DF,N,R4,P,Q,GF4,AVAR,STAT,Y,C,IC,WK,WK(0,4), |
|---|
| 220 | . WK(0,6),WK(0,7)) |
|---|
| 221 | 80 ERR = (R2-R1)/ (R1+R2) |
|---|
| 222 | IF (ERR*ERR+ONE.GT.ONE .AND. ERR.GT.1.0D-6) GO TO 60 |
|---|
| 223 | R1 = (R1+R2)*0.5D0 |
|---|
| 224 | C |
|---|
| 225 | C---CALCULATE SPLINE COEFFICIENTS--- |
|---|
| 226 | 90 CALL SPFIT1(X,AVH,DF,N,R1,P,Q,GF1,AVAR,STAT,Y,C,IC,WK,WK(0,4), |
|---|
| 227 | . WK(0,6),WK(0,7)) |
|---|
| 228 | 100 CALL SPCOF1(X,AVH,F,DF,N,P,Q,Y,C,IC,WK(0,6),WK(0,7)) |
|---|
| 229 | C |
|---|
| 230 | C---OPTIONALLY CALCULATE STANDARD ERROR ESTIMATES--- |
|---|
| 231 | IF (VAR.GE.ZERO) GO TO 110 |
|---|
| 232 | AVAR = STAT(6) |
|---|
| 233 | VAR = AVAR/ (AVDF*AVDF) |
|---|
| 234 | 110 IF (JOB.EQ.1) CALL SPERR1(X,AVH,DF,N,WK,P,AVAR,SE) |
|---|
| 235 | C |
|---|
| 236 | C---UNSCALE DF--- |
|---|
| 237 | DO 120 I = 1,N |
|---|
| 238 | DF(I) = DF(I)*AVDF |
|---|
| 239 | 120 CONTINUE |
|---|
| 240 | C |
|---|
| 241 | C--PUT STATISTICS IN WK--- |
|---|
| 242 | DO 130 I = 0,5 |
|---|
| 243 | WK(I,1) = STAT(I+1) |
|---|
| 244 | 130 CONTINUE |
|---|
| 245 | WK(5,1) = STAT(6)/ (AVDF*AVDF) |
|---|
| 246 | WK(6,1) = AVDF*AVDF |
|---|
| 247 | GO TO 150 |
|---|
| 248 | C |
|---|
| 249 | C---CHECK FOR ERROR CONDITION--- |
|---|
| 250 | 140 CONTINUE |
|---|
| 251 | C IF (IER.NE.0) CONTINUE |
|---|
| 252 | 150 RETURN |
|---|
| 253 | END |
|---|
| 254 | C ******************************************************* |
|---|
| 255 | SUBROUTINE SPINT1(X,AVH,Y,DY,AVDY,N,A,C,IC,R,T,IER) |
|---|
| 256 | C |
|---|
| 257 | C INITIALIZES THE ARRAYS C, R AND T FOR ONE DIMENSIONAL CUBIC |
|---|
| 258 | C SMOOTHING SPLINE FITTING BY SUBROUTINE SPFIT1. THE VALUES |
|---|
| 259 | C DF(I) ARE SCALED SO THAT THE SUM OF THEIR SQUARES IS N |
|---|
| 260 | C AND THE AVERAGE OF THE DIFFERENCES X(I+1) - X(I) IS CALCULATED |
|---|
| 261 | C IN AVH IN ORDER TO AVOID UNDERFLOW AND OVERFLOW PROBLEMS IN |
|---|
| 262 | C SPFIT1. |
|---|
| 263 | C |
|---|
| 264 | C SUBROUTINE SETS IER IF ELEMENTS OF X ARE NON-INCREASING, |
|---|
| 265 | C IF N IS LESS THAN 3, IF IC IS LESS THAN N-1 OR IF DY(I) IS |
|---|
| 266 | C NOT POSITIVE FOR SOME I. |
|---|
| 267 | C |
|---|
| 268 | C---SPECIFICATIONS FOR ARGUMENTS--- |
|---|
| 269 | INTEGER N,IC,IER |
|---|
| 270 | DOUBLE PRECISION X(N),Y(N),DY(N),A(N),C(IC,3),R(0:N+1,3), |
|---|
| 271 | . T(0:N+1,2),AVH,AVDY |
|---|
| 272 | C |
|---|
| 273 | C---SPECIFICATIONS FOR LOCAL VARIABLES--- |
|---|
| 274 | INTEGER I |
|---|
| 275 | DOUBLE PRECISION E,F,G,H,ZERO |
|---|
| 276 | DATA ZERO/0.0D0/ |
|---|
| 277 | C |
|---|
| 278 | C---INITIALIZATION AND INPUT CHECKING--- |
|---|
| 279 | IER = 0 |
|---|
| 280 | IF (N.LT.3) GO TO 60 |
|---|
| 281 | IF (IC.LT.N-1) GO TO 70 |
|---|
| 282 | C |
|---|
| 283 | C---GET AVERAGE X SPACING IN AVH--- |
|---|
| 284 | G = ZERO |
|---|
| 285 | DO 10 I = 1,N - 1 |
|---|
| 286 | H = X(I+1) - X(I) |
|---|
| 287 | IF (H.LE.ZERO) GO TO 80 |
|---|
| 288 | G = G + H |
|---|
| 289 | 10 CONTINUE |
|---|
| 290 | AVH = G/ (N-1) |
|---|
| 291 | C |
|---|
| 292 | C---SCALE RELATIVE WEIGHTS--- |
|---|
| 293 | G = ZERO |
|---|
| 294 | DO 20 I = 1,N |
|---|
| 295 | IF (DY(I).LE.ZERO) GO TO 90 |
|---|
| 296 | G = G + DY(I)*DY(I) |
|---|
| 297 | 20 CONTINUE |
|---|
| 298 | AVDY = DSQRT(G/N) |
|---|
| 299 | C |
|---|
| 300 | DO 30 I = 1,N |
|---|
| 301 | DY(I) = DY(I)/AVDY |
|---|
| 302 | 30 CONTINUE |
|---|
| 303 | C |
|---|
| 304 | C---INITIALIZE H,F--- |
|---|
| 305 | H = (X(2)-X(1))/AVH |
|---|
| 306 | F = (Y(2)-Y(1))/H |
|---|
| 307 | C |
|---|
| 308 | C---CALCULATE A,T,R--- |
|---|
| 309 | DO 40 I = 2,N - 1 |
|---|
| 310 | G = H |
|---|
| 311 | H = (X(I+1)-X(I))/AVH |
|---|
| 312 | E = F |
|---|
| 313 | F = (Y(I+1)-Y(I))/H |
|---|
| 314 | A(I) = F - E |
|---|
| 315 | T(I,1) = 2.0D0* (G+H)/3.0D0 |
|---|
| 316 | T(I,2) = H/3.0D0 |
|---|
| 317 | R(I,3) = DY(I-1)/G |
|---|
| 318 | R(I,1) = DY(I+1)/H |
|---|
| 319 | R(I,2) = -DY(I)/G - DY(I)/H |
|---|
| 320 | 40 CONTINUE |
|---|
| 321 | C |
|---|
| 322 | C---CALCULATE C = R'*R--- |
|---|
| 323 | R(N,2) = ZERO |
|---|
| 324 | R(N,3) = ZERO |
|---|
| 325 | R(N+1,3) = ZERO |
|---|
| 326 | DO 50 I = 2,N - 1 |
|---|
| 327 | C(I,1) = R(I,1)*R(I,1) + R(I,2)*R(I,2) + R(I,3)*R(I,3) |
|---|
| 328 | C(I,2) = R(I,1)*R(I+1,2) + R(I,2)*R(I+1,3) |
|---|
| 329 | C(I,3) = R(I,1)*R(I+2,3) |
|---|
| 330 | 50 CONTINUE |
|---|
| 331 | RETURN |
|---|
| 332 | C |
|---|
| 333 | C---ERROR CONDITIONS--- |
|---|
| 334 | 60 IER = 130 |
|---|
| 335 | RETURN |
|---|
| 336 | |
|---|
| 337 | 70 IER = 129 |
|---|
| 338 | RETURN |
|---|
| 339 | |
|---|
| 340 | 80 IER = 131 |
|---|
| 341 | RETURN |
|---|
| 342 | |
|---|
| 343 | 90 IER = 132 |
|---|
| 344 | RETURN |
|---|
| 345 | END |
|---|
| 346 | SUBROUTINE SPFIT1(X,AVH,DY,N,RHO,P,Q,FUN,VAR,STAT,A,C,IC,R,T,U,V) |
|---|
| 347 | C |
|---|
| 348 | C FITS A CUBIC SMOOTHING SPLINE TO DATA WITH RELATIVE |
|---|
| 349 | C WEIGHTING DY FOR A GIVEN VALUE OF THE SMOOTHING PARAMETER |
|---|
| 350 | C RHO USING AN ALGORITHM BASED ON THAT OF C.H. REINSCH (1967), |
|---|
| 351 | C NUMER. MATH. 10, 177-183. |
|---|
| 352 | C |
|---|
| 353 | C THE TRACE OF THE INFLUENCE MATRIX IS CALCULATED USING AN |
|---|
| 354 | C ALGORITHM DEVELOPED BY M.F.HUTCHINSON AND F.R.DE HOOG (NUMER. |
|---|
| 355 | C MATH., IN PRESS), ENABLING THE GENERALIZED CROSS VALIDATION |
|---|
| 356 | C AND RELATED STATISTICS TO BE CALCULATED IN ORDER N OPERATIONS. |
|---|
| 357 | C |
|---|
| 358 | C THE ARRAYS A, C, R AND T ARE ASSUMED TO HAVE BEEN INITIALIZED |
|---|
| 359 | C BY THE SUBROUTINE SPINT1. OVERFLOW AND UNDERFLOW PROBLEMS ARE |
|---|
| 360 | C AVOIDED BY USING P=RHO/(1 + RHO) AND Q=1/(1 + RHO) INSTEAD OF |
|---|
| 361 | C RHO AND BY SCALING THE DIFFERENCES X(I+1) - X(I) BY AVH. |
|---|
| 362 | C |
|---|
| 363 | C THE VALUES IN DF ARE ASSUMED TO HAVE BEEN SCALED SO THAT THE |
|---|
| 364 | C SUM OF THEIR SQUARED VALUES IS N. THE VALUE IN VAR, WHEN IT IS |
|---|
| 365 | C NON-NEGATIVE, IS ASSUMED TO HAVE BEEN SCALED TO COMPENSATE FOR |
|---|
| 366 | C THE SCALING OF THE VALUES IN DF. |
|---|
| 367 | C |
|---|
| 368 | C THE VALUE RETURNED IN FUN IS AN ESTIMATE OF THE TRUE MEAN SQUARE |
|---|
| 369 | C WHEN VAR IS NON-NEGATIVE, AND IS THE GENERALIZED CROSS VALIDATION |
|---|
| 370 | C WHEN VAR IS NEGATIVE. |
|---|
| 371 | C |
|---|
| 372 | C---SPECIFICATIONS FOR ARGUMENTS--- |
|---|
| 373 | INTEGER IC,N |
|---|
| 374 | DOUBLE PRECISION X(N),DY(N),RHO,STAT(6),A(N),C(IC,3),R(0:N+1,3), |
|---|
| 375 | . T(0:N+1,2),U(0:N+1),V(0:N+1),FUN,VAR,AVH,P,Q |
|---|
| 376 | C |
|---|
| 377 | C---LOCAL VARIABLES--- |
|---|
| 378 | INTEGER I |
|---|
| 379 | DOUBLE PRECISION E,F,G,H,ZERO,ONE,TWO,RHO1 |
|---|
| 380 | DATA ZERO,ONE,TWO/0.0D0,1.0D0,2.0D0/ |
|---|
| 381 | C |
|---|
| 382 | C---USE P AND Q INSTEAD OF RHO TO PREVENT OVERFLOW OR UNDERFLOW--- |
|---|
| 383 | RHO1 = ONE + RHO |
|---|
| 384 | P = RHO/RHO1 |
|---|
| 385 | Q = ONE/RHO1 |
|---|
| 386 | IF (RHO1.EQ.ONE) P = ZERO |
|---|
| 387 | IF (RHO1.EQ.RHO) Q = ZERO |
|---|
| 388 | C |
|---|
| 389 | C---RATIONAL CHOLESKY DECOMPOSITION OF P*C + Q*T--- |
|---|
| 390 | F = ZERO |
|---|
| 391 | G = ZERO |
|---|
| 392 | H = ZERO |
|---|
| 393 | DO 10 I = 0,1 |
|---|
| 394 | R(I,1) = ZERO |
|---|
| 395 | 10 CONTINUE |
|---|
| 396 | DO 20 I = 2,N - 1 |
|---|
| 397 | R(I-2,3) = G*R(I-2,1) |
|---|
| 398 | R(I-1,2) = F*R(I-1,1) |
|---|
| 399 | R(I,1) = ONE/ (P*C(I,1)+Q*T(I,1)-F*R(I-1,2)-G*R(I-2,3)) |
|---|
| 400 | F = P*C(I,2) + Q*T(I,2) - H*R(I-1,2) |
|---|
| 401 | G = H |
|---|
| 402 | H = P*C(I,3) |
|---|
| 403 | 20 CONTINUE |
|---|
| 404 | C |
|---|
| 405 | C---SOLVE FOR U--- |
|---|
| 406 | U(0) = ZERO |
|---|
| 407 | U(1) = ZERO |
|---|
| 408 | DO 30 I = 2,N - 1 |
|---|
| 409 | U(I) = A(I) - R(I-1,2)*U(I-1) - R(I-2,3)*U(I-2) |
|---|
| 410 | 30 CONTINUE |
|---|
| 411 | U(N) = ZERO |
|---|
| 412 | U(N+1) = ZERO |
|---|
| 413 | DO 40 I = N - 1,2,-1 |
|---|
| 414 | U(I) = R(I,1)*U(I) - R(I,2)*U(I+1) - R(I,3)*U(I+2) |
|---|
| 415 | 40 CONTINUE |
|---|
| 416 | C |
|---|
| 417 | C---CALCULATE RESIDUAL VECTOR V--- |
|---|
| 418 | E = ZERO |
|---|
| 419 | H = ZERO |
|---|
| 420 | DO 50 I = 1,N - 1 |
|---|
| 421 | G = H |
|---|
| 422 | H = (U(I+1)-U(I))/ ((X(I+1)-X(I))/AVH) |
|---|
| 423 | V(I) = DY(I)* (H-G) |
|---|
| 424 | E = E + V(I)*V(I) |
|---|
| 425 | 50 CONTINUE |
|---|
| 426 | V(N) = DY(N)* (-H) |
|---|
| 427 | E = E + V(N)*V(N) |
|---|
| 428 | C |
|---|
| 429 | C---CALCULATE UPPER THREE BANDS OF INVERSE MATRIX--- |
|---|
| 430 | R(N,1) = ZERO |
|---|
| 431 | R(N,2) = ZERO |
|---|
| 432 | R(N+1,1) = ZERO |
|---|
| 433 | DO 60 I = N - 1,2,-1 |
|---|
| 434 | G = R(I,2) |
|---|
| 435 | H = R(I,3) |
|---|
| 436 | R(I,2) = -G*R(I+1,1) - H*R(I+1,2) |
|---|
| 437 | R(I,3) = -G*R(I+1,2) - H*R(I+2,1) |
|---|
| 438 | R(I,1) = R(I,1) - G*R(I,2) - H*R(I,3) |
|---|
| 439 | 60 CONTINUE |
|---|
| 440 | C |
|---|
| 441 | C---CALCULATE TRACE--- |
|---|
| 442 | F = ZERO |
|---|
| 443 | G = ZERO |
|---|
| 444 | H = ZERO |
|---|
| 445 | DO 70 I = 2,N - 1 |
|---|
| 446 | F = F + R(I,1)*C(I,1) |
|---|
| 447 | G = G + R(I,2)*C(I,2) |
|---|
| 448 | H = H + R(I,3)*C(I,3) |
|---|
| 449 | 70 CONTINUE |
|---|
| 450 | F = F + TWO* (G+H) |
|---|
| 451 | C |
|---|
| 452 | C---CALCULATE STATISTICS--- |
|---|
| 453 | STAT(1) = P |
|---|
| 454 | STAT(2) = F*P |
|---|
| 455 | STAT(3) = N*E/ (F*F) |
|---|
| 456 | STAT(4) = E*P*P/N |
|---|
| 457 | STAT(6) = E*P/F |
|---|
| 458 | IF (VAR.GE.ZERO) GO TO 80 |
|---|
| 459 | STAT(5) = STAT(6) - STAT(4) |
|---|
| 460 | FUN = STAT(3) |
|---|
| 461 | GO TO 90 |
|---|
| 462 | |
|---|
| 463 | 80 STAT(5) = DMAX1(STAT(4)-TWO*VAR*STAT(2)/N+VAR,ZERO) |
|---|
| 464 | FUN = STAT(5) |
|---|
| 465 | 90 RETURN |
|---|
| 466 | END |
|---|
| 467 | SUBROUTINE SPERR1(X,AVH,DY,N,R,P,VAR,SE) |
|---|
| 468 | C |
|---|
| 469 | C CALCULATES BAYESIAN ESTIMATES OF THE STANDARD ERRORS OF THE FITTED |
|---|
| 470 | C VALUES OF A CUBIC SMOOTHING SPLINE BY CALCULATING THE DIAGONAL ELEMENTS |
|---|
| 471 | C OF THE INFLUENCE MATRIX. |
|---|
| 472 | C |
|---|
| 473 | C---SPECIFICATIONS FOR ARGUMENTS--- |
|---|
| 474 | INTEGER N |
|---|
| 475 | DOUBLE PRECISION X(N),DY(N),R(0:N+1,3),SE(N),AVH,P,VAR |
|---|
| 476 | C |
|---|
| 477 | C---SPECIFICATIONS FOR LOCAL VARIABLES--- |
|---|
| 478 | INTEGER I |
|---|
| 479 | DOUBLE PRECISION F,G,H,F1,G1,H1,ZERO,ONE |
|---|
| 480 | DATA ZERO,ONE/0.0D0,1.0D0/ |
|---|
| 481 | C |
|---|
| 482 | C---INITIALIZE--- |
|---|
| 483 | H = AVH/ (X(2)-X(1)) |
|---|
| 484 | SE(1) = ONE - P*DY(1)*DY(1)*H*H*R(2,1) |
|---|
| 485 | R(1,1) = ZERO |
|---|
| 486 | R(1,2) = ZERO |
|---|
| 487 | R(1,3) = ZERO |
|---|
| 488 | C |
|---|
| 489 | C---CALCULATE DIAGONAL ELEMENTS--- |
|---|
| 490 | DO 10 I = 2,N - 1 |
|---|
| 491 | F = H |
|---|
| 492 | H = AVH/ (X(I+1)-X(I)) |
|---|
| 493 | G = -F - H |
|---|
| 494 | F1 = F*R(I-1,1) + G*R(I-1,2) + H*R(I-1,3) |
|---|
| 495 | G1 = F*R(I-1,2) + G*R(I,1) + H*R(I,2) |
|---|
| 496 | H1 = F*R(I-1,3) + G*R(I,2) + H*R(I+1,1) |
|---|
| 497 | SE(I) = ONE - P*DY(I)*DY(I)* (F*F1+G*G1+H*H1) |
|---|
| 498 | 10 CONTINUE |
|---|
| 499 | SE(N) = ONE - P*DY(N)*DY(N)*H*H*R(N-1,1) |
|---|
| 500 | C |
|---|
| 501 | C---CALCULATE STANDARD ERROR ESTIMATES--- |
|---|
| 502 | DO 20 I = 1,N |
|---|
| 503 | SE(I) = DSQRT(DMAX1(SE(I)*VAR,ZERO))*DY(I) |
|---|
| 504 | 20 CONTINUE |
|---|
| 505 | RETURN |
|---|
| 506 | END |
|---|
| 507 | SUBROUTINE SPCOF1(X,AVH,Y,DY,N,P,Q,A,C,IC,U,V) |
|---|
| 508 | C |
|---|
| 509 | C CALCULATES COEFFICIENTS OF A CUBIC SMOOTHING SPLINE FROM |
|---|
| 510 | C PARAMETERS CALCULATED BY SUBROUTINE SPFIT1. |
|---|
| 511 | C |
|---|
| 512 | C---SPECIFICATIONS FOR ARGUMENTS--- |
|---|
| 513 | INTEGER IC,N |
|---|
| 514 | DOUBLE PRECISION X(N),Y(N),DY(N),P,Q,A(N),C(IC,3),U(0:N+1), |
|---|
| 515 | . V(0:N+1),AVH |
|---|
| 516 | C |
|---|
| 517 | C---SPECIFICATIONS FOR LOCAL VARIABLES--- |
|---|
| 518 | INTEGER I |
|---|
| 519 | DOUBLE PRECISION H,QH |
|---|
| 520 | C |
|---|
| 521 | C---CALCULATE A--- |
|---|
| 522 | QH = Q/ (AVH*AVH) |
|---|
| 523 | DO 10 I = 1,N |
|---|
| 524 | A(I) = Y(I) - P*DY(I)*V(I) |
|---|
| 525 | U(I) = QH*U(I) |
|---|
| 526 | 10 CONTINUE |
|---|
| 527 | C |
|---|
| 528 | C---CALCULATE C--- |
|---|
| 529 | DO 20 I = 1,N - 1 |
|---|
| 530 | H = X(I+1) - X(I) |
|---|
| 531 | C(I,3) = (U(I+1)-U(I))/ (3.0D0*H) |
|---|
| 532 | C(I,1) = (A(I+1)-A(I))/H - (H*C(I,3)+U(I))*H |
|---|
| 533 | C(I,2) = U(I) |
|---|
| 534 | 20 CONTINUE |
|---|
| 535 | RETURN |
|---|
| 536 | END |
|---|
| 537 | END MODULE |
|---|
