TreeMig Code
Loading...
Searching...
No Matches
GFT_do.f90
Go to the documentation of this file.
1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -*- Mode: F90 -*- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2!! GFT_do.f90 --- GFT forward/backward FFT routines
3!!
4!! Auteur : Jalel Chergui (CNRS/IDRIS) <Jalel.Chergui@idris.fr>
5!! Cr�� le : Tue Feb 19 10:25:04 2002
6!! Dern. mod. par : Jalel Chergui (CNRS/IDRIS) <Jalel.Chergui@idris.fr>
7!! Dern. mod. le : Thu Jul 11 17:07:50 2002
8!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
9!
10! Permission is granted to copy and distribute this file or modified
11! versions of this file for no fee, provided the copyright notice and
12! this permission notice are preserved on all copies.
13! Copyright (C) Februry 2002, CNRS/IDRIS, Jalel.Chergui@idris.fr.
14!
15MODULE gft_do
16 USE gft_common
17 USE gft_jmfft
18
19 IMPLICIT NONE
20
21 PRIVATE
22
23 PUBLIC :: gft_do_fft
24
25
26
27 INTERFACE gft_do_fft
28
29 MODULE PROCEDURE gft_do_cc_1d, gft_do_cc_2d, gft_do_cc_3d, &
30 gft_do_cr_1d, gft_do_cr_2d, gft_do_cr_3d, &
31 gft_do_rc_1d, gft_do_rc_2d, gft_do_rc_3d, &
32 gft_do_mcc_1d, gft_do_mcr_1d, gft_do_mrc_1d
33 END INTERFACE gft_do_fft
34
35 CONTAINS
36
37 SUBROUTINE gft_do_cc_1d(FFT, isign, scale, c_in, c_out, code)
38
39 IMPLICIT NONE
40
41 !... Input dummy Parameters
42 TYPE(gft_cc), INTENT(inout) :: FFT
43 INTEGER, INTENT(in) :: isign
44 REAL(kind=gft_prec), INTENT(in) :: scale
45 COMPLEX(kind=GFT_prec), INTENT(in), DIMENSION(0:) :: c_in
46
47 !... Output dummy Parameters
48 COMPLEX(kind=GFT_prec), INTENT(out), DIMENSION(0:) :: c_out
49 INTEGER, OPTIONAL, INTENT(out) :: code
50
51 !... Local variables
52 REAL(kind=gft_prec), DIMENSION(0:2*SIZE(c_in)-1) :: x
53 REAL(kind=gft_prec), DIMENSION(0:2*SIZE(c_out)-1) :: y
54 INTEGER :: i, ioff, nfact
55 INTEGER, DIMENSION(0:99) :: fact
56 CHARACTER(len=*), PARAMETER :: spname="GFT_do_cc_1d"
57! print *,spname
58
59
60 !... Check conditions
61 IF (.NOT.fft%Init) CALL gft_error(spname//": Error: GFT_init routine has not been called prior to this routine")
62 IF ( SIZE(c_in) < fft%Nx) CALL gft_error(spname//": Error: Size(C_IN) < Nx")
63 IF ( SIZE(c_out) < fft%Nx) CALL gft_error(spname//": Error: Size(C_OUT) < Nx")
64 IF (isign /=-1 .AND. isign /= 1) CALL gft_error(spname//": Error: sign has invalid value")
65
66 fft%Isign = isign
67 fft%Scale = scale
68
69! "dir$ ivdep"
70! ocl novrec
71! cdir nodep
72 DO i = 0,fft%Nx-1
73 x(2*i) = real(c_in(i),kind=gft_prec)
74 x(2*i+1) = aimag(c_in(i))
75 END DO
76 nfact = nint(fft%Table(0))
77 fact(0:nfact-1) = nint(fft%Table(0:nfact-1))
78
79 ! On copie le tableau d'entree dans le tableau de travail
80 ! On en profite pour premultiplier et pour tenir compte du signe
81 DO i = 0,fft%Nx-1
82 fft%Work(i) = scale* x(2*i)
83 fft%Work(fft%Nx+i) = isign*scale* x(2*i+1)
84 END DO
85 ioff = 0
86
87 ! On appelle le sous-programme principal
88 CALL jmccm1d(1,fft%Nx,fact,100,0,fft%Table,fft%TableSize,100+0,fft%Work,fft%WorkSize,ioff)
89
90 ! On recopie dans le tableau d'arrivee
91 DO i = 0,fft%Nx-1
92 y(2*i) = fft%Work(ioff +i)
93 y(2*i+1) = isign * fft%Work(ioff+fft%Nx+i)
94 END DO
95
96!dir$ ivdep
97!ocl novrec
98!cdir nodep
99 DO i = 0,fft%Nx-1
100 c_out(i) = cmplx(y(2*i),y(2*i+1),kind=gft_prec)
101 END DO
102
103 IF( PRESENT(code) ) code=0
104END SUBROUTINE gft_do_cc_1d
105!!!!!!
106SUBROUTINE gft_do_cc_2d(FFT, isign, scale, c_in, c_out, code)
107 IMPLICIT NONE
108
109 !... Input dummy Parameters
110 TYPE(gft_cc), INTENT(inout) :: FFT
111 INTEGER, INTENT(in) :: isign
112 REAL(kind=gft_prec), INTENT(in) :: scale
113 COMPLEX(kind=GFT_prec), INTENT(in), DIMENSION(0:,0:) :: c_in
114
115 !... Output dummy Parameters
116 COMPLEX(kind=GFT_prec), INTENT(out), DIMENSION(0:,0:) :: c_out
117 INTEGER, OPTIONAL, INTENT(out) :: code
118
119 !... Local variables
120 REAL(kind=gft_prec), DIMENSION(0:2*SIZE(c_in)-1) :: x
121 REAL(kind=gft_prec), DIMENSION(0:2*SIZE(c_out)-1) :: y
122 INTEGER :: i, j, ioff
123 INTEGER :: nfact, mfact
124 INTEGER, DIMENSION(0:99) :: fact
125 INTEGER :: ideb, ifin, jdeb, jfin, n_temp, m_temp, nwork_temp
126 LOGICAL :: debut, fin
127 CHARACTER(len=*), PARAMETER :: spname="GFT_do_cc_2d"
128 ! print *,spname
129
130 !... Check conditions
131 IF (.NOT.fft%Init) CALL gft_error(spname//": Error: GFT_init must be called prior to this routine")
132 fft%Ldx1 = SIZE(c_in ,dim=1)
133 fft%Ldy1 = SIZE(c_out,dim=1)
134 IF ( SIZE(c_in) < fft%Nx*fft%Ny ) CALL gft_error(spname//": Error: Size(C_IN) < Nx*Ny")
135 IF ( SIZE(c_out) < fft%Nx*fft%Ny ) CALL gft_error(spname//": Error: Size(C_OUT) < Nx*Ny")
136 IF ( fft%Ldx1 < fft%Nx ) CALL gft_error(spname//": Error: Size(C_IN,dim=1) < Nx")
137 IF ( fft%Ldy1 < fft%Nx ) CALL gft_error(spname//": Error: Size(C_OUT,dim=1) < Nx")
138 IF (isign /=-1 .AND. isign /= 1) CALL gft_error(spname//": Error: isign value is invalid")
139
140 fft%Isign = isign
141 fft%Scale = scale
142
143 DO j = 0, fft%Ny-1
144!dir$ ivdep
145!ocl novrec
146!cdir nodep
147 DO i = 0,fft%Nx-1
148 x(2*i +j*2*fft%Ldx1) = real(c_in(i,j),kind=gft_prec)
149 x(2*i+1+j*2*fft%Ldx1) = aimag(c_in(i,j))
150 END DO
151 END DO
152
153 nfact = nint(fft%Table(0))
154 mfact = nint(fft%Table(nfact)) + nfact
155 fact(0:mfact-1) = nint(fft%Table(0:mfact-1))
156
157 ! On fait les T.F. sur la premiere dimension en tronconnant sur la deuxieme
158 debut = .true.
159 DO
160
161 ! Tronconnage
162 ! Note : on met npair a .true. car il n'y a pas de restriction dans ce cas
163 CALL jmdecoup(fft%Ny,4*fft%Nx,fft%WorkSize,debut,.true.,m_temp,jdeb,jfin,nwork_temp,fin)
164
165 ! On copie le tableau d'entree dans le tableau de travail
166 ! On en profite pour premultiplier et pour tenir compte du signe
167 ! Note : On copie en transposant
168 DO i = 0,fft%Nx-1
169!dir$ ivdep
170!ocl novrec
171!cdir nodep
172 DO j = jdeb,jfin
173 fft%Work(j-jdeb+m_temp*i) = scale*x(2*i +j*2*fft%Ldx1)
174 fft%Work(j-jdeb+m_temp*(fft%Nx+i)) = isign*scale*x(2*i+1+j*2*fft%Ldx1)
175 END DO
176 END DO
177 ioff = 0
178
179 ! Attention : ioff1 est peut-etre modifie en sortie
180 CALL jmccm1d(m_temp,fft%Nx,fact,100,0 ,fft%Table,fft%TableSize,100+0 ,fft%Work,nwork_temp,ioff)
181
182 ! On recopie dans le tableau d'arrivee
183 DO i = 0,fft%Nx-1
184!dir$ ivdep
185!ocl novrec
186!cdir nodep
187 DO j = jdeb,jfin
188 y(2*i +j*2*fft%Ldy1) = fft%Work(ioff+j-jdeb+m_temp*i)
189 y(2*i+1+j*2*fft%Ldy1) = fft%Work(ioff+j-jdeb+m_temp*(fft%Nx+i))
190 END DO
191 END DO
192
193 ! A-t-on fini ?
194 IF (fin) THEN
195 EXIT
196 ELSE
197 debut = .false.
198 cycle
199 END IF
200
201 END DO
202
203 ! On fait les T.F. sur l'autre dimension
204 debut = .true.
205 DO
206
207 ! Tronconnage
208 CALL jmdecoup(fft%Nx,4*fft%Ny,fft%WorkSize,debut,.true.,n_temp,ideb,ifin,nwork_temp,fin)
209
210 ! On copie
211 DO j = 0,fft%Ny-1
212!dir$ ivdep
213!ocl novrec
214!cdir nodep
215 DO i = ideb,ifin
216 fft%Work(i-ideb+n_temp*j) = y(2*i +j*2*fft%Ldy1)
217 fft%Work(i-ideb+n_temp*(fft%Ny+j)) = y(2*i+1+j*2*fft%Ldy1)
218 END DO
219 END DO
220 ioff = 0
221
222 CALL jmccm1d(n_temp,fft%Ny,fact,100,nfact,fft%Table,fft%TableSize,100+2*fft%Nx,fft%Work,nwork_temp,ioff)
223
224 ! On recopie dans le tableau d'arrivee
225 DO j = 0,fft%Ny-1
226!dir$ ivdep
227!ocl novrec
228!cdir nodep
229 DO i = ideb,ifin
230 y(2*i +j*2*fft%Ldy1) = fft%Work(ioff+i-ideb+n_temp*j)
231 y(2*i+1+j*2*fft%Ldy1) = isign*fft%Work(ioff+i-ideb+n_temp*(fft%Ny+j))
232 END DO
233 END DO
234
235 ! A-t-on fini ?
236 IF (fin) THEN
237 EXIT
238 ELSE
239 debut = .false.
240 cycle
241 END IF
242
243 END DO
244
245 DO j = 0,fft%Ny-1
246!dir$ ivdep
247!ocl novrec
248!cdir nodep
249 DO i = 0,fft%Nx-1
250 c_out(i,j) = cmplx(y(2*i+j*2*fft%Ldy1), y(2*i+1+j*2*fft%Ldy1),kind=gft_prec)
251 END DO
252 END DO
253
254 IF( PRESENT(code) ) code=0
255END SUBROUTINE gft_do_cc_2d
256!!!!!!
257SUBROUTINE gft_do_cc_3d(FFT, isign, scale, c_in, c_out, code)
258 IMPLICIT NONE
259
260 !... Input dummy Parameters
261 TYPE(gft_cc), INTENT(inout) :: FFT
262 INTEGER, INTENT(in) :: isign
263 REAL(kind=gft_prec), INTENT(in) :: scale
264 COMPLEX(kind=GFT_prec), INTENT(in), DIMENSION(0:,0:,0:) :: c_in
265
266 !... Output dummy Parameters
267 COMPLEX(kind=GFT_prec), INTENT(out), DIMENSION(0:,0:,0:) :: c_out
268 INTEGER, OPTIONAL, INTENT(out) :: code
269
270 !... Local variables
271 INTEGER :: i, j, k, ioff
272 INTEGER :: nfact, mfact, lfact
273 INTEGER, DIMENSION(0:99) :: fact
274 INTEGER :: ideb, ifin, i1, i2, jdeb, jfin, j1, j2, kdeb, kfin
275 INTEGER :: nwork_temp, nmtemp, mltemp, nltemp, iwork
276 LOGICAL :: debut, fini
277 CHARACTER(len=*), PARAMETER :: spname="GFT_do_cc_3d"
278 REAL(kind=8), dimension(0:2*SIZE(c_in)-1) :: x
279 REAL(kind=8), dimension(0:2*SIZE(c_out)-1) :: y
280
281 ! print *,spname
282
283 !... Check conditions
284 IF (.NOT.fft%Init) CALL gft_error(spname//": Error: GFT_init must be called prior to this routine")
285 fft%Ldx1 = SIZE(c_in,dim=1)
286 fft%Ldx2 = SIZE(c_in,dim=2)
287 fft%Ldy1 = SIZE(c_out,dim=1)
288 fft%Ldy2 = SIZE(c_out,dim=2)
289 IF ( SIZE(c_in) < fft%Nx*fft%Ny*fft%Nz) CALL gft_error(spname//": Error: Size(C_IN) < Nx*Ny*Nz")
290 IF ( SIZE(c_out) < fft%Nx*fft%Ny*fft%Nz) CALL gft_error(spname//": Error: Size(C_OUT) < Nx*Ny*Nz")
291 IF ( fft%Ldx1 < fft%Nx ) CALL gft_error(spname//": Error: size(C_IN,dim=1) < Nx")
292 IF ( fft%Ldx2 < fft%Ny ) CALL gft_error(spname//": Error: size(C_IN,dim=2) < Ny")
293 IF ( fft%Ldy1 < fft%Nx ) CALL gft_error(spname//": Error: size(C_OUT,dim=1) < Nx")
294 IF ( fft%Ldy2 < fft%Ny ) CALL gft_error(spname//": Error: size(C_OUT,dim=2) < Ny")
295 IF (isign /=-1 .AND. isign /= 1) CALL gft_error(spname//": Error: sign has an invalid value")
296
297
298 fft%Isign = isign
299 fft%Scale = scale
300
301 DO k = 0, fft%Nz-1
302 DO j = 0, fft%Ny-1
303!dir$ ivdep
304!ocl novrec
305!cdir nodep
306 DO i = 0,fft%Nx-1
307 x(2*i +2*fft%Ldx1*j+2*fft%Ldx1*fft%Ldx2*k) = real(c_in(i,j,k),kind=gft_prec)
308 x(2*i+1+2*fft%Ldx1*j+2*fft%Ldx1*fft%Ldx2*k) = aimag(c_in(i,j,k))
309 END DO
310 END DO
311 END DO
312
313 nfact = nint(fft%Table(0))
314 mfact = nint(fft%Table(nfact)) + nfact
315 lfact = nint(fft%Table(mfact)) + mfact
316 fact(0:lfact-1) = nint(fft%Table(0:lfact-1))
317
318 ! On fait les T.F. sur la troisieme dimension en tronconnant sur la premiere
319 ! et la deuxieme
320 debut = .true.
321 fini = .false.
322 DO WHILE (.NOT.fini)
323
324 ! Tronconnage
325 ! Note : on met npair a .true. car il n'y a pas de restriction dans ce cas
326 CALL jmdecoup3(fft%Nx,fft%Ny,4*fft%Nz,fft%WorkSize,debut,.true.,ideb,ifin,jdeb,jfin,nmtemp,nwork_temp,fini)
327 debut = .false.
328
329 ! On copie le tableau d'entree dans le tableau de travail
330 ! On en profite pour premultiplier et pour tenir compte du signe
331 ! On prend garde a la gestion des extremites
332 DO k = 0,fft%Nz-1
333 iwork = 0
334 DO j = jdeb,jfin
335 i1 = 0
336 i2 = fft%Nx-1
337 IF (j == jdeb) i1 = ideb
338 IF (j == jfin) i2 = ifin
339!dir$ ivdep
340!ocl novrec
341!cdir nodep
342 DO i = i1,i2
343 fft%Work( iwork+k*nmtemp) = scale*x(2*i +2*fft%Ldx1*j+2*fft%Ldx1*fft%Ldx2*k)
344 fft%Work(nwork_temp/4+iwork+k*nmtemp) = isign*scale*x(2*i+1+2*fft%Ldx1*j+2*fft%Ldx1*fft%Ldx2*k)
345 iwork = iwork+1
346 END DO
347 END DO
348 END DO
349
350 ! On fait les T.F. sur la troisieme dimension
351 ioff = 0
352 CALL jmccm1d(nmtemp,fft%Nz,fact,100,mfact,fft%Table,fft%TableSize,100+2*(fft%Nx+fft%Ny),fft%Work,nwork_temp,ioff)
353
354 ! On recopie dans le tableau d'arrivee
355 DO k = 0,fft%Nz-1
356 iwork = 0
357 DO j = jdeb,jfin
358 i1 = 0
359 i2 = fft%Nx-1
360 IF (j == jdeb) i1 = ideb
361 IF (j == jfin) i2 = ifin
362!dir$ ivdep
363!ocl novrec
364!cdir nodep
365 DO i = i1,i2
366 y(2*i +2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k) = fft%Work(ioff +iwork+k*nmtemp)
367 y(2*i+1+2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k) = fft%Work(ioff+nwork_temp/4+iwork+k*nmtemp)
368 iwork = iwork+1
369 END DO
370 END DO
371 END DO
372
373 END DO
374
375 ! On fait les T.F. sur la deuxieme dimension en tronconnant sur la premiere
376 ! et la troisieme
377 debut = .true.
378 fini = .false.
379 DO WHILE (.NOT.fini)
380
381 ! Tronconnage
382 CALL jmdecoup3(fft%Nx,fft%Nz,4*fft%Ny,fft%WorkSize,debut,.true.,ideb,ifin,kdeb,kfin,nltemp,nwork_temp,fini)
383 debut = .false.
384
385 ! On copie le tableau d'entree dans le tableau de travail
386 ! On prend garde a la gestion des extremites
387 DO j = 0,fft%Ny-1
388 iwork = 0
389 DO k = kdeb,kfin
390 i1 = 0
391 i2 = fft%Nx-1
392 IF (k == kdeb) i1 = ideb
393 IF (k == kfin) i2 = ifin
394!dir$ ivdep
395!ocl novrec
396!cdir nodep
397 DO i = i1,i2
398 fft%Work( iwork+j*nltemp) = y(2*i +2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k)
399 fft%Work(nwork_temp/4+iwork+j*nltemp) = y(2*i+1+2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k)
400 iwork = iwork+1
401 END DO
402 END DO
403 END DO
404
405 ! On fait les T.F. sur la deuxieme dimension
406 ioff = 0
407 CALL jmccm1d(nltemp,fft%Ny,fact,100,nfact,fft%table,fft%TableSize,100+2*fft%Nx ,fft%Work,nwork_temp,ioff)
408
409 ! On recopie dans le tableau d'arrivee
410 DO j = 0,fft%Ny-1
411 iwork = 0
412 DO k = kdeb,kfin
413 i1 = 0
414 i2 = fft%Nx-1
415 IF (k == kdeb) i1 = ideb
416 IF (k == kfin) i2 = ifin
417!dir$ ivdep
418!ocl novrec
419!cdir nodep
420 DO i = i1,i2
421 y(2*i +2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k) = fft%Work(ioff +iwork+j*nltemp)
422 y(2*i+1+2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k) = fft%Work(ioff+nwork_temp/4+iwork+j*nltemp)
423 iwork = iwork+1
424 END DO
425 END DO
426 END DO
427
428 END DO
429
430 ! On fait les T.F. sur la premiere dimension en tronconnant sur la deuxieme
431 ! et la troisieme
432 debut = .true.
433 fini = .false.
434 DO WHILE (.NOT.fini)
435
436 ! Tronconnage
437 CALL jmdecoup3(fft%Ny,fft%Nz,4*fft%Nx,fft%WorkSize,debut,.true.,jdeb,jfin,kdeb,kfin,mltemp,nwork_temp,fini)
438 debut = .false.
439
440 ! On copie le tableau d'entree dans le tableau de travail
441 ! On prend garde a la gestion des extremites
442 DO i = 0,fft%Nx-1
443 iwork = 0
444 DO k = kdeb,kfin
445 j1 = 0
446 j2 = fft%Ny-1
447 IF (k == kdeb) j1 = jdeb
448 IF (k == kfin) j2 = jfin
449!dir$ ivdep
450!ocl novrec
451!cdir nodep
452 DO j = j1,j2
453 fft%Work( iwork+i*mltemp) = y(2*i +2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k)
454 fft%Work(nwork_temp/4+iwork+i*mltemp) = y(2*i+1+2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k)
455 iwork = iwork+1
456 END DO
457 END DO
458 END DO
459
460 ! On fait les T.F. sur la premiere dimension
461 ioff = 0
462 CALL jmccm1d(mltemp,fft%Nx,fact,100,0 ,fft%Table,fft%TableSize,100+0 ,fft%Work,nwork_temp,ioff)
463
464 ! On recopie dans le tableau d'arrivee en tenant compte du signe
465 DO i = 0,fft%Nx-1
466 iwork = 0
467 DO k = kdeb,kfin
468 j1 = 0
469 j2 = fft%Ny-1
470 IF (k == kdeb) j1 = jdeb
471 IF (k == kfin) j2 = jfin
472!dir$ ivdep
473!ocl novrec
474!cdir nodep
475 DO j = j1,j2
476 y(2*i +2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k) = fft%Work(ioff +iwork+i*mltemp)
477 y(2*i+1+2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k) = isign*fft%Work(ioff+nwork_temp/4+iwork+i*mltemp)
478 iwork = iwork+1
479 END DO
480 END DO
481 END DO
482
483 END DO
484
485
486 DO k = 0, fft%Nz-1
487 DO j = 0, fft%Ny-1
488!dir$ ivdep
489!ocl novrec
490!cdir nodep
491 DO i = 0, fft%Nx-1
492 c_out(i,j,k) = cmplx(y(2*i+2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k),y(2*i+1+2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k),kind=gft_prec)
493 END DO
494 END DO
495 END DO
496
497 IF( PRESENT(code) ) code = 0
498END SUBROUTINE gft_do_cc_3d
499!!!!!!
500SUBROUTINE gft_do_cr_1d(FFT, isign, scale, c_in, r_out, code)
501 IMPLICIT NONE
502
503 !... Input dummy Parameters
504 TYPE(gft_rcr), INTENT(inout) :: fft
505 INTEGER, INTENT(in) :: isign
506 REAL(kind=gft_prec), INTENT(in) :: scale
507 COMPLEX(kind=GFT_prec), INTENT(in), DIMENSION(0:) :: c_in
508
509 !... Output dummy Parameters
510 REAL(kind=gft_prec), INTENT(out), DIMENSION(0:) :: r_out
511 INTEGER, OPTIONAL, INTENT(out) :: code
512
513 ! Local variables
514 REAL(kind=gft_prec), DIMENSION(0:2*(SIZE(r_out)/2)+1) :: x
515 INTEGER :: i, nfact
516 INTEGER, DIMENSION(0:99) :: fact
517 INTEGER :: dimx, debx, incx, jumpx
518 INTEGER :: dimy, deby, incy, jumpy
519 CHARACTER(len=*), PARAMETER :: spname = "GFT_do_cr_1d"
520 ! print *,spname
521
522 !... Check conditions
523 IF (.NOT.fft%Init) CALL gft_error(spname//": Error: GFT_init must be called prior to this routine")
524 fft%Ldx1 = SIZE(c_in)
525 fft%Ldy1 = SIZE(r_out)
526 IF ( fft%Ldx1 < fft%Nx/2+1) CALL gft_error(spname//": Error: Size(C_IN) < Nx/2+1")
527 IF ( fft%Ldy1 < fft%Nx) CALL gft_error(spname//": Error: Size(R_OUT) < Nx")
528 IF (isign /=-1 .AND. isign /= 1) CALL gft_error(spname//": Error: sign has an invalid value")
529
530 fft%Isign = isign
531 fft%Scale = scale
532
533 DO i = 0,fft%Nx/2
534 x(2*i) = real(c_in(i),kind=gft_prec)
535 x(2*i+1) = aimag(c_in(i))
536 END DO
537 nfact = nint(fft%Table(0))
538 fact(0:nfact-1) = nint(fft%Table(0:nfact-1))
539
540 !... Do the job (don't forget that R_out is Y)
541 dimx = 2*(fft%Nx/2)+2 ; debx = 0 ; incx = 1 ; jumpx = 0
542 dimy = fft%Nx ; deby = 0 ; incy = 1 ; jumpy = 0
543 CALL jmcsm1dxy(1,fft%Nx,fact,100,0,fft%Table,fft%TableSize,100+0, &
544 fft%Work,fft%WorkSize,x,dimx,debx,incx,jumpy,r_out,dimy, &
545 deby,incy,jumpy,isign,scale)
546
547 IF( PRESENT(code) ) code = 0
548END SUBROUTINE gft_do_cr_1d
549!!!!!!
550SUBROUTINE gft_do_cr_2d(FFT, isign, scale, c_in, r_out, code)
551 IMPLICIT NONE
552
553 !... Input dummy Parameters
554 TYPE(gft_rcr), INTENT(inout) :: FFT
555 INTEGER, INTENT(in) :: isign
556 REAL(kind=gft_prec), INTENT(in) :: scale
557 COMPLEX(kind=GFT_prec), INTENT(in), DIMENSION(0:,0:) :: c_in
558
559 !... Output dummy Parameters
560 REAL(kind=gft_prec), INTENT(out), DIMENSION(0:,0:) :: r_out
561 INTEGER, OPTIONAL, INTENT(out) :: code
562
563 !... Local variables
564 REAL(kind=gft_prec), DIMENSION(0:2*SIZE(c_in)-1) :: x
565 REAL(kind=gft_prec), DIMENSION(0:SIZE(r_out)-1) :: y
566 INTEGER :: i, j, ioff, nfact, mfact
567 INTEGER, DIMENSION(0:99) :: fact
568 INTEGER :: ideb, ifin, jdeb, jfin
569 INTEGER :: n_temp, m_temp, nwork_temp
570 LOGICAL :: debut, fin
571 INTEGER :: dimy, deby, incy, jumpy
572 INTEGER :: signe
573 REAL(kind=gft_prec) :: scale_temp
574 CHARACTER(len=*), PARAMETER :: spname = "GFT_do_cr_2d"
575 ! print *,spname
576
577 !... Check conditions
578 IF (.NOT.fft%Init) CALL gft_error(spname//": Error: GFT_init must be called prior to this routine")
579 fft%Ldx1 = SIZE(c_in ,dim=1)
580 fft%Ldy1 = SIZE(r_out,dim=1)
581 IF ( fft%Ldx1 < fft%Nx/2+1 ) CALL gft_error(spname//": Error: Size(C_IN,dim=1) < Nx/2+1")
582 IF ( fft%Ldy1 < fft%Nx+2 ) CALL gft_error(spname//": Error: Size(R_OUT,dim=1) < Nx+2")
583 IF ( SIZE(r_out) < (fft%Nx+2)*fft%Ny ) CALL gft_error(spname//": Error: Size(R_OUT) < (Nx+2)*Ny")
584 IF ( SIZE(c_in) < (fft%Nx/2+1)*fft%Ny ) CALL gft_error(spname//": Error: Size(C_IN) < (Nx/2+1)*Ny")
585 IF (isign /=-1 .AND. isign /= 1) CALL gft_error(spname//": Error: isign has an invalid value")
586
587 fft%Isign = isign
588 fft%Scale = scale
589
590 DO j = 0, fft%Ny-1
591!dir$ ivdep
592!ocl novrec
593!cdir nodep
594 DO i = 0,fft%Nx/2
595 x(2*i +2*fft%Ldx1*j) = real(c_in(i,j),kind=gft_prec)
596 x(2*i+1+2*fft%Ldx1*j) = aimag(c_in(i,j))
597 END DO
598 END DO
599
600 nfact = nint(fft%Table(0))
601 mfact = nint(fft%Table(nfact)) + nfact
602 fact(0:mfact-1) = nint(fft%Table(0:mfact-1))
603
604 ! On fait les T.F. sur la premiere dimension en tronconnant sur la deuxieme
605 debut = .true.
606 DO
607
608 ! Tronconnage
609 CALL jmdecoup(fft%Nx/2+1,4*fft%Ny,fft%WorkSize,debut,fft%even_y,n_temp,ideb,ifin,nwork_temp,fin)
610
611 ! On copie le tableau d'entree dans le tableau de travail sans permuter
612 ! les dimensions (n en premier) pour faire d'abord la tf sur m
613 ! On en profite pour premultiplier et pour tenir compte du signe
614 DO j = 0,fft%Ny-1
615!dir$ ivdep
616!ocl novrec
617!cdir nodep
618 DO i = ideb,ifin
619 fft%Work( n_temp*j+i-ideb) = scale*x(2*i +2*fft%Ldx1*j)
620 fft%Work(nwork_temp/4+n_temp*j+i-ideb) = isign*scale*x(2*i+1+2*fft%Ldx1*j)
621 END DO
622 END DO
623 ioff = 0
624
625 ! On fait la FFT complexe -> complexe sur la deuxieme dimension (m)
626 CALL jmccm1d(n_temp,fft%Ny,fact,100,nfact,fft%Table,fft%TableSize,100+2*fft%Nx,fft%Work,nwork_temp,ioff)
627
628 ! On recopie dans le tableau d'arrivee
629 DO j = 0,fft%Ny-1
630!dir$ ivdep
631!ocl novrec
632!cdir nodep
633 DO i = ideb,ifin
634 y(2*i +fft%Ldy1*j) = fft%Work(ioff+ n_temp*j+i-ideb)
635 y(2*i+1+fft%Ldy1*j) = fft%Work(ioff+nwork_temp/4+n_temp*j+i-ideb)
636 END DO
637 END DO
638
639 ! A-t-on fini ?
640 IF (fin) THEN
641 EXIT
642 ELSE
643 debut = .false.
644 cycle
645 END IF
646
647 END DO
648
649 ! On fait les T.F. sur l'autre dimension
650 debut = .true.
651 DO
652
653 ! Tronconnage
654 CALL jmdecoup(fft%Ny,2*fft%Nx,fft%WorkSize,debut,fft%even_x,m_temp,jdeb,jfin,nwork_temp,fin)
655
656 ! On fait la FFT complexe -> reel sur le premiere dimension (n)
657 dimy = fft%Ldy1*fft%Ny ; deby = jdeb*fft%Ldy1 ; incy = 1 ; jumpy = fft%Ldy1
658 signe = 1
659 scale_temp = real(1,kind=gft_prec)
660 CALL jmcsm1dxy(m_temp,fft%Nx,fact,100,0,fft%Table,fft%TableSize,100+0, &
661 fft%Work,nwork_temp, y,dimy,deby,incy,jumpy,y,dimy,deby, &
662 incy,jumpy,signe,scale_temp)
663
664 ! A-t-on fini ?
665 IF (fin) THEN
666 EXIT
667 ELSE
668 debut = .false.
669 cycle
670 END IF
671
672 END DO
673
674 DO j = 0,fft%Ny-1
675!dir$ ivdep
676!ocl novrec
677!cdir nodep
678 DO i = 0, fft%Nx-1
679 r_out(i,j) = y(i+fft%Ldy1*j)
680 END DO
681 END DO
682
683 IF( PRESENT(code) ) code = 0
684END SUBROUTINE gft_do_cr_2d
685!!!!!!
686SUBROUTINE gft_do_cr_3d(FFT, isign, scale, c_in, r_out, code)
687 IMPLICIT NONE
688
689 !... Input dummy Parameters
690 TYPE(gft_rcr), INTENT(inout) :: FFT
691 INTEGER, INTENT(in) :: isign
692 REAL(kind=gft_prec), INTENT(in) :: scale
693 COMPLEX(kind=GFT_prec), INTENT(in), DIMENSION(0:,0:,0:) :: c_in
694
695 !... Output dummy Parameters
696 REAL(kind=gft_prec), INTENT(out), DIMENSION(0:,0:,0:) :: r_out
697 INTEGER, OPTIONAL, INTENT(out) :: code
698
699 ! Local Variables
700 REAL(kind=gft_prec), DIMENSION(0:2*SIZE(c_in)-1) :: x
701 REAL(kind=gft_prec), DIMENSION(0:SIZE(r_out)-1) :: y
702 INTEGER :: i, j, k, ioff
703 INTEGER :: nfact, mfact, lfact
704 INTEGER, DIMENSION(0:99) :: fact
705 INTEGER :: ideb, ifin, jdeb, jfin, kdeb, kfin, i1, i2, j1, j2
706 INTEGER :: nltemp, nmtemp, mltemp, nwork_temp, iwork
707 LOGICAL :: debut, fini
708 CHARACTER(len=*), PARAMETER :: spname = "GFT_do_cr_3d"
709
710! print *,spname
711
712 ! Check conditions
713 IF (.NOT.fft%Init) CALL gft_error(spname//": Error: GFT_init must be called prior to this routine")
714 fft%Ldx1 = SIZE(c_in ,dim=1)
715 fft%Ldx2 = SIZE(c_in ,dim=2)
716 fft%Ldy1 = SIZE(r_out,dim=1)
717 fft%Ldy2 = SIZE(r_out,dim=2)
718 IF (fft%Ldx1 < fft%Nx/2+1) CALL gft_error(spname//": Error: Size(C_IN,dim=1) < Nx/2+1")
719 IF (fft%Ldy1 < fft%Nx+2 ) CALL gft_error(spname//": Error: Size(R_OUT,dim=1) < Nx+2")
720 IF (fft%Ldx2 < fft%Ny ) CALL gft_error(spname//": Error: Size(C_IN,dim=2) < Ny")
721 IF (fft%Ldy2 < fft%Ny ) CALL gft_error(spname//": Error: Size(R_OUT,dim=2) < Ny")
722 IF ( SIZE(r_out) < (fft%Nx+2)*fft%Ny*fft%Nz ) &
723 CALL gft_error(spname//": Error: Size(R_OUT) < (Nx+2)*Ny*Nz")
724 IF ( SIZE(c_in) < (fft%Nx/2+1)*fft%Ny*fft%Nz ) &
725 CALL gft_error(spname//": Error: Size(C_IN) < (Nx/2+1)*Ny*Nz")
726 IF (isign /=-1 .AND. isign /= 1) &
727 CALL gft_error(spname//": Error: isign has an invalid value")
728
729 fft%Isign = isign
730 fft%Scale = scale
731
732 DO k = 0, fft%Nz-1
733 DO j = 0, fft%Ny-1
734!dir$ ivdep
735!ocl novrec
736!cdir nodep
737 DO i = 0,fft%Nx/2
738 x(2*i +2*fft%Ldx1*j+2*fft%Ldx1*fft%Ldx2*k) = real(c_in(i,j,k),kind=gft_prec)
739 x(2*i+1+2*fft%Ldx1*j+2*fft%Ldx1*fft%Ldx2*k) = aimag(c_in(i,j,k))
740 END DO
741 END DO
742 END DO
743
744 nfact = nint(fft%Table(0))
745 mfact = nint(fft%Table(nfact)) + nfact
746 lfact = nint(fft%Table(mfact)) + mfact
747 fact(0:lfact-1) = nint(fft%Table(0:lfact-1))
748
749 ! On fait les T.F. sur la troisieme dimension en tronconnant sur la premiere
750 ! et la deuxieme
751 debut = .true.
752 fini = .false.
753 DO WHILE (.NOT.fini)
754
755 ! Tronconnage
756 ! Note : on met npair a .true. car il n'y a pas de restriction dans ce cas
757 CALL jmdecoup3(fft%Nx/2+1,fft%Ny,4*fft%Nz,fft%WorkSize,debut,.true.,ideb,ifin,jdeb,jfin,nmtemp,nwork_temp,fini)
758 debut = .false.
759
760 ! On copie le tableau d'entree dans le tableau de travail
761 ! On en profite pour premultiplier et pour tenir compte du signe
762 ! On prend garde a la gestion des extremites
763 DO k = 0,fft%Nz-1
764 iwork = 0
765 DO j = jdeb,jfin
766 i1 = 0
767 i2 = fft%Nx/2
768 IF (j == jdeb) i1 = ideb
769 IF (j == jfin) i2 = ifin
770!dir$ ivdep
771!ocl novrec
772!cdir nodep
773 DO i = i1,i2
774 fft%Work( iwork+k*nmtemp) = &
775 & scale*x(2*i +2*fft%Ldx1*j+2*fft%Ldx1*fft%Ldx2*k)
776 fft%Work(nwork_temp/4+iwork+k*nmtemp) = &
777 & isign*scale*x(2*i+1+2*fft%Ldx1*j+2*fft%Ldx1*fft%Ldx2*k)
778 iwork = iwork+1
779 END DO
780 END DO
781 END DO
782
783 ! On fait les T.F. sur la troisieme dimension
784 ioff = 0
785 CALL jmccm1d(nmtemp,fft%Nz,fact,100,mfact,fft%Table,fft%TableSize,100+2*(fft%Nx+fft%Ny),fft%Work,nwork_temp,ioff)
786
787 ! On recopie dans le tableau d'arrivee
788 DO k = 0,fft%Nz-1
789 iwork = 0
790 DO j = jdeb,jfin
791 i1 = 0
792 i2 = fft%Nx/2
793 IF (j == jdeb) i1 = ideb
794 IF (j == jfin) i2 = ifin
795!dir$ ivdep
796!ocl novrec
797!cdir nodep
798 DO i = i1,i2
799 y(2*i +fft%Ldy1*j+fft%Ldy1*fft%Ldy2*k) = fft%Work(ioff+ iwork+k*nmtemp)
800 y(2*i+1+fft%Ldy1*j+fft%Ldy1*fft%Ldy2*k) = fft%Work(ioff+nwork_temp/4+iwork+k*nmtemp)
801 iwork = iwork+1
802 END DO
803 END DO
804 END DO
805
806 END DO
807
808 ! On fait les T.F. sur la deuxieme dimension en tronconnant sur la premiere
809 ! et la troisieme
810 debut = .true.
811 fini = .false.
812 DO WHILE (.NOT.fini)
813
814 ! Tronconnage
815 CALL jmdecoup3(fft%Nx/2+1,fft%Nz,4*fft%Ny,fft%WorkSize,debut,.true.,ideb,ifin,kdeb,kfin,nltemp,nwork_temp,fini)
816 debut = .false.
817
818 ! On copie le tableau d'entree dans le tableau de travail
819 ! On prend garde a la gestion des extremites
820 DO j = 0,fft%Ny-1
821 iwork = 0
822 DO k = kdeb,kfin
823 i1 = 0
824 i2 = fft%Nx/2
825 IF (k == kdeb) i1 = ideb
826 IF (k == kfin) i2 = ifin
827!dir$ ivdep
828!ocl novrec
829!cdir nodep
830 DO i = i1,i2
831 fft%Work( iwork+j*nltemp) = &
832 & y(2*i +fft%Ldy1*j+fft%Ldy1*fft%Ldy2*k)
833 fft%Work(nwork_temp/4+iwork+j*nltemp) = &
834 & y(2*i+1+fft%Ldy1*j+fft%Ldy1*fft%Ldy2*k)
835 iwork = iwork+1
836 END DO
837 END DO
838 END DO
839
840 ! On fait les T.F. sur la deuxieme dimension
841 ioff = 0
842 CALL jmccm1d(nltemp,fft%Ny,fact,100,nfact,fft%Table,fft%TableSize,100+2*fft%Nx ,fft%Work,nwork_temp,ioff)
843
844 ! On recopie dans le tableau d'arrivee
845 DO j = 0,fft%Ny-1
846 iwork = 0
847 DO k = kdeb,kfin
848 i1 = 0
849 i2 = fft%Nx/2
850 IF (k == kdeb) i1 = ideb
851 IF (k == kfin) i2 = ifin
852!dir$ ivdep
853!ocl novrec
854!cdir nodep
855 DO i = i1,i2
856 y(2*i +fft%Ldy1*j+fft%Ldy1*fft%Ldy2*k) = fft%Work(ioff +iwork+j*nltemp)
857 y(2*i+1+fft%Ldy1*j+fft%Ldy1*fft%Ldy2*k) = fft%Work(ioff+nwork_temp/4+iwork+j*nltemp)
858 iwork = iwork+1
859 END DO
860 END DO
861 END DO
862
863 END DO
864
865 ! On fait les T.F. sur la premiere dimension en tronconnant sur la deuxieme
866 ! et la troisieme
867 debut = .true.
868 fini = .false.
869 DO WHILE (.NOT.fini)
870
871 ! Tronconnage
872 CALL jmdecoup3(fft%Ny,fft%Nz,4*(fft%Nx/2+1),fft%WorkSize,debut,fft%even_x,jdeb,jfin,kdeb,kfin,mltemp,nwork_temp,fini)
873 debut = .false.
874
875 ! On copie le tableau d'entree dans le tableau de travail
876 ! On prend garde a la gestion des extremites
877 DO i = 0,fft%Nx/2
878 iwork = 0
879 DO k = kdeb,kfin
880 j1 = 0
881 j2 = fft%Ny-1
882 IF (k == kdeb) j1 = jdeb
883 IF (k == kfin) j2 = jfin
884!dir$ ivdep
885!ocl novrec
886!cdir nodep
887 DO j = j1,j2
888 fft%Work( iwork+i*mltemp) = y(2*i +fft%Ldy1*j+fft%Ldy1*fft%Ldy2*k)
889 fft%Work(nwork_temp/4+iwork+i*mltemp) = y(2*i+1+fft%Ldy1*j+fft%Ldy1*fft%Ldy2*k)
890 iwork = iwork+1
891 END DO
892 END DO
893 END DO
894
895 ! On fait les T.F. sur la premiere dimension
896 ioff = 0
897 CALL jmcsm1d(mltemp,fft%Nx,fact,100,0 ,fft%Table,fft%TableSize,100+0 ,fft%Work,nwork_temp,ioff)
898
899 ! On recopie dans le tableau d'arrivee
900 DO i = 0,fft%Nx-1
901 iwork = 0
902 DO k = kdeb,kfin
903 j1 = 0
904 j2 = fft%Ny-1
905 IF (k == kdeb) j1 = jdeb
906 IF (k == kfin) j2 = jfin
907!dir$ ivdep
908!ocl novrec
909!cdir nodep
910 DO j = j1,j2
911 y(i+fft%Ldy1*j+fft%Ldy1*fft%Ldy2*k) = fft%Work(ioff+iwork+i*mltemp)
912 iwork = iwork+1
913 END DO
914 END DO
915 END DO
916
917 END DO
918
919 DO k = 0, fft%Nz-1
920 DO j = 0, fft%Ny-1
921!dir$ ivdep
922!ocl novrec
923!cdir nodep
924 DO i = 0, fft%Nx-1
925 r_out(i,j,k) = y(i+fft%Ldy1*j+fft%Ldy1*fft%Ldy2*k)
926 END DO
927 END DO
928 END DO
929 IF( PRESENT(code) ) code = 0
930END SUBROUTINE gft_do_cr_3d
931!!!!!!
932SUBROUTINE gft_do_rc_1d(FFT, isign, scale, r_in, c_out, code)
933 IMPLICIT NONE
934
935 !... Input dummy Parameters
936 TYPE(gft_rcr), INTENT(inout) :: FFT
937 INTEGER, INTENT(in) :: isign
938 REAL(kind=gft_prec), INTENT(in) :: scale
939 REAL(kind=gft_prec), INTENT(in), DIMENSION(0:) :: r_in
940
941 !... Output dummy Parameters
942 COMPLEX(kind=GFT_prec), INTENT(out), DIMENSION(0:) :: c_out
943 INTEGER, OPTIONAL, INTENT(out) :: code
944
945 !... Local Variables
946 REAL(kind=gft_prec), DIMENSION(0:2*(SIZE(r_in)/2)+1) :: y
947 INTEGER :: i, nfact
948 INTEGER, DIMENSION(0:99) :: fact
949 INTEGER :: dimx, debx, incx, jumpx
950 INTEGER :: dimy, deby, incy, jumpy
951 CHARACTER(len=*), PARAMETER :: spname = "GFT_do_rc_1d"
952! print *,spname
953
954 !... Check conditions
955 IF (.NOT.fft%Init) CALL gft_error(spname//": Error: GFT_init must be called prior to this routine")
956 fft%Ldx1 = SIZE(r_in)
957 fft%Ldy1 = SIZE(c_out)
958 IF ( fft%Ldx1 < fft%Nx) CALL gft_error(spname//": Error: Size(R_IN) < Nx")
959 IF ( fft%Ldy1 < fft%Nx/2+1) CALL gft_error(spname//": Error: Size(C_OUT) < Nx/2+1")
960 IF (isign /=-1 .AND. isign /= 1) CALL gft_error(spname//": Error: sign has invalid value")
961
962 fft%Isign = isign
963 fft%Scale = scale
964
965 nfact = nint(fft%Table(0))
966 fact(0:nfact-1) = nint(fft%Table(0:nfact-1))
967
968 !... Do the job (remember that r_in is x)
969 dimx = fft%Nx ; debx = 0 ; incx = 1 ; jumpx = 0
970 dimy = 2*(fft%Nx/2)+2 ; deby = 0 ; incy = 1 ; jumpy = 0
971 CALL jmscm1dxy(1,fft%Nx,fact,100,0,fft%Table,fft%TableSize,100+0, &
972 fft%Work,fft%WorkSize,r_in,dimx,debx,incx,jumpx,y, &
973 dimy,deby,incy,jumpy,isign,scale)
974!dir$ ivdep
975!ocl novrec
976!cdir nodep
977 DO i = 0, fft%Nx/2
978 c_out(i) = cmplx(y(2*i), y(2*i+1),kind=gft_prec)
979 END DO
980
981 IF( PRESENT(code) ) code = 0
982END SUBROUTINE gft_do_rc_1d
983!!!!!!
984SUBROUTINE gft_do_rc_2d(FFT, isign, scale, r_in, c_out, code)
985 IMPLICIT NONE
986
987 !... Input dummy Parameters
988 TYPE(gft_rcr), INTENT(inout) :: FFT
989 INTEGER, INTENT(in) :: isign
990 REAL(kind=gft_prec), INTENT(in) :: scale
991 REAL(kind=gft_prec), INTENT(in), DIMENSION(0:,0:) :: r_in
992
993 !... Output dummy Parameters
994 COMPLEX(kind=GFT_prec), INTENT(out), DIMENSION(0:,0:) :: c_out
995 INTEGER, OPTIONAL, INTENT(out) :: code
996
997 !... Local Variables
998 REAL(kind=gft_prec), DIMENSION(0:SIZE(r_in)-1) :: x
999 REAL(kind=gft_prec), DIMENSION(0:2*SIZE(c_out)-1) :: y
1000 INTEGER :: i, j, ioff, nfact, mfact
1001 INTEGER, DIMENSION(0:99) :: fact
1002 INTEGER :: ideb, ifin, jdeb, jfin, n_temp, m_temp, nwork_temp
1003 LOGICAL :: debut, fin
1004 INTEGER :: dimx, debx, incx, jumpx
1005 INTEGER :: dimy, deby, incy, jumpy
1006 INTEGER :: signe
1007 CHARACTER(len=*), PARAMETER :: spname = "GFT_do_rc_2d"
1008
1009! print *,spname
1010
1011! print *, "in GFT_do_rc_2d"
1012
1013 !... Check conditions
1014 IF (.NOT.fft%Init) CALL gft_error(spname//": Error: GFT_init must be called prior to this routine")
1015 fft%Ldx1 = SIZE(r_in ,dim=1)
1016 fft%Ldy1 = SIZE(c_out,dim=1)
1017 IF (fft%Ldx1 < fft%Nx ) CALL gft_error(spname//": Error: Size(R_IN,dim=1) < Nx")
1018 IF (fft%Ldy1 < fft%Nx/2+1) CALL gft_error(spname//": Error: Size(C_OUT,dim=1) < Nx/2+1")
1019 IF (SIZE(r_in) < fft%Nx*fft%Ny) CALL gft_error(spname//": Error: Size(R_IN) < Nx*Ny")
1020 IF (SIZE(c_out) < (fft%Nx/2+1)*fft%Ny) CALL gft_error(spname//": Error: Size(C_OUT) < (Nx/2+1)*Ny")
1021
1022 nfact = nint(fft%Table(0))
1023 mfact = nint(fft%Table(nfact)) + nfact
1024 fact(0:mfact-1) = nint(fft%Table(0:mfact-1))
1025
1026 DO j=0, fft%Ny-1
1027!dir$ ivdep
1028!ocl novrec
1029!cdir nodep
1030 DO i = 0, fft%Nx-1
1031 x(i+fft%Ldx1*j) = r_in(i,j)
1032 END DO
1033 END DO
1034
1035 ! On fait les T.F. sur la premiere dimension en tronconnant sur la deuxieme
1036 debut = .true.
1037 DO
1038
1039 ! Tronconnage
1040 CALL jmdecoup(fft%Ny,2*fft%Nx,fft%WorkSize,debut,fft%even_x,m_temp,jdeb,jfin,nwork_temp,fin)
1041
1042 ! On fait les T.F. reelles sur la premiere dimension
1043 ! Tout se passe comme si on faisait une T.F. 1D multiple (d'ordre m)
1044 dimx = fft%Ldx1*fft%Ny ; debx = jdeb*fft%Ldx1 ; incx = 1 ; jumpx = fft%Ldx1
1045 dimy = 2*fft%Ldy1*fft%Ny ; deby = jdeb*2*fft%Ldy1 ; incy = 1 ; jumpy = 2*fft%Ldy1
1046 signe = 1
1047 CALL jmscm1dxy(m_temp,fft%Nx,fact,100,0,fft%Table,fft%TableSize,100+0, &
1048 fft%Work,nwork_temp, x,dimx,debx,incx,jumpx,y,dimy,deby,incy,jumpy, &
1049 signe,scale)
1050
1051 ! A-t-on fini ?
1052 IF (fin) THEN
1053 EXIT
1054 ELSE
1055 debut = .false.
1056 cycle
1057 END IF
1058
1059 END DO
1060
1061 ! On fait les T.F. sur l'autre dimension
1062 debut = .true.
1063 DO
1064
1065 ! Tronconnage
1066 CALL jmdecoup(fft%Nx/2+1,4*fft%Ny,fft%WorkSize,debut,fft%even_y,n_temp,ideb,ifin,nwork_temp,fin)
1067
1068 ! On copie
1069 DO j = 0,fft%Ny-1
1070!dir$ ivdep
1071!ocl novrec
1072!cdir nodep
1073 DO i = ideb,ifin
1074 fft%Work( n_temp*j+i-ideb) = y(2*i +2*fft%Ldy1*j)
1075 fft%Work(nwork_temp/4+n_temp*j+i-ideb) = y(2*i+1+2*fft%Ldy1*j)
1076 END DO
1077 END DO
1078 ioff = 0
1079
1080 ! On fait les T.F. sur l'autre dimension (divisee par deux bien sur)
1081 ! On va chercher l'autre table des cosinus
1082 CALL jmccm1d(n_temp,fft%Ny,fact,100,nfact,fft%Table,fft%TableSize,100+2*fft%Nx,fft%Work,nwork_temp,ioff)
1083
1084 ! On recopie dans le tableau d'arrivee
1085 DO j = 0,fft%Ny-1
1086!dir$ ivdep
1087!ocl novrec
1088!cdir nodep
1089 DO i = ideb,ifin
1090 y(2*i +2*fft%Ldy1*j) = fft%Work(ioff +n_temp*j+i-ideb)
1091 y(2*i+1+2*fft%Ldy1*j) = isign*fft%Work(ioff+nwork_temp/4+n_temp*j+i-ideb)
1092 END DO
1093 END DO
1094
1095 ! A-t-on fini ?
1096 IF (fin) THEN
1097 EXIT
1098 ELSE
1099 debut = .false.
1100 cycle
1101 END IF
1102
1103 END DO
1104
1105 DO j=0, fft%Ny-1
1106!dir$ ivdep
1107!ocl novrec
1108!cdir nodep
1109 DO i = 0, fft%Nx/2
1110 c_out(i,j) = cmplx(y(2*i+2*fft%Ldy1*j), y(2*i+1+2*fft%Ldy1*j),kind=gft_prec)
1111 END DO
1112 END DO
1113
1114 IF( PRESENT(code) ) code = 0
1115END SUBROUTINE gft_do_rc_2d
1116!!!!!!
1117SUBROUTINE gft_do_rc_3d(FFT, isign, scale, r_in, c_out, code)
1118 IMPLICIT NONE
1119
1120 !... Input dummy Parameters
1121 TYPE(gft_rcr), INTENT(inout) :: FFT
1122 INTEGER, INTENT(in) :: isign
1123 REAL(kind=gft_prec), INTENT(in) :: scale
1124 REAL(kind=gft_prec), INTENT(in), DIMENSION(0:,0:,0:) :: r_in
1125
1126 !... Output dummy Parameters
1127 COMPLEX(kind=GFT_prec), INTENT(out), DIMENSION(0:,0:,0:) :: c_out
1128 INTEGER, OPTIONAL, INTENT(out) :: code
1129
1130 !... Local Variables
1131 REAL(kind=gft_prec), DIMENSION(0:SIZE(r_in)-1) :: x
1132 REAL(kind=gft_prec), DIMENSION(0:2*SIZE(c_out)-1) :: y
1133 INTEGER :: i, j, k, ioff, nfact, mfact, lfact
1134 INTEGER, DIMENSION(0:99) :: fact
1135 INTEGER :: ideb, ifin, jdeb, jfin, kdeb, kfin, i1, i2, j1, j2
1136 INTEGER :: nltemp, nmtemp, mltemp, nwork_temp, iwork
1137 LOGICAL :: debut, fini
1138 CHARACTER(len=*), PARAMETER :: spname = "GFT_do_rc_3d"
1139 ! print *,spname
1140
1141 !... Check conditions
1142 IF (.NOT.fft%Init) &
1143 CALL gft_error(spname//": Error: GFT_init must be called prior to this routine")
1144 fft%Ldx1 = SIZE(r_in ,dim=1)
1145 fft%Ldx2 = SIZE(r_in ,dim=2)
1146 fft%Ldy1 = SIZE(c_out,dim=1)
1147 fft%Ldy2 = SIZE(c_out,dim=2)
1148 IF (fft%Ldx1 < fft%Nx ) CALL gft_error(spname//": Error: Size(R_IN,dim=1) < Nx")
1149 IF (fft%Ldx2 < fft%Ny ) CALL gft_error(spname//": Error: Size(R_IN,dim=2) < Ny")
1150 IF (fft%Ldy1 < fft%Nx/2+1) CALL gft_error(spname//": Error: Size(C_OUT,dim=1) < Nx/2+1")
1151 IF (fft%Ldy2 < fft%Ny) CALL gft_error(spname//": Error: Size(C_OUT,dim=2) < Ny")
1152 IF ( SIZE(r_in) < fft%Nx*fft%Ny*fft%Nz ) &
1153 CALL gft_error(spname//": Error: Size(R_IN) < Nx*Ny*Nz")
1154 IF ( SIZE(c_out) < (fft%Nx/2+1)*fft%Ny*fft%Nz ) &
1155 CALL gft_error(spname//": Error: Size(C_OUT) < (Nx/2+1)*Ny*Nz")
1156 IF (isign /=-1 .AND. isign /= 1) &
1157 CALL gft_error(spname//": Error: isign has an invalid value")
1158
1159 fft%Isign = isign
1160 fft%Scale = scale
1161
1162 DO k = 0, fft%Nz-1
1163 DO j = 0, fft%Ny-1
1164!dir$ ivdep
1165!ocl novrec
1166!cdir nodep
1167 DO i = 0,fft%Nx-1
1168 x(i+fft%Ldx1*j+fft%Ldx1*fft%Ldx2*k) = r_in(i,j,k)
1169 END DO
1170 END DO
1171 END DO
1172
1173 nfact = nint(fft%Table(0))
1174 mfact = nint(fft%Table(nfact)) + nfact
1175 lfact = nint(fft%Table(mfact)) + mfact
1176 fact(0:lfact-1) = nint(fft%Table(0:lfact-1))
1177
1178 ! On fait les T.F. sur la premiere dimension en tronconnant sur la deuxieme
1179 ! et la troisieme
1180 debut = .true.
1181 fini = .false.
1182 DO WHILE (.NOT.fini)
1183
1184 ! Tronconnage
1185 CALL jmdecoup3(fft%Ny,fft%Nz,4*(fft%Nx/2+1),fft%WorkSize,debut,fft%even_x,jdeb,jfin,kdeb,kfin,mltemp,nwork_temp,fini)
1186 debut = .false.
1187
1188 ! On copie le tableau d'entree dans le tableau de travail
1189 ! On prend garde a la gestion des extremites
1190 DO i = 0,fft%Nx-1
1191 iwork = 0
1192 DO k = kdeb,kfin
1193 j1 = 0
1194 j2 = fft%Ny-1
1195 IF (k == kdeb) j1 = jdeb
1196 IF (k == kfin) j2 = jfin
1197!dir$ ivdep
1198!ocl novrec
1199!cdir nodep
1200 DO j = j1,j2
1201 fft%Work(iwork+i*mltemp) = scale*x(i+fft%Ldx1*j+fft%Ldx1*fft%Ldx2*k)
1202 iwork = iwork+1
1203 END DO
1204 END DO
1205 END DO
1206
1207 ! On fait les T.F. sur la premiere dimension
1208 ioff = 0
1209 CALL jmscm1d(mltemp,fft%Nx,fact,100,0 ,fft%Table,fft%TableSize,100+0 ,fft%Work,nwork_temp,ioff)
1210
1211 ! On recopie dans le tableau d'arrivee
1212 DO i = 0,fft%Nx/2
1213 iwork = 0
1214 DO k = kdeb,kfin
1215 j1 = 0
1216 j2 = fft%Ny-1
1217 IF (k == kdeb) j1 = jdeb
1218 IF (k == kfin) j2 = jfin
1219!dir$ ivdep
1220!ocl novrec
1221!cdir nodep
1222 DO j = j1,j2
1223 y(2*i +2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k) = fft%Work(ioff+ iwork+i*mltemp)
1224 y(2*i+1+2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k) = fft%Work(ioff+nwork_temp/4+iwork+i*mltemp)
1225 iwork = iwork+1
1226 END DO
1227 END DO
1228 END DO
1229
1230 END DO
1231
1232 ! On fait les T.F. sur la troisieme dimension en tronconnant sur la premiere
1233 ! et la deuxieme
1234 debut = .true.
1235 fini = .false.
1236 DO WHILE (.NOT.fini)
1237
1238 ! Tronconnage
1239 ! Note : on met npair a .true. car il n'y a pas de restriction dans ce cas
1240 CALL jmdecoup3(fft%Nx/2+1,fft%Ny,4*fft%Nz,fft%WorkSize,debut,.true.,ideb,ifin,jdeb,jfin,nmtemp,nwork_temp,fini)
1241 debut = .false.
1242
1243 ! On copie le tableau d'entree dans le tableau de travail
1244 ! On en profite pour premultiplier et pour tenir compte du signe
1245 ! On prend garde a la gestion des extremites
1246 DO k = 0,fft%Nz-1
1247 iwork = 0
1248 DO j = jdeb,jfin
1249 i1 = 0
1250 i2 = fft%Nx/2
1251 IF (j == jdeb) i1 = ideb
1252 IF (j == jfin) i2 = ifin
1253!dir$ ivdep
1254!ocl novrec
1255!cdir nodep
1256 DO i = i1,i2
1257 fft%work( iwork+k*nmtemp) = y(2*i +2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k)
1258 fft%Work(nwork_temp/4+iwork+k*nmtemp) = y(2*i+1+2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k)
1259 iwork = iwork+1
1260 END DO
1261 END DO
1262 END DO
1263
1264 ! On fait les T.F. sur la troisieme dimension
1265 ioff = 0
1266 CALL jmccm1d(nmtemp,fft%Nz,fact,100,mfact,fft%Table,fft%TableSize,100+2*(fft%Nx+fft%Ny),fft%Work,nwork_temp,ioff)
1267
1268 ! On recopie dans le tableau d'arrivee
1269 DO k = 0,fft%Nz-1
1270 iwork = 0
1271 DO j = jdeb,jfin
1272 i1 = 0
1273 i2 = fft%Nx/2
1274 IF (j == jdeb) i1 = ideb
1275 IF (j == jfin) i2 = ifin
1276!dir$ ivdep
1277!ocl novrec
1278!cdir nodep
1279 DO i = i1,i2
1280 y(2*i +2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k) = fft%work(ioff+ iwork+k*nmtemp)
1281 y(2*i+1+2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k) = fft%work(ioff+nwork_temp/4+iwork+k*nmtemp)
1282 iwork = iwork+1
1283 END DO
1284 END DO
1285 END DO
1286
1287 END DO
1288
1289 ! On fait les T.F. sur la deuxieme dimension en tronconnant sur la premiere
1290 ! et la troisieme
1291 debut = .true.
1292 fini = .false.
1293 DO WHILE (.NOT.fini)
1294
1295 ! Tronconnage
1296 CALL jmdecoup3(fft%Nx/2+1,fft%Nz,4*fft%Ny,fft%WorkSize,debut,.true.,ideb,ifin,kdeb,kfin,nltemp,nwork_temp,fini)
1297 debut = .false.
1298
1299 ! On copie le tableau d'entree dans le tableau de travail
1300 ! On prend garde a la gestion des extremites
1301 DO j = 0,fft%Ny-1
1302 iwork = 0
1303 DO k = kdeb,kfin
1304 i1 = 0
1305 i2 = fft%Nx/2
1306 IF (k == kdeb) i1 = ideb
1307 IF (k == kfin) i2 = ifin
1308!dir$ ivdep
1309!ocl novrec
1310!cdir nodep
1311 DO i = i1,i2
1312 fft%work( iwork+j*nltemp) = y(2*i +2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k)
1313 fft%work(nwork_temp/4+iwork+j*nltemp) = y(2*i+1+2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k)
1314 iwork = iwork+1
1315 END DO
1316 END DO
1317 END DO
1318
1319 ! On fait les T.F. sur la deuxieme dimension
1320 ioff = 0
1321 CALL jmccm1d(nltemp,fft%Ny,fact,100,nfact,fft%Table,fft%TableSize,100+2*fft%Nx ,fft%Work,nwork_temp,ioff)
1322
1323 ! On recopie dans le tableau d'arrivee
1324 DO j = 0,fft%Ny-1
1325 iwork = 0
1326 DO k = kdeb,kfin
1327 i1 = 0
1328 i2 = fft%Nx/2
1329 IF (k == kdeb) i1 = ideb
1330 IF (k == kfin) i2 = ifin
1331!dir$ ivdep
1332!ocl novrec
1333!cdir nodep
1334 DO i = i1,i2
1335 y(2*i +2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k) = &
1336 & fft%Work(ioff +iwork+j*nltemp)
1337 y(2*i+1+2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k) = &
1338 & isign*fft%Work(ioff+nwork_temp/4+iwork+j*nltemp)
1339 iwork = iwork+1
1340 END DO
1341 END DO
1342 END DO
1343
1344 END DO
1345
1346 DO k = 0, fft%Nz-1
1347 DO j = 0, fft%Ny-1
1348!dir$ ivdep
1349!ocl novrec
1350!cdir nodep
1351 DO i = 0, fft%Nx/2
1352 c_out(i,j,k) = cmplx(y(2*i+2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k), &
1353 y(2*i+1+2*fft%Ldy1*j+2*fft%Ldy1*fft%Ldy2*k),kind=gft_prec)
1354 END DO
1355 END DO
1356 END DO
1357
1358 IF( PRESENT(code) ) code = 0
1359END SUBROUTINE gft_do_rc_3d
1360!!!!!!
1361SUBROUTINE gft_do_mcc_1d(FFT, isign, scale, c_in, c_out, code)
1362 IMPLICIT NONE
1363
1364 !... Input dummy Parameters
1365 TYPE(gft_mcc), INTENT(inout) :: FFT
1366 INTEGER, INTENT(in) :: isign
1367 REAL(kind=gft_prec), INTENT(in) :: scale
1368 COMPLEX(kind=GFT_prec), INTENT(in), DIMENSION(0:,0:) :: c_in
1369
1370 !... Output dummy Parameters
1371 COMPLEX(kind=GFT_prec), INTENT(out), DIMENSION(0:,0:) :: c_out
1372 INTEGER, OPTIONAL, INTENT(out) :: code
1373
1374
1375 !... Local variables
1376 REAL(kind=gft_prec), DIMENSION(0:2*SIZE(c_in)-1) :: x
1377 REAL(kind=gft_prec), DIMENSION(0:2*SIZE(c_out)-1) :: y
1378 INTEGER :: i, j, ioff, nfact
1379 INTEGER, DIMENSION(0:99) :: fact
1380 CHARACTER(len=*), PARAMETER :: spname = 'GFT_do_mcc_1d'
1381 ! print *,spname
1382
1383 !... Check conditions
1384 IF (.NOT.fft%Init) CALL gft_error(spname//": Error: GFT_init must be called prior to this routine")
1385 fft%Ldx1 = SIZE(c_in ,dim=1)
1386 fft%Ldy1 = SIZE(c_out,dim=1)
1387 IF ( SIZE(c_in) < fft%Nx*fft%Ny ) CALL gft_error(spname//": Error: Size(C_IN) < Nx*Ny")
1388 IF ( SIZE(c_out) < fft%Nx*fft%Ny ) CALL gft_error(spname//": Error: Size(C_OUT) < Nx*Ny")
1389 IF ( fft%Ldx1 < fft%Nx ) CALL gft_error(spname//": Error: Size(C_IN,dim=1) < Nx")
1390 IF ( fft%Ldy1 < fft%Nx ) CALL gft_error(spname//": Error: Size(C_OUT,dim=1) < Nx")
1391 IF (isign /=-1 .AND. isign /= 1) CALL gft_error(spname//": Error: isign value is invalid")
1392
1393 fft%Isign = isign
1394 fft%Scale = scale
1395
1396 DO j = 0, fft%Ny-1
1397!dir$ ivdep
1398!ocl novrec
1399!cdir nodep
1400 DO i = 0,fft%Nx-1
1401 x(2*i +j*2*fft%Ldx1) = real(c_in(i,j),kind=gft_prec)
1402 x(2*i+1+j*2*fft%Ldx1) = aimag(c_in(i,j))
1403 END DO
1404 END DO
1405
1406 nfact = nint(fft%table(0))
1407 fact(0:nfact-1) = nint(fft%table(0:nfact-1))
1408
1409 ! On copie le tableau d'entree dans le tableau de travail
1410 ! On en profite pour premultiplier et pour tenir compte du signe
1411 DO i = 0,fft%Nx-1
1412!dir$ ivdep
1413!ocl novrec
1414!cdir nodep
1415 DO j = 0,fft%Ny-1
1416 fft%Work(j+fft%Ny*i) = scale*x(2*i +2*fft%Ldx1*j)
1417 fft%Work(j+fft%Ny*(fft%Nx+i)) = isign*scale*x(2*i+1+2*fft%Ldx1*j)
1418 END DO
1419 END DO
1420
1421 ! On appelle le sous-programme principal
1422 ioff = 0
1423 CALL jmccm1d(fft%Ny, fft%Nx,fact,100,0,fft%Table,fft%TableSize,100+0,fft%Work,fft%WorkSize,ioff)
1424
1425 ! On recopie le tableau de travail dans le tableau de sortie
1426 DO i = 0,fft%Nx-1
1427!dir$ ivdep
1428!ocl novrec
1429!cdir nodep
1430 DO j = 0,fft%Ny-1
1431 y(2*i +2*fft%Ldy1*j) = fft%work(ioff+j+fft%Ny*i)
1432 y(2*i+1+2*fft%Ldy1*j) = isign*fft%work(ioff+j+fft%Ny*(fft%Nx+i))
1433 END DO
1434 END DO
1435
1436 DO j = 0, fft%Ny-1
1437 DO i = 0, fft%Nx-1
1438 c_out(i,j) = cmplx(y(2*i+2*fft%Ldy1*j), y(2*i+1+2*fft%Ldy1*j),kind=gft_prec)
1439 END DO
1440 END DO
1441
1442 IF( PRESENT(code) ) code = 0
1443END SUBROUTINE gft_do_mcc_1d
1444!!!!!!
1445SUBROUTINE gft_do_mcr_1d(FFT, isign, scale, c_in, r_out, code)
1446 IMPLICIT NONE
1447
1448 !... Input dummy Parameters
1449 TYPE(gft_mrcr), INTENT(inout) :: FFT
1450 INTEGER, INTENT(in) :: isign
1451 REAL(kind=gft_prec), INTENT(in) :: scale
1452 COMPLEX(kind=GFT_prec), INTENT(in), DIMENSION(0:,0:) :: c_in
1453
1454 !... Output dummy Parameters
1455 REAL(kind=gft_prec), INTENT(out), DIMENSION(0:,0:) :: r_out
1456 INTEGER, OPTIONAL, INTENT(out) :: code
1457
1458 !... Local variables
1459 REAL(kind=gft_prec), DIMENSION(0:2*SIZE(c_in)-1) :: x
1460 REAL(kind=gft_prec), DIMENSION(0:SIZE(r_out)-1) :: y
1461 INTEGER :: i, j, nfact
1462 INTEGER, DIMENSION(0:99) :: fact
1463 INTEGER :: dimx, debx, incx, jumpx
1464 INTEGER :: dimy, deby, incy, jumpy
1465 CHARACTER(len=*), PARAMETER :: spname = "GFT_do_mcr_1d"
1466! print *,spname
1467
1468 !... Check conditions
1469 IF (.NOT.fft%Init) CALL gft_error(spname//": Error: GFT_init must be called prior to this routine")
1470 fft%Ldx1 = SIZE(c_in ,dim=1)
1471 fft%Ldy1 = SIZE(r_out,dim=1)
1472 IF ( fft%Ldx1 < fft%Nx/2+1 ) CALL gft_error(spname//": Error: Size(C_IN,dim=1) < Nx/2+1")
1473 IF ( fft%Ldy1 < fft%Nx ) CALL gft_error(spname//": Error: Size(R_OUT,dim=1) < Nx")
1474 IF ( SIZE(r_out) < fft%Nx*fft%Ny ) CALL gft_error(spname//": Error: Size(R_OUT) < Nx*Ny")
1475 IF ( SIZE(c_in) < (fft%Nx/2+1)*fft%Ny ) CALL gft_error(spname//": Error: Size(C_IN) < (Nx/2+1)*Ny")
1476 IF (isign /=-1 .AND. isign /= 1) CALL gft_error(spname//": Error: isign has an invalid value")
1477
1478 fft%Isign = isign
1479 fft%Scale = scale
1480
1481 DO j = 0, fft%Ny-1
1482!dir$ ivdep
1483!ocl novrec
1484!cdir nodep
1485 DO i = 0,fft%Nx/2
1486 x(2*i +2*fft%Ldx1*j) = real(c_in(i,j),kind=gft_prec)
1487 x(2*i+1+2*fft%Ldx1*j) = aimag(c_in(i,j))
1488 END DO
1489 END DO
1490
1491 nfact = nint(fft%Table(0))
1492 fact(0:nfact-1) = nint(fft%Table(0:nfact-1))
1493
1494 ! On appelle le sous-programme principal
1495 dimx = 2*fft%Ldx1*fft%Ny ; debx = 0 ; incx = 1 ; jumpx = 2*fft%Ldx1
1496 dimy = fft%Ldy1*fft%Ny ; deby = 0 ; incy = 1 ; jumpy = fft%Ldy1
1497 CALL jmcsm1dxy(fft%Ny,fft%Nx,fact,100,0,fft%Table,fft%TableSize,100+0, &
1498 fft%Work,fft%WorkSize,x,dimx,debx,incx,jumpx,y,dimy,deby,&
1499 incy,jumpy,isign,scale)
1500
1501 DO j = 0, fft%Ny-1
1502 DO i = 0, fft%Nx-1
1503 r_out(i,j) = y(i + fft%Ldy1*j)
1504 END DO
1505 END DO
1506
1507 IF( PRESENT(code) ) code = 0
1508END SUBROUTINE gft_do_mcr_1d
1509
1510SUBROUTINE gft_do_mrc_1d(FFT, isign, scale, r_in, c_out, code)
1511 IMPLICIT NONE
1512
1513 !... Input dummy Parameters
1514 TYPE(gft_mrcr), INTENT(inout) :: FFT
1515 INTEGER, INTENT(in) :: isign
1516 REAL(kind=gft_prec), INTENT(in) :: scale
1517 REAL(kind=gft_prec), INTENT(in), DIMENSION(0:,0:) :: r_in
1518
1519 !... Output dummy Parameters
1520 COMPLEX(kind=GFT_prec), INTENT(out), DIMENSION(0:,0:) :: c_out
1521 INTEGER, OPTIONAL, INTENT(out) :: code
1522
1523 !... Local Variables
1524 REAL(kind=gft_prec), DIMENSION(0:SIZE(r_in)-1) :: x
1525 REAL(kind=gft_prec), DIMENSION(0:2*SIZE(c_out)-1) :: y
1526 INTEGER :: i, j, nfact
1527 INTEGER, DIMENSION(0:99) :: fact
1528 INTEGER :: dimx, debx, incx, jumpx
1529 INTEGER :: dimy, deby, incy, jumpy
1530 CHARACTER(len=*), PARAMETER :: spname = "GFT_do_mrc_1d"
1531! print *,spname
1532
1533 !... Check conditions
1534 IF (.NOT.fft%Init) CALL gft_error(spname//": Error: GFT_init must be called prior to this routine")
1535 fft%Ldx1 = SIZE(r_in ,dim=1)
1536 fft%Ldy1 = SIZE(c_out,dim=1)
1537 IF (fft%Ldx1 < fft%Nx ) CALL gft_error(spname//": Error: Size(R_IN,dim=1) < Nx")
1538 IF (fft%Ldy1 < fft%Nx/2+1) CALL gft_error(spname//": Error: Size(C_OUT,dim=1) < Nx/2+1")
1539 IF (SIZE(r_in) < fft%Nx*fft%Ny) CALL gft_error(spname//": Error: Size(R_IN) < Nx*Ny")
1540 IF (SIZE(c_out) < (fft%Nx/2+1)*fft%Ny) CALL gft_error(spname//": Error: Size(C_OUT) < (Nx/2+1)*Ny")
1541
1542 DO j=0, fft%Ny-1
1543!dir$ ivdep
1544!ocl novrec
1545!cdir nodep
1546 DO i = 0, fft%Nx-1
1547 x(i+fft%Ldx1*j) = r_in(i,j)
1548 END DO
1549 END DO
1550
1551 nfact = nint(fft%Table(0))
1552 fact(0:nfact-1) = nint(fft%Table(0:nfact-1))
1553
1554 ! On appelle le sous-programme principal
1555 dimx = fft%Ldx1*fft%Ny ; debx = 0 ; incx = 1 ; jumpx = fft%Ldx1
1556 dimy = 2*fft%Ldy1*fft%Ny ; deby = 0 ; incy = 1 ; jumpy = 2*fft%Ldy1
1557 CALL jmscm1dxy(fft%Ny,fft%Nx,fact,100,0,fft%Table,fft%TableSize,100+0, &
1558 fft%Work,fft%WorkSize,x,dimx,debx,incx,jumpx,y,dimy,deby,&
1559 incy,jumpy,isign,scale)
1560
1561 DO j = 0, fft%Ny-1
1562 DO i = 0, fft%Nx/2
1563 c_out(i,j) = cmplx(y(2*i+2*fft%Ldy1*j), y(2*i+1+2*fft%Ldy1*j),kind=gft_prec)
1564 END DO
1565 END DO
1566
1567 IF( PRESENT(code) ) code = 0
1568END SUBROUTINE gft_do_mrc_1d
1569
1570END MODULE gft_do
gft_do::gft_do_fft
Definition GFT_do.f90:27
gft_common
Definition GFT_common.f90:15
gft_common::gft_error
subroutine gft_error(message)
Definition GFT_common.f90:108
gft_common::gft_prec
integer, parameter gft_prec
Definition GFT_common.f90:24
gft_do
Definition GFT_do.f90:15
gft_jmfft
Definition GFT_jmfft.f90:15
gft_jmfft::jmcsm1d
subroutine jmcsm1d(m, n, fact, nfact, ifact, table, ntable, itable, work, nwork, ioff)
Definition GFT_jmfft.f90:1174
gft_jmfft::jmccm1d
subroutine jmccm1d(m, n, fact, nfact, ifact, table, ntable, itable, work, nwork, ioff)
Definition GFT_jmfft.f90:335
gft_jmfft::jmscm1d
subroutine jmscm1d(m, n, fact, nfact, ifact, table, ntable, itable, work, nwork, ioff)
Definition GFT_jmfft.f90:1582
gft_jmfft::jmscm1dxy
subroutine jmscm1dxy(m, n, fact, nfact, ifact, table, ntable, itable, work, nwork, x, dimx, debx, incx, jumpx, y, dimy, deby, incy, jumpy, isign, scale)
Definition GFT_jmfft.f90:1783
gft_jmfft::jmdecoup3
subroutine jmdecoup3(n, m, nmr, nwork, debut, lpair, ideb, ifin, jdeb, jfin, nmtemp, nwork_temp, fini)
Definition GFT_jmfft.f90:274
gft_jmfft::jmdecoup
subroutine jmdecoup(n, nr, nwork, debut, mpair, n_temp, ideb, ifin, nwork_temp, fin)
Definition GFT_jmfft.f90:227
gft_jmfft::jmcsm1dxy
subroutine jmcsm1dxy(m, n, fact, nfact, ifact, table, ntable, itable, work, nwork, x, dimx, debx, incx, jumpx, y, dimy, deby, incy, jumpy, isign, scale)
Definition GFT_jmfft.f90:1379
gft_common::gft_cc
Definition GFT_common.f90:30
gft_common::gft_mcc
Definition GFT_common.f90:72
gft_common::gft_mrcr
Definition GFT_common.f90:88
gft_common::gft_rcr
Definition GFT_common.f90:49