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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
|
C
C Copyright (c) 1997 Silvano Bonazzola
C
C This file is part of LORENE.
C
C LORENE is free software; you can redistribute it and/or modify
C it under the terms of the GNU General Public License as published by
C the Free Software Foundation; either version 2 of the License, or
C (at your option) any later version.
C
C LORENE is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public License
C along with LORENE; if not, write to the Free Software
C Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
C
C
subroutine chexms(n,n64,f,cc,cs)
implicit double precision(a-h,o-z)
C
C version du 04.12.1996
c
c## routine modifiee le 08.10.1993: ajout des save
c et suppression des variables declarees
c mais non utilisees
c
c
c *** routine pour le calcul simultane des coeff. de tche-
c *** bytchef de n64 fonctions. la routine est completment
c craytinisee.
c **** f doit etre echantillonee dans n+1 points.. dimension minimal
c *** des tableaux: =(n64+1)*(n+3), ou n+1 est le nombre des points
c *** d'echantillonage.
c
c cette routine est legerement plus rapide que la chems.
c l'output est dans cc au lieu de cs
c
c
c arguments de la routine:
c
c n = nombre des degres de liberte-1.
c n64 = nombre des fonctions d'ont on veut calculer la trans-
c formation de tchebytchev.
c f =tableau imput. les valeurs echantillonnes des n64
c fonctions sont stockees en 'parallel': dans f(1),f(2)..
c f(n64) on stocke le valeures de n64 fonctions en x=x0
c dans f(1+n64),f(2+n64),...f(n64+n64), les valeures
c des n64 fonctions dans le point x1 et ainsi de suite.
c si n64 est un multiple de 8, pour des raisons de
c craytinisation les donnes sont stockees dan la fa-
c con suivante: f(1),f(2),...f(n64), f(n64+1+1),
c f(n64+1+2),...f(n64+1+n64).
c cc =output. les coeff. de tchebytchev sont sto-
c en 'parallel'.
c cs =tableau de travail
c
c
c tous les tests de routine ont ete effectues les 20/11/1986
C
C $Id: chexms.f,v 1.2 2012/03/30 12:12:42 j_novak Exp $
C $Log: chexms.f,v $
C Revision 1.2 2012/03/30 12:12:42 j_novak
C Cleaning of fortran files
C
C Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
C LORENE
C
c Revision 1.1 1997/10/23 08:15:39 eric
c Initial revision
c
C
C $Header: /cvsroot/Lorene/F77/Source/Poisson2d/chexms.f,v 1.2 2012/03/30 12:12:42 j_novak Exp $
C
C
character*120 header
data header/'$Header: /cvsroot/Lorene/F77/Source/Poisson2d/chexms.f,v 1.2 2012/03/30 12:12:42 j_novak Exp $'/
dimension cc(*),cs(*),f(*)
dimension sen(513),fn11(513),y(513)
data ndim/0/
data nfon/0/
c##
save ndim, nfon, n1, n2, n21, an2, x0, pi, sen, n65
save n63, n65n65, nm65, nm651, nm652, n6519, jinit
c
n513=513
if(n.lt.n513) go to 12
200 format(10x,'dimensions insuffisantes dans chem64,n=',i4)
call exit
c
c preparation des quantites necessaires pour le calcul.
c ces quantites sont calculees la
c premiere fois qu'on appelle la routine.
c elles sont recalculees toutes les fois qu'on change
c n.
c
12 continue
c if(n.eq.ndim) go to 10
n1=n+1
n2=n/2
n21=n2+1
an2=1.d+00/n2
x0=0
pi=2*acos(x0)/n
do 14 l=1,n1
sen(l)=sin((l-1)*pi)
14 continue
10 continue
c if(nfon.eq.n64.and.n.eq.ndim) go to 18
c
c preparation des quantites necessaires pour le calcul.
c ces quantites sont calculees la premiere fois qu'on
c appelle la routine. elles sont recalculees toutes
c les fois qu'on change n64.
c
n65=n64
ndim=n
if((n64/8)*8.eq.n64) n65=n64+1
n63=n64-1
n65n65=n65+n65
nm65=n65*n
nm651=nm65+n65
nm652=n65*n21
n6519=n65*(n21-2)+1
jinit=3*n65+1
nfon=n64
18 continue
c
c
c pondesration des la partie antisimmetrique des fonc-
c tions.
c
c
do 3 l=1,n65
fn11(l)= (f(l)-f(l+nm65))*.5
3 continue
c
n21l=nm652
n20l=n6519
do 1 l=1,n2
n2lx=n21l-n20l+1
sen20l=sen(n21-l)+.5
lx64=n20l+n63
do 6 lx=n20l,lx64
cs(lx)=(f(lx)-f(lx+n2lx))*sen20l
6 continue
do 20 lx=n20l,lx64
f(lx)=f(lx)-cs(lx)
f(lx+n2lx)=f(lx+n2lx)+cs(lx)
20 continue
n21l=n21l+n65
n20l=n20l-n65
1 continue
c
c calcul de la t.f. fonction precedente.
c
call tfm1s(n,n64,f,cc)
c
do 4 m=1,n65
y(m)=0
cs(m+n65)=0
4 continue
c
jl=n65n65
do 9 l=jinit,nm65,n65n65
jl2=l-1
n65l=l+n63
do 5 m=l,n65l
f(m)=cs(m-n65n65)-cc(m)
5 continue
do 7 m=l,n65l
cs(m)=f(m)
7 continue
do 8 m=1,n65
y(m)=y(m)+cs(m+jl2)
8 continue
9 continue
c
do 13 m=1,n65
f(m)=(fn11(m)-y(m))*an2
13 continue
c
jm0=n65
do 11 l=1,n2
jm=jm0+1
jm65=jm0+n65
do 15 m=jm,jm65
cc(m)=cs(m)+f(m-jm0)
15 continue
jm0=jm0+n65n65
11 continue
110 format(1x,20i5)
100 format(1x,10d12.4)
101 format(1x,'chem64')
return
end
c
|
