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
208
209
210
211
212
213
214
215
216
217
218
|
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 CD2MRS(N,N64,IND,CC,CS,Y)
IMPLICIT double PRECISION (A-H,O-Z)
C
C CALCUL DE LA DERIVEE 2EME DE N64 FONCTION A ECHANTILLONNAGE
C RAREFIE A L'ORIGINE.
C
C ARGUMENTS DE LA SUBROUTINE:
C
C N =NOMBRE DES DEGRES DE LIBERTEE.
C IND =PARAMETRE: SI IND=0 ON TRAITE LES FONCTIONS
C PAIRES, SI IND=1 LES FONCTION IMPAIRES.
C N64 =NOMBRE DES FONCTION DONT ON VEUT CALCULER LA
C TRANSFORMEE.
C CC =TABLEAU D'ENTREE: IL CONTIENT LES COEFFICIENTS
C DE TCHEBYTCHEV DE LA FONCTION D'ONT ON VEUT CAL-
C CULER LA DERIVEE.
C CS =TABLEAU DE TRAVAIL ET TABLEAU OUTPUT SI
C INDEX=1.
C INDEX =PARAMETRE: SI INDEX=1 ON A EN SORTIE DANS CS LES
C COEFFICIENTS DE LA DERIVEE DE LA FONCTION, SI
C INDEX=2 ON EN SORTIE DANS Y LA FONCTION DERIVEE.
C Y =TABLEAU OUTPUT
C
C LES DIMENSIONS DES TABLEAUX DOIVENT ETRE GT.EQ. (N64+1)*(N+2)
C
C SUBROUTINE TESTEE LE 15/3/1987
C
C
C
C $Id: cd2mrs.f,v 1.2 2012/03/30 12:12:42 j_novak Exp $
C $Log: cd2mrs.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:31 e_gourgoulhon
C LORENE
C
c Revision 1.2 1997/05/23 11:34:44 hyc
c *** empty log message ***
c
C Revision 1.1 1997/03/17 20:34:13 hyc
C Initial revision
C
C
C $Header: /cvsroot/Lorene/F77/Source/Poisson2d/cd2mrs.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/cd2mrs.f,v 1.2 2012/03/30 12:12:42 j_novak Exp $'/
DIMENSION CC(*),CS(*),Y(*),SOMM(1025)
C
DATA NDIM,NEQ/0,0/
save N1,N0,N2,NN,N65,N63,N66,N166,N065,NM66
save N6565,NM65,N1M65,NALP3,NALP4,NALI3,NALI4
save NDIM,NEQ
C
C INITIALISATION
C
IF(NDIM.NE.N) THEN
N1=N+1
N0=N-1
N2=N-2
NN=N+N
ENDIF
C
IF(NDIM.EQ.N.AND.NEQ.EQ.N64) GO TO 800
C
N1025=1025
IF(N64.GT.1025) THEN
PRINT*,'DIMENSIONS INSUFF. DANS CD2MRS, N64=',N64
ENDIF
C
N65=N64
IF((N64/8)*8.EQ.N64) N65=N64+1
N63=N64-1
N66=N65+1
N166=N0*N65+1
N065=N166-1
NM66=N166+N65
N6565=N65+N65
NM65=N*N65
N1M65=NM65+N65
NALP3=4*N0**2
NALP4=4*N**3
NALI3=(N0+N2)**2
NALI4=(N0+N)**3
NDIM=N
NEQ=N64
800 CONTINUE
C
C CALCUL DE LA DERIVEE 2ME D'UNE FONCTION PAIRE.
C
IF(IND.EQ.0) THEN
C
AL3=NALP3
AL4=NALP4
L1=NN
JM1=NM65
JM2=JM1-N65
C
DO 1 M=1,N64
Y(M)=CC(M+JM1)
1 CONTINUE
C
DO 2 M=1,N64
CS(M)=Y(M)*N
SOMM(M)=Y(M)*AL4
2 CONTINUE
C
DO 3 M=1,N64
Y(M+JM2)=SOMM(M)-CS(M)*AL3
3 CONTINUE
C
JM1=JM2
JM2=JM2-N65
L1=N0+N0
DO 4 L=2,N
AL3=(L1-2)**2
AL4=L1**3
DO 5 M=1,N64
Y(M)=CC(M+JM1)
5 CONTINUE
DO 6 M=1,N64
CS(M)=CS(M)+Y(M)*L1
SOMM(M)=SOMM(M)+Y(M)*AL4
6 CONTINUE
DO 7 M=1,N64
Y(M+JM2)=SOMM(M)-CS(M)*AL3
7 CONTINUE
JM1=JM2
JM2=JM2-N65
L1=L1-2
4 CONTINUE
C
DO 8 M=NM66,N1M65
Y(M)=0
8 CONTINUE
RETURN
ENDIF
C
C CALCUL DERIVEE DANS LE CAS D'UNE FONCTION IMPAIRE
C
IF(IND.EQ.1) THEN
C
JM1=N065
JM2=JM1-N65
L1=N0+N
AL3=NALI3
AL4=NALI4
C
DO 9 M=1,N64
Y(M)=CC(M+JM1)
9 CONTINUE
DO 10 M=1,N64
CS(M)=Y(M)*L1
SOMM(M)=Y(M)*AL4
10 CONTINUE
DO 11 M=1,N64
Y(M+JM2)=SOMM(M)-CS(M)*AL3
11 CONTINUE
C
JM1=JM2
JM2=JM2-N65
L1=L1-2
DO 12 L=3,N
AL3=(L1-2)**2
AL4=L1**3
DO 13 M=1,N64
Y(M)=CC(M+JM1)
13 CONTINUE
DO 14 M=1,N64
CS(M)=CS(M)+Y(M)*L1
SOMM(M)=SOMM(M)+Y(M)*AL4
14 CONTINUE
DO 16 M=1,N64
Y(M+JM2)=SOMM(M)-CS(M)*AL3
16 CONTINUE
JM1=JM2
JM2=JM2-N65
L1=L1-2
12 CONTINUE
C
DO 17 M=N166,N1M65
Y(M)=0
17 JM1=JM2
JM2=JM2-N65
L1=L1-2
ENDIF
RETURN
C
100 FORMAT(1X,10E12.4)
101 FORMAT(1X,' ')
END
|
