Free cookie consent management tool by TermsFeed Policy Generator

source: branches/2789_MathNetNumerics-Exploration/HeuristicLab.Algorithms.DataAnalysis.Experimental/642.f @ 16946

Last change on this file since 16946 was 15449, checked in by gkronber, 7 years ago

#2789 created x64 build of CUBGCV algorithm, fixed some bugs

File size: 18.0 KB
Line 
1       MODULE MODCUBGCV
2       CONTAINS
3C     ALGORITHM 642 COLLECTED ALGORITHMS FROM ACM.
4C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.12, NO. 2,
5C     JUN., 1986, P. 150.
6C   SUBROUTINE NAME     - CUBGCV
7C
8C--------------------------------------------------------------------------
9C
10C   COMPUTER            - VAX/DOUBLE
11C
12C   AUTHOR              - M.F.HUTCHINSON
13C                         CSIRO DIVISION OF MATHEMATICS AND STATISTICS
14C                         P.O. BOX 1965
15C                         CANBERRA, ACT 2601
16C                         AUSTRALIA
17C
18C   LATEST REVISION     - 15 AUGUST 1985
19C
20C   PURPOSE             - CUBIC SPLINE DATA SMOOTHER
21C
22C   USAGE               - CALL CUBGCV (X,F,DF,N,Y,C,IC,VAR,JOB,SE,WK,IER)
23C
24C   ARGUMENTS    X      - VECTOR OF LENGTH N CONTAINING THE
25C                           ABSCISSAE OF THE N DATA POINTS
26C                           (X(I),F(I)) I=1..N. (INPUT) X
27C                           MUST BE ORDERED SO THAT
28C                           X(I) .LT. X(I+1).
29C                F      - VECTOR OF LENGTH N CONTAINING THE
30C                           ORDINATES (OR FUNCTION VALUES)
31C                           OF THE N DATA POINTS (INPUT).
32C                DF     - VECTOR OF LENGTH N. (INPUT/OUTPUT)
33C                           DF(I) IS THE RELATIVE STANDARD DEVIATION
34C                           OF THE ERROR ASSOCIATED WITH DATA POINT I.
35C                           EACH DF(I) MUST BE POSITIVE.  THE VALUES IN
36C                           DF ARE SCALED BY THE SUBROUTINE SO THAT
37C                           THEIR MEAN SQUARE VALUE IS 1, AND UNSCALED
38C                           AGAIN ON NORMAL EXIT.
39C                           THE MEAN SQUARE VALUE OF THE DF(I) IS RETURNED
40C                           IN WK(7) ON NORMAL EXIT.
41C                           IF THE ABSOLUTE STANDARD DEVIATIONS ARE KNOWN,
42C                           THESE SHOULD BE PROVIDED IN DF AND THE ERROR
43C                           VARIANCE PARAMETER VAR (SEE BELOW) SHOULD THEN
44C                           BE SET TO 1.
45C                           IF THE RELATIVE STANDARD DEVIATIONS ARE UNKNOWN,
46C                           SET EACH DF(I)=1.
47C                N      - NUMBER OF DATA POINTS (INPUT).
48C                           N MUST BE .GE. 3.
49C                Y,C    - SPLINE COEFFICIENTS. (OUTPUT) Y
50C                           IS A VECTOR OF LENGTH N. C IS
51C                           AN N-1 BY 3 MATRIX. THE VALUE
52C                           OF THE SPLINE APPROXIMATION AT T IS
53C                           S(T)=((C(I,3)*D+C(I,2))*D+C(I,1))*D+Y(I)
54C                           WHERE X(I).LE.T.LT.X(I+1) AND
55C                           D = T-X(I).
56C                IC     - ROW DIMENSION OF MATRIX C EXACTLY
57C                           AS SPECIFIED IN THE DIMENSION
58C                           STATEMENT IN THE CALLING PROGRAM. (INPUT)
59C                VAR    - ERROR VARIANCE. (INPUT/OUTPUT)
60C                           IF VAR IS NEGATIVE (I.E. UNKNOWN) THEN
61C                           THE SMOOTHING PARAMETER IS DETERMINED
62C                           BY MINIMIZING THE GENERALIZED CROSS VALIDATION
63C                           AND AN ESTIMATE OF THE ERROR VARIANCE IS
64C                           RETURNED IN VAR.
65C                           IF VAR IS NON-NEGATIVE (I.E. KNOWN) THEN THE
66C                           SMOOTHING PARAMETER IS DETERMINED TO MINIMIZE
67C                           AN ESTIMATE, WHICH DEPENDS ON VAR, OF THE TRUE
68C                           MEAN SQUARE ERROR, AND VAR IS UNCHANGED.
69C                           IN PARTICULAR, IF VAR IS ZERO, THEN AN
70C                           INTERPOLATING NATURAL CUBIC SPLINE IS CALCULATED.
71C                           VAR SHOULD BE SET TO 1 IF ABSOLUTE STANDARD
72C                           DEVIATIONS HAVE BEEN PROVIDED IN DF (SEE ABOVE).
73C                JOB    - JOB SELECTION PARAMETER. (INPUT)
74C                         JOB = 0 SHOULD BE SELECTED IF POINT STANDARD ERROR
75C                           ESTIMATES ARE NOT REQUIRED IN SE.
76C                         JOB = 1 SHOULD BE SELECTED IF POINT STANDARD ERROR
77C                           ESTIMATES ARE REQUIRED IN SE.
78C                SE     - VECTOR OF LENGTH N CONTAINING BAYESIAN STANDARD
79C                           ERROR ESTIMATES OF THE FITTED SPLINE VALUES IN Y.
80C                           SE IS NOT REFERENCED IF JOB=0. (OUTPUT)
81C                WK     - WORK VECTOR OF LENGTH 7*(N + 2). ON NORMAL EXIT THE
82C                           FIRST 7 VALUES OF WK ARE ASSIGNED AS FOLLOWS:-
83C
84C                           WK(1) = SMOOTHING PARAMETER (= RHO/(RHO + 1))
85C                           WK(2) = ESTIMATE OF THE NUMBER OF DEGREES OF
86C                                   FREEDOM OF THE RESIDUAL SUM OF SQUARES
87C                           WK(3) = GENERALIZED CROSS VALIDATION
88C                           WK(4) = MEAN SQUARE RESIDUAL
89C                           WK(5) = ESTIMATE OF THE TRUE MEAN SQUARE ERROR
90C                                   AT THE DATA POINTS
91C                           WK(6) = ESTIMATE OF THE ERROR VARIANCE
92C                           WK(7) = MEAN SQUARE VALUE OF THE DF(I)
93C
94C                           IF WK(1)=0 (RHO=0) AN INTERPOLATING NATURAL CUBIC
95C                           SPLINE HAS BEEN CALCULATED.
96C                           IF WK(1)=1 (RHO=INFINITE) A LEAST SQUARES
97C                           REGRESSION LINE HAS BEEN CALCULATED.
98C                           WK(2) IS AN ESTIMATE OF THE NUMBER OF DEGREES OF
99C                           FREEDOM OF THE RESIDUAL WHICH REDUCES TO THE
100C                           USUAL VALUE OF N-2 WHEN A LEAST SQUARES REGRESSION
101C                           LINE IS CALCULATED.
102C                           WK(3),WK(4),WK(5) ARE CALCULATED WITH THE DF(I)
103C                           SCALED TO HAVE MEAN SQUARE VALUE 1.  THE
104C                           UNSCALED VALUES OF WK(3),WK(4),WK(5) MAY BE
105C                           CALCULATED BY DIVIDING BY WK(7).
106C                           WK(6) COINCIDES WITH THE OUTPUT VALUE OF VAR IF
107C                           VAR IS NEGATIVE ON INPUT.  IT IS CALCULATED WITH
108C                           THE UNSCALED VALUES OF THE DF(I) TO FACILITATE
109C                           COMPARISONS WITH A PRIORI VARIANCE ESTIMATES.
110C
111C                IER    - ERROR PARAMETER. (OUTPUT)
112C                         TERMINAL ERROR
113C                           IER = 129, IC IS LESS THAN N-1.
114C                           IER = 130, N IS LESS THAN 3.
115C                           IER = 131, INPUT ABSCISSAE ARE NOT
116C                             ORDERED SO THAT X(I).LT.X(I+1).
117C                           IER = 132, DF(I) IS NOT POSITIVE FOR SOME I.
118C                           IER = 133, JOB IS NOT 0 OR 1.
119C
120C   PRECISION/HARDWARE  - DOUBLE
121C
122C   REQUIRED ROUTINES   - SPINT1,SPFIT1,SPCOF1,SPERR1
123C
124C   REMARKS      THE NUMBER OF ARITHMETIC OPERATIONS REQUIRED BY THE
125C                SUBROUTINE IS PROPORTIONAL TO N.  THE SUBROUTINE
126C                USES AN ALGORITHM DEVELOPED BY M.F. HUTCHINSON AND
127C                F.R. DE HOOG, 'SMOOTHING NOISY DATA WITH SPLINE
128C                FUNCTIONS', NUMER. MATH. (IN PRESS)
129C
130C-----------------------------------------------------------------------
131C
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
136C
137C---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)
141C
142C---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
145C
146      DATA RATIO/2.0D0/
147      DATA TAU/1.618033989D0/
148      DATA ZERO,ONE/0.0D0,1.0D0/
149C
150C---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
157C
158C---CHECK FOR ZERO VARIANCE---
159      IF (VAR.NE.ZERO) GO TO 10
160      R1 = ZERO
161      GO TO 90
162C
163C---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
171C
172C---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
183C
184C---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))
200C
201C---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
224C
225C---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))
229C
230C---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)
235C
236C---UNSCALE DF---
237      DO 120 I = 1,N
238         DF(I) = DF(I)*AVDF
239  120 CONTINUE
240C
241C--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
248C
249C---CHECK FOR ERROR CONDITION---
250  140 CONTINUE
251C     IF (IER.NE.0) CONTINUE
252  150 RETURN
253      END
254C     *******************************************************
255      SUBROUTINE SPINT1(X,AVH,Y,DY,AVDY,N,A,C,IC,R,T,IER)
256C
257C INITIALIZES THE ARRAYS C, R AND T FOR ONE DIMENSIONAL CUBIC
258C SMOOTHING SPLINE FITTING BY SUBROUTINE SPFIT1.  THE VALUES
259C DF(I) ARE SCALED SO THAT THE SUM OF THEIR SQUARES IS N
260C AND THE AVERAGE OF THE DIFFERENCES X(I+1) - X(I) IS CALCULATED
261C IN AVH IN ORDER TO AVOID UNDERFLOW AND OVERFLOW PROBLEMS IN
262C SPFIT1.
263C
264C SUBROUTINE SETS IER IF ELEMENTS OF X ARE NON-INCREASING,
265C IF N IS LESS THAN 3, IF IC IS LESS THAN N-1 OR IF DY(I) IS
266C NOT POSITIVE FOR SOME I.
267C
268C---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
272C
273C---SPECIFICATIONS FOR LOCAL VARIABLES---
274      INTEGER I
275      DOUBLE PRECISION E,F,G,H,ZERO
276      DATA ZERO/0.0D0/
277C
278C---INITIALIZATION AND INPUT CHECKING---
279      IER = 0
280      IF (N.LT.3) GO TO 60
281      IF (IC.LT.N-1) GO TO 70
282C
283C---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)
291C
292C---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)
299C
300      DO 30 I = 1,N
301         DY(I) = DY(I)/AVDY
302   30 CONTINUE
303C
304C---INITIALIZE H,F---
305      H = (X(2)-X(1))/AVH
306      F = (Y(2)-Y(1))/H
307C
308C---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
321C
322C---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
332C
333C---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)
347C
348C FITS A CUBIC SMOOTHING SPLINE TO DATA WITH RELATIVE
349C WEIGHTING DY FOR A GIVEN VALUE OF THE SMOOTHING PARAMETER
350C RHO USING AN ALGORITHM BASED ON THAT OF C.H. REINSCH (1967),
351C NUMER. MATH. 10, 177-183.
352C
353C THE TRACE OF THE INFLUENCE MATRIX IS CALCULATED USING AN
354C ALGORITHM DEVELOPED BY M.F.HUTCHINSON AND F.R.DE HOOG (NUMER.
355C MATH., IN PRESS), ENABLING THE GENERALIZED CROSS VALIDATION
356C AND RELATED STATISTICS TO BE CALCULATED IN ORDER N OPERATIONS.
357C
358C THE ARRAYS A, C, R AND T ARE ASSUMED TO HAVE BEEN INITIALIZED
359C BY THE SUBROUTINE SPINT1.  OVERFLOW AND UNDERFLOW PROBLEMS ARE
360C AVOIDED BY USING P=RHO/(1 + RHO) AND Q=1/(1 + RHO) INSTEAD OF
361C RHO AND BY SCALING THE DIFFERENCES X(I+1) - X(I) BY AVH.
362C
363C THE VALUES IN DF ARE ASSUMED TO HAVE BEEN SCALED SO THAT THE
364C SUM OF THEIR SQUARED VALUES IS N.  THE VALUE IN VAR, WHEN IT IS
365C NON-NEGATIVE, IS ASSUMED TO HAVE BEEN SCALED TO COMPENSATE FOR
366C THE SCALING OF THE VALUES IN DF.
367C
368C THE VALUE RETURNED IN FUN IS AN ESTIMATE OF THE TRUE MEAN SQUARE
369C WHEN VAR IS NON-NEGATIVE, AND IS THE GENERALIZED CROSS VALIDATION
370C WHEN VAR IS NEGATIVE.
371C
372C---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
376C
377C---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/
381C
382C---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
388C
389C---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
404C
405C---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
416C
417C---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)
428C
429C---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
440C
441C---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)
451C
452C---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)
468C
469C CALCULATES BAYESIAN ESTIMATES OF THE STANDARD ERRORS OF THE FITTED
470C VALUES OF A CUBIC SMOOTHING SPLINE BY CALCULATING THE DIAGONAL ELEMENTS
471C OF THE INFLUENCE MATRIX.
472C
473C---SPECIFICATIONS FOR ARGUMENTS---
474      INTEGER N
475      DOUBLE PRECISION X(N),DY(N),R(0:N+1,3),SE(N),AVH,P,VAR
476C
477C---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/
481C
482C---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
488C
489C---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)
500C
501C---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)
508C
509C CALCULATES COEFFICIENTS OF A CUBIC SMOOTHING SPLINE FROM
510C PARAMETERS CALCULATED BY SUBROUTINE SPFIT1.
511C
512C---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
516C
517C---SPECIFICATIONS FOR LOCAL VARIABLES---
518      INTEGER I
519      DOUBLE PRECISION H,QH
520C
521C---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
527C
528C---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
Note: See TracBrowser for help on using the repository browser.