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
149
150
151
|
*DECK DP1VLU
SUBROUTINE DP1VLU (L, NDER, X, YFIT, YP, A)
C***BEGIN PROLOGUE DP1VLU
C***PURPOSE Use the coefficients generated by DPOLFT 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 DOUBLE 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 Abstract
C
C The subroutine DP1VLU uses the coefficients generated by DPOLFT
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 DP1VLU are
C
C Input -- ALL TYPE REAL variables are DOUBLE PRECISION
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 DPOLFT .
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 DPOLFT .
C
C Output -- ALL TYPE REAL variables are DOUBLE PRECISION
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 890911 Removed unnecessary intrinsics. (WRB)
C 891006 Cosmetic changes to prologue. (WRB)
C 891006 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 DP1VLU
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
INTEGER I,IC,ILO,IN,INP1,IUP,K1,K1I,K2,K3,K3P1,K3PN,K4,K4P1,K4PN,
* KC,L,LM1,LP1,MAXORD,N,NDER,NDO,NDP1,NORD
DOUBLE PRECISION A(*),CC,DIF,VAL,X,YFIT,YP(*)
CHARACTER*8 XERN1, XERN2
C***FIRST EXECUTABLE STATEMENT DP1VLU
IF (L .LT. 0) GO TO 12
NDO = MAX(NDER,0)
NDO = MIN(NDO,L)
MAXORD = A(1) + 0.5D0
K1 = MAXORD + 1
K2 = K1 + MAXORD
K3 = K2 + MAXORD + 2
NORD = A(K3) + 0.5D0
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.0D0
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.0D0
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', 'DP1VLU',
* 'THE ORDER OF POLYNOMIAL EVALUATION, L = ' // XERN1 //
* ' REQUESTED EXCEEDS THE HIGHEST ORDER FIT, NORD = ' // XERN2 //
* ', COMPUTED BY DPOLFT -- EXECUTION TERMINATED.', 8, 2)
RETURN
C
12 CALL XERMSG ('SLATEC', 'DP1VLU',
+ 'INVALID INPUT PARAMETER. ORDER OF POLYNOMIAL EVALUATION ' //
+ 'REQUESTED IS NEGATIVE.', 2, 2)
RETURN
END
|
