1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
|
*DECK PVALUE
SUBROUTINE PVALUE (L, NDER, X, YFIT, YP, A)
C***BEGIN PROLOGUE PVALUE
C***PURPOSE Use the coefficients generated by POLFIT to evaluate the
C polynomial fit of degree L, along with the first NDER of
C its derivatives, at a specified point.
C***LIBRARY SLATEC
C***CATEGORY K6
C***TYPE SINGLE PRECISION (PVALUE-S, DP1VLU-D)
C***KEYWORDS CURVE FITTING, LEAST SQUARES, POLYNOMIAL APPROXIMATION
C***AUTHOR Shampine, L. F., (SNLA)
C Davenport, S. M., (SNLA)
C***DESCRIPTION
C
C Written by L. F. Shampine and S. M. Davenport.
C
C Abstract
C
C The subroutine PVALUE uses the coefficients generated by POLFIT
C to evaluate the polynomial fit of degree L , along with the first
C NDER of its derivatives, at a specified point. Computationally
C stable recurrence relations are used to perform this task.
C
C The parameters for PVALUE are
C
C Input --
C L - the degree of polynomial to be evaluated. L may be
C any non-negative integer which is less than or equal
C to NDEG , the highest degree polynomial provided
C by POLFIT .
C NDER - the number of derivatives to be evaluated. NDER
C may be 0 or any positive value. If NDER is less
C than 0, it will be treated as 0.
C X - the argument at which the polynomial and its
C derivatives are to be evaluated.
C A - work and output array containing values from last
C call to POLFIT .
C
C Output --
C YFIT - value of the fitting polynomial of degree L at X
C YP - array containing the first through NDER derivatives
C of the polynomial of degree L . YP must be
C dimensioned at least NDER in the calling program.
C
C***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston,
C Curve fitting by polynomials in one variable, Report
C SLA-74-0270, Sandia Laboratories, June 1974.
C***ROUTINES CALLED XERMSG
C***REVISION HISTORY (YYMMDD)
C 740601 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890531 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
C 900510 Convert XERRWV calls to XERMSG calls. (RWC)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE PVALUE
DIMENSION YP(*),A(*)
CHARACTER*8 XERN1, XERN2
C***FIRST EXECUTABLE STATEMENT PVALUE
IF (L .LT. 0) GO TO 12
NDO = MAX(NDER,0)
NDO = MIN(NDO,L)
MAXORD = A(1) + 0.5
K1 = MAXORD + 1
K2 = K1 + MAXORD
K3 = K2 + MAXORD + 2
NORD = A(K3) + 0.5
IF (L .GT. NORD) GO TO 11
K4 = K3 + L + 1
IF (NDER .LT. 1) GO TO 2
DO 1 I = 1,NDER
1 YP(I) = 0.0
2 IF (L .GE. 2) GO TO 4
IF (L .EQ. 1) GO TO 3
C
C L IS 0
C
VAL = A(K2+1)
GO TO 10
C
C L IS 1
C
3 CC = A(K2+2)
VAL = A(K2+1) + (X-A(2))*CC
IF (NDER .GE. 1) YP(1) = CC
GO TO 10
C
C L IS GREATER THAN 1
C
4 NDP1 = NDO + 1
K3P1 = K3 + 1
K4P1 = K4 + 1
LP1 = L + 1
LM1 = L - 1
ILO = K3 + 3
IUP = K4 + NDP1
DO 5 I = ILO,IUP
5 A(I) = 0.0
DIF = X - A(LP1)
KC = K2 + LP1
A(K4P1) = A(KC)
A(K3P1) = A(KC-1) + DIF*A(K4P1)
A(K3+2) = A(K4P1)
C
C EVALUATE RECURRENCE RELATIONS FOR FUNCTION VALUE AND DERIVATIVES
C
DO 9 I = 1,LM1
IN = L - I
INP1 = IN + 1
K1I = K1 + INP1
IC = K2 + IN
DIF = X - A(INP1)
VAL = A(IC) + DIF*A(K3P1) - A(K1I)*A(K4P1)
IF (NDO .LE. 0) GO TO 8
DO 6 N = 1,NDO
K3PN = K3P1 + N
K4PN = K4P1 + N
6 YP(N) = DIF*A(K3PN) + N*A(K3PN-1) - A(K1I)*A(K4PN)
C
C SAVE VALUES NEEDED FOR NEXT EVALUATION OF RECURRENCE RELATIONS
C
DO 7 N = 1,NDO
K3PN = K3P1 + N
K4PN = K4P1 + N
A(K4PN) = A(K3PN)
7 A(K3PN) = YP(N)
8 A(K4P1) = A(K3P1)
9 A(K3P1) = VAL
C
C NORMAL RETURN OR ABORT DUE TO ERROR
C
10 YFIT = VAL
RETURN
C
11 WRITE (XERN1, '(I8)') L
WRITE (XERN2, '(I8)') NORD
CALL XERMSG ('SLATEC', 'PVALUE',
* 'THE ORDER OF POLYNOMIAL EVALUATION, L = ' // XERN1 //
* ' REQUESTED EXCEEDS THE HIGHEST ORDER FIT, NORD = ' // XERN2 //
* ', COMPUTED BY POLFIT -- EXECUTION TERMINATED.', 8, 2)
RETURN
C
12 CALL XERMSG ('SLATEC', 'PVALUE',
+ 'INVALID INPUT PARAMETER. ORDER OF POLYNOMIAL EVALUATION ' //
+ 'REQUESTED IS NEGATIVE -- EXECUTION TERMINATED.', 2, 2)
RETURN
END
|
