TreeMig Code
Loading...
Searching...
No Matches
GFT_jmfft.f90
Go to the documentation of this file.
1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -*- Mode: F90 -*- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2!! GFT_jmfft.f90 --- J.-M. FFT routines
3!!
4!! Auteur : Jean-Marie Teuler, CNRS-IDRIS
5!! Créé le : Tue Feb 19 10:22:47 2002
6!! Dern. mod. par : Jalel Chergui (CNRS/IDRIS) <Jalel.Chergui@idris.fr>
7!! Dern. mod. le : Thu Jun 13 12:01:48 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, Jean-Marie Teuler.
14!
15MODULE gft_jmfft
16USE gft_common
17
18CONTAINS
19
20SUBROUTINE jmtable(table,ntable,itable,n)
21
22 IMPLICIT NONE
23
24 ! Arguments
25 INTEGER, INTENT(in) :: ntable, itable, n
26 REAL(kind=gft_prec), INTENT(out), DIMENSION(0:ntable-1) :: table
27
28 ! Variables locales
29 REAL(kind=gft_prec), SAVE :: twopi
30 REAL(kind=gft_prec) :: temp, temp1, temp2
31
32 INTEGER :: i
33
34 twopi = 2.0_gft_prec * acos(-1.0_gft_prec)
35
36 ! Calcul des sinus et cosinus
37
38 ! Si n est multiple de 4, astuces en serie
39 IF (mod(n,4) == 0) THEN
40 ! On se debarrasse des cas limite
41 table(itable+ 0) = 1
42 table(itable+n+ 0) = 0
43 table(itable+ n/4) = 0
44 table(itable+n+ n/4) = 1
45 table(itable+ n/2) = -1
46 table(itable+n+ n/2) = 0
47 table(itable+ 3*n/4) = 0
48 table(itable+n+3*n/4) = -1
49 ! Cas general
50!dir$ ivdep
51!ocl novrec
52!cdir nodep
53 DO i = 1,n/4-1
54 temp = cos(twopi*real(i,kind=gft_prec)/real(n,kind=gft_prec))
55 table(itable+ i) = temp
56 table(itable+ n/2-i) = -temp
57 table(itable+ n/2+i) = -temp
58 table(itable+ n-i) = temp
59 table(itable+n+ n/4+i) = temp
60 table(itable+n+ n/4-i) = temp
61 table(itable+n+3*n/4+i) = -temp
62 table(itable+n+3*n/4-i) = -temp
63 END DO
64
65 ! Si n est simplement multiple de 2 (sans etre multiple de 4)
66 ELSE IF (mod(n,2) == 0) THEN
67 ! On se debarrasse des cas limite
68 table(itable+ 0) = 1
69 table(itable+n+ 0) = 0
70 table(itable+ n/2) = -1
71 table(itable+n+n/2) = 0
72 ! Cas general
73!dir$ ivdep
74!ocl novrec
75!cdir nodep
76 DO i = 1,n/2-1
77 temp1 = cos(twopi*real(i,kind=gft_prec)/real(n,kind=gft_prec))
78 table(itable+ i) = temp1
79 table(itable+ n/2+i) = -temp1
80 temp2 = sin(twopi*real(i,kind=gft_prec)/real(n,kind=gft_prec))
81 table(itable+n+ i) = temp2
82 table(itable+n+n/2+i) = -temp2
83 END DO
84
85 ! Si n est impair
86 ELSE
87 ! On se debarrasse des cas limite
88 table(itable+ 0) = 1
89 table(itable+n+0) = 0
90!dir$ ivdep
91!ocl novrec
92!cdir nodep
93 DO i = 1,n/2
94 temp1 = cos(twopi*real(i,kind=gft_prec)/real(n,kind=gft_prec))
95 table(itable+ i) = temp1
96 table(itable+ n-i) = temp1
97 temp2 = sin(twopi*real(i,kind=gft_prec)/real(n,kind=gft_prec))
98 table(itable+n+ i) = temp2
99 table(itable+n+n-i) = -temp2
100 END DO
101
102 END IF
103
104END SUBROUTINE jmtable
105
106SUBROUTINE jmfact(n,fact,nfact,ideb,ifin)
107
108 IMPLICIT NONE
109
110 ! Arguments
111 INTEGER, INTENT(in) :: n, nfact, ideb
112 INTEGER, INTENT(out) :: ifin
113 INTEGER, INTENT(inout), DIMENSION(0:nfact-1) :: fact
114
115 ! Variables locales
116 INTEGER :: m
117 INTEGER :: n2, p2, n3, p3, n5, p5
118 CHARACTER(len=*), PARAMETER :: nomsp = 'JMFACT'
119 ! Nombres premiers
120 INTEGER, PARAMETER :: npremiers = 7
121 INTEGER, DIMENSION(0:npremiers-1) :: premiers = (/7,11,13,17,19,23,29/)
122 INTEGER :: ip, premier, pp, np
123
124 m = n
125
126 ! Etude des puissances de deux
127 p2 = 0
128 n2 = 1
129 DO
130 IF (mod(m,2) == 0) THEN
131 p2 = p2+1
132 n2 = n2*2
133 m = m/2
134 ELSE
135 EXIT
136 END IF
137 END DO
138 ifin = ideb+3
139 IF (ifin > nfact) &
140 & CALL jmerreur2(nomsp,7,nfact,ifin)
141 fact(ideb+1) = n2
142 fact(ideb+2) = p2
143
144 ! Etude des puissances de trois
145 p3 = 0
146 n3 = 1
147 DO
148 IF (mod(m,3) == 0) THEN
149 p3 = p3+1
150 n3 = n3*3
151 m = m/3
152 ELSE
153 EXIT
154 END IF
155 END DO
156 ifin = ifin+2
157 IF (ifin > nfact) &
158 & CALL jmerreur2(nomsp,7,nfact,ifin)
159 fact(ideb+3) = n3
160 fact(ideb+4) = p3
161
162 ! Etude des puissances de cinq
163 p5 = 0
164 n5 = 1
165 DO
166 IF (mod(m,5) == 0) THEN
167 p5 = p5+1
168 n5 = n5*5
169 m = m/5
170 ELSE
171 EXIT
172 END IF
173 END DO
174 ifin = ifin+2
175 IF (ifin > nfact) &
176 & CALL jmerreur2(nomsp,7,nfact,ifin)
177 fact(ideb+5) = n5
178 fact(ideb+6) = p5
179
180 ! On met a jour le nombre de termes
181 fact(ideb) = 7
182
183 ! Si on a fini
184 IF (n2*n3*n5 == n) RETURN
185
186 ! Il reste maintenant des facteurs premiers bizarres
187 ! On va boucler tant qu'on n'a pas fini ou tant qu'on n'a pas epuise la liste
188
189 DO ip = 0,npremiers-1
190
191 premier = premiers(ip)
192
193 pp = 0
194 np = 1
195 DO
196 IF (mod(m,premier) == 0) THEN
197 pp = pp+1
198 np = np*premier
199 m = m/premier
200 ELSE
201 EXIT
202 END IF
203 END DO
204 ifin = ifin+2
205 IF (ifin > nfact) &
206 & CALL jmerreur2(nomsp,7,nfact,ifin)
207 fact(ifin-2) = pp
208 fact(ifin-1) = premier
209 fact(ideb) = fact(ideb) + 2
210
211 ! Si le nombre est completement factorise, inutile de continuer
212 IF (m == 1) EXIT
213
214 END DO
215
216 ! On regarde si la factorisation est terminee
217 IF (m == 1) THEN
218 RETURN
219 ELSE
220 CALL jmerreur1(nomsp,3,n)
221 END IF
222
223END SUBROUTINE jmfact
224
225! Pour tronconner de facon a tenir dans le nwork disponible
226SUBROUTINE jmdecoup(n,nr,nwork,debut,mpair,n_temp,ideb,ifin,nwork_temp,fin)
227
228 IMPLICIT NONE
229
230 ! Arguments
231 INTEGER, INTENT(in) :: n, nr, nwork
232 LOGICAL, INTENT(in) :: debut, mpair
233 INTEGER, INTENT(out) :: n_temp, ideb, nwork_temp
234 INTEGER, INTENT(inout) :: ifin
235 LOGICAL, INTENT(out) :: fin
236
237 ! Variables locales
238 CHARACTER(len=*), PARAMETER :: nomsp = 'JMDECOUP'
239
240 ! n*nr est l'espace total qu'il faudrait pour work.
241 ! Malheureusement, on n'a que nwork au plus
242 ! On va donc decouper n en morceaux pour tenir
243
244 ! Gestion de debut
245 IF (debut) THEN
246 ideb = 0
247 ELSE
248 ideb = ifin+1
249 END IF
250
251 ! Gestion de n_temp et ifin
252 n_temp = nwork/nr
253 ! Si m impair, on doit eviter que n_temp soit impair (routine cs et sc)
254 IF (.NOT.mpair .AND. mod(n_temp,2) /= 0) n_temp = n_temp-1
255 ifin = min(ideb+n_temp-1,n-1)
256 n_temp = ifin-ideb+1
257 ! On verifie que n_temp n'est pas nul
258 IF (n_temp <= 0) THEN
259 CALL jmerreur3(nomsp,6,n,nr,nwork)
260 END IF
261 nwork_temp = n_temp*nr
262
263 ! Gestion de fin
264 IF (ifin == n-1) THEN
265 fin = .true.
266 ELSE
267 fin = .false.
268 END IF
269
270END SUBROUTINE jmdecoup
271
272! Pour tronconner en dimension 3 de facon a tenir dans le nwork disponible
273SUBROUTINE jmdecoup3(n,m,nmr,nwork,debut,lpair,ideb,ifin,jdeb,jfin,nmtemp,nwork_temp,fini)
274
275 IMPLICIT NONE
276
277 ! Arguments
278 INTEGER, INTENT(in) :: n, m, nmr, nwork
279 LOGICAL, INTENT(in) :: debut, lpair
280 INTEGER, INTENT(out) :: nmtemp, nwork_temp
281 INTEGER, INTENT(out) :: ideb, jdeb
282 INTEGER, INTENT(inout) :: ifin, jfin
283 LOGICAL, INTENT(out) :: fini
284
285 ! Variables locales
286 INTEGER :: ijdeb, ijfin
287 CHARACTER(len=*), PARAMETER :: nomsp = 'JMDECOUP3'
288
289 ! n*m*nr est l'espace total qu'il faudrait pour work.
290 ! Malheureusement, on n'a que nwork au plus
291 ! On va donc decouper n et m en morceaux pour tenir
292
293 ! Gestion de debut
294 IF (debut) THEN
295 ideb = 0
296 jdeb = 0
297 ELSE
298 IF (ifin < n-1) THEN
299 ideb = ifin+1
300 jdeb = jfin
301 ELSE
302 ideb = 0
303 jdeb = jfin+1
304 END IF
305 END IF
306
307 ! Gestion de nmtemp
308 nmtemp = nwork/nmr
309 ! Si l impair, on doit eviter que nmtemp soit impair (routine cs et sc)
310 IF (.NOT.lpair .AND. mod(nmtemp,2) /= 0) nmtemp = nmtemp-1
311 ! Pour simplifier, on passe par des indices 2d
312 ijdeb = ideb+jdeb*n
313 ijfin = min(ijdeb+nmtemp-1,n*m-1)
314 nmtemp = ijfin-ijdeb+1
315 ! On verifie que nmtemp n'est pas nul
316 IF (nmtemp <= 0) THEN
317 CALL jmerreur4(nomsp,6,n,m,nmr,nwork)
318 END IF
319 nwork_temp = nmtemp*nmr
320
321 ! On deduit ifin et jfin de ijfin
322 jfin = ijfin/n
323 ifin = ijfin-n*jfin
324
325 ! Gestion de fin
326 IF (ifin == n-1 .AND. jfin == m-1) THEN
327 fini = .true.
328 ELSE
329 fini = .false.
330 END IF
331
332END SUBROUTINE jmdecoup3
333
334SUBROUTINE jmccm1d(m,n,fact,nfact,ifact,table,ntable,itable,work,nwork,ioff)
335
336 IMPLICIT NONE
337
338 ! Arguments
339 INTEGER, INTENT(in) :: m, n
340 INTEGER, INTENT(in) :: nfact, ifact
341 INTEGER, INTENT(in), DIMENSION(0:nfact-1) :: fact
342 INTEGER, INTENT(in) :: ntable,itable
343 REAL(kind=gft_prec), INTENT(in), DIMENSION(0:ntable-1) :: table
344 INTEGER, INTENT(in) :: nwork
345 REAL(kind=gft_prec), INTENT(inout), DIMENSION(0:nwork-1) :: work
346 INTEGER, INTENT(inout) :: ioff
347
348 ! Variables locales
349 INTEGER :: nterms
350 INTEGER :: np, pp, lastnp, premier
351 INTEGER :: nprod, nprod1, nprod2
352 INTEGER :: n2, p2, n3, p3, n5, p5
353 INTEGER :: i
354
355 ! On recupere les facteurs
356 nterms = fact(ifact)
357 n2 = fact(ifact+1)
358 p2 = fact(ifact+2)
359 n3 = fact(ifact+3)
360 p3 = fact(ifact+4)
361 n5 = fact(ifact+5)
362 p5 = fact(ifact+6)
363 nprod = n2*n3*n5
364 DO i = 7,nterms-1,2
365 nprod = nprod*fact(ifact+i+1)**fact(ifact+i)
366 END DO
367
368 ! On fait n3*n5 T.F. de n2 (qui est en puissances de 2)
369 IF (n2 /= 1) THEN
370 CALL jmccm1d2(p2,n2,m*(nprod/n2),table,ntable,itable,n,n/n2,work,nwork,ioff)
371 END IF
372
373 ! On transpose (on tient compte de ioff) en permutant les deux parties
374 ! On en profite pour multiplier par le bon wij
375 IF (n2 /= 1 .AND. nprod /= n2) THEN
376 CALL jmcctranspcs(m,n,n2,nprod/n2,table,ntable,itable,work,nwork,ioff)
377 END IF
378
379 ! On fait n5*n2 T.F. de n3 (en puissances de 3)
380 IF (n3 /= 1) THEN
381 CALL jmccm1d3(p3,n3,m*(nprod/n3),table,ntable,itable,n,n/n3,work,nwork,ioff)
382 END IF
383
384 ! On transpose (on tient compte de ioff) en permutant les deux parties
385 ! On en profite pour multiplier par le bon wij
386 IF (n3 /= 1 .AND. nprod /= n3) THEN
387 CALL jmcctranspcs(m*n2,n,n3,nprod/(n2*n3), &
388 & table,ntable,itable,work,nwork,ioff)
389 END IF
390
391 ! On fait n2*n3 T.F. de n5 (en puissances de 5)
392 IF (n5 /= 1) THEN
393 CALL jmccm1d5(p5,n5,m*(nprod/n5),table,ntable,itable,n,n/n5,work,nwork,ioff)
394 END IF
395
396 ! On transpose s'il y a lieu (si on a fait quelque chose et s'il reste des
397 ! termes a traiter
398 IF (n5 /= 1 .AND. nprod /= n5 .AND. nterms > 7) THEN
399 CALL jmcctranspcs(m*n2*n3,n,n5,nprod/(n2*n3*n5), &
400 & table,ntable,itable,work,nwork,ioff)
401 END IF
402 nprod1 = m*n2*n3
403 nprod2 = n2*n3*n5
404 lastnp = n5
405
406 ! On passe aux nombres premiers autres que 2, 3 et 5
407 DO i = 7,nterms-1,2
408
409 pp = fact(ifact+i)
410 premier = fact(ifact+i+1)
411 np = premier**pp
412
413 CALL jmccm1dp(premier,pp,m*(nprod/np), &
414 & table,ntable,itable,n,n/np,work,nwork,ioff)
415
416 nprod1 = nprod1 * lastnp
417 nprod2 = nprod2 * np
418 IF (np /= 1 .AND. nprod /= np .AND. nterms > i+1) THEN
419 CALL jmcctranspcs(nprod1,n,np,nprod/nprod2, &
420 & table,ntable,itable,work,nwork,ioff)
421 END IF
422 lastnp = np
423
424 END DO
425
426END SUBROUTINE jmccm1d
427
428SUBROUTINE jmccm1d2(p,n,m,table,ntable,itable,ntable2,mtable,work,nwork,ioff)
429
430 IMPLICIT NONE
431
432 ! Arguments
433 INTEGER, INTENT(in) :: p, n, m
434 INTEGER, INTENT(in) :: ntable,itable,ntable2,mtable
435 REAL(kind=gft_prec), INTENT(in), DIMENSION(0:ntable-1) :: table
436 INTEGER, INTENT(in) :: nwork
437 REAL(kind=gft_prec), INTENT(inout), DIMENSION(0:nwork-1) :: work
438 INTEGER, INTENT(inout) :: ioff
439
440 ! Variables locales
441 INTEGER :: k, jl
442 INTEGER :: it1,iu1,it2,iu2
443 INTEGER :: jt1,ju1,jt2,ju2
444 REAL(kind=gft_prec) :: x1, x2, y1, y2
445 REAL(kind=gft_prec) :: c, s
446 INTEGER :: ioff1, ioff2
447
448 ! On joue sur ioff pour alterner entre le haut et le bas de work
449 ioff1 = ioff
450 ioff2 = nwork/2-ioff1
451
452 IF (mod(p,2)==0) THEN
453
454 ! Si p est pair, on peut travailler entierement en base 4
455 CALL jmccm1d4(p/2,n,m,table,ntable,itable,ntable2,mtable,work,nwork,ioff1)
456 ioff = ioff1
457
458 ELSE
459
460 ! On fait les premieres etapes en base 4
461 CALL jmccm1d4(p/2,n,2*m,table,ntable,itable,ntable2,mtable*2,work,nwork,ioff1)
462 ioff2 = nwork/2-ioff1
463 ! On fait la derniere etape en base 2
464 IF (m >= 16 .OR. 2**(p-1) < 8) THEN
465 DO k = 0,2**(p-1)-1
466
467 ! Les sinus et cosinus
468 c = table(itable+ mtable*k)
469 s = table(itable+ntable2+mtable*k)
470
471 ! Les indices
472 it1 = ioff1 +m*(k*2 )
473 iu1 = ioff1+nwork/4+m*(k*2 )
474 it2 = ioff1 +m*(k*2+1)
475 iu2 = ioff1+nwork/4+m*(k*2+1)
476 jt1 = ioff2 +m*( k )
477 ju1 = ioff2+nwork/4+m*( k )
478 jt2 = ioff2 +m*((k+2**(p-1)))
479 ju2 = ioff2+nwork/4+m*((k+2**(p-1)))
480
481!dir$ ivdep
482!ocl novrec
483!cdir nodep
484 DO jl = 0,m-1
485 x1 = work(it1+jl)
486 y1 = work(iu1+jl)
487 x2 = work(it2+jl)
488 y2 = work(iu2+jl)
489 work(jt1+jl) = x1 + ( x2*c - y2*s )
490 work(ju1+jl) = y1 + ( x2*s + y2*c )
491 work(jt2+jl) = x1 - ( x2*c - y2*s )
492 work(ju2+jl) = y1 - ( x2*s + y2*c )
493 END DO
494 END DO
495 ELSE
496 DO jl = 0,m-1
497!dir$ ivdep
498!ocl novrec
499!cdir nodep
500 DO k = 0,2**(p-1)-1
501 x1 = work(ioff1+jl +m*(k*2 ))
502 y1 = work(ioff1+jl+nwork/4+m*(k*2 ))
503 x2 = work(ioff1+jl +m*(k*2+1))
504 y2 = work(ioff1+jl+nwork/4+m*(k*2+1))
505 ! Les sinus et cosinus
506 c = table(itable+ mtable*k)
507 s = table(itable+ntable2+mtable*k)
508 work(ioff2+jl +m*( k )) = x1 + ( x2*c - y2*s )
509 work(ioff2+jl+nwork/4+m*( k )) = y1 + ( x2*s + y2*c )
510 work(ioff2+jl +m*((k+2**(p-1)))) = x1 - ( x2*c - y2*s )
511 work(ioff2+jl+nwork/4+m*((k+2**(p-1)))) = y1 - ( x2*s + y2*c )
512 END DO
513 END DO
514 END IF
515
516 ioff = ioff2
517
518 END IF
519
520END SUBROUTINE jmccm1d2
521
522SUBROUTINE jmccm1d3(p,n,m,table,ntable,itable,ntable2,mtable,work,nwork,ioff)
523
524 IMPLICIT NONE
525
526 ! Arguments
527 INTEGER, INTENT(in) :: p, n, m
528 INTEGER, INTENT(in) :: ntable,itable,ntable2, mtable
529 REAL(kind=gft_prec), INTENT(in), DIMENSION(0:ntable-1) :: table
530 INTEGER, INTENT(in) :: nwork
531 REAL(kind=gft_prec), INTENT(inout), DIMENSION(0:nwork-1) :: work
532 INTEGER, INTENT(inout) :: ioff
533
534 ! Variables locales
535 INTEGER :: i, k, jl
536 REAL(kind=gft_prec) :: x1, x2, x3, y1, y2, y3, t1, t2, t3, u1, u2, u3
537 REAL(kind=gft_prec) :: c2, s2, c3, s3
538 INTEGER :: it1,iu1,it2,iu2,it3,iu3
539 INTEGER :: jt1,ju1,jt2,ju2,jt3,ju3
540 REAL(kind=gft_prec) :: r,s,t,u
541
542 ! Gestion des constantes cosinus
543 REAL(kind=gft_prec), SAVE :: ctwopi3, stwopi3
544 LOGICAL, SAVE :: first = .true.
545
546 INTEGER :: ioff1, ioff2
547
548 ! On recupere cos et sin de 2*pi/3
549 IF (first) THEN
550 first = .false.
551 ctwopi3 = -1.0_gft_prec/2.0_gft_prec
552 stwopi3 = sqrt(3.0_gft_prec)/2.0_gft_prec
553 END IF
554
555 ! On joue sur ioff pour alterner entre le haut et le bas de work
556 ioff1 = ioff
557 ioff2 = nwork/2-ioff1
558
559 ! Boucle sur les etapes
560 DO i = 0, p-1
561
562 IF (m*3**(p-i-1) >= 16 .OR. 3**i < 8) THEN
563
564 DO k = 0,3**i-1
565
566 ! Les sinus et cosinus
567 c2 = table(itable+ mtable* 3**(p-i-1)*k)
568 s2 = table(itable+ntable2+mtable* 3**(p-i-1)*k)
569 c3 = table(itable+ mtable*2*3**(p-i-1)*k)
570 s3 = table(itable+ntable2+mtable*2*3**(p-i-1)*k)
571
572 ! Les indices
573 it1 = ioff1 +m*(k*3**(p-i) )
574 iu1 = ioff1+nwork/4+m*(k*3**(p-i) )
575 it2 = ioff1 +m*(k*3**(p-i)+ 3**(p-i-1))
576 iu2 = ioff1+nwork/4+m*(k*3**(p-i)+ 3**(p-i-1))
577 it3 = ioff1 +m*(k*3**(p-i)+2*3**(p-i-1))
578 iu3 = ioff1+nwork/4+m*(k*3**(p-i)+2*3**(p-i-1))
579 jt1 = ioff2 +m*( k *3**(p-i-1))
580 ju1 = ioff2+nwork/4+m*( k *3**(p-i-1))
581 jt2 = ioff2 +m*((k+ 3**i)*3**(p-i-1))
582 ju2 = ioff2+nwork/4+m*((k+ 3**i)*3**(p-i-1))
583 jt3 = ioff2 +m*((k+2*3**i)*3**(p-i-1))
584 ju3 = ioff2+nwork/4+m*((k+2*3**i)*3**(p-i-1))
585
586!dir$ ivdep
587!ocl novrec
588!cdir nodep
589 DO jl = 0,m*3**(p-i-1)-1
590
591 r = (c2*work(it2+jl))-(s2*work(iu2+jl))
592 s = (c2*work(iu2+jl))+(s2*work(it2+jl))
593 t = (c3*work(it3+jl))-(s3*work(iu3+jl))
594 u = (c3*work(iu3+jl))+(s3*work(it3+jl))
595 x1 = work(it1+jl)
596 y1 = work(iu1+jl)
597 work(jt1+jl) = x1 + r + t
598 work(ju1+jl) = y1 + s + u
599 work(jt2+jl) = x1 + ctwopi3*(r+t) - stwopi3*(s-u)
600 work(ju2+jl) = y1 + ctwopi3*(s+u) + stwopi3*(r-t)
601 work(jt3+jl) = x1 + ctwopi3*(r+t) + stwopi3*(s-u)
602 work(ju3+jl) = y1 + ctwopi3*(s+u) - stwopi3*(r-t)
603
604 END DO
605
606 END DO
607
608 ELSE
609
610 DO jl = 0,m*3**(p-i-1)-1
611
612!dir$ ivdep
613!ocl novrec
614!cdir nodep
615 DO k = 0,3**i-1
616
617 t1 = work(ioff1+jl +m*(k*3**(p-i) ))
618 u1 = work(ioff1+jl+nwork/4+m*(k*3**(p-i) ))
619 t2 = work(ioff1+jl +m*(k*3**(p-i)+ 3**(p-i-1)))
620 u2 = work(ioff1+jl+nwork/4+m*(k*3**(p-i)+ 3**(p-i-1)))
621 t3 = work(ioff1+jl +m*(k*3**(p-i)+2*3**(p-i-1)))
622 u3 = work(ioff1+jl+nwork/4+m*(k*3**(p-i)+2*3**(p-i-1)))
623
624 ! Les sinus et cosinus
625 c2 = table(itable+ mtable* 3**(p-i-1)*k)
626 s2 = table(itable+ntable2+mtable* 3**(p-i-1)*k)
627 c3 = table(itable+ mtable*2*3**(p-i-1)*k)
628 s3 = table(itable+ntable2+mtable*2*3**(p-i-1)*k)
629
630 ! On premultiplie
631 x1 = t1
632 y1 = u1
633 x2 = c2*t2-s2*u2
634 y2 = c2*u2+s2*t2
635 x3 = c3*t3-s3*u3
636 y3 = c3*u3+s3*t3
637
638 ! Il reste a multiplier par les twopi3
639 work(ioff2+jl +m*( k *3**(p-i-1))) = &
640 & x1 + x2 + x3
641 work(ioff2+jl+nwork/4+m*( k *3**(p-i-1))) = &
642 & y1 + y2 + y3
643 work(ioff2+jl +m*((k+ 3**i)*3**(p-i-1))) = &
644 & x1 + ctwopi3*x2-stwopi3*y2 + ctwopi3*x3+stwopi3*y3
645 work(ioff2+jl+nwork/4+m*((k+ 3**i)*3**(p-i-1))) = &
646 & y1 + ctwopi3*y2+stwopi3*x2 + ctwopi3*y3-stwopi3*x3
647 work(ioff2+jl +m*((k+2*3**i)*3**(p-i-1))) = &
648 & x1 + ctwopi3*x2+stwopi3*y2 + ctwopi3*x3-stwopi3*y3
649 work(ioff2+jl+nwork/4+m*((k+2*3**i)*3**(p-i-1))) = &
650 & y1 + ctwopi3*y2-stwopi3*x2 + ctwopi3*y3+stwopi3*x3
651
652 END DO
653
654 END DO
655
656 END IF
657
658 ! On alterne les offsets
659 ioff1 = nwork/2-ioff1
660 ioff2 = nwork/2-ioff2
661
662 ! Fin boucle sur le nombre d'etapes
663 END DO
664
665 ioff = ioff1
666
667END SUBROUTINE jmccm1d3
668
669SUBROUTINE jmccm1d4(p,n,m,table,ntable,itable,ntable2,mtable,work,nwork,ioff)
670
671 IMPLICIT NONE
672
673 ! Arguments
674 INTEGER, INTENT(in) :: p, n, m
675 INTEGER, INTENT(in) :: ntable,itable,ntable2, mtable
676 REAL(kind=gft_prec), INTENT(in), DIMENSION(0:ntable-1) :: table
677 INTEGER, INTENT(in) :: nwork
678 REAL(kind=gft_prec), INTENT(inout), DIMENSION(0:nwork-1) :: work
679 INTEGER, INTENT(inout) :: ioff
680
681 ! Variables locales
682 INTEGER :: i, k, jl
683 REAL(kind=gft_prec) :: x0,x1,x2,x3,y0,y1,y2,y3,t0,t1,t2,t3,u0,u1,u2,u3
684 REAL(kind=gft_prec) :: x0px2,x0mx2,x1px3,x1mx3
685 REAL(kind=gft_prec) :: y0py2,y0my2,y1py3,y1my3
686 REAL(kind=gft_prec) :: c1, s1, c2, s2, c3, s3
687 INTEGER :: ioff1, ioff2
688 INTEGER :: it0,iu0,it1,iu1,it2,iu2,it3,iu3
689 INTEGER :: jt0,ju0,jt1,ju1,jt2,ju2,jt3,ju3
690
691 ! On joue sur ioff pour alterner entre le haut et le bas de work
692 ioff1 = ioff
693 ioff2 = nwork/2-ioff1
694
695 ! Boucle sur les etapes
696 DO i = 0, p-1
697
698 IF (m*4**(p-i-1) >= 16 .OR. 4**i < 8) THEN
699
700 DO k = 0,4**i-1
701
702 ! Les sinus et cosinus
703 c1 = table(itable+ mtable* 4**(p-i-1)*k)
704 s1 = table(itable+ntable2+mtable* 4**(p-i-1)*k)
705 c2 = table(itable+ mtable*2*4**(p-i-1)*k)
706 s2 = table(itable+ntable2+mtable*2*4**(p-i-1)*k)
707 c3 = table(itable+ mtable*3*4**(p-i-1)*k)
708 s3 = table(itable+ntable2+mtable*3*4**(p-i-1)*k)
709
710 ! Les indices
711 it0 = ioff1 +m*(k*4**(p-i) )
712 iu0 = ioff1+nwork/4+m*(k*4**(p-i) )
713 it1 = ioff1 +m*(k*4**(p-i)+ 4**(p-i-1))
714 iu1 = ioff1+nwork/4+m*(k*4**(p-i)+ 4**(p-i-1))
715 it2 = ioff1 +m*(k*4**(p-i)+2*4**(p-i-1))
716 iu2 = ioff1+nwork/4+m*(k*4**(p-i)+2*4**(p-i-1))
717 it3 = ioff1 +m*(k*4**(p-i)+3*4**(p-i-1))
718 iu3 = ioff1+nwork/4+m*(k*4**(p-i)+3*4**(p-i-1))
719 jt0 = ioff2 +m*( k *4**(p-i-1))
720 ju0 = ioff2+nwork/4+m*( k *4**(p-i-1))
721 jt1 = ioff2 +m*((k+ 4**i)*4**(p-i-1))
722 ju1 = ioff2+nwork/4+m*((k+ 4**i)*4**(p-i-1))
723 jt2 = ioff2 +m*((k+2*4**i)*4**(p-i-1))
724 ju2 = ioff2+nwork/4+m*((k+2*4**i)*4**(p-i-1))
725 jt3 = ioff2 +m*((k+3*4**i)*4**(p-i-1))
726 ju3 = ioff2+nwork/4+m*((k+3*4**i)*4**(p-i-1))
727
728!dir$ ivdep
729!ocl novrec
730!cdir nodep
731 DO jl = 0,m*4**(p-i-1)-1
732
733 x0px2 = work(it0+jl) + (c2*work(it2+jl)-s2*work(iu2+jl))
734 x0mx2 = work(it0+jl) - (c2*work(it2+jl)-s2*work(iu2+jl))
735 y0py2 = work(iu0+jl) + (c2*work(iu2+jl)+s2*work(it2+jl))
736 y0my2 = work(iu0+jl) - (c2*work(iu2+jl)+s2*work(it2+jl))
737 x1px3 = (c1*work(it1+jl)-s1*work(iu1+jl))+(c3*work(it3+jl)-s3*work(iu3+jl))
738 x1mx3 = (c1*work(it1+jl)-s1*work(iu1+jl))-(c3*work(it3+jl)-s3*work(iu3+jl))
739 y1py3 = (c1*work(iu1+jl)+s1*work(it1+jl))+(c3*work(iu3+jl)+s3*work(it3+jl))
740 y1my3 = (c1*work(iu1+jl)+s1*work(it1+jl))-(c3*work(iu3+jl)+s3*work(it3+jl))
741
742 ! Il reste a multiplier par les twopi4
743 work(jt0+jl) = (x0px2)+(x1px3)
744 work(jt2+jl) = (x0px2)-(x1px3)
745 work(ju0+jl) = (y0py2)+(y1py3)
746 work(ju2+jl) = (y0py2)-(y1py3)
747 work(jt1+jl) = (x0mx2)-(y1my3)
748 work(jt3+jl) = (x0mx2)+(y1my3)
749 work(ju1+jl) = (y0my2)+(x1mx3)
750 work(ju3+jl) = (y0my2)-(x1mx3)
751
752 END DO
753
754 END DO
755
756 ELSE
757
758 DO jl = 0,m*4**(p-i-1)-1
759
760!dir$ ivdep
761!ocl novrec
762!cdir nodep
763 DO k = 0,4**i-1
764
765 t0 = work(ioff1+jl +m*(k*4**(p-i) ))
766 u0 = work(ioff1+jl+nwork/4+m*(k*4**(p-i) ))
767 t1 = work(ioff1+jl +m*(k*4**(p-i)+ 4**(p-i-1)))
768 u1 = work(ioff1+jl+nwork/4+m*(k*4**(p-i)+ 4**(p-i-1)))
769 t2 = work(ioff1+jl +m*(k*4**(p-i)+2*4**(p-i-1)))
770 u2 = work(ioff1+jl+nwork/4+m*(k*4**(p-i)+2*4**(p-i-1)))
771 t3 = work(ioff1+jl +m*(k*4**(p-i)+3*4**(p-i-1)))
772 u3 = work(ioff1+jl+nwork/4+m*(k*4**(p-i)+3*4**(p-i-1)))
773
774 ! Les sinus et cosinus
775 c1 = table(itable+ mtable* 4**(p-i-1)*k)
776 s1 = table(itable+ntable2+mtable* 4**(p-i-1)*k)
777 c2 = table(itable+ mtable*2*4**(p-i-1)*k)
778 s2 = table(itable+ntable2+mtable*2*4**(p-i-1)*k)
779 c3 = table(itable+ mtable*3*4**(p-i-1)*k)
780 s3 = table(itable+ntable2+mtable*3*4**(p-i-1)*k)
781
782 ! On premultiplie
783 x0 = t0
784 y0 = u0
785 x1 = c1*t1-s1*u1
786 y1 = c1*u1+s1*t1
787 x2 = c2*t2-s2*u2
788 y2 = c2*u2+s2*t2
789 x3 = c3*t3-s3*u3
790 y3 = c3*u3+s3*t3
791
792 ! Il reste a multiplier par les twopi4
793 work(ioff2+jl +m*( k *4**(p-i-1))) = x0+x1+x2+x3
794 work(ioff2+jl+nwork/4+m*( k *4**(p-i-1))) = y0+y1+y2+y3
795 work(ioff2+jl +m*((k+ 4**i)*4**(p-i-1))) = x0-y1-x2+y3
796 work(ioff2+jl+nwork/4+m*((k+ 4**i)*4**(p-i-1))) = y0+x1-y2-x3
797 work(ioff2+jl +m*((k+2*4**i)*4**(p-i-1))) = x0-x1+x2-x3
798 work(ioff2+jl+nwork/4+m*((k+2*4**i)*4**(p-i-1))) = y0-y1+y2-y3
799 work(ioff2+jl +m*((k+3*4**i)*4**(p-i-1))) = x0+y1-x2-y3
800 work(ioff2+jl+nwork/4+m*((k+3*4**i)*4**(p-i-1))) = y0-x1-y2+x3
801
802 END DO
803
804 END DO
805
806 END IF
807
808 ! On alterne les offsets
809 ioff1 = nwork/2-ioff1
810 ioff2 = nwork/2-ioff2
811
812 ! Fin boucle sur le nombre d'etapes
813 END DO
814
815 ioff = ioff1
816
817END SUBROUTINE jmccm1d4
818
819SUBROUTINE jmccm1d5(p,n,m,table,ntable,itable,ntable2,mtable,work,nwork,ioff)
820
821 IMPLICIT NONE
822
823 ! Arguments
824 INTEGER, INTENT(in) :: p, n, m
825 INTEGER, INTENT(in) :: ntable,itable,ntable2, mtable
826 REAL(kind=gft_prec), INTENT(in), DIMENSION(0:ntable-1) :: table
827 INTEGER, INTENT(in) :: nwork
828 REAL(kind=gft_prec), INTENT(inout), DIMENSION(0:nwork-1) :: work
829 INTEGER, INTENT(inout) :: ioff
830
831 ! Variables locales
832 INTEGER :: i, k, jl
833 REAL(kind=gft_prec) :: x0,x1,x2,x3,x4,y0,y1,y2,y3,y4,t0,t1,t2,t3,t4,u0,u1,u2,u3,u4
834 REAL(kind=gft_prec) :: c1, s1, c2, s2, c3, s3, c4, s4
835 INTEGER :: ioff1, ioff2
836
837 ! Gestion des constantes cosinus
838 REAL(kind=gft_prec), SAVE :: twopi5
839 REAL(kind=gft_prec), SAVE :: ctwopi51, ctwopi52, ctwopi53, ctwopi54
840 REAL(kind=gft_prec), SAVE :: stwopi51, stwopi52, stwopi53, stwopi54
841 LOGICAL, SAVE :: first = .true.
842
843 ! On recupere cos et sin de 2*pi/5
844 IF (first) THEN
845 first = .false.
846 twopi5 = 2.0_gft_prec*acos(-1.0_gft_prec)/5.0_gft_prec
847 ctwopi51 = cos( twopi5)
848 stwopi51 = sin( twopi5)
849 ctwopi52 = cos(2.0_gft_prec*twopi5)
850 stwopi52 = sin(2.0_gft_prec*twopi5)
851 ctwopi53 = cos(3.0_gft_prec*twopi5)
852 stwopi53 = sin(3.0_gft_prec*twopi5)
853 ctwopi54 = cos(4.0_gft_prec*twopi5)
854 stwopi54 = sin(4.0_gft_prec*twopi5)
855 END IF
856
857 ! On joue sur ioff pour alterner entre le haut et le bas de work
858 ioff1 = ioff
859 ioff2 = nwork/2-ioff1
860
861 ! Boucle sur les etapes
862 DO i = 0, p-1
863
864 IF (m*5**(p-i-1) >= 16 .OR. 5**i < 8) THEN
865
866 DO k = 0,5**i-1
867
868!dir$ ivdep
869!ocl novrec
870!cdir nodep
871 DO jl = 0,m*5**(p-i-1)-1
872
873 t0 = work(ioff1+jl +m*(k*5**(p-i) ))
874 u0 = work(ioff1+jl+nwork/4+m*(k*5**(p-i) ))
875 t1 = work(ioff1+jl +m*(k*5**(p-i)+ 5**(p-i-1)))
876 u1 = work(ioff1+jl+nwork/4+m*(k*5**(p-i)+ 5**(p-i-1)))
877 t2 = work(ioff1+jl +m*(k*5**(p-i)+2*5**(p-i-1)))
878 u2 = work(ioff1+jl+nwork/4+m*(k*5**(p-i)+2*5**(p-i-1)))
879 t3 = work(ioff1+jl +m*(k*5**(p-i)+3*5**(p-i-1)))
880 u3 = work(ioff1+jl+nwork/4+m*(k*5**(p-i)+3*5**(p-i-1)))
881 t4 = work(ioff1+jl +m*(k*5**(p-i)+4*5**(p-i-1)))
882 u4 = work(ioff1+jl+nwork/4+m*(k*5**(p-i)+4*5**(p-i-1)))
883
884 ! Les sinus et cosinus
885 c1 = table(itable+ mtable* 5**(p-i-1)*k)
886 s1 = table(itable+ntable2+mtable* 5**(p-i-1)*k)
887 c2 = table(itable+ mtable*2*5**(p-i-1)*k)
888 s2 = table(itable+ntable2+mtable*2*5**(p-i-1)*k)
889 c3 = table(itable+ mtable*3*5**(p-i-1)*k)
890 s3 = table(itable+ntable2+mtable*3*5**(p-i-1)*k)
891 c4 = table(itable+ mtable*4*5**(p-i-1)*k)
892 s4 = table(itable+ntable2+mtable*4*5**(p-i-1)*k)
893
894 ! On premultiplie
895 x0 = t0
896 y0 = u0
897 x1 = c1*t1-s1*u1
898 y1 = c1*u1+s1*t1
899 x2 = c2*t2-s2*u2
900 y2 = c2*u2+s2*t2
901 x3 = c3*t3-s3*u3
902 y3 = c3*u3+s3*t3
903 x4 = c4*t4-s4*u4
904 y4 = c4*u4+s4*t4
905
906 ! Il reste a multiplier par les twopi5
907 work(ioff2+jl +m*( k *5**(p-i-1))) = &
908 & x0 + x1 + x2 &
909 & + x3 + x4
910 work(ioff2+jl+nwork/4+m*( k *5**(p-i-1))) = &
911 & y0 + y1 + y2 &
912 & + y3 + y4
913 work(ioff2+jl +m*((k+ 5**i)*5**(p-i-1))) = &
914 & x0 + ctwopi51*x1 - stwopi51*y1 &
915 & + ctwopi52*x2 - stwopi52*y2 &
916 & + ctwopi53*x3 - stwopi53*y3 &
917 & + ctwopi54*x4 - stwopi54*y4
918 work(ioff2+jl+nwork/4+m*((k+ 5**i)*5**(p-i-1))) = &
919 & y0 + ctwopi51*y1 + stwopi51*x1 &
920 & + ctwopi52*y2 + stwopi52*x2 &
921 & + ctwopi53*y3 + stwopi53*x3 &
922 & + ctwopi54*y4 + stwopi54*x4
923 work(ioff2+jl +m*((k+2*5**i)*5**(p-i-1))) = &
924 & x0 + ctwopi52*x1 - stwopi52*y1 &
925 & + ctwopi54*x2 - stwopi54*y2 &
926 & + ctwopi51*x3 - stwopi51*y3 &
927 & + ctwopi53*x4 - stwopi53*y4
928 work(ioff2+jl+nwork/4+m*((k+2*5**i)*5**(p-i-1))) = &
929 & y0 + ctwopi52*y1 + stwopi52*x1 &
930 & + ctwopi54*y2 + stwopi54*x2 &
931 & + ctwopi51*y3 + stwopi51*x3 &
932 & + ctwopi53*y4 + stwopi53*x4
933 work(ioff2+jl +m*((k+3*5**i)*5**(p-i-1))) = &
934 & x0 + ctwopi53*x1 - stwopi53*y1 &
935 & + ctwopi51*x2 - stwopi51*y2 &
936 & + ctwopi54*x3 - stwopi54*y3 &
937 & + ctwopi52*x4 - stwopi52*y4
938 work(ioff2+jl+nwork/4+m*((k+3*5**i)*5**(p-i-1))) = &
939 & y0 + ctwopi53*y1 + stwopi53*x1 &
940 & + ctwopi51*y2 + stwopi51*x2 &
941 & + ctwopi54*y3 + stwopi54*x3 &
942 & + ctwopi52*y4 + stwopi52*x4
943 work(ioff2+jl +m*((k+4*5**i)*5**(p-i-1))) = &
944 & x0 + ctwopi54*x1 - stwopi54*y1 &
945 & + ctwopi53*x2 - stwopi53*y2 &
946 & + ctwopi52*x3 - stwopi52*y3 &
947 & + ctwopi51*x4 - stwopi51*y4
948 work(ioff2+jl+nwork/4+m*((k+4*5**i)*5**(p-i-1))) = &
949 & y0 + ctwopi54*y1 + stwopi54*x1 &
950 & + ctwopi53*y2 + stwopi53*x2 &
951 & + ctwopi52*y3 + stwopi52*x3 &
952 & + ctwopi51*y4 + stwopi51*x4
953
954 END DO
955
956 END DO
957
958 ELSE
959
960 DO jl = 0,m*5**(p-i-1)-1
961
962!dir$ ivdep
963!ocl novrec
964!cdir nodep
965 DO k = 0,5**i-1
966
967 t0 = work(ioff1+jl +m*(k*5**(p-i) ))
968 u0 = work(ioff1+jl+nwork/4+m*(k*5**(p-i) ))
969 t1 = work(ioff1+jl +m*(k*5**(p-i)+ 5**(p-i-1)))
970 u1 = work(ioff1+jl+nwork/4+m*(k*5**(p-i)+ 5**(p-i-1)))
971 t2 = work(ioff1+jl +m*(k*5**(p-i)+2*5**(p-i-1)))
972 u2 = work(ioff1+jl+nwork/4+m*(k*5**(p-i)+2*5**(p-i-1)))
973 t3 = work(ioff1+jl +m*(k*5**(p-i)+3*5**(p-i-1)))
974 u3 = work(ioff1+jl+nwork/4+m*(k*5**(p-i)+3*5**(p-i-1)))
975 t4 = work(ioff1+jl +m*(k*5**(p-i)+4*5**(p-i-1)))
976 u4 = work(ioff1+jl+nwork/4+m*(k*5**(p-i)+4*5**(p-i-1)))
977
978 ! Les sinus et cosinus
979 c1 = table(itable+ mtable* 5**(p-i-1)*k)
980 s1 = table(itable+ntable2+mtable* 5**(p-i-1)*k)
981 c2 = table(itable+ mtable*2*5**(p-i-1)*k)
982 s2 = table(itable+ntable2+mtable*2*5**(p-i-1)*k)
983 c3 = table(itable+ mtable*3*5**(p-i-1)*k)
984 s3 = table(itable+ntable2+mtable*3*5**(p-i-1)*k)
985 c4 = table(itable+ mtable*4*5**(p-i-1)*k)
986 s4 = table(itable+ntable2+mtable*4*5**(p-i-1)*k)
987
988 ! On premultiplie
989 x0 = t0
990 y0 = u0
991 x1 = c1*t1-s1*u1
992 y1 = c1*u1+s1*t1
993 x2 = c2*t2-s2*u2
994 y2 = c2*u2+s2*t2
995 x3 = c3*t3-s3*u3
996 y3 = c3*u3+s3*t3
997 x4 = c4*t4-s4*u4
998 y4 = c4*u4+s4*t4
999
1000 ! Il reste a multiplier par les twopi5
1001 work(ioff2+jl +m*( k *5**(p-i-1))) = &
1002 & x0 + x1 + x2 &
1003 & + x3 + x4
1004 work(ioff2+jl+nwork/4+m*( k *5**(p-i-1))) = &
1005 & y0 + y1 + y2 &
1006 & + y3 + y4
1007 work(ioff2+jl +m*((k+ 5**i)*5**(p-i-1))) = &
1008 & x0 + ctwopi51*x1 - stwopi51*y1 &
1009 & + ctwopi52*x2 - stwopi52*y2 &
1010 & + ctwopi53*x3 - stwopi53*y3 &
1011 & + ctwopi54*x4 - stwopi54*y4
1012 work(ioff2+jl+nwork/4+m*((k+ 5**i)*5**(p-i-1))) = &
1013 & y0 + ctwopi51*y1 + stwopi51*x1 &
1014 & + ctwopi52*y2 + stwopi52*x2 &
1015 & + ctwopi53*y3 + stwopi53*x3 &
1016 & + ctwopi54*y4 + stwopi54*x4
1017 work(ioff2+jl +m*((k+2*5**i)*5**(p-i-1))) = &
1018 & x0 + ctwopi52*x1 - stwopi52*y1 &
1019 & + ctwopi54*x2 - stwopi54*y2 &
1020 & + ctwopi51*x3 - stwopi51*y3 &
1021 & + ctwopi53*x4 - stwopi53*y4
1022 work(ioff2+jl+nwork/4+m*((k+2*5**i)*5**(p-i-1))) = &
1023 & y0 + ctwopi52*y1 + stwopi52*x1 &
1024 & + ctwopi54*y2 + stwopi54*x2 &
1025 & + ctwopi51*y3 + stwopi51*x3 &
1026 & + ctwopi53*y4 + stwopi53*x4
1027 work(ioff2+jl +m*((k+3*5**i)*5**(p-i-1))) = &
1028 & x0 + ctwopi53*x1 - stwopi53*y1 &
1029 & + ctwopi51*x2 - stwopi51*y2 &
1030 & + ctwopi54*x3 - stwopi54*y3 &
1031 & + ctwopi52*x4 - stwopi52*y4
1032 work(ioff2+jl+nwork/4+m*((k+3*5**i)*5**(p-i-1))) = &
1033 & y0 + ctwopi53*y1 + stwopi53*x1 &
1034 & + ctwopi51*y2 + stwopi51*x2 &
1035 & + ctwopi54*y3 + stwopi54*x3 &
1036 & + ctwopi52*y4 + stwopi52*x4
1037 work(ioff2+jl +m*((k+4*5**i)*5**(p-i-1))) = &
1038 & x0 + ctwopi54*x1 - stwopi54*y1 &
1039 & + ctwopi53*x2 - stwopi53*y2 &
1040 & + ctwopi52*x3 - stwopi52*y3 &
1041 & + ctwopi51*x4 - stwopi51*y4
1042 work(ioff2+jl+nwork/4+m*((k+4*5**i)*5**(p-i-1))) = &
1043 & y0 + ctwopi54*y1 + stwopi54*x1 &
1044 & + ctwopi53*y2 + stwopi53*x2 &
1045 & + ctwopi52*y3 + stwopi52*x3 &
1046 & + ctwopi51*y4 + stwopi51*x4
1047
1048 END DO
1049
1050 END DO
1051
1052 END IF
1053
1054 ! On alterne les offsets
1055 ioff1 = nwork/2-ioff1
1056 ioff2 = nwork/2-ioff2
1057
1058 ! Fin boucle sur le nombre d'etapes
1059 END DO
1060
1061 ioff = ioff1
1062
1063END SUBROUTINE jmccm1d5
1064
1065SUBROUTINE jmccm1dp(p,q,m,table,ntable,itable,ntable2,mtable,work,nwork,ioff)
1066
1067 ! On fait m t.f. d'ordre q en base p (p**q)
1068 ! Note : n n'est pas utilise ici
1069
1070 IMPLICIT NONE
1071
1072 ! Arguments
1073 INTEGER, INTENT(in) :: p, q, m
1074 INTEGER, INTENT(in) :: ntable,itable,ntable2, mtable
1075 REAL(kind=gft_prec), INTENT(in), DIMENSION(0:ntable-1) :: table
1076 INTEGER, INTENT(in) :: nwork
1077 REAL(kind=gft_prec), INTENT(inout), DIMENSION(0:nwork-1) :: work
1078 INTEGER, INTENT(inout) :: ioff
1079
1080 ! Variables locales
1081 INTEGER :: i, k, jl, jp, kp
1082 REAL(kind=gft_prec) :: ck, sk, tk, uk, cpjk, spjk
1083 INTEGER :: pqq, pi, pqi, pqii
1084 INTEGER :: ikpr, ikpi, ijpr, ijpi
1085 INTEGER :: itr, iti, jtr, jti
1086 INTEGER :: ioff1, ioff2
1087 REAL(kind=gft_prec) :: c11, c12, c21, c22
1088
1089 ! On joue sur ioff pour alterner entre le haut et le bas de work
1090 ioff1 = ioff
1091 ioff2 = nwork/2-ioff1
1092
1093 ! Pour le calcul des cos(2*pi/p)
1094 pqq = p**(q-1)
1095
1096 ! Boucle sur les etapes
1097 DO i = 0, q-1
1098
1099 pi = p**i
1100 pqi = p**(q-i)
1101 pqii = p**(q-i-1)
1102
1103 DO k = 0,pi-1
1104
1105 DO jp = 0,p-1
1106
1107 DO jl = 0,m*pqii-1
1108
1109 ijpr = ioff2 + jl + m*((k+jp*pi)*pqii)
1110 ijpi = ijpr + nwork/4
1111
1112 work(ijpr) = 0
1113 work(ijpi) = 0
1114
1115 END DO
1116
1117 END DO
1118
1119 DO kp = 0,p-1
1120
1121 itr = itable+mtable*kp*pqii*k
1122 iti = itr + ntable2
1123 ck = table(itr)
1124 sk = table(iti)
1125
1126 DO jp = 0,p-1
1127
1128 ! Gymanstique infernale pour recuperer cos(2*pi/p) etc
1129 jtr = itable+mtable*pqq*mod(jp*kp,p)
1130 jti = jtr + ntable2
1131 cpjk = table(jtr)
1132 spjk = table(jti)
1133 c11 = (cpjk*ck-spjk*sk)
1134 c12 = (cpjk*sk+spjk*ck)
1135 c21 = (cpjk*sk+spjk*ck)
1136 c22 = (cpjk*ck-spjk*sk)
1137
1138!dir$ ivdep
1139!ocl novrec
1140!cdir nodep
1141 DO jl = 0,m*pqii-1
1142
1143 ikpr = ioff1+jl+m*(k*pqi+kp*pqii)
1144 ikpi = ikpr + nwork/4
1145 tk = work(ikpr)
1146 uk = work(ikpi)
1147
1148 ijpr = ioff2+jl+m*((k+jp*pi)*pqii)
1149 ijpi = ijpr + nwork/4
1150
1151 work(ijpr) = work(ijpr) + tk*c11-uk*c12
1152 work(ijpi) = work(ijpi) + tk*c21+uk*c22
1153
1154 END DO
1155
1156 END DO
1157
1158 END DO
1159
1160 END DO
1161
1162 ! On alterne les offsets
1163 ioff1 = nwork/2 - ioff1
1164 ioff2 = nwork/2 - ioff2
1165
1166 ! Fin boucle sur le nombre d'etapes
1167 END DO
1168
1169 ioff = ioff1
1170
1171END SUBROUTINE jmccm1dp
1172
1173SUBROUTINE jmcsm1d(m,n,fact,nfact,ifact,table,ntable,itable,work,nwork,ioff)
1174
1175 IMPLICIT NONE
1176
1177 ! Arguments
1178 INTEGER, INTENT(in) :: m, n
1179 INTEGER, INTENT(in) :: nfact, ifact
1180 INTEGER, INTENT(inout), DIMENSION(0:nfact-1) :: fact
1181 INTEGER, INTENT(in) :: ntable,itable
1182 REAL(kind=gft_prec), INTENT(inout), DIMENSION(0:ntable-1) :: table
1183 INTEGER, INTENT(in) :: nwork
1184 REAL(kind=gft_prec), INTENT(inout), DIMENSION(0:nwork-1) :: work
1185 INTEGER, INTENT(inout) :: ioff
1186
1187 ! Variables locales
1188 INTEGER :: ioff1, ioff2
1189 INTEGER :: i, j
1190 REAL(kind=gft_prec) :: t, u, v, w, tt, uu, vv, ww
1191 REAL(kind=gft_prec) :: c, s
1192 INTEGER :: it
1193
1194 ! Gestion de work
1195 ioff1 = ioff
1196 ioff2 = nwork/2 - ioff1
1197
1198 ! On doit faire m T.F. complexes -> reelles de longueur n
1199 ! Si m est pair
1200 IF (mod(m,2) == 0) THEN
1201
1202 ! On distribue
1203
1204 IF (m/2 >= 16 .OR. n/2 < 8) THEN
1205
1206 DO i = 0,n/2
1207!dir$ ivdep
1208!ocl novrec
1209!cdir nodep
1210 DO j = 0,m/2-1
1211 it = n-i
1212 IF (i == 0) it = 0
1213 t = work(ioff1 +i*m+j )
1214 u = work(ioff1+nwork/4+i*m+j )
1215 v = work(ioff1 +i*m+j+m/2)
1216 w = work(ioff1+nwork/4+i*m+j+m/2)
1217 work(ioff2 + i*m/2+j) = (t-w)
1218 work(ioff2+nwork/4+ i*m/2+j) = (u+v)
1219 work(ioff2 +it*m/2+j) = (t+w)
1220 work(ioff2+nwork/4+it*m/2+j) = (v-u)
1221 END DO
1222 END DO
1223
1224 ELSE
1225
1226 DO j = 0,m/2-1
1227!dir$ ivdep
1228!ocl novrec
1229!cdir nodep
1230 DO i = 0,n/2
1231 it = n-i
1232 IF (i == 0) it = 0
1233 t = work(ioff1 +i*m+j )
1234 u = work(ioff1+nwork/4+i*m+j )
1235 v = work(ioff1 +i*m+j+m/2)
1236 w = work(ioff1+nwork/4+i*m+j+m/2)
1237 work(ioff2 + i*m/2+j) = (t-w)
1238 work(ioff2+nwork/4+ i*m/2+j) = (u+v)
1239 work(ioff2 +it*m/2+j) = (t+w)
1240 work(ioff2+nwork/4+it*m/2+j) = (v-u)
1241 END DO
1242 END DO
1243
1244 END IF
1245
1246 ! On fait m/2 t.f. complexes -> complexes de longueur n
1247 CALL jmccm1d(m/2,n,fact,nfact,ifact,table,ntable,itable,work,nwork,ioff2)
1248 ioff1 = nwork/2 - ioff2
1249
1250 ! On reconstitue
1251
1252 IF (m/2 >= 16 .OR. n < 8) THEN
1253
1254 DO i = 0,n-1
1255!dir$ ivdep
1256!ocl novrec
1257!cdir nodep
1258 DO j = 0,m/2-1
1259 work(ioff1+i*m+j ) = work(ioff2 +i*m/2+j)
1260 work(ioff1+i*m+j+m/2) = work(ioff2+nwork/4+i*m/2+j)
1261 END DO
1262 END DO
1263
1264 ELSE
1265
1266 DO j = 0,m/2-1
1267!dir$ ivdep
1268!ocl novrec
1269!cdir nodep
1270 DO i = 0,n-1
1271 work(ioff1+i*m+j ) = work(ioff2 +i*m/2+j)
1272 work(ioff1+i*m+j+m/2) = work(ioff2+nwork/4+i*m/2+j)
1273 END DO
1274 END DO
1275
1276 END IF
1277
1278 ! Si m n'est pas pair mais que n l'est
1279 ELSE IF (mod(n,2) == 0) THEN
1280
1281 ! On distribue
1282
1283 IF (m >= 16 .OR. n/2 < 8) THEN
1284
1285 DO i = 0,n/2-1
1286!dir$ ivdep
1287!ocl novrec
1288!cdir nodep
1289 DO j = 0,m-1
1290 ! Note : Signe - sur les parties imaginaires pour inversion
1291 t = work(ioff1 + i*m+j)
1292 u = -work(ioff1+nwork/4+ i*m+j)
1293 v = work(ioff1 +(n/2-i)*m+j)
1294 w = -work(ioff1+nwork/4+(n/2-i)*m+j)
1295 c = table(itable+i)
1296 s = table(itable+i+n)
1297 tt = (t+v)/2.0_gft_prec
1298 uu = (u-w)/2.0_gft_prec
1299 vv = (c*(t-v)+s*(u+w))/2.0_gft_prec
1300 ww = (c*(u+w)-s*(t-v))/2.0_gft_prec
1301 ! Note : le facteur 2 et le signe - viennent de l'inversion Fourier
1302 work(ioff2 +m*i+j) = 2.0_gft_prec*(tt-ww)
1303 work(ioff2+nwork/4+m*i+j) = -2.0_gft_prec*(uu+vv)
1304 END DO
1305 END DO
1306
1307 ELSE
1308
1309 DO j = 0,m-1
1310!dir$ ivdep
1311!ocl novrec
1312!cdir nodep
1313 DO i = 0,n/2-1
1314 ! Note : Signe - sur les parties imaginaires pour inversion
1315 t = work(ioff1 + i*m+j)
1316 u = -work(ioff1+nwork/4+ i*m+j)
1317 v = work(ioff1 +(n/2-i)*m+j)
1318 w = -work(ioff1+nwork/4+(n/2-i)*m+j)
1319 c = table(itable+i)
1320 s = table(itable+i+n)
1321 tt = (t+v)/2.0_gft_prec
1322 uu = (u-w)/2.0_gft_prec
1323 vv = (c*(t-v)+s*(u+w))/2.0_gft_prec
1324 ww = (c*(u+w)-s*(t-v))/2.0_gft_prec
1325 ! Note : le facteur 2 et le signe - viennent de l'inversion Fourier
1326 work(ioff2 +m*i+j) = 2.0_gft_prec*(tt-ww)
1327 work(ioff2+nwork/4+m*i+j) = -2.0_gft_prec*(uu+vv)
1328 END DO
1329 END DO
1330
1331 END IF
1332
1333 ! On fait m t.f. complexes de taille n/2
1334 fact(ifact+1) = fact(ifact+1)/2 ! Revient a remplacer n2 par n2/2
1335 fact(ifact+2) = fact(ifact+2)-1 ! Revient a remplacer p2 par p2-1
1336 CALL jmccm1d(m,n,fact,nfact,ifact,table,ntable,itable,work,nwork,ioff2)
1337 fact(ifact+1) = fact(ifact+1)*2 ! On retablit les valeurs initiales
1338 fact(ifact+2) = fact(ifact+2)+1
1339 ioff1 = nwork/2 - ioff2
1340
1341 ! On reconstitue
1342
1343 IF (m >= 16 .OR. n/2 < 8) THEN
1344
1345 DO i = 0, n/2-1
1346!dir$ ivdep
1347!ocl novrec
1348!cdir nodep
1349 DO j = 0, m-1
1350 ! Note : le signe - vient de l'inversion
1351 work(ioff1+m*(2*i )+j) = work(ioff2 +m*i+j)
1352 work(ioff1+m*(2*i+1)+j) = -work(ioff2+nwork/4+m*i+j)
1353 END DO
1354 END DO
1355
1356 ELSE
1357
1358 DO j = 0, m-1
1359!dir$ ivdep
1360!ocl novrec
1361!cdir nodep
1362 DO i = 0, n/2-1
1363 ! Note : le signe - vient de l'inversion
1364 work(ioff1+m*(2*i )+j) = work(ioff2 +m*i+j)
1365 work(ioff1+m*(2*i+1)+j) = -work(ioff2+nwork/4+m*i+j)
1366 END DO
1367 END DO
1368
1369 END IF
1370
1371 END IF
1372
1373 ioff = ioff1
1374
1375END SUBROUTINE jmcsm1d
1376
1377! Variante de jmcsm1d ou on fournit x en entree et en sortie
1378SUBROUTINE jmcsm1dxy(m,n,fact,nfact,ifact,table,ntable,itable,work,nwork,x,dimx,debx,incx,jumpx,y,dimy,deby,incy,jumpy,isign,scale)
1379
1380 IMPLICIT NONE
1381
1382 ! Arguments
1383 INTEGER, INTENT(in) :: m, n
1384 INTEGER, INTENT(in) :: nfact, ifact
1385 INTEGER, INTENT(inout), DIMENSION(0:nfact-1) :: fact
1386 INTEGER, INTENT(in) :: ntable,itable
1387 REAL(kind=gft_prec), INTENT(in), DIMENSION(0:ntable-1) :: table
1388 INTEGER, INTENT(in) :: nwork
1389 REAL(kind=gft_prec), INTENT(inout), DIMENSION(0:nwork-1) :: work
1390 INTEGER, INTENT(in) :: dimx, debx, incx, jumpx
1391 INTEGER, INTENT(in) :: dimy, deby, incy, jumpy
1392 REAL(kind=gft_prec), INTENT(in), DIMENSION(0:dimx-1) :: x
1393 REAL(kind=gft_prec), INTENT(out), DIMENSION(0:dimy-1) :: y
1394 INTEGER, INTENT(in) :: isign
1395 REAL(kind=gft_prec), INTENT(in) :: scale
1396
1397 ! Variables locales
1398 INTEGER :: i, j
1399 REAL(kind=gft_prec) :: t, u, v, w, tt, uu, vv, ww
1400 REAL(kind=gft_prec) :: c, s
1401 INTEGER :: it
1402 INTEGER :: ioff
1403
1404 ! On doit faire m T.F. complexes -> reelles de longueur n
1405 ! Si m est pair
1406 IF (mod(m,2) == 0) THEN
1407
1408 ! On distribue
1409
1410 IF (m/2 >= 16 .OR. n/2 < 8) THEN
1411
1412 DO i = 0,n/2
1413!dir$ ivdep
1414!ocl novrec
1415!cdir nodep
1416 DO j = 0,m/2-1
1417 it = n-i
1418 IF (i == 0) it = 0
1419 t = scale*x(debx+incx*(2*i )+jumpx*(j ))
1420 u = isign*scale*x(debx+incx*(2*i+1)+jumpx*(j ))
1421 v = scale*x(debx+incx*(2*i )+jumpx*(j+m/2))
1422 w = isign*scale*x(debx+incx*(2*i+1)+jumpx*(j+m/2))
1423 work( i*m/2+j) = (t-w)
1424 work(nwork/4+ i*m/2+j) = (u+v)
1425 work( it*m/2+j) = (t+w)
1426 work(nwork/4+it*m/2+j) = (v-u)
1427 END DO
1428 END DO
1429
1430 ELSE
1431
1432 DO j = 0,m/2-1
1433!dir$ ivdep
1434!ocl novrec
1435!cdir nodep
1436 DO i = 0,n/2
1437 it = n-i
1438 IF (i == 0) it = 0
1439 t = scale*x(debx+incx*(2*i )+jumpx*(j ))
1440 u = isign*scale*x(debx+incx*(2*i+1)+jumpx*(j ))
1441 v = scale*x(debx+incx*(2*i )+jumpx*(j+m/2))
1442 w = isign*scale*x(debx+incx*(2*i+1)+jumpx*(j+m/2))
1443 work( i*m/2+j) = (t-w)
1444 work(nwork/4+ i*m/2+j) = (u+v)
1445 work( it*m/2+j) = (t+w)
1446 work(nwork/4+it*m/2+j) = (v-u)
1447 END DO
1448 END DO
1449
1450 END IF
1451
1452 ! On fait m/2 t.f. complexes -> complexes de longueur n
1453 ioff = 0
1454 CALL jmccm1d(m/2,n,fact,nfact,ifact,table,ntable,itable,work,nwork,ioff)
1455
1456 ! On reconstitue
1457
1458 IF (m/2 >= 16 .OR. n < 8) THEN
1459
1460 DO i = 0,n-1
1461!dir$ ivdep
1462!ocl novrec
1463!cdir nodep
1464 DO j = 0,m/2-1
1465 y(deby+jumpy*(j )+incy*i) = work(ioff +i*m/2+j)
1466 y(deby+jumpy*(j+m/2)+incy*i) = work(ioff+nwork/4+i*m/2+j)
1467 END DO
1468 END DO
1469
1470 ELSE
1471
1472 DO j = 0,m/2-1
1473!dir$ ivdep
1474!ocl novrec
1475!cdir nodep
1476 DO i = 0,n-1
1477 y(deby+jumpy*(j )+incy*i) = work(ioff +i*m/2+j)
1478 y(deby+jumpy*(j+m/2)+incy*i) = work(ioff+nwork/4+i*m/2+j)
1479 END DO
1480 END DO
1481
1482 END IF
1483
1484 ! Si m n'est pas pair mais que n l'est
1485 ELSE IF (mod(n,2) == 0) THEN
1486
1487 ! On distribue
1488
1489 IF (m >= 16 .OR. n/2 < 8) THEN
1490
1491 DO i = 0,n/2-1
1492!dir$ ivdep
1493!ocl novrec
1494!cdir nodep
1495 DO j = 0,m-1
1496 ! Note : Signe - sur les parties imaginaires pour inversion
1497 t = scale*x(debx+(2*i) *incx+j*jumpx)
1498 u = -isign*scale*x(debx+(2*i+1)*incx+j*jumpx)
1499 v = scale*x(debx+(2*(n/2-i) )*incx+j*jumpx)
1500 w = -isign*scale*x(debx+(2*(n/2-i)+1)*incx+j*jumpx)
1501 c = table(itable+i)
1502 s = table(itable+i+n)
1503 tt = (t+v)/2.0_gft_prec
1504 uu = (u-w)/2.0_gft_prec
1505 vv = (c*(t-v)+s*(u+w))/2.0_gft_prec
1506 ww = (c*(u+w)-s*(t-v))/2.0_gft_prec
1507 ! Note : le facteur 2 et le signe - viennent de l'inversion Fourier
1508 work( m*i+j) = 2.0_gft_prec*(tt-ww)
1509 work(nwork/4+m*i+j) = -2.0_gft_prec*(uu+vv)
1510 END DO
1511 END DO
1512
1513 ELSE
1514
1515 DO j = 0,m-1
1516!dir$ ivdep
1517!ocl novrec
1518!cdir nodep
1519 DO i = 0,n/2-1
1520 ! Note : Signe - sur les parties imaginaires pour inversion
1521 t = scale*x(debx+(2*i) *incx+j*jumpx)
1522 u = -isign*scale*x(debx+(2*i+1)*incx+j*jumpx)
1523 v = scale*x(debx+(2*(n/2-i) )*incx+j*jumpx)
1524 w = -isign*scale*x(debx+(2*(n/2-i)+1)*incx+j*jumpx)
1525 c = table(itable+i)
1526 s = table(itable+i+n)
1527 tt = (t+v)/2.0_gft_prec
1528 uu = (u-w)/2.0_gft_prec
1529 vv = (c*(t-v)+s*(u+w))/2.0_gft_prec
1530 ww = (c*(u+w)-s*(t-v))/2.0_gft_prec
1531 ! Note : le facteur 2 et le signe - viennent de l'inversion Fourier
1532 work( m*i+j) = 2.0_gft_prec*(tt-ww)
1533 work(nwork/4+m*i+j) = -2.0_gft_prec*(uu+vv)
1534 END DO
1535 END DO
1536
1537 END IF
1538
1539 ! On fait m t.f. complexes de taille n/2
1540 ioff = 0
1541 fact(ifact+1) = fact(ifact+1)/2 ! Revient a remplacer n2 par n2/2
1542 fact(ifact+2) = fact(ifact+2)-1 ! Revient a remplacer p2 par p2-1
1543 CALL jmccm1d(m,n,fact,nfact,ifact,table,ntable,itable,work,nwork,ioff)
1544 fact(ifact+1) = fact(ifact+1)*2 ! On retablit les valeurs initiales
1545 fact(ifact+2) = fact(ifact+2)+1
1546
1547 ! On reconstitue
1548
1549 IF (m >= 16 .OR. n/2 < 8) THEN
1550
1551 DO i = 0, n/2-1
1552!dir$ ivdep
1553!ocl novrec
1554!cdir nodep
1555 DO j = 0, m-1
1556 ! Note : le signe - vient de l'inversion
1557 y(deby+incy*(2*i )+jumpy*j) = work(ioff +m*i+j)
1558 y(deby+incy*(2*i+1)+jumpy*j) = -work(ioff+nwork/4+m*i+j)
1559 END DO
1560 END DO
1561
1562 ELSE
1563
1564!dir$ ivdep
1565 DO j = 0, m-1
1566!ocl novrec
1567!cdir nodep
1568 DO i = 0, n/2-1
1569 ! Note : le signe - vient de l'inversion
1570 y(deby+incy*(2*i )+jumpy*j) = work(ioff +m*i+j)
1571 y(deby+incy*(2*i+1)+jumpy*j) = -work(ioff+nwork/4+m*i+j)
1572 END DO
1573 END DO
1574
1575 END IF
1576
1577 END IF
1578
1579END SUBROUTINE jmcsm1dxy
1580
1581SUBROUTINE jmscm1d(m,n,fact,nfact,ifact,table,ntable,itable,work,nwork,ioff)
1582
1583 IMPLICIT NONE
1584
1585 ! Arguments
1586 INTEGER, INTENT(in) :: m, n
1587 INTEGER, INTENT(in) :: nfact, ifact
1588 INTEGER, INTENT(inout), DIMENSION(0:nfact-1) :: fact
1589 INTEGER, INTENT(in) :: ntable,itable
1590 REAL(kind=gft_prec), INTENT(in), DIMENSION(0:ntable-1) :: table
1591 INTEGER, INTENT(in) :: nwork
1592 REAL(kind=gft_prec), INTENT(inout), DIMENSION(0:nwork-1) :: work
1593 INTEGER, INTENT(inout) :: ioff
1594
1595 ! Variables locales
1596 INTEGER :: ioff1, ioff2
1597 INTEGER :: i, j
1598 REAL(kind=gft_prec) :: t, u, v, w
1599 REAL(kind=gft_prec) :: c, s
1600 INTEGER :: is, it
1601
1602 ! Gestion de work
1603 ioff1 = ioff
1604 ioff2 = nwork/2 - ioff1
1605
1606 ! On doit faire m T.F. reelles de longueur n
1607 ! Si m est pair
1608 IF (mod(m,2) == 0) THEN
1609
1610 ! On distribue
1611 IF (m/2 >= 16 .OR. n < 8) THEN
1612
1613 DO i = 0,n-1
1614!dir$ ivdep
1615!ocl novrec
1616!cdir nodep
1617 DO j = 0,m/2-1
1618 work(ioff2 +i*m/2+j) = work(ioff1+i*m+j )
1619 work(ioff2+nwork/4+i*m/2+j) = work(ioff1+i*m+j+m/2)
1620 END DO
1621 END DO
1622
1623 ELSE
1624
1625 DO j = 0,m/2-1
1626!dir$ ivdep
1627!ocl novrec
1628!cdir nodep
1629 DO i = 0,n-1
1630 work(ioff2 +i*m/2+j) = work(ioff1+i*m+j )
1631 work(ioff2+nwork/4+i*m/2+j) = work(ioff1+i*m+j+m/2)
1632 END DO
1633 END DO
1634
1635 END IF
1636
1637 ! On fait m/2 t.f. complexes de longueur n
1638 CALL jmccm1d(m/2,n,fact,nfact,ifact,table,ntable,itable,work,nwork,ioff2)
1639 ioff1 = nwork/2 - ioff2
1640
1641 ! On regenere le resultat
1642 IF (m/2 >= 16 .OR. n/2 < 8) THEN
1643
1644 DO i = 0,n/2
1645!dir$ ivdep
1646!ocl novrec
1647!cdir nodep
1648 DO j = 0,m/2-1
1649 it = n-i
1650 IF (i == 0) it = 0
1651 t = work(ioff2 + i*m/2+j)
1652 u = work(ioff2+nwork/4+ i*m/2+j)
1653 v = work(ioff2 +it*m/2+j)
1654 w = work(ioff2+nwork/4+it*m/2+j)
1655 work(ioff1 +i*m+j ) = (t+v)/2.0_gft_prec
1656 work(ioff1+nwork/4+i*m+j ) = (u-w)/2.0_gft_prec
1657 work(ioff1 +i*m+j+m/2) = (u+w)/2.0_gft_prec
1658 work(ioff1+nwork/4+i*m+j+m/2) = (v-t)/2.0_gft_prec
1659 END DO
1660 END DO
1661
1662 ELSE
1663
1664 DO j = 0,m/2-1
1665!dir$ ivdep
1666!ocl novrec
1667!cdir nodep
1668 DO i = 0,n/2
1669 it = n-i
1670 IF (i == 0) it = 0
1671 t = work(ioff2 + i*m/2+j)
1672 u = work(ioff2+nwork/4+ i*m/2+j)
1673 v = work(ioff2 +it*m/2+j)
1674 w = work(ioff2+nwork/4+it*m/2+j)
1675 work(ioff1 +i*m+j ) = (t+v)/2.0_gft_prec
1676 work(ioff1+nwork/4+i*m+j ) = (u-w)/2.0_gft_prec
1677 work(ioff1 +i*m+j+m/2) = (u+w)/2.0_gft_prec
1678 work(ioff1+nwork/4+i*m+j+m/2) = (v-t)/2.0_gft_prec
1679 END DO
1680 END DO
1681
1682 END IF
1683
1684 ! Si m n'est pas pair mais que n l'est
1685 ELSE IF (mod(n,2) == 0) THEN
1686
1687 ! On distribue les indices pairs et impairs selon n
1688
1689 IF (m >= 16 .OR. n/2 < 8) THEN
1690
1691 DO i = 0, n/2-1
1692!dir$ ivdep
1693!ocl novrec
1694!cdir nodep
1695 DO j = 0, m-1
1696 work(ioff2+ m*i+j) = work(ioff1+m*(2*i )+j)
1697 work(ioff2+nwork/4+m*i+j) = work(ioff1+m*(2*i+1)+j)
1698 END DO
1699 END DO
1700
1701 ELSE
1702
1703 DO j = 0, m-1
1704!dir$ ivdep
1705!ocl novrec
1706!cdir nodep
1707 DO i = 0, n/2-1
1708 work(ioff2+ m*i+j) = work(ioff1+m*(2*i )+j)
1709 work(ioff2+nwork/4+m*i+j) = work(ioff1+m*(2*i+1)+j)
1710 END DO
1711 END DO
1712
1713 END IF
1714
1715 ! On fait m t.f. complexes de taille n/2
1716 fact(ifact+1) = fact(ifact+1)/2 ! Revient a remplacer n2 par n2/2
1717 fact(ifact+2) = fact(ifact+2)-1 ! Revient a remplacer p2 par p2-1
1718 CALL jmccm1d(m,n,fact,nfact,ifact,table,ntable,itable,work,nwork,ioff2)
1719 fact(ifact+1) = fact(ifact+1)*2 ! On retablit les valeurs initiales
1720 fact(ifact+2) = fact(ifact+2)+1
1721 ioff1 = nwork/2 - ioff2
1722
1723 ! Maintenant, il faut reconstituer la t.f. reelle
1724
1725 IF (m >= 16 .OR. n/2 < 8) THEN
1726
1727 DO i = 0,n/2
1728!dir$ ivdep
1729!ocl novrec
1730!cdir nodep
1731 DO j = 0,m-1
1732 is = i
1733 it = n/2-i
1734 IF (i == 0 .OR. i == n/2) THEN
1735 is = 0
1736 it = 0
1737 END IF
1738 t = work(ioff2 +is*m+j)
1739 u = work(ioff2+nwork/4+is*m+j)
1740 v = work(ioff2 +it*m+j)
1741 w = work(ioff2+nwork/4+it*m+j)
1742 c = table(itable+i)
1743 s = table(itable+i+n)
1744 work(ioff1 +i*m+j) = (t+v)/2.0_gft_prec + c*(u+w)/2.0_gft_prec - s*(v-t)/2.0_gft_prec
1745 work(ioff1+nwork/4+i*m+j) = (u-w)/2.0_gft_prec + c*(v-t)/2.0_gft_prec + s*(u+w)/2.0_gft_prec
1746 END DO
1747 END DO
1748
1749 ELSE
1750
1751 DO j = 0,m-1
1752!dir$ ivdep
1753!ocl novrec
1754!cdir nodep
1755 DO i = 0,n/2
1756 is = i
1757 it = n/2-i
1758 IF (i == 0 .OR. i == n/2) THEN
1759 is = 0
1760 it = 0
1761 END IF
1762 t = work(ioff2 +is*m+j)
1763 u = work(ioff2+nwork/4+is*m+j)
1764 v = work(ioff2 +it*m+j)
1765 w = work(ioff2+nwork/4+it*m+j)
1766 c = table(itable+i)
1767 s = table(itable+i+n)
1768 work(ioff1 +i*m+j) = (t+v)/2.0_gft_prec + c*(u+w)/2.0_gft_prec - s*(v-t)/2.0_gft_prec
1769 work(ioff1+nwork/4+i*m+j) = (u-w)/2.0_gft_prec + c*(v-t)/2.0_gft_prec + s*(u+w)/2.0_gft_prec
1770 END DO
1771 END DO
1772
1773 END IF
1774
1775 END IF
1776
1777 ioff = ioff1
1778
1779END SUBROUTINE jmscm1d
1780
1781! Variante de jmscm1d ou on fournit x en entree et en sortie
1782SUBROUTINE jmscm1dxy(m,n,fact,nfact,ifact,table,ntable,itable,work,nwork,x,dimx,debx,incx,jumpx,y,dimy,deby,incy,jumpy,isign,scale)
1783
1784 IMPLICIT NONE
1785
1786 ! Arguments
1787 INTEGER, INTENT(in) :: m, n
1788 INTEGER, INTENT(in) :: nfact, ifact
1789 INTEGER, INTENT(inout), DIMENSION(0:nfact-1) :: fact
1790 INTEGER, INTENT(in) :: ntable,itable
1791 REAL(kind=gft_prec), INTENT(in), DIMENSION(0:ntable-1) :: table
1792 INTEGER, INTENT(in) :: nwork
1793 REAL(kind=gft_prec), INTENT(inout), DIMENSION(0:nwork-1) :: work
1794 INTEGER, INTENT(in) :: dimx, debx, incx, jumpx
1795 INTEGER, INTENT(in) :: dimy, deby, incy, jumpy
1796 REAL(kind=gft_prec), INTENT(in) :: scale
1797 REAL(kind=gft_prec), INTENT(in), DIMENSION(0:dimx-1) :: x
1798 REAL(kind=gft_prec), INTENT(out), DIMENSION(0:dimy-1) :: y
1799 INTEGER, INTENT(in) :: isign
1800
1801 ! Variables locales
1802 INTEGER :: i, j
1803 REAL(kind=gft_prec) :: t, u, v, w
1804 REAL(kind=gft_prec) :: c, s
1805 INTEGER :: is, it
1806 INTEGER :: ioff
1807
1808 ! On doit faire m T.F. reelles de longueur n
1809 ! Si m est pair
1810 IF (mod(m,2) == 0) THEN
1811
1812 ! On distribue
1813 IF (m/2 >= 16 .OR. n < 8) THEN
1814
1815 DO i = 0,n-1
1816!dir$ ivdep
1817!ocl novrec
1818!cdir nodep
1819 DO j = 0,m/2-1
1820 work( i*m/2+j) = x(debx+i*incx+(j) *jumpx)
1821 work(nwork/4+i*m/2+j) = x(debx+i*incx+(j+m/2)*jumpx)
1822 END DO
1823 END DO
1824
1825 ELSE
1826
1827 DO j = 0,m/2-1
1828!dir$ ivdep
1829!ocl novrec
1830!cdir nodep
1831 DO i = 0,n-1
1832 work( i*m/2+j) = x(debx+i*incx+(j) *jumpx)
1833 work(nwork/4+i*m/2+j) = x(debx+i*incx+(j+m/2)*jumpx)
1834 END DO
1835 END DO
1836
1837 END IF
1838
1839 ! On fait m/2 t.f. complexes de longueur n
1840 ioff = 0
1841 CALL jmccm1d(m/2,n,fact,nfact,ifact,table,ntable,itable,work,nwork,ioff)
1842
1843 ! On regenere le resultat
1844 IF (m/2 >= 16 .OR. n/2 < 8) THEN
1845
1846 DO i = 0,n/2
1847!dir$ ivdep
1848!ocl novrec
1849!cdir nodep
1850 DO j = 0,m/2-1
1851 it = n-i
1852 IF (i == 0) it = 0
1853 t = work(ioff + i*m/2+j)
1854 u = work(ioff+nwork/4+ i*m/2+j)
1855 v = work(ioff +it*m/2+j)
1856 w = work(ioff+nwork/4+it*m/2+j)
1857 y(deby+(2*i) *incy+(j) *jumpy) = scale*(t+v)/2.0_gft_prec
1858 y(deby+(2*i+1)*incy+(j) *jumpy) = isign*scale*(u-w)/2.0_gft_prec
1859 y(deby+(2*i) *incy+(j+m/2)*jumpy) = scale*(u+w)/2.0_gft_prec
1860 y(deby+(2*i+1)*incy+(j+m/2)*jumpy) = isign*scale*(v-t)/2.0_gft_prec
1861 END DO
1862 END DO
1863
1864 ELSE
1865
1866 DO j = 0,m/2-1
1867!dir$ ivdep
1868!ocl novrec
1869!cdir nodep
1870 DO i = 0,n/2
1871 it = n-i
1872 IF (i == 0) it = 0
1873 t = work(ioff + i*m/2+j)
1874 u = work(ioff+nwork/4+ i*m/2+j)
1875 v = work(ioff +it*m/2+j)
1876 w = work(ioff+nwork/4+it*m/2+j)
1877 y(deby+(2*i) *incy+(j) *jumpy) = scale*(t+v)/2.0_gft_prec
1878 y(deby+(2*i+1)*incy+(j) *jumpy) = isign*scale*(u-w)/2.0_gft_prec
1879 y(deby+(2*i) *incy+(j+m/2)*jumpy) = scale*(u+w)/2.0_gft_prec
1880 y(deby+(2*i+1)*incy+(j+m/2)*jumpy) = isign*scale*(v-t)/2.0_gft_prec
1881 END DO
1882 END DO
1883
1884 END IF
1885
1886 ! Si m n'est pas pair mais que n l'est
1887 ELSE IF (mod(n,2) == 0) THEN
1888
1889 ! On distribue les indices pairs et impairs selon n
1890
1891 IF (m >= 16 .OR. n/2 < 8) THEN
1892
1893 DO i = 0, n/2-1
1894!dir$ ivdep
1895!ocl novrec
1896!cdir nodep
1897 DO j = 0, m-1
1898 work( m*i+j) = x(debx+incx*(2*i )+jumpx*j)
1899 work(nwork/4+m*i+j) = x(debx+incx*(2*i+1)+jumpx*j)
1900 END DO
1901 END DO
1902
1903 ELSE
1904
1905 DO j = 0, m-1
1906!dir$ ivdep
1907!ocl novrec
1908!cdir nodep
1909 DO i = 0, n/2-1
1910 work( m*i+j) = x(debx+incx*(2*i )+jumpx*j)
1911 work(nwork/4+m*i+j) = x(debx+incx*(2*i+1)+jumpx*j)
1912 END DO
1913 END DO
1914
1915 END IF
1916
1917 ! On fait m t.f. complexes de taille n/2
1918 ioff = 0
1919 fact(ifact+1) = fact(ifact+1)/2 ! Revient a remplacer n2 par n2/2
1920 fact(ifact+2) = fact(ifact+2)-1 ! Revient a remplacer p2 par p2-1
1921 CALL jmccm1d(m,n,fact,nfact,ifact,table,ntable,itable,work,nwork,ioff)
1922 fact(ifact+1) = fact(ifact+1)*2 ! On retablit les valeurs initiales
1923 fact(ifact+2) = fact(ifact+2)+1
1924
1925 ! Maintenant, il faut reconstituer la t.f. reelle
1926
1927 IF (m >= 16 .OR. n/2 < 8) THEN
1928
1929 DO i = 0,n/2
1930!dir$ ivdep
1931!ocl novrec
1932!cdir nodep
1933 DO j = 0,m-1
1934 is = i
1935 it = n/2-i
1936 IF (i == 0 .OR. i == n/2) THEN
1937 is = 0
1938 it = 0
1939 END IF
1940 t = work(ioff +is*m+j)
1941 u = work(ioff+nwork/4+is*m+j)
1942 v = work(ioff +it*m+j)
1943 w = work(ioff+nwork/4+it*m+j)
1944 c = table(itable+i)
1945 s = table(itable+i+n)
1946 y(deby+(2*i )*incy+j*jumpy) = &
1947 & scale*((t+v)/2.0_gft_prec + c*(u+w)/2.0_gft_prec - s*(v-t)/2.0_gft_prec)
1948 y(deby+(2*i+1)*incy+j*jumpy) = &
1949 & isign*scale*((u-w)/2.0_gft_prec + c*(v-t)/2.0_gft_prec + s*(u+w)/2.0_gft_prec)
1950 END DO
1951 END DO
1952
1953 ELSE
1954
1955 DO j = 0,m-1
1956!dir$ ivdep
1957!ocl novrec
1958!cdir nodep
1959 DO i = 0,n/2
1960 is = i
1961 it = n/2-i
1962 IF (i == 0 .OR. i == n/2) THEN
1963 is = 0
1964 it = 0
1965 END IF
1966 t = work(ioff +is*m+j)
1967 u = work(ioff+nwork/4+is*m+j)
1968 v = work(ioff +it*m+j)
1969 w = work(ioff+nwork/4+it*m+j)
1970 c = table(itable+i)
1971 s = table(itable+i+n)
1972 y(deby+(2*i )*incy+j*jumpy) = &
1973 & scale*((t+v)/2.0_gft_prec + c*(u+w)/2.0_gft_prec - s*(v-t)/2.0_gft_prec)
1974 y(deby+(2*i+1)*incy+j*jumpy) = &
1975 & isign*scale*((u-w)/2.0_gft_prec + c*(v-t)/2.0_gft_prec + s*(u+w)/2.0_gft_prec)
1976 END DO
1977 END DO
1978
1979 END IF
1980
1981 END IF
1982
1983END SUBROUTINE jmscm1dxy
1984
1985SUBROUTINE jmtransp(n,m,l,work,nwork,ioff)
1986
1987 IMPLICIT NONE
1988
1989 ! Arguments
1990 INTEGER, INTENT(in) :: n, m, l
1991 INTEGER, INTENT(in) :: nwork
1992 REAL(kind=gft_prec), INTENT(inout), DIMENSION(0:nwork-1) :: work
1993 INTEGER, INTENT(inout) :: ioff
1994
1995 ! Variables locales
1996 INTEGER :: ioff1, ioff2
1997 INTEGER :: ij, k
1998
1999 ioff1 = ioff
2000 ioff2 = nwork/2-ioff1
2001
2002 ! On transpose (nm)(l) en (l)(nm) en distinguant les parties reelles et im.
2003 IF (m*n >= 16 .OR. l < 8) THEN
2004
2005 DO k = 0,l-1
2006!dir$ ivdep
2007!ocl novrec
2008!cdir nodep
2009 DO ij = 0,m*n-1
2010 work(ioff2+ ij*l+k) = work(ioff1+ k*n*m+ij)
2011 work(ioff2+n*m*l+ij*l+k) = work(ioff1+n*m*l+k*n*m+ij)
2012 END DO
2013 END DO
2014
2015 ELSE
2016
2017 DO ij = 0,m*n-1
2018!dir$ ivdep
2019!ocl novrec
2020!cdir nodep
2021 DO k = 0,l-1
2022 work(ioff2+ ij*l+k) = work(ioff1+ k*n*m+ij)
2023 work(ioff2+n*m*l+ij*l+k) = work(ioff1+n*m*l+k*n*m+ij)
2024 END DO
2025 END DO
2026
2027 END IF
2028
2029 ioff = ioff2
2030
2031END SUBROUTINE jmtransp
2032
2033SUBROUTINE jmcctranspcs(m,n,n2,n3,table,ntable,itable,work,nwork,ioff)
2034
2035 ! Cette subroutine permute le contenu du tableau work de la facon suivante
2036 ! On considere qu'a l'origine ce tableau est en (m,n3,n2)
2037 ! On doit transposer le terme d'ordre (k,j,i) en (k,i,j)
2038 ! On en profite pour faire les multiplications par wij
2039 ! Le role de n est seulement d'attaquer les bonnes valeurs du tableau table
2040 ! (il y a un ecart de n entre les cos et les sin, et le stride entre
2041 ! les cos est de n/(n2*n3)
2042 ! Note : le sens de n2 et n3 ici n'a rien a voir avec celui de jmccm1d
2043
2044 IMPLICIT NONE
2045
2046 ! Arguments
2047 INTEGER, INTENT(in) :: m, n
2048 INTEGER, INTENT(in) :: n2, n3
2049 INTEGER, INTENT(in) :: ntable,itable
2050 REAL(kind=gft_prec), INTENT(in), DIMENSION(0:ntable-1) :: table
2051 INTEGER, INTENT(in) :: nwork
2052 REAL(kind=gft_prec), INTENT(inout), DIMENSION(0:nwork-1) :: work
2053 INTEGER, INTENT(inout) :: ioff
2054
2055 ! Variables locales
2056 INTEGER :: i, j, k
2057 REAL(kind=gft_prec) :: t, u, c, s
2058 INTEGER :: ioff1, ioff2
2059 INTEGER :: is
2060
2061 ! Gestion des offsets
2062 ioff1 = ioff
2063 ioff2 = nwork/2-ioff
2064
2065 ! Gestion du stride
2066 is = n/(n2*n3)
2067
2068 IF ( m >= 16 .OR. (n2 < 8 .AND. n3 < 8) ) THEN
2069
2070 DO i = 0,n2-1
2071 DO j = 0,n3-1
2072!dir$ ivdep
2073!ocl novrec
2074!cdir nodep
2075 DO k = 0,m-1
2076 t = work(ioff1 +k+m*(j+n3*i))
2077 u = work(ioff1+nwork/4+k+m*(j+n3*i))
2078 c = table(itable+ is*i*j)
2079 s = table(itable+n+is*i*j)
2080 work(ioff2 +k+m*(i+n2*j)) = c*t-s*u
2081 work(ioff2+nwork/4+k+m*(i+n2*j)) = c*u+s*t
2082 END DO
2083 END DO
2084 END DO
2085
2086 ELSE IF ( n2 >= 16 .OR. n3 < 8 ) THEN
2087
2088 DO j = 0,n3-1
2089 DO k = 0,m-1
2090!dir$ ivdep
2091!ocl novrec
2092!cdir nodep
2093 DO i = 0,n2-1
2094 t = work(ioff1 +k+m*(j+n3*i))
2095 u = work(ioff1+nwork/4+k+m*(j+n3*i))
2096 c = table(itable+ is*i*j)
2097 s = table(itable+n+is*i*j)
2098 work(ioff2 +k+m*(i+n2*j)) = c*t-s*u
2099 work(ioff2+nwork/4+k+m*(i+n2*j)) = c*u+s*t
2100 END DO
2101 END DO
2102 END DO
2103
2104 ELSE
2105
2106 DO i = 0,n2-1
2107 DO k = 0,m-1
2108!dir$ ivdep
2109!ocl novrec
2110!cdir nodep
2111 DO j = 0,n3-1
2112 t = work(ioff1 +k+m*(j+n3*i))
2113 u = work(ioff1+nwork/4+k+m*(j+n3*i))
2114 c = table(itable+ is*i*j)
2115 s = table(itable+n+is*i*j)
2116 work(ioff2 +k+m*(i+n2*j)) = c*t-s*u
2117 work(ioff2+nwork/4+k+m*(i+n2*j)) = c*u+s*t
2118 END DO
2119 END DO
2120 END DO
2121
2122 END IF
2123
2124 ioff = ioff2
2125
2126END SUBROUTINE jmcctranspcs
2127
2128SUBROUTINE jmsetnwork(nwork)
2129
2130 ! Subroutine appelee par l'utilisateur pour augmenter le nwork
2131 ! des routines 2d et 3d
2132
2133 IMPLICIT NONE
2134
2135 ! Arguments
2136 INTEGER, INTENT(in) :: nwork
2137
2138 ! Variables locales
2139 CHARACTER(len=*), PARAMETER :: nomsp = 'JMSETNWORK'
2140 INTEGER :: nwork2
2141
2142 IF (nwork <= 0) THEN
2143 CALL jmerreur1(nomsp,4,nwork)
2144 END IF
2145
2146 nwork2 = nwork
2147 CALL jmgetsetnwork(nwork2,'s')
2148
2149END SUBROUTINE jmsetnwork
2150
2151SUBROUTINE jmgetnwork(nwork,nwork_def,nwork_min)
2152
2153 ! On recupere la valeur de nwork si elle a ete augmentee par l'utilisateur
2154 ! Sinon on prend la valeur par defaut
2155 ! Il s'agit du nwork des routines 2d et 3d
2156
2157 IMPLICIT NONE
2158
2159 ! Arguments
2160 INTEGER, INTENT(out) :: nwork
2161 INTEGER, INTENT(in) :: nwork_def, nwork_min
2162
2163 ! Variables locales
2164 INTEGER :: nwork_loc
2165 CHARACTER(len=*), PARAMETER :: nomsp = 'JMGETNWORK'
2166
2167 CALL jmgetsetnwork(nwork_loc,'g')
2168
2169 ! Valeur par defaut
2170 IF (nwork_loc == -1) THEN
2171 nwork = nwork_def
2172 ! Valeur invalide (trop petite)
2173 ELSE IF (nwork_loc < nwork_min) THEN
2174 CALL jmerreur2(nomsp,5,nwork_loc,nwork_min)
2175 ! Valeur correcte
2176 ELSE
2177 nwork = nwork_loc
2178 END IF
2179
2180END SUBROUTINE jmgetnwork
2181
2182
2183SUBROUTINE jmgetsetnwork(nwork,TYPE)
2184
2185 ! Subroutine qui permet de stocker une valeur statique
2186 ! Ceci evite de recourir a un common ...
2187
2188 IMPLICIT NONE
2189
2190 ! Arguments
2191 INTEGER, INTENT(inout) :: nwork
2192 CHARACTER(len=1), INTENT(in) :: TYPE
2193
2194 ! Variables locales
2195
2196 ! Variable statique
2197 INTEGER, SAVE :: nwork_last = -1
2198
2199 IF (TYPE == 's') then
2200 nwork_last = nwork
2201 ELSE IF (TYPE == 'g') then
2202 nwork = nwork_last
2203 END IF
2204
2205END SUBROUTINE jmgetsetnwork
2206
2207SUBROUTINE jmgeterreur(arret)
2208
2209 IMPLICIT NONE
2210
2211 ! Arguments
2212 LOGICAL, INTENT(out) :: arret
2213
2214 ! Variables locales
2215
2216 CALL jmgetseterreur(arret,'g')
2217
2218END SUBROUTINE jmgeterreur
2219
2220SUBROUTINE jmerreur1(nomsp,code,var1)
2221
2222 IMPLICIT NONE
2223
2224 ! Arguments
2225 CHARACTER(len=*), INTENT(in) :: nomsp
2226 INTEGER, INTENT(in) :: code
2227 INTEGER, INTENT(in) :: var1
2228
2229 ! Variables locales
2230 INTEGER :: arret
2231 CHARACTER(len=80) :: message
2232
2233 CALL jmgetstop(arret)
2234 IF (arret == 1) THEN
2235 CALL jmgetmessage(code,message)
2236 print *,'JMFFT error in ',trim(nomsp),' : ',trim(message), &
2237 & ' (',var1,')'
2238 stop 1
2239 ELSE
2240 CALL jmsetcode(code)
2241 END IF
2242
2243END SUBROUTINE jmerreur1
2244
2245
2246SUBROUTINE jmerreur2(nomsp,code,var1,var2)
2247
2248 IMPLICIT NONE
2249
2250 ! Arguments
2251 CHARACTER(len=*), INTENT(in) :: nomsp
2252 INTEGER, INTENT(in) :: code
2253 INTEGER, INTENT(in) :: var1, var2
2254
2255 ! Variables locales
2256 INTEGER :: arret
2257 CHARACTER(len=80) :: message
2258
2259 CALL jmgetstop(arret)
2260 IF (arret == 1) THEN
2261 CALL jmgetmessage(code,message)
2262 print *,'JMFFT error in ',trim(nomsp),' : ',trim(message), &
2263 & ' (',var1,var2,')'
2264 stop 1
2265 ELSE
2266 CALL jmsetcode(code)
2267 END IF
2268
2269END SUBROUTINE jmerreur2
2270
2271SUBROUTINE jmerreur3(nomsp,code,var1,var2,var3)
2272
2273 IMPLICIT NONE
2274
2275 ! Arguments
2276 CHARACTER(len=*), INTENT(in) :: nomsp
2277 INTEGER, INTENT(in) :: code
2278 INTEGER, INTENT(in) :: var1, var2, var3
2279
2280 ! Variables locales
2281 INTEGER :: arret
2282 CHARACTER(len=80) :: message
2283
2284 CALL jmgetstop(arret)
2285 IF (arret == 1) THEN
2286 CALL jmgetmessage(code,message)
2287 print *,'JMFFT error in ',trim(nomsp),' : ',trim(message), &
2288 & ' (',var1,var2,var3,')'
2289 stop 1
2290 ELSE
2291 CALL jmsetcode(code)
2292 END IF
2293
2294END SUBROUTINE jmerreur3
2295
2296SUBROUTINE jmerreur4(nomsp,code,var1,var2,var3,var4)
2297
2298 IMPLICIT NONE
2299
2300 ! Arguments
2301 CHARACTER(len=*), INTENT(in) :: nomsp
2302 INTEGER, INTENT(in) :: code
2303 INTEGER, INTENT(in) :: var1, var2, var3,var4
2304
2305 ! Variables locales
2306 INTEGER :: arret
2307 CHARACTER(len=80) :: message
2308
2309 CALL jmgetstop(arret)
2310 IF (arret == 1) THEN
2311 CALL jmgetmessage(code,message)
2312 print *,'JMFFT error in ',trim(nomsp),' : ',trim(message), &
2313 & ' (',var1,var2,var3,var4,')'
2314 stop 1
2315 ELSE
2316 CALL jmsetcode(code)
2317 END IF
2318
2319END SUBROUTINE jmerreur4
2320
2321SUBROUTINE jmgetstop(arret)
2322
2323 IMPLICIT NONE
2324
2325 ! Arguments
2326 INTEGER, INTENT(out) :: arret
2327
2328 ! Variables locales
2329
2330 CALL jmgetsetstop(arret,'g')
2331
2332END SUBROUTINE jmgetstop
2333
2334SUBROUTINE jmgetsetstop(arret,TYPE)
2335
2336 ! Subroutine qui permet de stocker une valeur statique
2337 ! Ceci evite de recourir a un common ...
2338
2339 IMPLICIT NONE
2340
2341 ! Arguments
2342 INTEGER, INTENT(inout) :: arret
2343 CHARACTER(len=1), INTENT(in) :: TYPE
2344
2345 ! Variables locales
2346
2347 ! Variable statique
2348 INTEGER, SAVE :: arret_last = 1
2349
2350 IF (TYPE == 's') then
2351 arret_last = arret
2352 ELSE IF (TYPE == 'g') then
2353 arret = arret_last
2354 END IF
2355
2356END SUBROUTINE jmgetsetstop
2357
2358SUBROUTINE jmgetmessage(code,message)
2359
2360 IMPLICIT NONE
2361
2362 ! Arguments
2363 INTEGER, INTENT(in) :: code
2364 CHARACTER(len=*), INTENT(out) :: message
2365
2366 ! Variables locales
2367 INTEGER, PARAMETER :: mm = 26
2368 CHARACTER(len=34), DIMENSION(0:mm-1) :: messages = (/ &
2369 & "No error ", &
2370 & "Isign must be equal -1 or 1 ", &
2371 & "Isign must be equal 0, -1 or 1 ", &
2372 & "Prime numbers to large ", &
2373 & "Size of work is less or equal 0 ", &
2374 & "Size of work is to small ", &
2375 & "Impossible to decompose ", &
2376 & "Too many prime factors ", &
2377 & "Nz must be >= 1 ", &
2378 & "ldx doit etre >= n ", &
2379 & "ldx doit etre >= n/2+1 ", &
2380 & "ldx1 doit etre >= n ", &
2381 & "ldx1 doit etre >= n/2+1 ", &
2382 & "ldx2 doit etre >= m ", &
2383 & "ldy doit etre >= n ", &
2384 & "ldy doit etre >= n+2 ", &
2385 & "ldy doit etre >= n/2+1 ", &
2386 & "ldy1 doit etre >= n ", &
2387 & "ldy1 doit etre >= n+2 ", &
2388 & "ldy1 doit etre >= n/2+1 ", &
2389 & "ldy2 doit etre >= m ", &
2390 & "Nz must be >= 1 ", &
2391 & "Ny or Nx must be even ", &
2392 & "Nx must be >= 1 ", &
2393 & "Nx must be even ", &
2394 & "Nx or Ny or Nz must be even " &
2395 & /)
2396
2397 IF (code < 0 .OR. code >= mm) THEN
2398 print *,'JMFFT GETMESSAGE invalid code : ',code
2399 END IF
2400
2401 message = messages(code)
2402
2403END SUBROUTINE jmgetmessage
2404
2405SUBROUTINE jmsetcode(code)
2406
2407 IMPLICIT NONE
2408
2409 ! Arguments
2410 INTEGER, INTENT(in) :: code
2411
2412 ! Variables locales
2413 INTEGER :: errcode
2414
2415 errcode = code
2416 CALL jmgetsetcode(errcode,'s')
2417
2418END SUBROUTINE jmsetcode
2419
2420SUBROUTINE jmgetsetcode(code,TYPE)
2421
2422 ! Subroutine qui permet de stocker le dernier code de retour obtenu
2423 ! Ceci evite de recourir a un common ...
2424
2425 IMPLICIT NONE
2426
2427 ! Arguments
2428 INTEGER, INTENT(inout) :: code
2429 CHARACTER(len=1), INTENT(in) :: TYPE
2430
2431 ! Variables locales
2432
2433 ! Variable statique
2434 INTEGER, SAVE :: code_last = 0
2435
2436 IF (TYPE == 's') then
2437 code_last = code
2438 ELSE IF (TYPE == 'g') then
2439 code = code_last
2440 END IF
2441
2442END SUBROUTINE jmgetsetcode
2443
2444SUBROUTINE jmgetseterreur(arret,TYPE)
2445
2446 ! Subroutine qui permet de stocker une valeur statique
2447 ! Ceci evite de recourir a un common ...
2448
2449 IMPLICIT NONE
2450
2451 ! Arguments
2452 LOGICAL, INTENT(inout) :: arret
2453 CHARACTER(len=1), INTENT(in) :: TYPE
2454
2455 ! Variables locales
2456
2457 ! Variable statique
2458 LOGICAL, SAVE :: arret_last = .true.
2459
2460 IF (TYPE == 's') then
2461 arret_last = arret
2462 ELSE IF (TYPE == 'g') then
2463 arret = arret_last
2464 END IF
2465
2466END SUBROUTINE jmgetseterreur
2467
2468END MODULE gft_jmfft
gft_common
Definition GFT_common.f90:15
gft_common::gft_prec
integer, parameter gft_prec
Definition GFT_common.f90:24
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::jmccm1d5
subroutine jmccm1d5(p, n, m, table, ntable, itable, ntable2, mtable, work, nwork, ioff)
Definition GFT_jmfft.f90:820
gft_jmfft::jmccm1d
subroutine jmccm1d(m, n, fact, nfact, ifact, table, ntable, itable, work, nwork, ioff)
Definition GFT_jmfft.f90:335
gft_jmfft::jmccm1d4
subroutine jmccm1d4(p, n, m, table, ntable, itable, ntable2, mtable, work, nwork, ioff)
Definition GFT_jmfft.f90:670
gft_jmfft::jmgetsetcode
subroutine jmgetsetcode(code, type)
Definition GFT_jmfft.f90:2421
gft_jmfft::jmgetmessage
subroutine jmgetmessage(code, message)
Definition GFT_jmfft.f90:2359
gft_jmfft::jmccm1dp
subroutine jmccm1dp(p, q, m, table, ntable, itable, ntable2, mtable, work, nwork, ioff)
Definition GFT_jmfft.f90:1066
gft_jmfft::jmscm1d
subroutine jmscm1d(m, n, fact, nfact, ifact, table, ntable, itable, work, nwork, ioff)
Definition GFT_jmfft.f90:1582
gft_jmfft::jmgetseterreur
subroutine jmgetseterreur(arret, type)
Definition GFT_jmfft.f90:2445
gft_jmfft::jmccm1d3
subroutine jmccm1d3(p, n, m, table, ntable, itable, ntable2, mtable, work, nwork, ioff)
Definition GFT_jmfft.f90:523
gft_jmfft::jmgeterreur
subroutine jmgeterreur(arret)
Definition GFT_jmfft.f90:2208
gft_jmfft::jmgetnwork
subroutine jmgetnwork(nwork, nwork_def, nwork_min)
Definition GFT_jmfft.f90:2152
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::jmcctranspcs
subroutine jmcctranspcs(m, n, n2, n3, table, ntable, itable, work, nwork, ioff)
Definition GFT_jmfft.f90:2034
gft_jmfft::jmgetsetnwork
subroutine jmgetsetnwork(nwork, type)
Definition GFT_jmfft.f90:2184
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::jmgetstop
subroutine jmgetstop(arret)
Definition GFT_jmfft.f90:2322
gft_jmfft::jmtable
subroutine jmtable(table, ntable, itable, n)
Definition GFT_jmfft.f90:21
gft_jmfft::jmtransp
subroutine jmtransp(n, m, l, work, nwork, ioff)
Definition GFT_jmfft.f90:1986
gft_jmfft::jmsetcode
subroutine jmsetcode(code)
Definition GFT_jmfft.f90:2406
gft_jmfft::jmfact
subroutine jmfact(n, fact, nfact, ideb, ifin)
Definition GFT_jmfft.f90:107
gft_jmfft::jmerreur2
subroutine jmerreur2(nomsp, code, var1, var2)
Definition GFT_jmfft.f90:2247
gft_jmfft::jmsetnwork
subroutine jmsetnwork(nwork)
Definition GFT_jmfft.f90:2129
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_jmfft::jmerreur3
subroutine jmerreur3(nomsp, code, var1, var2, var3)
Definition GFT_jmfft.f90:2272
gft_jmfft::jmerreur1
subroutine jmerreur1(nomsp, code, var1)
Definition GFT_jmfft.f90:2221
gft_jmfft::jmgetsetstop
subroutine jmgetsetstop(arret, type)
Definition GFT_jmfft.f90:2335
gft_jmfft::jmccm1d2
subroutine jmccm1d2(p, n, m, table, ntable, itable, ntable2, mtable, work, nwork, ioff)
Definition GFT_jmfft.f90:429
gft_jmfft::jmerreur4
subroutine jmerreur4(nomsp, code, var1, var2, var3, var4)
Definition GFT_jmfft.f90:2297