30 USE iso_c_binding,
ONLY: c_f_pointer,&
53#include "../base/base_uses.f90"
57 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'fft_tools'
63 INTEGER :: nx = 0, ny = 0, nz = 0
64 INTEGER :: lmax = 0, mmax = 0, nmax = 0
65 INTEGER :: mx1 = 0, mx2 = 0, mx3 = 0
66 INTEGER :: my1 = 0, my2 = 0, my3 = 0
67 INTEGER :: mz1 = 0, mz2 = 0, mz3 = 0
68 INTEGER :: mcz1 = 0, mcz2 = 0, mcy3 = 0, mcx2 = 0
69 INTEGER :: lg = 0, mg = 0
70 INTEGER :: nbx = 0, nbz = 0
71 INTEGER :: nmray = 0, nyzray = 0
73 INTEGER,
DIMENSION(2) :: g_pos = 0, r_pos = 0, r_dim = 0
74 INTEGER :: numtask = 0
78 INTEGER :: fft_scratch_id = 0
79 INTEGER :: tf_type = -1
80 LOGICAL :: in_use = .true.
82 INTEGER,
DIMENSION(3) :: nfft = -1
85 INTEGER,
DIMENSION(2) :: dim = -1, pos = -1
87 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
POINTER,
CONTIGUOUS &
88 :: ziptr => null(), zoptr => null()
90 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER &
91 :: p1buf => null(), p2buf => null(), p3buf => null(), p4buf => null(), &
92 p5buf => null(), p6buf => null(), p7buf => null()
94 COMPLEX(KIND=dp),
DIMENSION(:, :),
POINTER,
CONTIGUOUS &
95 :: r1buf => null(), r2buf => null()
96 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
POINTER,
CONTIGUOUS &
99 COMPLEX(KIND=dp),
DIMENSION(:, :),
POINTER,
CONTIGUOUS &
100 :: a1buf => null(), a2buf => null(), a3buf => null(), &
101 a4buf => null(), a5buf => null(), a6buf => null()
103 INTEGER,
DIMENSION(:),
CONTIGUOUS,
POINTER :: scount => null(), rcount => null(), sdispl => null(), rdispl => null()
104 INTEGER,
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: pgcube => null()
105 INTEGER,
DIMENSION(:),
CONTIGUOUS,
POINTER :: xzcount => null(), yzcount => null(), xzdispl => null(), yzdispl => null()
106 INTEGER :: in = 0, mip = -1
107 REAL(kind=
dp) :: rsratio = 1.0_dp
108 COMPLEX(KIND=dp),
DIMENSION(:),
POINTER,
CONTIGUOUS &
109 :: xzbuf => null(), yzbuf => null()
110 COMPLEX(KIND=sp),
DIMENSION(:),
POINTER,
CONTIGUOUS &
111 :: xzbuf_sgl => null(), yzbuf_sgl => null()
112 COMPLEX(KIND=dp),
DIMENSION(:, :),
POINTER,
CONTIGUOUS &
113 :: rbuf1 => null(), rbuf2 => null(), rbuf3 => null(), rbuf4 => null(), &
114 rbuf5 => null(), rbuf6 => null(), rr => null()
115 COMPLEX(KIND=sp),
DIMENSION(:, :),
POINTER,
CONTIGUOUS &
116 :: ss => null(), tt => null()
117 INTEGER,
DIMENSION(:, :),
POINTER,
CONTIGUOUS :: pgrid => null()
118 INTEGER,
DIMENSION(:),
POINTER,
CONTIGUOUS :: xcor => null(), zcor => null(), pzcoord => null()
121 INTEGER :: last_tick = -1
139 PUBLIC :: fft_alloc, fft_dealloc
148 INTEGER,
PARAMETER :: fft_radix_allowed = 495, fft_radix_disallowed = 496
151 REAL(kind=
dp),
PARAMETER :: ratio_sparse_alltoall = 0.5_dp
154 LOGICAL,
SAVE :: alltoall_sgl = .false.
163 MODULE PROCEDURE fft3d_s, fft3d_ps, fft3d_pb
179 SUBROUTINE init_fft(fftlib, alltoall, pool_limit, wisdom_file, &
182 CHARACTER(LEN=*),
INTENT(IN) :: fftlib
183 LOGICAL,
INTENT(IN) :: alltoall
184 INTEGER,
INTENT(IN) :: pool_limit
185 CHARACTER(LEN=*),
INTENT(IN) :: wisdom_file
186 INTEGER,
INTENT(IN) :: plan_style
188 alltoall_sgl = alltoall
191 CALL release_fft_scratch_pool()
210 CHARACTER(LEN=*),
INTENT(IN) :: wisdom_file
211 LOGICAL,
INTENT(IN),
OPTIONAL :: export_wisdom
215 CALL release_fft_scratch_pool()
219 CALL fft_do_cleanup(wisdom_file, para_env%is_source(), export_wisdom)
234 INTEGER,
INTENT(IN) :: radix_in
235 INTEGER,
INTENT(OUT) :: radix_out
236 INTEGER,
INTENT(IN) :: operation
238 INTEGER,
PARAMETER :: number_of_small_radices = 81
239 INTEGER,
DIMENSION(number_of_small_radices),
PARAMETER :: small_radices = [2, 4, 6, 8, 9, 12,&
240 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 54, 60, 64, 72, 75, 80, 81, 90, 96, &
241 100, 108, 120, 125, 128, 135, 144, 150, 160, 162, 180, 192, 200, 216, 225, 240, 243, 256, &
242 270, 288, 300, 320, 324, 360, 375, 384, 400, 405, 432, 450, 480, 486, 500, 512, 540, 576, &
243 600, 625, 640, 648, 675, 720, 729, 750, 768, 800, 810, 864, 900, 960, 972, 1000, 1024]
245 INTEGER :: i, iloc, ldata, trial_number
246 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: data
249 ALLOCATE (
DATA(ldata))
252 DATA(1:number_of_small_radices) = small_radices(1:number_of_small_radices)
256 DO trial_number = 1, 2
260 IF (
DATA(i) == radix_in)
THEN
264 IF (operation == fft_radix_allowed)
THEN
266 ELSE IF (
DATA(i) > radix_in)
THEN
274 IF (trial_number == 1)
THEN
281 radix_out = maxval(data)
283 CALL cp_abort(__location__, &
284 "No suitable FFT lengths found! Either enable FFTW or reduce the grid size!")
291 IF (operation == fft_radix_allowed)
THEN
292 IF (
DATA(iloc) == radix_in)
THEN
293 radix_out = fft_radix_allowed
295 radix_out = fft_radix_disallowed
299 IF (
DATA(iloc) == radix_in)
THEN
300 radix_out =
DATA(iloc)
302 IF (abs(
DATA(iloc - 1) - radix_in) <= &
303 abs(
DATA(iloc) - radix_in))
THEN
304 radix_out =
DATA(iloc - 1)
306 radix_out =
DATA(iloc)
311 radix_out =
DATA(iloc)
315 IF (mod(
DATA(i), 2) == 1)
THEN
320 IF (mod(radix_out, 2) == 0)
THEN
322 IF (trial_number == 1)
THEN
325 radix_out = radix_in + 1 - mod(radix_in, 2)
351 SUBROUTINE fft_1d_many(sign, n, m, trans_in, trans_out, ldx_in, ldx_out, zin, zout, scale, stat)
352 INTEGER,
INTENT(in) :: sign, n, m
353 LOGICAL,
INTENT(in) :: trans_in, trans_out
354 INTEGER,
INTENT(in) :: ldx_in, ldx_out
355 COMPLEX(kind=dp),
DIMENSION(*),
INTENT(inout) :: zin, zout
356 REAL(kind=
dp),
INTENT(in) :: scale
357 INTEGER,
INTENT(out) :: stat
359 CHARACTER(len=*),
PARAMETER :: routinen =
'fft_1d_many'
361 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:), &
366 CALL timeset(routinen, handle)
369 cpabort(
"Unknown FFT-direction!")
374 CALL fft_alloc(ztmp, [ldx_in*n])
376 CALL fft_alloc(ztmp, [ldx_in*m])
379 CALL fft_dealloc(ztmp)
383 CALL timestop(handle)
398 SUBROUTINE fft3d_s(fsign, n, zin, zout, status, debug)
400 INTEGER,
INTENT(IN) :: fsign
401 INTEGER,
DIMENSION(:),
INTENT(INOUT) :: n
402 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :, :), &
404 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :, :), &
405 INTENT(INOUT),
OPTIONAL,
TARGET :: zout
406 INTEGER,
INTENT(OUT),
OPTIONAL :: status
407 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
409 CHARACTER(len=*),
PARAMETER :: routineN =
'fft3d_s'
411 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :, :), &
413 COMPLEX(KIND=dp),
DIMENSION(1, 1, 1),
TARGET :: zdum
414 INTEGER :: handle, ld(3), lo(3), output_unit, sign, &
416 LOGICAL :: fft_in_place, test
417 REAL(KIND=
dp) :: in_sum, norm, out_sum
420 CALL timeset(routinen, handle)
423 IF (fsign ==
fwfft)
THEN
424 norm = 1.0_dp/real(product(n), kind=
dp)
425 ELSE IF (fsign ==
bwfft)
THEN
428 cpabort(
"Unknown FFT direction!")
431 IF (
PRESENT(debug))
THEN
437 IF (
PRESENT(zout))
THEN
438 fft_in_place = .false.
440 fft_in_place = .true.
444 in_sum = sum(abs(zin))
451 IF (n(1) /= ld(1) .OR. n(2) /= ld(2) .OR. n(3) /= ld(3))
THEN
452 cpabort(
"Size and dimension (zin) have to be the same.")
458 IF (fft_in_place)
THEN
460 IF (fsign ==
fwfft)
THEN
461 CALL fft_3d(fft_scratch%fft_plan(1), norm, zin, zoptr, stat)
463 CALL fft_3d(fft_scratch%fft_plan(2), norm, zin, zoptr, stat)
466 IF (fsign ==
fwfft)
THEN
467 CALL fft_3d(fft_scratch%fft_plan(3), norm, zin, zout, stat)
469 CALL fft_3d(fft_scratch%fft_plan(4), norm, zin, zout, stat)
475 IF (
PRESENT(zout))
THEN
476 lo(1) =
SIZE(zout, 1)
477 lo(2) =
SIZE(zout, 2)
478 lo(3) =
SIZE(zout, 3)
479 IF (n(1) /= lo(1) .OR. n(2) /= lo(2) .OR. n(3) /= lo(3))
THEN
480 cpabort(
"Size and dimension (zout) have to be the same.")
484 IF (
PRESENT(status))
THEN
488 IF (test .AND. output_unit > 0)
THEN
489 IF (
PRESENT(zout))
THEN
490 out_sum = sum(abs(zout))
491 WRITE (output_unit,
'(A)')
" Out of place 3D FFT (local) : fft3d_s"
492 WRITE (output_unit,
'(A,T60,3I7)')
" Transform lengths ", n
493 WRITE (output_unit,
'(A,T60,3I7)')
" Input array dimensions ", ld
494 WRITE (output_unit,
'(A,T60,3I7)')
" Output array dimensions ", lo
495 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of input data ", in_sum
496 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of output data ", out_sum
498 out_sum = sum(abs(zin))
499 WRITE (output_unit,
'(A)')
" In place 3D FFT (local) : fft3d_s"
500 WRITE (output_unit,
'(A,T60,3I7)')
" Transform lengths ", n
501 WRITE (output_unit,
'(A,T60,3I7)')
" Input/output array dimensions ", ld
502 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of input data ", in_sum
503 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of output data ", out_sum
507 CALL timestop(handle)
509 END SUBROUTINE fft3d_s
524 SUBROUTINE fft3d_ps(fsign, n, cin, gin, rs_group, yzp, nyzray, &
527 INTEGER,
INTENT(IN) :: fsign
528 INTEGER,
DIMENSION(:),
INTENT(IN) :: n
529 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :, :), &
531 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
534 INTEGER,
CONTIGUOUS,
DIMENSION(:, :, 0:), &
536 INTEGER,
CONTIGUOUS,
DIMENSION(0:),
INTENT(IN) :: nyzray
537 INTEGER,
CONTIGUOUS,
DIMENSION(:, :, 0:, :), &
539 INTEGER,
INTENT(OUT),
OPTIONAL :: status
540 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
542 CHARACTER(len=*),
PARAMETER :: routineN =
'fft3d_ps'
544 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
545 POINTER :: pbuf, qbuf, rbuf, sbuf
546 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :, :), &
548 INTEGER :: g_pos, handle, lg, lmax, mcx2, mcz1, mcz2, mg, mmax, mx1, mx2, my1, mz2, n1, n2, &
549 nmax, numtask, nx, ny, nz, output_unit, r_dim(2), r_pos(2), rp, sign, stat
550 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: p2p
552 REAL(KIND=
dp) :: norm, sum_data
556 CALL timeset(routinen, handle)
559 IF (
PRESENT(debug))
THEN
565 g_pos = rs_group%mepos
566 numtask = rs_group%num_pe
567 r_dim = rs_group%num_pe_cart
568 r_pos = rs_group%mepos_cart
570 IF (fsign ==
fwfft)
THEN
571 norm = 1.0_dp/real(product(n), kind=
dp)
572 ELSE IF (fsign ==
bwfft)
THEN
575 cpabort(
"Unknown FFT direction!")
592 lmax = max(lg, (nx*ny*nz)/mmax + 1)
594 ALLOCATE (p2p(0:numtask - 1))
596 CALL rs_group%rank_compare(rs_group, p2p)
599 mx1 = bo(2, 1, rp, 1) - bo(1, 1, rp, 1) + 1
600 my1 = bo(2, 2, rp, 1) - bo(1, 2, rp, 1) + 1
601 mx2 = bo(2, 1, rp, 2) - bo(1, 1, rp, 2) + 1
602 mz2 = bo(2, 3, rp, 2) - bo(1, 3, rp, 2) + 1
604 n1 = maxval(bo(2, 1, :, 1) - bo(1, 1, :, 1) + 1)
605 n2 = maxval(bo(2, 2, :, 1) - bo(1, 2, :, 1) + 1)
606 nmax = max((2*n2)/numtask, 2)*mx2*mz2
607 nmax = max(nmax, n1*maxval(nyzray))
608 n1 = maxval(bo(2, 1, :, 2))
609 n2 = maxval(bo(2, 3, :, 2))
611 fft_scratch_size%nx = nx
612 fft_scratch_size%ny = ny
613 fft_scratch_size%nz = nz
614 fft_scratch_size%lmax = lmax
615 fft_scratch_size%mmax = mmax
616 fft_scratch_size%mx1 = mx1
617 fft_scratch_size%mx2 = mx2
618 fft_scratch_size%my1 = my1
619 fft_scratch_size%mz2 = mz2
620 fft_scratch_size%lg = lg
621 fft_scratch_size%mg = mg
622 fft_scratch_size%nbx = n1
623 fft_scratch_size%nbz = n2
624 mcz1 = maxval(bo(2, 3, :, 1) - bo(1, 3, :, 1) + 1)
625 mcx2 = maxval(bo(2, 1, :, 2) - bo(1, 1, :, 2) + 1)
626 mcz2 = maxval(bo(2, 3, :, 2) - bo(1, 3, :, 2) + 1)
627 fft_scratch_size%mcz1 = mcz1
628 fft_scratch_size%mcx2 = mcx2
629 fft_scratch_size%mcz2 = mcz2
630 fft_scratch_size%nmax = nmax
631 fft_scratch_size%nmray = maxval(nyzray)
632 fft_scratch_size%nyzray = nyzray(g_pos)
633 fft_scratch_size%rs_group = rs_group
634 fft_scratch_size%g_pos = g_pos
635 fft_scratch_size%r_pos = r_pos
636 fft_scratch_size%r_dim = r_dim
637 fft_scratch_size%numtask = numtask
640 IF (g_pos == 0 .AND. output_unit > 0)
THEN
641 WRITE (output_unit,
'(A)')
" Parallel 3D FFT : fft3d_ps"
642 WRITE (output_unit,
'(A,T60,3I7)')
" Transform lengths ", n
643 WRITE (output_unit,
'(A,T67,2I7)')
" Array dimensions (gin) ", lg, mg
644 WRITE (output_unit,
'(A,T60,3I7)')
" Array dimensions (cin) ", nx, ny, nz
648 IF (r_dim(2) > 1)
THEN
655 IF (r_dim(1) == 1)
THEN
656 cpabort(
"This processor distribution is not supported.")
658 CALL get_fft_scratch(fft_scratch, tf_type=300, n=n, fft_sizes=fft_scratch_size)
660 IF (sign ==
fwfft)
THEN
664 sum_data = abs(sum(cin))
665 CALL rs_group%sum(sum_data)
666 IF (g_pos == 0 .AND. output_unit > 0)
THEN
667 WRITE (output_unit,
'(A)')
" Two step communication algorithm "
668 WRITE (output_unit,
'(A,T60,3I7)')
" Transform Z ", n(3), mx1*my1
669 WRITE (output_unit,
'(A,T60,3I7)')
" Transform Y ", n(2), mx2*mz2
670 WRITE (output_unit,
'(A,T67,2I7)')
" Transform X ", n(1), nyzray(g_pos)
671 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(1) ", sum_data
675 pbuf => fft_scratch%p1buf
676 qbuf => fft_scratch%p2buf
679 CALL fft_1d(fft_scratch%fft_plan(1), cin, qbuf, norm, stat)
681 rbuf => fft_scratch%p3buf
684 sum_data = abs(sum(qbuf))
685 CALL rs_group%sum(sum_data)
686 IF (g_pos == 0 .AND. output_unit > 0)
THEN
687 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(2) T", sum_data
692 CALL cube_transpose_2(qbuf, bo(:, :, :, 1), bo(:, :, :, 2), rbuf, fft_scratch)
695 sum_data = abs(sum(rbuf))
696 CALL rs_group%sum(sum_data)
697 IF (g_pos == 0 .AND. output_unit > 0)
THEN
698 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(3) T", sum_data
702 pbuf => fft_scratch%p4buf
705 CALL fft_1d(fft_scratch%fft_plan(2), rbuf, pbuf, 1.0_dp, stat)
707 qbuf => fft_scratch%p5buf
710 sum_data = abs(sum(pbuf))
711 CALL rs_group%sum(sum_data)
712 IF (g_pos == 0 .AND. output_unit > 0)
THEN
713 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(4) TS", sum_data
718 CALL xz_to_yz(pbuf, rs_group, r_dim, g_pos, p2p, yzp, nyzray, &
719 bo(:, :, :, 2), qbuf, fft_scratch)
722 sum_data = abs(sum(qbuf))
723 CALL rs_group%sum(sum_data)
724 IF (g_pos == 0 .AND. output_unit > 0)
THEN
725 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(5) TS", sum_data
730 CALL fft_1d(fft_scratch%fft_plan(3), qbuf, gin, 1.0_dp, stat)
733 sum_data = abs(sum(gin))
734 CALL rs_group%sum(sum_data)
735 IF (g_pos == 0 .AND. output_unit > 0)
THEN
736 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(6) ", sum_data
740 ELSE IF (sign ==
bwfft)
THEN
744 sum_data = abs(sum(gin))
745 CALL rs_group%sum(sum_data)
746 IF (g_pos == 0 .AND. output_unit > 0)
THEN
747 WRITE (output_unit,
'(A)')
" Two step communication algorithm "
748 WRITE (output_unit,
'(A,T67,2I7)')
" Transform X ", n(1), nyzray(g_pos)
749 WRITE (output_unit,
'(A,T60,3I7)')
" Transform Y ", n(2), mx2*mz2
750 WRITE (output_unit,
'(A,T60,3I7)')
" Transform Z ", n(3), mx1*my1
751 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(1) ", sum_data
755 pbuf => fft_scratch%p7buf
758 CALL fft_1d(fft_scratch%fft_plan(4), gin, pbuf, norm, stat)
760 qbuf => fft_scratch%p4buf
763 sum_data = abs(sum(pbuf))
764 CALL rs_group%sum(sum_data)
765 IF (g_pos == 0 .AND. output_unit > 0)
THEN
766 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(2) TS", sum_data
771 CALL yz_to_xz(pbuf, rs_group, r_dim, g_pos, p2p, yzp, nyzray, &
772 bo(:, :, :, 2), qbuf, fft_scratch)
775 sum_data = abs(sum(qbuf))
776 CALL rs_group%sum(sum_data)
777 IF (g_pos == 0 .AND. output_unit > 0)
THEN
778 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(3) TS", sum_data
782 rbuf => fft_scratch%p3buf
785 CALL fft_1d(fft_scratch%fft_plan(5), qbuf, rbuf, 1.0_dp, stat)
787 pbuf => fft_scratch%p2buf
790 sum_data = abs(sum(rbuf))
791 CALL rs_group%sum(sum_data)
792 IF (g_pos == 0 .AND. output_unit > 0)
THEN
793 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(4) T", sum_data
798 CALL cube_transpose_1(rbuf, bo(:, :, :, 2), bo(:, :, :, 1), pbuf, fft_scratch)
801 sum_data = abs(sum(pbuf))
802 CALL rs_group%sum(sum_data)
803 IF (g_pos == 0 .AND. output_unit > 0)
THEN
804 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(5) T", sum_data
808 qbuf => fft_scratch%p1buf
811 CALL fft_1d(fft_scratch%fft_plan(6), pbuf, cin, 1.0_dp, stat)
814 sum_data = abs(sum(cin))
815 CALL rs_group%sum(sum_data)
816 IF (g_pos == 0 .AND. output_unit > 0)
THEN
817 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(6) ", sum_data
823 cpabort(
"Illegal fsign parameter.")
837 CALL get_fft_scratch(fft_scratch, tf_type=200, n=n, fft_sizes=fft_scratch_size)
839 sbuf => fft_scratch%r1buf
840 tbuf => fft_scratch%tbuf
845 IF (sign ==
fwfft)
THEN
849 sum_data = abs(sum(cin))
850 CALL rs_group%sum(sum_data)
851 IF (g_pos == 0 .AND. output_unit > 0)
THEN
852 WRITE (output_unit,
'(A)')
" One step communication algorithm "
853 WRITE (output_unit,
'(A,T60,3I7)')
" Transform YZ ", n(2), n(3), nx
854 WRITE (output_unit,
'(A,T67,2I7)')
" Transform X ", n(1), nyzray(g_pos)
855 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(1) ", sum_data
860 CALL fft_1d(fft_scratch%fft_plan(1), cin, sbuf, 1._dp, stat)
861 CALL fft_1d(fft_scratch%fft_plan(2), sbuf, tbuf, 1._dp, stat)
864 sum_data = abs(sum(tbuf))
865 CALL rs_group%sum(sum_data)
866 IF (g_pos == 0 .AND. output_unit > 0)
THEN
867 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(2) TS", sum_data
872 CALL yz_to_x(tbuf, rs_group, g_pos, p2p, yzp, nyzray, &
873 bo(:, :, :, 2), sbuf, fft_scratch)
876 sum_data = abs(sum(sbuf))
877 CALL rs_group%sum(sum_data)
878 IF (g_pos == 0 .AND. output_unit > 0)
THEN
879 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(3) TS", sum_data
883 CALL fft_1d(fft_scratch%fft_plan(3), sbuf, gin, norm, stat)
886 sum_data = abs(sum(gin))
887 CALL rs_group%sum(sum_data)
888 IF (g_pos == 0 .AND. output_unit > 0)
THEN
889 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(4) ", sum_data
893 ELSE IF (sign ==
bwfft)
THEN
897 sum_data = abs(sum(gin))
898 CALL rs_group%sum(sum_data)
899 IF (g_pos == 0 .AND. output_unit > 0)
THEN
900 WRITE (output_unit,
'(A)')
" One step communication algorithm "
901 WRITE (output_unit,
'(A,T67,2I7)')
" Transform X ", n(1), nyzray(g_pos)
902 WRITE (output_unit,
'(A,T60,3I7)')
" Transform YZ ", n(2), n(3), nx
903 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(1) ", sum_data
908 CALL fft_1d(fft_scratch%fft_plan(4), gin, sbuf, norm, stat)
911 sum_data = abs(sum(sbuf))
912 CALL rs_group%sum(sum_data)
913 IF (g_pos == 0 .AND. output_unit > 0)
THEN
914 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(2) TS", sum_data
919 CALL x_to_yz(sbuf, rs_group, g_pos, p2p, yzp, nyzray, &
920 bo(:, :, :, 2), tbuf, fft_scratch)
923 sum_data = abs(sum(tbuf))
924 CALL rs_group%sum(sum_data)
925 IF (g_pos == 0 .AND. output_unit > 0)
THEN
926 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(3) TS", sum_data
931 CALL fft_1d(fft_scratch%fft_plan(5), tbuf, sbuf, 1._dp, stat)
932 CALL fft_1d(fft_scratch%fft_plan(6), sbuf, cin, 1._dp, stat)
935 sum_data = abs(sum(cin))
936 CALL rs_group%sum(sum_data)
937 IF (g_pos == 0 .AND. output_unit > 0)
THEN
938 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(4) ", sum_data
942 cpabort(
"Illegal fsign parameter.")
951 IF (
PRESENT(status))
THEN
954 CALL timestop(handle)
956 END SUBROUTINE fft3d_ps
969 SUBROUTINE fft3d_pb(fsign, n, zin, gin, group, bo, status, debug)
971 INTEGER,
INTENT(IN) :: fsign
972 INTEGER,
DIMENSION(3),
INTENT(IN) :: n
973 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :, :), &
975 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
978 INTEGER,
CONTIGUOUS,
DIMENSION(:, :, 0:, :), &
980 INTEGER,
INTENT(OUT),
OPTIONAL :: status
981 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
983 CHARACTER(len=*),
PARAMETER :: routineN =
'fft3d_pb'
985 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
986 POINTER :: abuf, bbuf
987 INTEGER :: handle, lg(2), lz(3), mcx2, mcy3, mcz1, &
988 mcz2, mx1, mx2, mx3, my1, my2, my3, &
989 my_pos, mz1, mz2, mz3, output_unit, &
991 INTEGER,
DIMENSION(2) :: dim
993 REAL(KIND=
dp) :: norm, sum_data
1017 CALL timeset(routinen, handle)
1020 dim = group%num_pe_cart
1021 my_pos = group%mepos
1023 IF (
PRESENT(debug))
THEN
1029 IF (fsign ==
fwfft)
THEN
1030 norm = 1.0_dp/real(product(n), kind=
dp)
1031 ELSE IF (fsign ==
bwfft)
THEN
1034 cpabort(
"Unknown FFT direction!")
1040 lg(1) =
SIZE(gin, 1)
1041 lg(2) =
SIZE(gin, 2)
1042 lz(1) =
SIZE(zin, 1)
1043 lz(2) =
SIZE(zin, 2)
1044 lz(3) =
SIZE(zin, 3)
1045 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1046 WRITE (output_unit,
'(A)')
" Parallel 3D FFT : fft3d_pb"
1047 WRITE (output_unit,
'(A,T60,3I7)')
" Transform lengths ", n
1048 WRITE (output_unit,
'(A,T67,2I7)')
" Array dimensions (gin) ", lg
1049 WRITE (output_unit,
'(A,T60,3I7)')
" Array dimensions (cin) ", lz
1053 mx1 = bo(2, 1, my_pos, 1) - bo(1, 1, my_pos, 1) + 1
1054 my1 = bo(2, 2, my_pos, 1) - bo(1, 2, my_pos, 1) + 1
1055 mz1 = bo(2, 3, my_pos, 1) - bo(1, 3, my_pos, 1) + 1
1056 mx2 = bo(2, 1, my_pos, 2) - bo(1, 1, my_pos, 2) + 1
1057 my2 = bo(2, 2, my_pos, 2) - bo(1, 2, my_pos, 2) + 1
1058 mz2 = bo(2, 3, my_pos, 2) - bo(1, 3, my_pos, 2) + 1
1059 mx3 = bo(2, 1, my_pos, 3) - bo(1, 1, my_pos, 3) + 1
1060 my3 = bo(2, 2, my_pos, 3) - bo(1, 2, my_pos, 3) + 1
1061 mz3 = bo(2, 3, my_pos, 3) - bo(1, 3, my_pos, 3) + 1
1062 fft_scratch_size%mx1 = mx1
1063 fft_scratch_size%mx2 = mx2
1064 fft_scratch_size%mx3 = mx3
1065 fft_scratch_size%my1 = my1
1066 fft_scratch_size%my2 = my2
1067 fft_scratch_size%my3 = my3
1068 fft_scratch_size%mz1 = mz1
1069 fft_scratch_size%mz2 = mz2
1070 fft_scratch_size%mz3 = mz3
1071 mcz1 = maxval(bo(2, 3, :, 1) - bo(1, 3, :, 1) + 1)
1072 mcx2 = maxval(bo(2, 1, :, 2) - bo(1, 1, :, 2) + 1)
1073 mcz2 = maxval(bo(2, 3, :, 2) - bo(1, 3, :, 2) + 1)
1074 mcy3 = maxval(bo(2, 2, :, 3) - bo(1, 2, :, 3) + 1)
1075 fft_scratch_size%mcz1 = mcz1
1076 fft_scratch_size%mcx2 = mcx2
1077 fft_scratch_size%mcz2 = mcz2
1078 fft_scratch_size%mcy3 = mcy3
1079 fft_scratch_size%rs_group = group
1080 fft_scratch_size%g_pos = my_pos
1081 fft_scratch_size%numtask = dim(1)*dim(2)
1083 IF (dim(1) > 1 .AND. dim(2) > 1)
THEN
1089 CALL get_fft_scratch(fft_scratch, tf_type=100, n=n, fft_sizes=fft_scratch_size)
1091 IF (sign ==
fwfft)
THEN
1094 bbuf => fft_scratch%a2buf
1097 sum_data = abs(sum(zin))
1098 CALL group%sum(sum_data)
1099 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1100 WRITE (output_unit,
'(A)')
" Two step communication algorithm "
1101 WRITE (output_unit,
'(A,T67,2I7)')
" Transform Z ", n(3), mx1*my1
1102 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(1) ", sum_data
1107 CALL fft_1d(fft_scratch%fft_plan(1), zin, bbuf, norm, stat)
1109 abuf => fft_scratch%a3buf
1112 sum_data = abs(sum(bbuf))
1113 CALL group%sum(sum_data)
1114 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1115 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(2) T", sum_data
1119 CALL cube_transpose_2(bbuf, bo(:, :, :, 1), bo(:, :, :, 2), abuf, fft_scratch)
1121 bbuf => fft_scratch%a4buf
1124 sum_data = abs(sum(abuf))
1125 CALL group%sum(sum_data)
1126 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1127 WRITE (output_unit,
'(A,T67,2I7)')
" Transform Y ", n(2), mx2*mz2
1128 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(3) ", sum_data
1133 CALL fft_1d(fft_scratch%fft_plan(2), abuf, bbuf, 1.0_dp, stat)
1135 abuf => fft_scratch%a5buf
1138 sum_data = abs(sum(bbuf))
1139 CALL group%sum(sum_data)
1140 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1141 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(4) T", sum_data
1145 CALL cube_transpose_4(bbuf, bo(:, :, :, 2), bo(:, :, :, 3), abuf, fft_scratch)
1148 sum_data = abs(sum(abuf))
1149 CALL group%sum(sum_data)
1150 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1151 WRITE (output_unit,
'(A,T67,2I7)')
" Transform X ", n(1), my3*mz3
1152 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(5) ", sum_data
1157 CALL fft_1d(fft_scratch%fft_plan(3), abuf, gin, 1.0_dp, stat)
1160 sum_data = abs(sum(gin))
1161 CALL group%sum(sum_data)
1162 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1163 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(6) ", sum_data
1167 ELSE IF (sign ==
bwfft)
THEN
1170 bbuf => fft_scratch%a5buf
1173 sum_data = abs(sum(gin))
1174 CALL group%sum(sum_data)
1175 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1176 WRITE (output_unit,
'(A)')
" Two step communication algorithm "
1177 WRITE (output_unit,
'(A,T67,2I7)')
" Transform X ", n(1), my3*mz3
1178 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(1) ", sum_data
1183 CALL fft_1d(fft_scratch%fft_plan(4), gin, bbuf, 1.0_dp, stat)
1185 abuf => fft_scratch%a4buf
1188 sum_data = abs(sum(bbuf))
1189 CALL group%sum(sum_data)
1190 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1191 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(2) T", sum_data
1195 CALL cube_transpose_3(bbuf, bo(:, :, :, 3), bo(:, :, :, 2), abuf, fft_scratch)
1197 bbuf => fft_scratch%a3buf
1200 sum_data = abs(sum(abuf))
1201 CALL group%sum(sum_data)
1202 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1203 WRITE (output_unit,
'(A,T67,2I7)')
" Transform Y ", n(2), mx2*mz2
1204 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(3) ", sum_data
1209 CALL fft_1d(fft_scratch%fft_plan(5), abuf, bbuf, 1.0_dp, stat)
1211 abuf => fft_scratch%a2buf
1214 sum_data = abs(sum(bbuf))
1215 CALL group%sum(sum_data)
1216 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1217 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(4) T", sum_data
1221 CALL cube_transpose_1(bbuf, bo(:, :, :, 2), bo(:, :, :, 1), abuf, fft_scratch)
1224 sum_data = abs(sum(abuf))
1225 CALL group%sum(sum_data)
1226 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1227 WRITE (output_unit,
'(A,T67,2I7)')
" Transform Z ", n(3), mx1*my1
1228 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(5) ", sum_data
1233 CALL fft_1d(fft_scratch%fft_plan(6), abuf, zin, norm, stat)
1236 sum_data = abs(sum(zin))
1237 CALL group%sum(sum_data)
1238 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1239 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(6) ", sum_data
1244 cpabort(
"Illegal fsign parameter.")
1249 ELSE IF (dim(2) == 1)
THEN
1255 CALL get_fft_scratch(fft_scratch, tf_type=101, n=n, fft_sizes=fft_scratch_size)
1257 IF (sign ==
fwfft)
THEN
1261 sum_data = abs(sum(zin))
1262 CALL group%sum(sum_data)
1263 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1264 WRITE (output_unit,
'(A)')
" one step communication algorithm "
1265 WRITE (output_unit,
'(A,T67,2I7)')
" Transform Z ", n(3), mx1*my1
1266 WRITE (output_unit,
'(A,T67,2I7)')
" Transform Y ", n(2), mx1*mz1
1267 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(1) ", sum_data
1271 abuf => fft_scratch%a3buf
1272 bbuf => fft_scratch%a4buf
1274 CALL fft_1d(fft_scratch%fft_plan(1), zin, abuf, norm, stat)
1275 CALL fft_1d(fft_scratch%fft_plan(2), abuf, bbuf, 1.0_dp, stat)
1277 abuf => fft_scratch%a5buf
1280 sum_data = abs(sum(bbuf))
1281 CALL group%sum(sum_data)
1282 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1283 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(2) T", sum_data
1287 CALL cube_transpose_6(bbuf, group, bo(:, :, :, 1), bo(:, :, :, 3), abuf, fft_scratch)
1290 sum_data = abs(sum(abuf))
1291 CALL group%sum(sum_data)
1292 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1293 WRITE (output_unit,
'(A,T67,2I7)')
" Transform X ", n(1), my3*mz3
1294 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(3) ", sum_data
1299 CALL fft_1d(fft_scratch%fft_plan(3), abuf, gin, 1.0_dp, stat)
1302 sum_data = abs(sum(gin))
1303 CALL group%sum(sum_data)
1304 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1305 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(4) ", sum_data
1309 ELSE IF (sign ==
bwfft)
THEN
1313 sum_data = abs(sum(gin))
1314 CALL group%sum(sum_data)
1315 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1316 WRITE (output_unit,
'(A)')
" one step communication algorithm "
1317 WRITE (output_unit,
'(A,T67,2I7)')
" Transform X ", n(1), my3*mz3
1318 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(1) ", sum_data
1322 bbuf => fft_scratch%a5buf
1325 CALL fft_1d(fft_scratch%fft_plan(4), gin, bbuf, 1.0_dp, stat)
1327 abuf => fft_scratch%a4buf
1330 sum_data = abs(sum(bbuf))
1331 CALL group%sum(sum_data)
1332 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1333 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(2) T", sum_data
1337 CALL cube_transpose_5(bbuf, group, bo(:, :, :, 3), bo(:, :, :, 1), abuf, fft_scratch)
1339 bbuf => fft_scratch%a3buf
1342 sum_data = abs(sum(abuf))
1343 CALL group%sum(sum_data)
1344 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1345 WRITE (output_unit,
'(A,T67,2I7)')
" Transform Y ", n(2), mx1*mz1
1346 WRITE (output_unit,
'(A,T67,2I7)')
" Transform Z ", n(3), mx1*my1
1347 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(3) ", sum_data
1352 CALL fft_1d(fft_scratch%fft_plan(5), abuf, bbuf, 1.0_dp, stat)
1355 CALL fft_1d(fft_scratch%fft_plan(6), bbuf, zin, norm, stat)
1358 sum_data = abs(sum(zin))
1359 CALL group%sum(sum_data)
1360 IF (my_pos == 0 .AND. output_unit > 0)
THEN
1361 WRITE (output_unit,
'(A,T61,E20.14)')
" Sum of data(4) ", sum_data
1366 cpabort(
"Illegal fsign parameter.")
1373 cpabort(
"Partition not implemented.")
1377 IF (
PRESENT(status))
THEN
1381 CALL timestop(handle)
1383 END SUBROUTINE fft3d_pb
1400 SUBROUTINE x_to_yz(sb, group, my_pos, p2p, yzp, nray, bo, tb, fft_scratch)
1402 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
1406 INTEGER,
INTENT(IN) :: my_pos
1407 INTEGER,
CONTIGUOUS,
DIMENSION(0:),
INTENT(IN) :: p2p
1408 INTEGER,
CONTIGUOUS,
DIMENSION(:, :, 0:), &
1410 INTEGER,
CONTIGUOUS,
DIMENSION(0:),
INTENT(IN) :: nray
1411 INTEGER,
CONTIGUOUS,
DIMENSION(:, :, 0:), &
1413 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :, :), &
1417 CHARACTER(len=*),
PARAMETER :: routinen =
'x_to_yz'
1419 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
1421 COMPLEX(KIND=sp),
CONTIGUOUS,
DIMENSION(:, :), &
1423 INTEGER :: handle, ip, ir, ix, ixx, iy, iz, mpr, &
1425 INTEGER,
CONTIGUOUS,
DIMENSION(:),
POINTER :: rcount, rdispl, scount, sdispl
1427 CALL timeset(routinen, handle)
1430 scount => fft_scratch%scount
1431 rcount => fft_scratch%rcount
1432 sdispl => fft_scratch%sdispl
1433 rdispl => fft_scratch%rdispl
1435 IF (alltoall_sgl)
THEN
1436 ss => fft_scratch%ss
1437 tt => fft_scratch%tt
1438 ss(:, :) = cmplx(sb(:, :), kind=
sp)
1441 rr => fft_scratch%rr
1445 nm = maxval(nray(0:np - 1))
1452 nx = bo(2, 1, ix) - bo(1, 1, ix) + 1
1454 sdispl(ip) = nr*(bo(1, 1, ix) - 1)
1457 nx = bo(2, 1, mpr) - bo(1, 1, mpr) + 1
1464 rdispl(ip) = nm*nx*ip
1467 IF (alltoall_sgl)
THEN
1468 CALL group%alltoall(ss, scount, sdispl, tt, rcount, rdispl)
1470 CALL group%alltoall(sb, scount, sdispl, rr, rcount, rdispl)
1473 nx = bo(2, 1, mpr) - bo(1, 1, mpr) + 1
1479 ixx = nray(ip)*(ix - 1)
1480 IF (alltoall_sgl)
THEN
1484 tb(iy, iz, ix) = tt(ir + ixx, ip)
1490 tb(iy, iz, ix) = rr(ir + ixx, ip)
1497 CALL timestop(handle)
1516 SUBROUTINE yz_to_x(tb, group, my_pos, p2p, yzp, nray, bo, sb, fft_scratch)
1518 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :, :), &
1522 INTEGER,
INTENT(IN) :: my_pos
1523 INTEGER,
CONTIGUOUS,
DIMENSION(0:),
INTENT(IN) :: p2p
1524 INTEGER,
CONTIGUOUS,
DIMENSION(:, :, 0:), &
1526 INTEGER,
CONTIGUOUS,
DIMENSION(0:),
INTENT(IN) :: nray
1527 INTEGER,
DIMENSION(:, :, 0:),
INTENT(IN) :: bo
1528 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
1532 CHARACTER(len=*),
PARAMETER :: routinen =
'yz_to_x'
1534 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
1536 COMPLEX(KIND=sp),
CONTIGUOUS,
DIMENSION(:, :), &
1538 INTEGER :: handle, ip, ir, ix, ixx, iy, iz, mpr, &
1540 INTEGER,
CONTIGUOUS,
DIMENSION(:),
POINTER :: rcount, rdispl, scount, sdispl
1542 CALL timeset(routinen, handle)
1546 scount => fft_scratch%scount
1547 rcount => fft_scratch%rcount
1548 sdispl => fft_scratch%sdispl
1549 rdispl => fft_scratch%rdispl
1551 IF (alltoall_sgl)
THEN
1552 ss => fft_scratch%ss
1553 tt => fft_scratch%tt
1557 rr => fft_scratch%rr
1560 nx = bo(2, 1, mpr) - bo(1, 1, mpr) + 1
1566 ixx = nray(ip)*(ix - 1)
1567 IF (alltoall_sgl)
THEN
1571 tt(ir + ixx, ip) = cmplx(tb(iy, iz, ix), kind=
sp)
1577 rr(ir + ixx, ip) = tb(iy, iz, ix)
1583 nm = maxval(nray(0:np - 1))
1590 nx = bo(2, 1, ix) - bo(1, 1, ix) + 1
1592 rdispl(ip) = nr*(bo(1, 1, ix) - 1)
1595 nx = bo(2, 1, mpr) - bo(1, 1, mpr) + 1
1602 sdispl(ip) = nm*nx*ip
1606 IF (alltoall_sgl)
THEN
1607 CALL group%alltoall(tt, scount, sdispl, ss, rcount, rdispl)
1610 CALL group%alltoall(rr, scount, sdispl, sb, rcount, rdispl)
1613 CALL timestop(handle)
1633 SUBROUTINE yz_to_xz(sb, group, dims, my_pos, p2p, yzp, nray, bo, tb, fft_scratch)
1635 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
1639 INTEGER,
DIMENSION(2),
INTENT(IN) :: dims
1640 INTEGER,
INTENT(IN) :: my_pos
1641 INTEGER,
CONTIGUOUS,
DIMENSION(0:),
INTENT(IN) :: p2p
1642 INTEGER,
CONTIGUOUS,
DIMENSION(:, :, 0:),
INTENT(IN) :: yzp
1643 INTEGER,
CONTIGUOUS,
DIMENSION(0:),
INTENT(IN) :: nray
1644 INTEGER,
CONTIGUOUS,
DIMENSION(:, :, 0:),
INTENT(IN) :: bo
1645 COMPLEX(KIND=dp),
DIMENSION(:, :),
INTENT(INOUT),
CONTIGUOUS :: tb
1648 CHARACTER(len=*),
PARAMETER :: routinen =
'yz_to_xz'
1650 COMPLEX(KIND=dp),
DIMENSION(:),
POINTER,
CONTIGUOUS :: xzbuf, yzbuf
1651 COMPLEX(KIND=sp),
DIMENSION(:),
POINTER,
CONTIGUOUS :: xzbuf_sgl, yzbuf_sgl
1652 INTEGER :: handle, icrs, ip, ipl, ipr, ir, ix, iz, &
1653 jj, jx, jy, jz, myx, myz, np, npx, &
1655 INTEGER,
DIMENSION(:),
POINTER,
CONTIGUOUS :: pzcoord, rcount, rdispl, scount, sdispl, &
1657 INTEGER,
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: pgrid
1659 CALL timeset(routinen, handle)
1663 rs_pos = p2p(my_pos)
1665 IF (alltoall_sgl)
THEN
1666 yzbuf_sgl => fft_scratch%yzbuf_sgl
1667 xzbuf_sgl => fft_scratch%xzbuf_sgl
1669 yzbuf => fft_scratch%yzbuf
1670 xzbuf => fft_scratch%xzbuf
1674 pgrid => fft_scratch%pgrid
1675 xcor => fft_scratch%xcor
1676 zcor => fft_scratch%zcor
1677 pzcoord => fft_scratch%pzcoord
1678 scount => fft_scratch%scount
1679 rcount => fft_scratch%rcount
1680 sdispl => fft_scratch%sdispl
1681 rdispl => fft_scratch%rdispl
1687 IF (fft_scratch%in == 0)
THEN
1693 xcor(bo(1, 1, ip):bo(2, 1, ip)) = ix
1697 zcor(bo(1, 3, ip):bo(2, 3, ip)) = iz
1700 IF (alltoall_sgl)
THEN
1701 DO ir = 1, nray(my_pos)
1702 jy = yzp(1, ir, my_pos)
1703 jz = yzp(2, ir, my_pos)
1704 ip = pgrid(xcor(jx), zcor(jz))
1705 scount(ip) = scount(ip) + 1
1708 DO ir = 1, nray(my_pos)
1709 jy = yzp(1, ir, my_pos)
1710 jz = yzp(2, ir, my_pos)
1711 ip = pgrid(xcor(jx), zcor(jz))
1712 scount(ip) = scount(ip) + 1
1717 CALL group%alltoall(scount, rcount, 1)
1718 fft_scratch%yzcount = scount
1719 fft_scratch%xzcount = rcount
1725 sdispl(ip) = sdispl(ip - 1) + scount(ip - 1)
1726 rdispl(ip) = rdispl(ip - 1) + rcount(ip - 1)
1729 fft_scratch%yzdispl = sdispl
1730 fft_scratch%xzdispl = rdispl
1734 IF (scount(ip) /= 0) icrs = icrs + 1
1735 IF (rcount(ip) /= 0) icrs = icrs + 1
1737 CALL group%sum(icrs)
1738 fft_scratch%rsratio = real(icrs, kind=
dp)/(real(2*np, kind=
dp)*real(np, kind=
dp))
1742 scount = fft_scratch%yzcount
1743 rcount = fft_scratch%xzcount
1744 sdispl = fft_scratch%yzdispl
1745 rdispl = fft_scratch%xzdispl
1755 IF (scount(ip) == 0) cycle
1758 nx = bo(2, 1, ipl) - bo(1, 1, ipl) + 1
1759 DO ir = 1, nray(my_pos)
1760 jz = yzp(2, ir, my_pos)
1761 IF (zcor(jz) == pzcoord(ipl))
THEN
1763 jy = yzp(1, ir, my_pos)
1764 IF (alltoall_sgl)
THEN
1766 yzbuf_sgl(sdispl(ip) + jj + jx*scount(ip)/nx) = cmplx(sb(ir, jx + bo(1, 1, ipl)), kind=
sp)
1770 yzbuf(sdispl(ip) + jj + jx*scount(ip)/nx) = sb(ir, jx + bo(1, 1, ipl))
1778 IF (alltoall_sgl)
THEN
1779 CALL group%alltoall(yzbuf_sgl, scount, sdispl, xzbuf_sgl, rcount, rdispl)
1781 IF (fft_scratch%rsratio < ratio_sparse_alltoall)
THEN
1782 CALL sparse_alltoall(yzbuf, scount, sdispl, xzbuf, rcount, rdispl, group)
1784 CALL group%alltoall(yzbuf, scount, sdispl, xzbuf, rcount, rdispl)
1788 myx = fft_scratch%sizes%r_pos(1)
1789 myz = fft_scratch%sizes%r_pos(2)
1790 nz = bo(2, 3, rs_pos) - bo(1, 3, rs_pos) + 1
1800 DO jx = 0, bo(2, 1, rs_pos) - bo(1, 1, rs_pos)
1803 IF (alltoall_sgl)
THEN
1804 IF (zcor(jz) == myz)
THEN
1807 jz = jz - bo(1, 3, rs_pos) + 1
1808 tb(jy, jz + jx*nz) = xzbuf_sgl(jj + rdispl(ipr))
1811 IF (zcor(jz) == myz)
THEN
1814 jz = jz - bo(1, 3, rs_pos) + 1
1815 tb(jy, jz + jx*nz) = xzbuf(jj + rdispl(ipr))
1823 CALL timestop(handle)
1843 SUBROUTINE xz_to_yz(sb, group, dims, my_pos, p2p, yzp, nray, bo, tb, fft_scratch)
1845 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
1849 INTEGER,
DIMENSION(2),
INTENT(IN) :: dims
1850 INTEGER,
INTENT(IN) :: my_pos
1851 INTEGER,
CONTIGUOUS,
DIMENSION(0:),
INTENT(IN) :: p2p
1852 INTEGER,
CONTIGUOUS,
DIMENSION(:, :, 0:),
INTENT(IN) :: yzp
1853 INTEGER,
CONTIGUOUS,
DIMENSION(0:),
INTENT(IN) :: nray
1854 INTEGER,
CONTIGUOUS,
DIMENSION(:, :, 0:),
INTENT(IN) :: bo
1855 COMPLEX(KIND=dp),
DIMENSION(:, :),
INTENT(INOUT),
CONTIGUOUS :: tb
1858 CHARACTER(len=*),
PARAMETER :: routinen =
'xz_to_yz'
1860 COMPLEX(KIND=dp),
DIMENSION(:),
POINTER,
CONTIGUOUS :: xzbuf, yzbuf
1861 COMPLEX(KIND=sp),
DIMENSION(:),
POINTER,
CONTIGUOUS :: xzbuf_sgl, yzbuf_sgl
1862 INTEGER :: handle, icrs, ip, ipl, ir, ix, ixx, iz, &
1863 jj, jx, jy, jz, mp, myx, myz, np, npx, &
1865 INTEGER,
DIMENSION(:),
POINTER,
CONTIGUOUS :: pzcoord, rcount, rdispl, scount, sdispl, &
1867 INTEGER,
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: pgrid
1869 CALL timeset(routinen, handle)
1873 IF (alltoall_sgl)
THEN
1874 yzbuf_sgl => fft_scratch%yzbuf_sgl
1875 xzbuf_sgl => fft_scratch%xzbuf_sgl
1877 yzbuf => fft_scratch%yzbuf
1878 xzbuf => fft_scratch%xzbuf
1882 pgrid => fft_scratch%pgrid
1883 xcor => fft_scratch%xcor
1884 zcor => fft_scratch%zcor
1885 pzcoord => fft_scratch%pzcoord
1886 scount => fft_scratch%scount
1887 rcount => fft_scratch%rcount
1888 sdispl => fft_scratch%sdispl
1889 rdispl => fft_scratch%rdispl
1893 IF (fft_scratch%in == 0)
THEN
1896 nx = maxval(bo(2, 1, :))
1900 xcor(bo(1, 1, ip):bo(2, 1, ip)) = ix
1904 zcor(bo(1, 3, ip):bo(2, 3, ip)) = iz
1907 DO ir = 1, nray(my_pos)
1908 jy = yzp(1, ir, my_pos)
1909 jz = yzp(2, ir, my_pos)
1910 ip = pgrid(xcor(jx), zcor(jz))
1911 rcount(ip) = rcount(ip) + 1
1915 CALL group%alltoall(rcount, scount, 1)
1916 fft_scratch%xzcount = scount
1917 fft_scratch%yzcount = rcount
1923 sdispl(ip) = sdispl(ip - 1) + scount(ip - 1)
1924 rdispl(ip) = rdispl(ip - 1) + rcount(ip - 1)
1927 fft_scratch%xzdispl = sdispl
1928 fft_scratch%yzdispl = rdispl
1932 IF (scount(ip) /= 0) icrs = icrs + 1
1933 IF (rcount(ip) /= 0) icrs = icrs + 1
1935 CALL group%sum(icrs)
1936 fft_scratch%rsratio = real(icrs, kind=
dp)/(real(2*np, kind=
dp)*real(np, kind=
dp))
1940 scount = fft_scratch%xzcount
1941 rcount = fft_scratch%yzcount
1942 sdispl = fft_scratch%xzdispl
1943 rdispl = fft_scratch%yzdispl
1947 myx = fft_scratch%sizes%r_pos(1)
1948 myz = fft_scratch%sizes%r_pos(2)
1950 nz = bo(2, 3, mp) - bo(1, 3, mp) + 1
1951 nx = bo(2, 1, mp) - bo(1, 1, mp) + 1
1963 IF (zcor(jz) == myz)
THEN
1966 jz = yzp(2, ir, ip) - bo(1, 3, mp) + 1
1967 IF (alltoall_sgl)
THEN
1969 ixx = jj + jx*scount(ipl)/nx
1970 xzbuf_sgl(ixx + sdispl(ipl)) = cmplx(sb(jy, jz + jx*nz), kind=
sp)
1974 ixx = jj + jx*scount(ipl)/nx
1975 xzbuf(ixx + sdispl(ipl)) = sb(jy, jz + jx*nz)
1983 IF (alltoall_sgl)
THEN
1984 CALL group%alltoall(xzbuf_sgl, scount, sdispl, yzbuf_sgl, rcount, rdispl)
1986 IF (fft_scratch%rsratio < ratio_sparse_alltoall)
THEN
1987 CALL sparse_alltoall(xzbuf, scount, sdispl, yzbuf, rcount, rdispl, group)
1989 CALL group%alltoall(xzbuf, scount, sdispl, yzbuf, rcount, rdispl)
1999 IF (rcount(ip) == 0) cycle
2002 nx = bo(2, 1, ipl) - bo(1, 1, ipl) + 1
2003 DO ir = 1, nray(my_pos)
2004 jz = yzp(2, ir, my_pos)
2005 IF (zcor(jz) == pzcoord(ipl))
THEN
2007 jy = yzp(1, ir, my_pos)
2008 IF (alltoall_sgl)
THEN
2010 tb(ir, jx + bo(1, 1, ipl)) = yzbuf_sgl(rdispl(ip) + jj + jx*rcount(ip)/nx)
2014 tb(ir, jx + bo(1, 1, ipl)) = yzbuf(rdispl(ip) + jj + jx*rcount(ip)/nx)
2022 CALL timestop(handle)
2039 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
2041 INTEGER,
CONTIGUOUS,
DIMENSION(:, :, 0:), &
2042 INTENT(IN) :: boin, boout
2043 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
2047 CHARACTER(len=*),
PARAMETER :: routinen =
'cube_transpose_1'
2049 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
2051 INTEGER :: handle, ip, ipl, ir, is, ixy, iz, mip, &
2053 INTEGER,
CONTIGUOUS,
DIMENSION(:),
POINTER :: rcount, rdispl, scount, sdispl
2054 INTEGER,
CONTIGUOUS,
DIMENSION(:, :),
POINTER :: pgrid
2055 INTEGER,
DIMENSION(2) :: dim, pos
2057 CALL timeset(routinen, handle)
2059 mip = fft_scratch%mip
2060 dim = fft_scratch%dim
2061 pos = fft_scratch%pos
2062 scount => fft_scratch%scount
2063 rcount => fft_scratch%rcount
2064 sdispl => fft_scratch%sdispl
2065 rdispl => fft_scratch%rdispl
2066 pgrid => fft_scratch%pgcube
2069 nx = boin(2, 1, mip) - boin(1, 1, mip) + 1
2070 nz = boin(2, 3, mip) - boin(1, 3, mip) + 1
2077 ny = boout(2, 2, ipl) - boout(1, 2, ipl) + 1
2078 scount(ip) = nx*nz*ny
2079 sdispl(ip) = nx*nz*(boout(1, 2, ipl) - 1)
2082 ny = boout(2, 2, mip) - boout(1, 2, mip) + 1
2083 mz = maxval(boin(2, 3, :) - boin(1, 3, :) + 1)
2089 nz = boin(2, 3, ipl) - boin(1, 3, ipl) + 1
2090 rcount(ip) = nx*nz*ny
2091 rdispl(ip) = nx*ny*mz*ip
2095 rbuf => fft_scratch%rbuf1
2097 CALL fft_scratch%cart_sub_comm(2)%alltoall(cin, scount, sdispl, rbuf, rcount, rdispl)
2105 nz = boin(2, 3, ipl) - boin(1, 3, ipl) + 1
2107 is = boin(1, 3, ipl) + iz - 1
2108 ir = iz + nz*(ixy - 1)
2109 sout(is, ixy) = rbuf(ir, ip)
2115 CALL timestop(handle)
2129 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
2131 INTEGER,
CONTIGUOUS,
DIMENSION(:, :, 0:), &
2132 INTENT(IN) :: boin, boout
2133 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
2137 CHARACTER(len=*),
PARAMETER :: routinen =
'cube_transpose_2'
2139 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
2141 INTEGER :: handle, ip, ipl, ir, ixy, iz, mip, mz, &
2143 INTEGER,
CONTIGUOUS,
DIMENSION(:),
POINTER :: rcount, rdispl, scount, sdispl
2144 INTEGER,
CONTIGUOUS,
DIMENSION(:, :),
POINTER :: pgrid
2145 INTEGER,
DIMENSION(2) :: dim, pos
2148 CALL timeset(routinen, handle)
2150 sub_group = fft_scratch%cart_sub_comm(2)
2151 mip = fft_scratch%mip
2152 dim = fft_scratch%dim
2153 pos = fft_scratch%pos
2154 scount => fft_scratch%scount
2155 rcount => fft_scratch%rcount
2156 sdispl => fft_scratch%sdispl
2157 rdispl => fft_scratch%rdispl
2158 pgrid => fft_scratch%pgcube
2161 nx = boin(2, 1, mip) - boin(1, 1, mip) + 1
2162 ny = boin(2, 2, mip) - boin(1, 2, mip) + 1
2163 mz = maxval(boout(2, 3, :) - boout(1, 3, :) + 1)
2165 rbuf => fft_scratch%rbuf2
2174 nz = boout(2, 3, ipl) - boout(1, 3, ipl) + 1
2175 DO iz = boout(1, 3, ipl), boout(2, 3, ipl)
2176 ir = iz - boout(1, 3, ipl) + 1 + (ixy - 1)*nz
2177 rbuf(ir, ip) = cin(iz, ixy)
2185 nz = boout(2, 3, ipl) - boout(1, 3, ipl) + 1
2186 scount(ip) = nx*ny*nz
2187 sdispl(ip) = nx*ny*mz*ip
2191 nz = boout(2, 3, mip) - boout(1, 3, mip) + 1
2197 ny = boin(2, 2, ipl) - boin(1, 2, ipl) + 1
2198 rcount(ip) = nx*ny*nz
2199 rdispl(ip) = nx*nz*(boin(1, 2, ipl) - 1)
2203 CALL sub_group%alltoall(rbuf, scount, sdispl, sout, rcount, rdispl)
2205 CALL timestop(handle)
2217 SUBROUTINE cube_transpose_3(cin, boin, boout, sout, fft_scratch)
2219 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
2221 INTEGER,
CONTIGUOUS,
DIMENSION(:, :, 0:), &
2222 INTENT(IN) :: boin, boout
2223 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
2227 CHARACTER(len=*),
PARAMETER :: routinen =
'cube_transpose_3'
2229 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
2231 INTEGER :: handle, ip, ipl, ir, is, ixz, iy, lb, &
2232 mip, my, my_id, np, num_threads, nx, &
2234 INTEGER,
CONTIGUOUS,
DIMENSION(:),
POINTER :: rcount, rdispl, scount, sdispl
2235 INTEGER,
CONTIGUOUS,
DIMENSION(:, :),
POINTER :: pgrid
2236 INTEGER,
DIMENSION(2) :: dim, pos
2239 CALL timeset(routinen, handle)
2241 sub_group = fft_scratch%cart_sub_comm(1)
2242 mip = fft_scratch%mip
2243 dim = fft_scratch%dim
2244 pos = fft_scratch%pos
2246 scount => fft_scratch%scount
2247 rcount => fft_scratch%rcount
2248 sdispl => fft_scratch%sdispl
2249 rdispl => fft_scratch%rdispl
2250 pgrid => fft_scratch%pgcube
2252 ny = boin(2, 2, mip) - boin(1, 2, mip) + 1
2253 nz = boin(2, 3, mip) - boin(1, 3, mip) + 1
2259 nx = boout(2, 1, ipl) - boout(1, 1, ipl) + 1
2260 scount(ip) = nx*nz*ny
2261 sdispl(ip) = ny*nz*(boout(1, 1, ipl) - 1)
2264 nx = boout(2, 1, mip) - boout(1, 1, mip) + 1
2265 my = maxval(boin(2, 2, :) - boin(1, 2, :) + 1)
2271 ny = boin(2, 2, ipl) - boin(1, 2, ipl) + 1
2272 rcount(ip) = nx*nz*ny
2273 rdispl(ip) = nx*my*nz*ip
2277 rbuf => fft_scratch%rbuf3
2285 IF (my_id < num_threads)
THEN
2286 lb = (
SIZE(rbuf, 2)*my_id)/num_threads
2287 ub = (
SIZE(rbuf, 2)*(my_id + 1))/num_threads - 1
2288 rbuf(:, lb:ub) = 0.0_dp
2292 CALL sub_group%alltoall(cin, scount, sdispl, rbuf, rcount, rdispl)
2300 ny = boin(2, 2, ipl) - boin(1, 2, ipl) + 1
2302 is = boin(1, 2, ipl) + iy - 1
2303 ir = iy + ny*(ixz - 1)
2304 sout(is, ixz) = rbuf(ir, ip)
2310 CALL timestop(handle)
2312 END SUBROUTINE cube_transpose_3
2322 SUBROUTINE cube_transpose_4(cin, boin, boout, sout, fft_scratch)
2324 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
2326 INTEGER,
CONTIGUOUS,
DIMENSION(:, :, 0:), &
2327 INTENT(IN) :: boin, boout
2328 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
2332 CHARACTER(len=*),
PARAMETER :: routinen =
'cube_transpose_4'
2334 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
2336 INTEGER :: handle, ip, ipl, ir, iy, izx, lb, mip, &
2337 my, my_id, np, num_threads, nx, ny, &
2339 INTEGER,
CONTIGUOUS,
DIMENSION(:),
POINTER :: rcount, rdispl, scount, sdispl
2340 INTEGER,
CONTIGUOUS,
DIMENSION(:, :),
POINTER :: pgrid
2341 INTEGER,
DIMENSION(2) :: dim, pos
2344 CALL timeset(routinen, handle)
2346 sub_group = fft_scratch%cart_sub_comm(1)
2347 mip = fft_scratch%mip
2348 dim = fft_scratch%dim
2349 pos = fft_scratch%pos
2351 scount => fft_scratch%scount
2352 rcount => fft_scratch%rcount
2353 sdispl => fft_scratch%sdispl
2354 rdispl => fft_scratch%rdispl
2355 pgrid => fft_scratch%pgcube
2357 nx = boin(2, 1, mip) - boin(1, 1, mip) + 1
2358 nz = boin(2, 3, mip) - boin(1, 3, mip) + 1
2359 my = maxval(boout(2, 2, :) - boout(1, 2, :) + 1)
2361 rbuf => fft_scratch%rbuf4
2369 IF (my_id < num_threads)
THEN
2370 lb = (
SIZE(rbuf, 2)*my_id)/num_threads
2371 ub = (
SIZE(rbuf, 2)*(my_id + 1))/num_threads - 1
2372 rbuf(:, lb:ub) = 0.0_dp
2380 ny = boout(2, 2, ipl) - boout(1, 2, ipl) + 1
2381 DO iy = boout(1, 2, ipl), boout(2, 2, ipl)
2382 ir = iy - boout(1, 2, ipl) + 1 + (izx - 1)*ny
2383 rbuf(ir, ip) = cin(iy, izx)
2391 ny = boout(2, 2, ipl) - boout(1, 2, ipl) + 1
2392 scount(ip) = nx*ny*nz
2393 sdispl(ip) = nx*nz*my*ip
2397 ny = boout(2, 2, mip) - boout(1, 2, mip) + 1
2403 nx = boin(2, 1, ipl) - boin(1, 1, ipl) + 1
2404 rcount(ip) = nx*ny*nz
2405 rdispl(ip) = ny*nz*(boin(1, 1, ipl) - 1)
2409 CALL sub_group%alltoall(rbuf, scount, sdispl, sout, rcount, rdispl)
2411 CALL timestop(handle)
2413 END SUBROUTINE cube_transpose_4
2424 SUBROUTINE cube_transpose_5(cin, group, boin, boout, sout, fft_scratch)
2426 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
2430 INTEGER,
CONTIGUOUS,
DIMENSION(:, :, 0:),
INTENT(IN) :: boin, boout
2431 COMPLEX(KIND=dp),
DIMENSION(:, :),
INTENT(OUT),
CONTIGUOUS :: sout
2434 CHARACTER(len=*),
PARAMETER :: routinen =
'cube_transpose_5'
2436 COMPLEX(KIND=dp),
DIMENSION(:, :),
POINTER,
CONTIGUOUS :: rbuf
2437 INTEGER :: handle, ip, ir, is, ixz, iy, lb, mip, &
2438 my, my_id, np, num_threads, nx, ny, &
2440 INTEGER,
DIMENSION(:),
POINTER,
CONTIGUOUS :: rcount, rdispl, scount, sdispl
2442 CALL timeset(routinen, handle)
2444 np = fft_scratch%sizes%numtask
2445 mip = fft_scratch%mip
2446 scount => fft_scratch%scount
2447 rcount => fft_scratch%rcount
2448 sdispl => fft_scratch%sdispl
2449 rdispl => fft_scratch%rdispl
2451 ny = boin(2, 2, mip) - boin(1, 2, mip) + 1
2452 nz = boin(2, 3, mip) - boin(1, 3, mip) + 1
2457 nx = boout(2, 1, ip) - boout(1, 1, ip) + 1
2458 scount(ip) = nx*nz*ny
2459 sdispl(ip) = ny*nz*(boout(1, 1, ip) - 1)
2462 nx = boout(2, 1, mip) - boout(1, 1, mip) + 1
2463 my = maxval(boin(2, 2, :) - boin(1, 2, :) + 1)
2468 ny = boin(2, 2, ip) - boin(1, 2, ip) + 1
2469 rcount(ip) = nx*nz*ny
2470 rdispl(ip) = nx*my*nz*ip
2474 rbuf => fft_scratch%rbuf5
2482 IF (my_id < num_threads)
THEN
2483 lb = (
SIZE(rbuf, 2)*my_id)/num_threads
2484 ub = (
SIZE(rbuf, 2)*(my_id + 1))/num_threads - 1
2485 rbuf(:, lb:ub) = 0.0_dp
2489 CALL group%alltoall(cin, scount, sdispl, rbuf, rcount, rdispl)
2496 ny = boin(2, 2, ip) - boin(1, 2, ip) + 1
2498 is = boin(1, 2, ip) + iy - 1
2499 ir = iy + ny*(ixz - 1)
2500 sout(is, ixz) = rbuf(ir, ip)
2506 CALL timestop(handle)
2508 END SUBROUTINE cube_transpose_5
2519 SUBROUTINE cube_transpose_6(cin, group, boin, boout, sout, fft_scratch)
2521 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:, :), &
2525 INTEGER,
CONTIGUOUS,
DIMENSION(:, :, 0:),
INTENT(IN) :: boin, boout
2526 COMPLEX(KIND=dp),
DIMENSION(:, :),
INTENT(OUT),
CONTIGUOUS :: sout
2529 CHARACTER(len=*),
PARAMETER :: routinen =
'cube_transpose_6'
2531 COMPLEX(KIND=dp),
DIMENSION(:, :),
POINTER,
CONTIGUOUS :: rbuf
2532 INTEGER :: handle, ip, ir, iy, izx, lb, mip, my, &
2533 my_id, np, num_threads, nx, ny, nz, ub
2534 INTEGER,
DIMENSION(:),
POINTER,
CONTIGUOUS :: rcount, rdispl, scount, sdispl
2536 CALL timeset(routinen, handle)
2538 np = fft_scratch%sizes%numtask
2539 mip = fft_scratch%mip
2540 scount => fft_scratch%scount
2541 rcount => fft_scratch%rcount
2542 sdispl => fft_scratch%sdispl
2543 rdispl => fft_scratch%rdispl
2545 nx = boin(2, 1, mip) - boin(1, 1, mip) + 1
2546 nz = boin(2, 3, mip) - boin(1, 3, mip) + 1
2547 my = maxval(boout(2, 2, :) - boout(1, 2, :) + 1)
2549 rbuf => fft_scratch%rbuf5
2557 IF (my_id < num_threads)
THEN
2558 lb = (
SIZE(rbuf, 2)*my_id)/num_threads
2559 ub = (
SIZE(rbuf, 2)*(my_id + 1))/num_threads - 1
2560 rbuf(:, lb:ub) = 0.0_dp
2567 ny = boout(2, 2, ip) - boout(1, 2, ip) + 1
2568 DO iy = boout(1, 2, ip), boout(2, 2, ip)
2569 ir = iy - boout(1, 2, ip) + 1 + (izx - 1)*ny
2570 rbuf(ir, ip) = cin(iy, izx)
2577 ny = boout(2, 2, ip) - boout(1, 2, ip) + 1
2578 scount(ip) = nx*ny*nz
2579 sdispl(ip) = nx*nz*my*ip
2583 ny = boout(2, 2, mip) - boout(1, 2, mip) + 1
2588 nx = boin(2, 1, ip) - boin(1, 1, ip) + 1
2589 rcount(ip) = nx*ny*nz
2590 rdispl(ip) = ny*nz*(boin(1, 1, ip) - 1)
2594 CALL group%alltoall(rbuf, scount, sdispl, sout, rcount, rdispl)
2596 CALL timestop(handle)
2598 END SUBROUTINE cube_transpose_6
2605 CALL release_fft_scratch_pool()
2621 SUBROUTINE deallocate_fft_scratch_type(fft_scratch)
2624#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2626 COMPLEX(KIND=dp),
POINTER :: dummy_ptr_z
2630 IF (
ASSOCIATED(fft_scratch%ziptr))
THEN
2631 CALL fft_dealloc(fft_scratch%ziptr)
2633 IF (
ASSOCIATED(fft_scratch%zoptr))
THEN
2634 CALL fft_dealloc(fft_scratch%zoptr)
2636 IF (
ASSOCIATED(fft_scratch%p1buf))
THEN
2637 CALL fft_dealloc(fft_scratch%p1buf)
2639 IF (
ASSOCIATED(fft_scratch%p2buf))
THEN
2640#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2641 dummy_ptr_z => fft_scratch%p2buf(1, 1)
2644 CALL fft_dealloc(fft_scratch%p2buf)
2647 IF (
ASSOCIATED(fft_scratch%p3buf))
THEN
2648#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2649 dummy_ptr_z => fft_scratch%p3buf(1, 1)
2652 CALL fft_dealloc(fft_scratch%p3buf)
2655 IF (
ASSOCIATED(fft_scratch%p4buf))
THEN
2656#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2657 dummy_ptr_z => fft_scratch%p4buf(1, 1)
2660 CALL fft_dealloc(fft_scratch%p4buf)
2663 IF (
ASSOCIATED(fft_scratch%p5buf))
THEN
2664#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2665 dummy_ptr_z => fft_scratch%p5buf(1, 1)
2668 CALL fft_dealloc(fft_scratch%p5buf)
2671 IF (
ASSOCIATED(fft_scratch%p6buf))
THEN
2672 CALL fft_dealloc(fft_scratch%p6buf)
2674 IF (
ASSOCIATED(fft_scratch%p7buf))
THEN
2675#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2676 dummy_ptr_z => fft_scratch%p7buf(1, 1)
2679 CALL fft_dealloc(fft_scratch%p7buf)
2682 IF (
ASSOCIATED(fft_scratch%r1buf))
THEN
2683#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2684 dummy_ptr_z => fft_scratch%r1buf(1, 1)
2687 CALL fft_dealloc(fft_scratch%r1buf)
2690 IF (
ASSOCIATED(fft_scratch%r2buf))
THEN
2691 CALL fft_dealloc(fft_scratch%r2buf)
2693 IF (
ASSOCIATED(fft_scratch%tbuf))
THEN
2694#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2695 dummy_ptr_z => fft_scratch%tbuf(1, 1, 1)
2698 CALL fft_dealloc(fft_scratch%tbuf)
2701 IF (
ASSOCIATED(fft_scratch%a1buf))
THEN
2702 CALL fft_dealloc(fft_scratch%a1buf)
2704 IF (
ASSOCIATED(fft_scratch%a2buf))
THEN
2705 CALL fft_dealloc(fft_scratch%a2buf)
2707 IF (
ASSOCIATED(fft_scratch%a3buf))
THEN
2708 CALL fft_dealloc(fft_scratch%a3buf)
2710 IF (
ASSOCIATED(fft_scratch%a4buf))
THEN
2711 CALL fft_dealloc(fft_scratch%a4buf)
2713 IF (
ASSOCIATED(fft_scratch%a5buf))
THEN
2714 CALL fft_dealloc(fft_scratch%a5buf)
2716 IF (
ASSOCIATED(fft_scratch%a6buf))
THEN
2717 CALL fft_dealloc(fft_scratch%a6buf)
2719 IF (
ASSOCIATED(fft_scratch%scount))
THEN
2720 DEALLOCATE (fft_scratch%scount, fft_scratch%rcount, &
2721 fft_scratch%sdispl, fft_scratch%rdispl)
2723 IF (
ASSOCIATED(fft_scratch%rr))
THEN
2724 DEALLOCATE (fft_scratch%rr)
2726 IF (
ASSOCIATED(fft_scratch%xzbuf))
THEN
2727 DEALLOCATE (fft_scratch%xzbuf)
2729 IF (
ASSOCIATED(fft_scratch%yzbuf))
THEN
2730 DEALLOCATE (fft_scratch%yzbuf)
2732 IF (
ASSOCIATED(fft_scratch%xzbuf_sgl))
THEN
2733 DEALLOCATE (fft_scratch%xzbuf_sgl)
2735 IF (
ASSOCIATED(fft_scratch%yzbuf_sgl))
THEN
2736 DEALLOCATE (fft_scratch%yzbuf_sgl)
2738 IF (
ASSOCIATED(fft_scratch%ss))
THEN
2739 DEALLOCATE (fft_scratch%ss)
2741 IF (
ASSOCIATED(fft_scratch%tt))
THEN
2742 DEALLOCATE (fft_scratch%tt)
2744 IF (
ASSOCIATED(fft_scratch%pgrid))
THEN
2745 DEALLOCATE (fft_scratch%pgrid)
2747 IF (
ASSOCIATED(fft_scratch%pgcube))
THEN
2748 DEALLOCATE (fft_scratch%pgcube)
2750 IF (
ASSOCIATED(fft_scratch%xcor))
THEN
2751 DEALLOCATE (fft_scratch%xcor, fft_scratch%zcor)
2753 IF (
ASSOCIATED(fft_scratch%pzcoord))
THEN
2754 DEALLOCATE (fft_scratch%pzcoord)
2756 IF (
ASSOCIATED(fft_scratch%xzcount))
THEN
2757 DEALLOCATE (fft_scratch%xzcount, fft_scratch%yzcount)
2758 DEALLOCATE (fft_scratch%xzdispl, fft_scratch%yzdispl)
2760 fft_scratch%rsratio = 1._dp
2762 IF (
ASSOCIATED(fft_scratch%rbuf1))
THEN
2763 DEALLOCATE (fft_scratch%rbuf1)
2765 IF (
ASSOCIATED(fft_scratch%rbuf2))
THEN
2766 DEALLOCATE (fft_scratch%rbuf2)
2768 IF (
ASSOCIATED(fft_scratch%rbuf3))
THEN
2769 DEALLOCATE (fft_scratch%rbuf3)
2771 IF (
ASSOCIATED(fft_scratch%rbuf4))
THEN
2772 DEALLOCATE (fft_scratch%rbuf4)
2774 IF (
ASSOCIATED(fft_scratch%rbuf5))
THEN
2775 DEALLOCATE (fft_scratch%rbuf5)
2777 IF (
ASSOCIATED(fft_scratch%rbuf6))
THEN
2778 DEALLOCATE (fft_scratch%rbuf6)
2782 CALL fft_scratch%cart_sub_comm(1)%free()
2785 CALL fft_scratch%cart_sub_comm(2)%free()
2795 END SUBROUTINE deallocate_fft_scratch_type
2800 SUBROUTINE release_fft_scratch_pool()
2809 IF (
ASSOCIATED(fft_scratch))
THEN
2810 fft_scratch_current => fft_scratch
2811 fft_scratch => fft_scratch_current%fft_scratch_next
2812 NULLIFY (fft_scratch_current%fft_scratch_next)
2814 CALL deallocate_fft_scratch_type(fft_scratch_current%fft_scratch)
2816 DEALLOCATE (fft_scratch_current%fft_scratch)
2817 DEALLOCATE (fft_scratch_current)
2825 END SUBROUTINE release_fft_scratch_pool
2830 SUBROUTINE resize_fft_scratch_pool()
2832 INTEGER :: last_tick, nscratch
2837 last_tick = huge(last_tick)
2838 NULLIFY (fft_scratch_old)
2843 IF (
ASSOCIATED(fft_scratch_current))
THEN
2844 nscratch = nscratch + 1
2846 IF (.NOT. fft_scratch_current%fft_scratch%in_use)
THEN
2847 IF (fft_scratch_current%fft_scratch%last_tick < last_tick)
THEN
2848 last_tick = fft_scratch_current%fft_scratch%last_tick
2849 fft_scratch_old => fft_scratch_current
2852 fft_scratch_current => fft_scratch_current%fft_scratch_next
2861 IF (
ASSOCIATED(fft_scratch_old))
THEN
2864 IF (
ASSOCIATED(fft_scratch_current))
THEN
2866 IF (
ASSOCIATED(fft_scratch_current%fft_scratch_next, fft_scratch_old))
THEN
2868 fft_scratch_current%fft_scratch_next => fft_scratch_old%fft_scratch_next
2871 CALL deallocate_fft_scratch_type(fft_scratch_old%fft_scratch)
2872 DEALLOCATE (fft_scratch_old%fft_scratch)
2873 DEALLOCATE (fft_scratch_old)
2876 fft_scratch_current => fft_scratch_current%fft_scratch_next
2884 cpwarn(
"The number of the scratches exceeded the limit, but none could be deallocated")
2888 END SUBROUTINE resize_fft_scratch_pool
2899 INTEGER,
INTENT(IN) :: tf_type
2900 INTEGER,
DIMENSION(:),
INTENT(IN) :: n
2902 OPTIONAL :: fft_sizes
2904 CHARACTER(len=*),
PARAMETER :: routinen =
'get_fft_scratch'
2906 INTEGER :: coord(2), dim(2), handle, i, ix, iz, lg, lmax, m1, m2, &
2907 mcx2, mcy3, mcz1, mcz2, mg, mmax, mx1, mx2, my1, my3, mz1, mz2, mz3, &
2908 nbx, nbz, nm, nmax, nmray, np, nx, ny, nyzray, nz, pos(2)
2909 INTEGER,
DIMENSION(3) :: pcoord
2911 LOGICAL,
DIMENSION(2) :: dims
2915#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2917 INTEGER(KIND=C_SIZE_T) :: length
2918 TYPE(c_ptr) :: cptr_r1buf, cptr_tbuf, &
2919 cptr_p2buf, cptr_p3buf, cptr_p4buf, cptr_p5buf, cptr_p7buf
2921 CALL timeset(routinen, handle)
2925 CALL resize_fft_scratch_pool()
2932 IF (
ASSOCIATED(fft_scratch_current))
THEN
2933 IF (fft_scratch_current%fft_scratch%in_use)
THEN
2934 fft_scratch_last => fft_scratch_current
2935 fft_scratch_current => fft_scratch_current%fft_scratch_next
2938 IF (tf_type /= fft_scratch_current%fft_scratch%tf_type)
THEN
2939 fft_scratch_last => fft_scratch_current
2940 fft_scratch_current => fft_scratch_current%fft_scratch_next
2943 IF (.NOT. all(n == fft_scratch_current%fft_scratch%nfft))
THEN
2944 fft_scratch_last => fft_scratch_current
2945 fft_scratch_current => fft_scratch_current%fft_scratch_next
2948 IF (
PRESENT(fft_sizes))
THEN
2949 IF (fft_sizes%rs_group /= fft_scratch_current%fft_scratch%group)
THEN
2950 fft_scratch_last => fft_scratch_current
2951 fft_scratch_current => fft_scratch_current%fft_scratch_next
2954 CALL is_equal(fft_sizes, fft_scratch_current%fft_scratch%sizes, equal)
2955 IF (.NOT. equal)
THEN
2956 fft_scratch_last => fft_scratch_current
2957 fft_scratch_current => fft_scratch_current%fft_scratch_next
2962 fft_scratch => fft_scratch_current%fft_scratch
2963 fft_scratch_current%fft_scratch%in_use = .true.
2968 ALLOCATE (fft_scratch_new)
2969 ALLOCATE (fft_scratch_new%fft_scratch)
2971 IF (tf_type /= 400)
THEN
2972 fft_scratch_new%fft_scratch%sizes = fft_sizes
2973 np = fft_sizes%numtask
2974 ALLOCATE (fft_scratch_new%fft_scratch%scount(0:np - 1), fft_scratch_new%fft_scratch%rcount(0:np - 1), &
2975 fft_scratch_new%fft_scratch%sdispl(0:np - 1), fft_scratch_new%fft_scratch%rdispl(0:np - 1), &
2976 fft_scratch_new%fft_scratch%pgcube(0:np - 1, 2))
2979 SELECT CASE (tf_type)
2981 cpabort(
"Invalid scratch type.")
2983 cpassert(
PRESENT(fft_sizes))
2990 CALL fft_alloc(fft_scratch_new%fft_scratch%a1buf, [mx1*my1, n(3)])
2991 CALL fft_alloc(fft_scratch_new%fft_scratch%a2buf, [n(3), mx1*my1])
2992 CALL fft_alloc(fft_scratch_new%fft_scratch%a3buf, [mx2*mz2, n(2)])
2993 CALL fft_alloc(fft_scratch_new%fft_scratch%a4buf, [n(2), mx2*mz2])
2994 CALL fft_alloc(fft_scratch_new%fft_scratch%a5buf, [my3*mz3, n(1)])
2995 CALL fft_alloc(fft_scratch_new%fft_scratch%a6buf, [n(1), my3*mz3])
2996 fft_scratch_new%fft_scratch%group = fft_sizes%rs_group
2998 dim = fft_sizes%rs_group%num_pe_cart
2999 pos = fft_sizes%rs_group%mepos_cart
3000 fft_scratch_new%fft_scratch%mip = fft_sizes%rs_group%mepos
3001 fft_scratch_new%fft_scratch%dim = dim
3002 fft_scratch_new%fft_scratch%pos = pos
3003 mcz1 = fft_sizes%mcz1
3004 mcx2 = fft_sizes%mcx2
3005 mcz2 = fft_sizes%mcz2
3006 mcy3 = fft_sizes%mcy3
3007 ALLOCATE (fft_scratch_new%fft_scratch%rbuf1(mx2*my1*mcz2, 0:dim(2) - 1))
3008 ALLOCATE (fft_scratch_new%fft_scratch%rbuf2(mx1*my1*mcz2, 0:dim(2) - 1))
3009 ALLOCATE (fft_scratch_new%fft_scratch%rbuf3(mx2*mz3*mcy3, 0:dim(1) - 1))
3010 ALLOCATE (fft_scratch_new%fft_scratch%rbuf4(mx2*mz2*mcy3, 0:dim(1) - 1))
3012 dims = [.true., .false.]
3013 CALL fft_scratch_new%fft_scratch%cart_sub_comm(1)%from_sub(fft_sizes%rs_group, dims)
3014 dims = [.false., .true.]
3015 CALL fft_scratch_new%fft_scratch%cart_sub_comm(2)%from_sub(fft_sizes%rs_group, dims)
3018 DO i = 0, dim(1) - 1
3020 CALL fft_sizes%rs_group%rank_cart(coord, fft_scratch_new%fft_scratch%pgcube(i, 1))
3022 DO i = 0, dim(2) - 1
3024 CALL fft_sizes%rs_group%rank_cart(coord, fft_scratch_new%fft_scratch%pgcube(i, 2))
3029 mx1*my1, n(3), n(3), mx1*my1, &
3030 fft_scratch_new%fft_scratch%a1buf, fft_scratch_new%fft_scratch%a2buf)
3032 mx2*mz2, n(2), n(2), mx2*mz2, &
3033 fft_scratch_new%fft_scratch%a3buf, fft_scratch_new%fft_scratch%a4buf)
3035 my3*mz3, n(1), n(1), my3*mz3, &
3036 fft_scratch_new%fft_scratch%a5buf, fft_scratch_new%fft_scratch%a6buf)
3038 n(1), my3*mz3, n(1), my3*mz3, &
3039 fft_scratch_new%fft_scratch%a6buf, fft_scratch_new%fft_scratch%a5buf)
3041 n(2), mx2*mz2, n(2), mx2*mz2, &
3042 fft_scratch_new%fft_scratch%a4buf, fft_scratch_new%fft_scratch%a3buf)
3044 n(3), mx1*my1, n(3), mx1*my1, &
3045 fft_scratch_new%fft_scratch%a2buf, fft_scratch_new%fft_scratch%a1buf)
3048 cpassert(
PRESENT(fft_sizes))
3054 CALL fft_alloc(fft_scratch_new%fft_scratch%a1buf, [mx1*my1, n(3)])
3055 CALL fft_alloc(fft_scratch_new%fft_scratch%a2buf, [n(3), mx1*my1])
3056 fft_scratch_new%fft_scratch%group = fft_sizes%rs_group
3057 CALL fft_alloc(fft_scratch_new%fft_scratch%a3buf, [mx1*mz1, n(2)])
3058 CALL fft_alloc(fft_scratch_new%fft_scratch%a4buf, [n(2), mx1*mz1])
3059 CALL fft_alloc(fft_scratch_new%fft_scratch%a5buf, [my3*mz3, n(1)])
3060 CALL fft_alloc(fft_scratch_new%fft_scratch%a6buf, [n(1), my3*mz3])
3062 dim = fft_sizes%rs_group%num_pe_cart
3063 pos = fft_sizes%rs_group%mepos_cart
3064 fft_scratch_new%fft_scratch%mip = fft_sizes%rs_group%mepos
3065 fft_scratch_new%fft_scratch%dim = dim
3066 fft_scratch_new%fft_scratch%pos = pos
3067 mcy3 = fft_sizes%mcy3
3068 ALLOCATE (fft_scratch_new%fft_scratch%rbuf5(mx1*mz3*mcy3, 0:dim(1) - 1))
3069 ALLOCATE (fft_scratch_new%fft_scratch%rbuf6(mx1*mz1*mcy3, 0:dim(1) - 1))
3073 mx1*my1, n(3), n(3), mx1*my1, &
3074 fft_scratch_new%fft_scratch%a1buf, fft_scratch_new%fft_scratch%a3buf)
3076 mx1*mz1, n(2), n(2), mx1*mz1, &
3077 fft_scratch_new%fft_scratch%a3buf, fft_scratch_new%fft_scratch%a4buf)
3079 my3*mz3, n(1), n(1), my3*mz3, &
3080 fft_scratch_new%fft_scratch%a5buf, fft_scratch_new%fft_scratch%a6buf)
3082 n(1), my3*mz3, n(1), my3*mz3, &
3083 fft_scratch_new%fft_scratch%a6buf, fft_scratch_new%fft_scratch%a5buf)
3085 n(2), mx1*mz1, n(2), mx1*mz1, &
3086 fft_scratch_new%fft_scratch%a4buf, fft_scratch_new%fft_scratch%a3buf)
3088 n(3), mx1*my1, n(3), mx1*my1, &
3089 fft_scratch_new%fft_scratch%a3buf, fft_scratch_new%fft_scratch%a1buf)
3092 cpassert(
PRESENT(fft_sizes))
3097 lmax = fft_sizes%lmax
3098 mmax = fft_sizes%mmax
3101 np = fft_sizes%numtask
3102 nmray = fft_sizes%nmray
3103 nyzray = fft_sizes%nyzray
3104#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
3105 length = int(2*
dp_size*max(mmax, 1)*max(lmax, 1), kind=c_size_t)
3108 CALL c_f_pointer(cptr_r1buf, fft_scratch_new%fft_scratch%r1buf, [max(mmax, 1), max(lmax, 1)])
3109 length = int(2*
dp_size*max(ny, 1)*max(nz, 1)*max(nx, 1), kind=c_size_t)
3112 CALL c_f_pointer(cptr_tbuf, fft_scratch_new%fft_scratch%tbuf, [max(ny, 1), max(nz, 1), max(nx, 1)])
3114 CALL fft_alloc(fft_scratch_new%fft_scratch%r1buf, [mmax, lmax])
3115 CALL fft_alloc(fft_scratch_new%fft_scratch%tbuf, [ny, nz, nx])
3117 fft_scratch_new%fft_scratch%group = fft_sizes%rs_group
3118 CALL fft_alloc(fft_scratch_new%fft_scratch%r2buf, [lg, mg])
3120 IF (alltoall_sgl)
THEN
3121 ALLOCATE (fft_scratch_new%fft_scratch%ss(mmax, lmax))
3122 ALLOCATE (fft_scratch_new%fft_scratch%tt(nm, 0:np - 1))
3124 ALLOCATE (fft_scratch_new%fft_scratch%rr(nm, 0:np - 1))
3129 nx*ny, nz, nz, nx*ny, &
3130 fft_scratch_new%fft_scratch%tbuf, fft_scratch_new%fft_scratch%r1buf)
3132 nx*nz, ny, ny, nx*nz, &
3133 fft_scratch_new%fft_scratch%r1buf, fft_scratch_new%fft_scratch%tbuf)
3135 nyzray, n(1), n(1), nyzray, &
3136 fft_scratch_new%fft_scratch%r1buf, fft_scratch_new%fft_scratch%r2buf)
3138 n(1), nyzray, n(1), nyzray, &
3139 fft_scratch_new%fft_scratch%r2buf, fft_scratch_new%fft_scratch%r1buf)
3141 ny, nx*nz, ny, nx*nz, &
3142 fft_scratch_new%fft_scratch%tbuf, fft_scratch_new%fft_scratch%r1buf)
3144 nz, nx*ny, nz, nx*ny, &
3145 fft_scratch_new%fft_scratch%r1buf, fft_scratch_new%fft_scratch%tbuf)
3148 cpassert(
PRESENT(fft_sizes))
3153 mcx2 = fft_sizes%mcx2
3156 nmax = fft_sizes%nmax
3157 nmray = fft_sizes%nmray
3158 nyzray = fft_sizes%nyzray
3159 m1 = fft_sizes%r_dim(1)
3160 m2 = fft_sizes%r_dim(2)
3163 CALL fft_alloc(fft_scratch_new%fft_scratch%p1buf, [mx1*my1, n(3)])
3164 CALL fft_alloc(fft_scratch_new%fft_scratch%p6buf, [lg, mg])
3165#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
3166 length = int(2*
dp_size*max(n(3), 1)*max(mx1*my1, 1), kind=c_size_t)
3169 CALL c_f_pointer(cptr_p2buf, fft_scratch_new%fft_scratch%p2buf, [max(n(3), 1), max(mx1*my1, 1)])
3170 length = int(2*
dp_size*max(mx2*mz2, 1)*max(n(2), 1), kind=c_size_t)
3173 CALL c_f_pointer(cptr_p3buf, fft_scratch_new%fft_scratch%p3buf, [max(mx2*mz2, 1), max(n(2), 1)])
3174 length = int(2*
dp_size*max(n(2), 1)*max(mx2*mz2, 1), kind=c_size_t)
3177 CALL c_f_pointer(cptr_p4buf, fft_scratch_new%fft_scratch%p4buf, [max(n(2), 1), max(mx2*mz2, 1)])
3178 length = int(2*
dp_size*max(nyzray, 1)*max(n(1), 1), kind=c_size_t)
3181 CALL c_f_pointer(cptr_p5buf, fft_scratch_new%fft_scratch%p5buf, [max(nyzray, 1), max(n(1), 1)])
3182 length = int(2*
dp_size*max(mg, 1)*max(lg, 1), kind=c_size_t)
3185 CALL c_f_pointer(cptr_p7buf, fft_scratch_new%fft_scratch%p7buf, [max(mg, 1), max(lg, 1)])
3187 CALL fft_alloc(fft_scratch_new%fft_scratch%p2buf, [n(3), mx1*my1])
3188 CALL fft_alloc(fft_scratch_new%fft_scratch%p3buf, [mx2*mz2, n(2)])
3189 CALL fft_alloc(fft_scratch_new%fft_scratch%p4buf, [n(2), mx2*mz2])
3190 CALL fft_alloc(fft_scratch_new%fft_scratch%p5buf, [nyzray, n(1)])
3191 CALL fft_alloc(fft_scratch_new%fft_scratch%p7buf, [mg, lg])
3193 IF (alltoall_sgl)
THEN
3194 ALLOCATE (fft_scratch_new%fft_scratch%yzbuf_sgl(mg*lg))
3195 ALLOCATE (fft_scratch_new%fft_scratch%xzbuf_sgl(n(2)*mx2*mz2))
3197 ALLOCATE (fft_scratch_new%fft_scratch%yzbuf(mg*lg))
3198 ALLOCATE (fft_scratch_new%fft_scratch%xzbuf(n(2)*mx2*mz2))
3200 ALLOCATE (fft_scratch_new%fft_scratch%pgrid(0:m1 - 1, 0:m2 - 1))
3201 ALLOCATE (fft_scratch_new%fft_scratch%xcor(nbx))
3202 ALLOCATE (fft_scratch_new%fft_scratch%zcor(nbz))
3203 ALLOCATE (fft_scratch_new%fft_scratch%pzcoord(0:np - 1))
3204 ALLOCATE (fft_scratch_new%fft_scratch%xzcount(0:np - 1), &
3205 fft_scratch_new%fft_scratch%yzcount(0:np - 1))
3206 ALLOCATE (fft_scratch_new%fft_scratch%xzdispl(0:np - 1), &
3207 fft_scratch_new%fft_scratch%yzdispl(0:np - 1))
3208 fft_scratch_new%fft_scratch%group = fft_sizes%rs_group
3210 dim = fft_sizes%rs_group%num_pe_cart
3211 pos = fft_sizes%rs_group%mepos_cart
3212 fft_scratch_new%fft_scratch%mip = fft_sizes%rs_group%mepos
3213 fft_scratch_new%fft_scratch%dim = dim
3214 fft_scratch_new%fft_scratch%pos = pos
3215 mcz1 = fft_sizes%mcz1
3216 mcz2 = fft_sizes%mcz2
3217 ALLOCATE (fft_scratch_new%fft_scratch%rbuf1(mx2*my1*mcz2, 0:dim(2) - 1))
3218 ALLOCATE (fft_scratch_new%fft_scratch%rbuf2(mx1*my1*mcz2, 0:dim(2) - 1))
3220 dims = [.false., .true.]
3221 CALL fft_scratch_new%fft_scratch%cart_sub_comm(2)%from_sub(fft_sizes%rs_group, dims)
3224 DO i = 0, dim(2) - 1
3226 CALL fft_sizes%rs_group%rank_cart(coord, fft_scratch_new%fft_scratch%pgcube(i, 2))
3233 CALL fft_sizes%rs_group%rank_cart(coord, fft_scratch_new%fft_scratch%pgrid(ix, iz))
3239 CALL fft_sizes%rs_group%coords(i, pcoord)
3240 fft_scratch_new%fft_scratch%pzcoord(i) = pcoord(2)
3248 cpassert(
SIZE(fft_scratch_new%fft_scratch%p1buf) >= mx1*my1*n(3))
3250 mx1*my1, n(3), n(3), mx1*my1, &
3251 fft_scratch_new%fft_scratch%p1buf, fft_scratch_new%fft_scratch%p2buf)
3252 cpassert(
SIZE(fft_scratch_new%fft_scratch%p3buf) >= mx2*mz2*n(2))
3254 mx2*mz2, n(2), n(2), mx2*mz2, &
3255 fft_scratch_new%fft_scratch%p3buf, fft_scratch_new%fft_scratch%p4buf)
3256 cpassert(
SIZE(fft_scratch_new%fft_scratch%p5buf) >= nyzray*n(1))
3258 nyzray, n(1), n(1), nyzray, &
3259 fft_scratch_new%fft_scratch%p5buf, fft_scratch_new%fft_scratch%p6buf)
3260 cpassert(
SIZE(fft_scratch_new%fft_scratch%p6buf) >= n(1)*nyzray)
3262 n(1), nyzray, n(1), nyzray, &
3263 fft_scratch_new%fft_scratch%p6buf, fft_scratch_new%fft_scratch%p7buf)
3264 cpassert(
SIZE(fft_scratch_new%fft_scratch%p4buf) >= n(2)*mx2*mz2)
3266 n(2), mx2*mz2, n(2), mx2*mz2, &
3267 fft_scratch_new%fft_scratch%p4buf, fft_scratch_new%fft_scratch%p3buf)
3268 cpassert(
SIZE(fft_scratch_new%fft_scratch%p2buf) >= n(3)*mx1*my1)
3270 n(3), mx1*my1, n(3), mx1*my1, &
3271 fft_scratch_new%fft_scratch%p2buf, fft_scratch_new%fft_scratch%p1buf)
3275 CALL fft_alloc(fft_scratch_new%fft_scratch%ziptr, n)
3276 CALL fft_alloc(fft_scratch_new%fft_scratch%zoptr, n)
3280 fft_scratch_new%fft_scratch%ziptr, fft_scratch_new%fft_scratch%ziptr)
3282 fft_scratch_new%fft_scratch%ziptr, fft_scratch_new%fft_scratch%ziptr)
3285 fft_scratch_new%fft_scratch%ziptr, fft_scratch_new%fft_scratch%zoptr)
3287 fft_scratch_new%fft_scratch%ziptr, fft_scratch_new%fft_scratch%zoptr)
3291 NULLIFY (fft_scratch_new%fft_scratch_next)
3292 fft_scratch_new%fft_scratch%fft_scratch_id = &
3293 fft_scratch_last%fft_scratch%fft_scratch_id + 1
3294 fft_scratch_new%fft_scratch%in_use = .true.
3295 fft_scratch_new%fft_scratch%nfft = n
3296 fft_scratch_last%fft_scratch_next => fft_scratch_new
3297 fft_scratch_new%fft_scratch%tf_type = tf_type
3298 fft_scratch => fft_scratch_new%fft_scratch
3307 CALL timestop(handle)
3319 INTEGER :: scratch_id
3322 scratch_id = fft_scratch%fft_scratch_id
3326 IF (
ASSOCIATED(fft_scratch_current))
THEN
3327 IF (scratch_id == fft_scratch_current%fft_scratch%fft_scratch_id)
THEN
3328 fft_scratch%in_use = .false.
3329 NULLIFY (fft_scratch)
3332 fft_scratch_current => fft_scratch_current%fft_scratch_next
3335 cpabort(
"Invalid scratch type.")
3352 SUBROUTINE sparse_alltoall(rs, scount, sdispl, rq, rcount, rdispl, group)
3353 COMPLEX(KIND=dp),
DIMENSION(:),
POINTER :: rs
3354 INTEGER,
DIMENSION(:),
POINTER :: scount, sdispl
3355 COMPLEX(KIND=dp),
DIMENSION(:),
POINTER :: rq
3356 INTEGER,
DIMENSION(:),
POINTER :: rcount, rdispl
3360 COMPLEX(KIND=dp),
DIMENSION(:),
POINTER :: msgin, msgout
3361 INTEGER :: ip, n, nr, ns, pos
3367 ALLOCATE (sreq(0:n - 1))
3368 ALLOCATE (rreq(0:n - 1))
3371 IF (rcount(ip) == 0) cycle
3372 IF (ip == pos) cycle
3373 msgout => rq(rdispl(ip) + 1:rdispl(ip) + rcount(ip))
3374 CALL group%irecv(msgout, ip, rreq(nr))
3379 IF (scount(ip) == 0) cycle
3380 IF (ip == pos) cycle
3381 msgin => rs(sdispl(ip) + 1:sdispl(ip) + scount(ip))
3382 CALL group%isend(msgin, ip, sreq(ns))
3385 IF (rcount(pos) /= 0)
THEN
3386 IF (rcount(pos) /= scount(pos)) cpabort(
"Invalid count.")
3387 rq(rdispl(pos) + 1:rdispl(pos) + rcount(pos)) = rs(sdispl(pos) + 1:sdispl(pos) + scount(pos))
3395 END SUBROUTINE sparse_alltoall
3404 SUBROUTINE is_equal(fft_size_1, fft_size_2, equal)
3410 equal = equal .AND. fft_size_1%nx == fft_size_2%nx
3411 equal = equal .AND. fft_size_1%ny == fft_size_2%ny
3412 equal = equal .AND. fft_size_1%nz == fft_size_2%nz
3414 equal = equal .AND. fft_size_1%lmax == fft_size_2%lmax
3415 equal = equal .AND. fft_size_1%mmax == fft_size_2%mmax
3416 equal = equal .AND. fft_size_1%nmax == fft_size_2%nmax
3418 equal = equal .AND. fft_size_1%mx1 == fft_size_2%mx1
3419 equal = equal .AND. fft_size_1%mx2 == fft_size_2%mx2
3420 equal = equal .AND. fft_size_1%mx3 == fft_size_2%mx3
3422 equal = equal .AND. fft_size_1%my1 == fft_size_2%my1
3423 equal = equal .AND. fft_size_1%my2 == fft_size_2%my2
3424 equal = equal .AND. fft_size_1%my3 == fft_size_2%my3
3426 equal = equal .AND. fft_size_1%mcz1 == fft_size_2%mcz1
3427 equal = equal .AND. fft_size_1%mcx2 == fft_size_2%mcx2
3428 equal = equal .AND. fft_size_1%mcz2 == fft_size_2%mcz2
3429 equal = equal .AND. fft_size_1%mcy3 == fft_size_2%mcy3
3431 equal = equal .AND. fft_size_1%lg == fft_size_2%lg
3432 equal = equal .AND. fft_size_1%mg == fft_size_2%mg
3434 equal = equal .AND. fft_size_1%nbx == fft_size_2%nbx
3435 equal = equal .AND. fft_size_1%nbz == fft_size_2%nbz
3437 equal = equal .AND. fft_size_1%nmray == fft_size_2%nmray
3438 equal = equal .AND. fft_size_1%nyzray == fft_size_2%nyzray
3440 equal = equal .AND. fft_size_1%rs_group == fft_size_2%rs_group
3442 equal = equal .AND. all(fft_size_1%g_pos == fft_size_2%g_pos)
3443 equal = equal .AND. all(fft_size_1%r_pos == fft_size_2%r_pos)
3444 equal = equal .AND. all(fft_size_1%r_dim == fft_size_2%r_dim)
3446 equal = equal .AND. fft_size_1%numtask == fft_size_2%numtask
3448 END SUBROUTINE is_equal
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
subroutine, public fft_3d(plan, scale, zin, zout, stat)
...
subroutine, public fft_destroy_plan(plan)
...
subroutine, public fft_1d(plan, zin, zout, scale, stat)
...
subroutine, public fft_do_cleanup(wisdom_file, ionode, export_wisdom)
...
subroutine, public fft_create_plan_1d(plan, fsign, trans_in, trans_out, ldx_in, ldx_out, n, m, zin, zout)
...
subroutine, public fft_do_init(fftlib, plan_style, wisdom_file)
...
subroutine, public fft_get_lengths(data, max_length)
...
subroutine, public fft_create_plan_3d(plan, fft_in_place, fsign, n, zin, zout)
...
Type to store data about a (1D or 3D) FFT, including FFTW plan.
Defines the basic variable types.
integer, parameter, public dp_size
integer, parameter, public dp
integer, parameter, public sp
Definition of mathematical constants and functions.
complex(kind=dp), parameter, public z_zero
Interface to the message passing library MPI.
type(mp_comm_type), parameter, public mp_comm_null
Fortran API for the offload package, which is written in C.
integer function, public offload_free_pinned_mem(buffer)
free pinned memory
integer function, public offload_malloc_pinned_mem(buffer, length)
allocate pinned memory.