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
205
206
207
|
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 CHEYMS(N,N64,F,CC,CS)
IMPLICIT double PRECISION (A-H,O-Z)
C
C version du 04.12.1996
C
C
C *** ROUTINE POUR LE CALCUL SIMULTANE DES COEFFICIENTS DE TCHE-
C *** BYTCHEF DE N64 FONCTIONS. LA ROUTINE EST COMPLETEMENT
C *** CRAYTINISEE.
C *** F DOIT ETRE ECHANTILLONNEE DANS N+1 POINTS. DIMENSION MINIMALE
C *** DES TABLEAUX: =(N64+1)*(N+3), OU N+1 EST LE NOMBRE DES POINTS
C *** D'ECHANTILLONNAGE.
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 FORMEE DE TCHEBYTCHEV.
C F =TABLEAU INPUT. LES VALEURS ECHANTILLONNEES DES N64
C FONCTIONS SONT STOCKEES EN 'PARALLEL': DANS F(1),F(2)..
C F(N64) ON STOCKE LE VALEURS DE N64 FONCTIONS EN X=X0
C DANS F(1+N64),F(2+N64),...F(N64+N64), LES VALEURS
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 DONNEES 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 COEFFICIENTS DE TCHEBYTCHEV SONT STO-
C CKES EN 'PARALLEL'.
C CS =TABLEAU DE TRAVAIL
C
C
C Tous les tests de routine ont ete effectues les 20/11/1986
C
C
C $Id: cheyms.f,v 1.2 2012/03/30 12:12:42 j_novak Exp $
C $Log: cheyms.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:28:14 eric
c Initial revision
c
C
C $Header: /cvsroot/Lorene/F77/Source/Poisson2d/cheyms.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/cheyms.f,v 1.2 2012/03/30 12:12:42 j_novak Exp $'/
PARAMETER (NLOC = 513)
DIMENSION CC(*),CS(*),F(*),SEN(NLOC),FN11(NLOC),Y(NLOC)
DATA NDIM/0/
DATA NFON/0/
C
save NDIM,NFON,SEN,AN2,N2,N21,N65,N63,N65N65,NM65,NM652,N6519,JINIT
C
IF(N.LT.NLOC) GO TO 12
200 FORMAT(10X,'DIMENSIONS INSUFFISANTES DANS CHEYMS, N=',I4)
CALL EXIT
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 LES
C FOIS QU'ON CHANGE N.
C
12 CONTINUE
C
C Initialisation
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
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 PONDERATION DES LA PARTIE ANTISYMETRIQUE DES FONC-
C TIONS.
C
C
DO 19 L=1,NM65
CS(L)=0
19 CONTINUE
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 TFM2S(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
C
DO 11 L=1,N2
JM=JM0+1
JM65=JM0+N65
call titi
if (jm.lt.1) stop 'CHEYMS'
DO 15 M=JM,JM65
CC(M)=CS(M)+F(M-JM0)
15 CONTINUE
JM0=JM0+N65N65
11 CONTINUE
C
110 FORMAT(1X,20I5)
100 FORMAT(1X,10D12.4)
101 FORMAT(1X,'CHEM64')
RETURN
END
|
