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
|
subroutine dqcheb(x,fval,cheb12,cheb24)
c***begin prologue dqcheb
c***refer to dqc25c,dqc25f,dqc25s
c***routines called (none)
c***revision date 830518 (yymmdd)
c***keywords chebyshev series expansion, fast fourier transform
c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
c de doncker,elise,appl. math. & progr. div. - k.u.leuven
c***purpose this routine computes the chebyshev series expansion
c of degrees 12 and 24 of a function using a
c fast fourier transform method
c f(x) = sum(k=1,..,13) (cheb12(k)*t(k-1,x)),
c f(x) = sum(k=1,..,25) (cheb24(k)*t(k-1,x)),
c where t(k,x) is the chebyshev polynomial of degree k.
c***description
c
c chebyshev series expansion
c standard fortran subroutine
c double precision version
c
c parameters
c on entry
c x - double precision
c vector of dimension 11 containing the
c values cos(k*pi/24), k = 1, ..., 11
c
c fval - double precision
c vector of dimension 25 containing the
c function values at the points
c (b+a+(b-a)*cos(k*pi/24))/2, k = 0, ...,24,
c where (a,b) is the approximation interval.
c fval(1) and fval(25) are divided by two
c (these values are destroyed at output).
c
c on return
c cheb12 - double precision
c vector of dimension 13 containing the
c chebyshev coefficients for degree 12
c
c cheb24 - double precision
c vector of dimension 25 containing the
c chebyshev coefficients for degree 24
c
c***end prologue dqcheb
c
double precision alam,alam1,alam2,cheb12,cheb24,fval,part1,part2,
* part3,v,x
integer i,j
c
dimension cheb12(13),cheb24(25),fval(25),v(12),x(11)
c
c***first executable statement dqcheb
do 10 i=1,12
j = 26-i
v(i) = fval(i)-fval(j)
fval(i) = fval(i)+fval(j)
10 continue
alam1 = v(1)-v(9)
alam2 = x(6)*(v(3)-v(7)-v(11))
cheb12(4) = alam1+alam2
cheb12(10) = alam1-alam2
alam1 = v(2)-v(8)-v(10)
alam2 = v(4)-v(6)-v(12)
alam = x(3)*alam1+x(9)*alam2
cheb24(4) = cheb12(4)+alam
cheb24(22) = cheb12(4)-alam
alam = x(9)*alam1-x(3)*alam2
cheb24(10) = cheb12(10)+alam
cheb24(16) = cheb12(10)-alam
part1 = x(4)*v(5)
part2 = x(8)*v(9)
part3 = x(6)*v(7)
alam1 = v(1)+part1+part2
alam2 = x(2)*v(3)+part3+x(10)*v(11)
cheb12(2) = alam1+alam2
cheb12(12) = alam1-alam2
alam = x(1)*v(2)+x(3)*v(4)+x(5)*v(6)+x(7)*v(8)
* +x(9)*v(10)+x(11)*v(12)
cheb24(2) = cheb12(2)+alam
cheb24(24) = cheb12(2)-alam
alam = x(11)*v(2)-x(9)*v(4)+x(7)*v(6)-x(5)*v(8)
* +x(3)*v(10)-x(1)*v(12)
cheb24(12) = cheb12(12)+alam
cheb24(14) = cheb12(12)-alam
alam1 = v(1)-part1+part2
alam2 = x(10)*v(3)-part3+x(2)*v(11)
cheb12(6) = alam1+alam2
cheb12(8) = alam1-alam2
alam = x(5)*v(2)-x(9)*v(4)-x(1)*v(6)
* -x(11)*v(8)+x(3)*v(10)+x(7)*v(12)
cheb24(6) = cheb12(6)+alam
cheb24(20) = cheb12(6)-alam
alam = x(7)*v(2)-x(3)*v(4)-x(11)*v(6)+x(1)*v(8)
* -x(9)*v(10)-x(5)*v(12)
cheb24(8) = cheb12(8)+alam
cheb24(18) = cheb12(8)-alam
do 20 i=1,6
j = 14-i
v(i) = fval(i)-fval(j)
fval(i) = fval(i)+fval(j)
20 continue
alam1 = v(1)+x(8)*v(5)
alam2 = x(4)*v(3)
cheb12(3) = alam1+alam2
cheb12(11) = alam1-alam2
cheb12(7) = v(1)-v(5)
alam = x(2)*v(2)+x(6)*v(4)+x(10)*v(6)
cheb24(3) = cheb12(3)+alam
cheb24(23) = cheb12(3)-alam
alam = x(6)*(v(2)-v(4)-v(6))
cheb24(7) = cheb12(7)+alam
cheb24(19) = cheb12(7)-alam
alam = x(10)*v(2)-x(6)*v(4)+x(2)*v(6)
cheb24(11) = cheb12(11)+alam
cheb24(15) = cheb12(11)-alam
do 30 i=1,3
j = 8-i
v(i) = fval(i)-fval(j)
fval(i) = fval(i)+fval(j)
30 continue
cheb12(5) = v(1)+x(8)*v(3)
cheb12(9) = fval(1)-x(8)*fval(3)
alam = x(4)*v(2)
cheb24(5) = cheb12(5)+alam
cheb24(21) = cheb12(5)-alam
alam = x(8)*fval(2)-fval(4)
cheb24(9) = cheb12(9)+alam
cheb24(17) = cheb12(9)-alam
cheb12(1) = fval(1)+fval(3)
alam = fval(2)+fval(4)
cheb24(1) = cheb12(1)+alam
cheb24(25) = cheb12(1)-alam
cheb12(13) = v(1)-v(3)
cheb24(13) = cheb12(13)
alam = 0.1d+01/0.6d+01
do 40 i=2,12
cheb12(i) = cheb12(i)*alam
40 continue
alam = 0.5d+00*alam
cheb12(1) = cheb12(1)*alam
cheb12(13) = cheb12(13)*alam
do 50 i=2,24
cheb24(i) = cheb24(i)*alam
50 continue
cheb24(1) = 0.5d+00*alam*cheb24(1)
cheb24(25) = 0.5d+00*alam*cheb24(25)
return
end
|
