1!! Copyright (C) 2011 J. Alberdi-Rodriguez, P. Garcia Risueño, M. Oliveira
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
21!> The includes for the PFFT
22module pfft_params_oct_m
23  use, intrinsic :: iso_c_binding
24  use fftw_params_oct_m
25  implicit none
26
27  private
28
29#ifdef HAVE_PFFT
30
31  ! make public some symbols from the pfft library
32  public :: pfft_cleanup,               &
33            pfft_create_procmesh_2d,    &
34            pfft_destroy_plan,          &
35            pfft_execute,               &
36            pfft_init,                  &
37            pfft_local_size_dft_r2c_3d, &
38            pfft_local_size_dft_3d,     &
39            pfft_plan_dft_c2r_3d,       &
40            pfft_plan_dft_r2c_3d,       &
41            pfft_plan_dft_3d,           &
42            PFFT_INT,                   &
43            PFFT_PTRDIFF_T,             &
44            PFFT_FLOAT,                 &
45            PFFT_DOUBLE,                &
46            PFFT_LDOUBLE,               &
47            PFFT_UNSIGNED,              &
48            PFFT_TRANSPOSED_IN,         &
49            PFFT_TRANSPOSED_OUT,        &
50            PFFT_ESTIMATE,              &
51            PFFT_TUNE,                  &
52            PFFT_DESTROY_INPUT,         &
53            PFFT_MEASURE,               &
54            PFFT_FORWARD,               &
55            PFFT_BACKWARD
56
57  include "pfft.f03"
58#endif
59
60end module pfft_params_oct_m
61
62!> The low level module to work with the PFFT library.
63!! http://www-user.tu-chemnitz.de/~mpip/software.php?lang=en
64module pfft_oct_m
65  use global_oct_m
66  use, intrinsic :: iso_c_binding
67  use math_oct_m
68  use messages_oct_m
69  use pfft_params_oct_m
70  use profiling_oct_m
71
72  implicit none
73
74  private
75
76  public ::                    &
77    pfft_decompose,            &
78    pfft_prepare_plan_r2c,     &
79    pfft_prepare_plan_c2r,     &
80    pfft_prepare_plan_c2c,     &
81    pfft_get_dims
82
83contains
84
85  ! ---------------------------------------------------------
86  !> Decompose all available processors in 2D processor grid, most equally possible
87  subroutine pfft_decompose(n_proc, dim1, dim2)
88    integer, intent(in)  :: n_proc !< Number of processors
89    integer, intent(out) :: dim1   !< First out dimension
90    integer, intent(out) :: dim2   !< Second out dimension
91
92    integer :: np, i
93
94    PUSH_SUB(pfft_decompose)
95
96    ASSERT(n_proc > 0)
97
98    dim1 = 1
99    dim2 = 1
100    np = n_proc
101    i = n_proc-1
102
103    if (is_prime(n_proc)) then
104      dim1 = n_proc
105    else
106      do
107        if (i <= 1) exit
108        if(mod(np,i) /= 0.or.(.not. is_prime(i))) then
109          i=i-1
110          cycle
111        end if
112        np=np/i
113        if (dim1 <= dim2) then
114          dim1 = dim1*i
115        else
116          dim2 = dim2*i
117        end if
118      end do
119    end if
120
121    ASSERT(dim1*dim2 == n_proc)
122
123    POP_SUB(pfft_decompose)
124  end subroutine pfft_decompose
125
126  ! ---------------------------------------------------------
127  !> Octopus subroutine to prepare a PFFT plan real to complex
128  subroutine pfft_prepare_plan_r2c(plan, n, in, out, sign, flags, mpi_comm)
129    type(C_PTR),      intent(out)   :: plan       !< The plan that is created by PFFT
130    integer,          intent(in)    :: n(:)       !< The size of the global matrix
131    FLOAT,   pointer, intent(inout) :: in(:,:,:)  !< The input matrix that is going to be used to do the transform
132    CMPLX,   pointer, intent(inout) :: out(:,:,:) !< The output matrix that is going to be used to do the transform
133    integer,          intent(in)    :: sign       !< Sign flag to decide FFT direction. Has to be FFTW_FORWARD
134    integer,          intent(in)    :: flags      !< Flags for FFT library. Could be changed with the input variable
135                                                  !! FFTPreparePlan. Default value is FFTW_MEASURE
136    integer,          intent(in)    :: mpi_comm   !< MPI communicator
137
138    integer(C_INTPTR_T) :: tmp_n(3)
139    type(profile_t), save   :: prof
140
141    PUSH_SUB(pfft_prepare_plan_r2c)
142    call profiling_in(prof,"PFFT_PLAN_R2C")
143
144#ifdef HAVE_PFFT
145    ASSERT(sign == PFFT_FORWARD)
146#endif
147
148    tmp_n(1:3) = n(3:1:-1)
149
150#ifdef HAVE_PFFT
151    plan = pfft_plan_dft_r2c_3d(tmp_n,in,out,mpi_comm,sign,&
152      PFFT_TRANSPOSED_OUT + PFFT_MEASURE + PFFT_DESTROY_INPUT + PFFT_TUNE)
153#else
154    plan = C_NULL_PTR
155#endif
156
157    call profiling_out(prof)
158    POP_SUB(pfft_prepare_plan_r2c)
159  end subroutine pfft_prepare_plan_r2c
160
161  ! ---------------------------------------------------------
162  !> Octopus subroutine to prepare a PFFT plan real to complex
163  subroutine pfft_prepare_plan_c2r(plan, n, in, out, sign, flags, mpi_comm)
164    type(C_PTR),      intent(out)   :: plan       !< The plan that is created by PFFT
165    integer,          intent(in)    :: n(:)       !< The size of the global matrix
166    CMPLX,   pointer, intent(inout) :: in(:,:,:)  !< The input matrix that is going to be used to do the transform
167    FLOAT,   pointer, intent(inout) :: out(:,:,:) !< The output matrix that is going to be used to do the transform
168    integer,          intent(in)    :: sign       !< Sign flag to decide FFT direction. Has to be FFTW_BACKWARD
169    integer,          intent(in)    :: flags      !< Flags for FFT library. Could be changed with the input variable
170                                                  !! FFTPreparePlan. Default value is FFTW_MEASURE
171    integer,          intent(in)    :: mpi_comm   !< MPI communicator
172
173    integer(C_INTPTR_T) :: tmp_n(3)
174    type(profile_t), save   :: prof
175    PUSH_SUB(pfft_prepare_plan_c2r)
176    call profiling_in(prof,"PFFT_PLAN_C2R")
177
178#ifdef HAVE_PFFT
179    ASSERT(sign == PFFT_BACKWARD)
180#endif
181
182    tmp_n(1:3) = n(3:1:-1)
183
184#ifdef HAVE_PFFT
185    plan = pfft_plan_dft_c2r_3d(tmp_n,in,out,mpi_comm, sign, &
186      PFFT_TRANSPOSED_IN + PFFT_MEASURE + PFFT_DESTROY_INPUT + PFFT_TUNE)
187    !call PDFFT(plan_dft_c2r_3d) (plan, tmp_n(1), in(1,1,1), out(1,1,1), mpi_comm, sign, &
188    !     PFFT_TRANSPOSED_IN + PFFT_MEASURE + PFFT_DESTROY_INPUT + PFFT_TUNE)
189#else
190    plan = C_NULL_PTR
191#endif
192
193    call profiling_out(prof)
194    POP_SUB(pfft_prepare_plan_c2r)
195  end subroutine pfft_prepare_plan_c2r
196
197  ! ---------------------------------------------------------
198  !> Octopus subroutine to prepare a PFFT plan real to complex
199  subroutine pfft_prepare_plan_c2c(plan, n, in, out, sign, flags, mpi_comm)
200    type(C_PTR),      intent(out)   :: plan       !< The plan that is created by PFFT
201    integer,          intent(in)    :: n(:)       !< The size of the global matrix
202    CMPLX,   pointer, intent(inout) :: in(:,:,:)  !< The input matrix that is going to be used to do the transform
203    CMPLX,   pointer, intent(inout) :: out(:,:,:) !< The output matrix that is going to be used to do the transform
204    integer,          intent(in)    :: sign       !< Sign flag to decide FFT direction.
205                                                         !! Has to be FFTW_FORWARD or FFTW_BACKWARD
206    integer,          intent(in)    :: flags      !< Flags for FFT library. Could be changed with the input variable
207                                                  !! FFTPreparePlan. Default value is FFTW_MEASURE
208    integer,          intent(in)    :: mpi_comm   !< MPI communicator
209
210#ifdef HAVE_PFFT
211    integer(C_INT) :: pfft_flags
212#endif
213    integer(C_INTPTR_T) :: tmp_n(3)
214
215    PUSH_SUB(pfft_prepare_plan_c2c)
216
217    tmp_n(1:3) = n(3:1:-1)
218
219#ifdef HAVE_PFFT
220    if (sign == PFFT_FORWARD) then
221      pfft_flags = PFFT_TRANSPOSED_OUT
222    else
223      pfft_flags = PFFT_TRANSPOSED_IN
224    end if
225
226    plan = pfft_plan_dft_3d(tmp_n,in,out,mpi_comm, sign, pfft_flags)
227    !call PDFFT(plan_dft_3d) (plan, tmp_n(1), in(1,1,1), out(1,1,1), mpi_comm, sign, pfft_flags)
228#else
229    plan = C_NULL_PTR
230#endif
231
232    POP_SUB(pfft_prepare_plan_c2c)
233  end subroutine pfft_prepare_plan_c2c
234
235  ! ---------------------------------------------------------
236  subroutine pfft_get_dims(rs_n_global, mpi_comm, is_real, alloc_size, fs_n_global, rs_n, fs_n, rs_istart, fs_istart)
237    integer, intent(in)  :: rs_n_global(1:3) !< The general number of elements in each dimension in real space
238    integer, intent(in)  :: mpi_comm         !< MPI comunicator
239    logical, intent(in)  :: is_real          !< The input is real or complex
240    integer, intent(out) :: alloc_size       !< Number of elements that has to be allocated locally
241    integer, intent(out) :: fs_n_global(1:3) !< The general number of elements in each dimension in Fourier space
242    integer, intent(out) :: rs_n(1:3)        !< Local number of elements in each direction in real space
243    integer, intent(out) :: fs_n(1:3)        !< Local number of elements in each direction in Fourier space
244    integer, intent(out) :: rs_istart(1:3)   !< Where does the local portion of the function start in real space
245    integer, intent(out) :: fs_istart(1:3)   !< Where does the local portion of the function start in Fourier space
246
247    integer(C_INTPTR_T) :: tmp_alloc_size, tmp_n(3)
248    integer(C_INTPTR_T) :: tmp_rs_n(3), tmp_rs_istart(3)
249    integer(C_INTPTR_T) :: tmp_fs_n(3), tmp_fs_istart(3)
250
251    PUSH_SUB(pfft_get_dims)
252
253    fs_n_global(1:3) = rs_n_global(1:3)
254    ! if calling the ISO_C_BINDING routines from PFFT, one must take into account that
255    ! the array order is the C array order
256    tmp_n(1:3) = rs_n_global(3:1:-1)
257
258    tmp_alloc_size = 0
259    if (is_real) then
260#ifdef HAVE_PFFT
261       tmp_alloc_size = pfft_local_size_dft_r2c_3d(tmp_n, mpi_comm, PFFT_TRANSPOSED_OUT, &
262         tmp_rs_n, tmp_rs_istart, tmp_fs_n, tmp_fs_istart)
263#endif
264       !call PDFFT(local_size_dft_r2c_3d)(tmp_alloc_size, tmp_n(1), mpi_comm, PFFT_TRANSPOSED_OUT, &
265       !       tmp_rs_n(1), tmp_rs_istart(1), tmp_fs_n(1), tmp_fs_istart(1))
266       fs_n_global(1) = rs_n_global(1)/2 + 1
267     else
268#ifdef HAVE_PFFT
269       tmp_alloc_size = pfft_local_size_dft_3d(tmp_n,mpi_comm, PFFT_TRANSPOSED_OUT, &
270         tmp_rs_n,tmp_rs_istart,tmp_fs_n,tmp_fs_istart)
271#endif
272       !call PDFFT(local_size_dft_3d)(tmp_alloc_size, tmp_n(1), mpi_comm, PFFT_TRANSPOSED_OUT, &
273       !             tmp_rs_n(1), tmp_rs_istart(1), tmp_fs_n(1), tmp_fs_istart(1))
274    end if
275
276    alloc_size       = tmp_alloc_size
277    rs_n(1:3)        = tmp_rs_n(3:1:-1)
278    fs_n(1:3)        = tmp_fs_n(3:1:-1)
279    rs_istart(1:3)   = tmp_rs_istart(3:1:-1)+1
280    fs_istart(1:3)   = tmp_fs_istart(3:1:-1)+1
281
282    POP_SUB(pfft_get_dims)
283  end subroutine pfft_get_dims
284
285end module pfft_oct_m
286
287!! Local Variables:
288!! mode: f90
289!! coding: utf-8
290!! End:
291