1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
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
21module fourier_space_oct_m
22  use accel_oct_m
23  use cube_oct_m
24  use cube_function_oct_m
25  use global_oct_m
26  use loct_pointer_oct_m
27  use math_oct_m
28  use messages_oct_m
29  use fft_oct_m
30#ifdef HAVE_OPENMP
31  use omp_lib
32#endif
33#ifdef HAVE_PFFT
34  use pfft_oct_m
35#endif
36  use profiling_oct_m
37  use types_oct_m
38
39  implicit none
40  private
41  public ::                     &
42    cube_function_alloc_fs,     &
43    cube_function_free_fs,      &
44    dcube_function_rs2fs,       &
45    zcube_function_rs2fs,       &
46    dcube_function_fs2rs,       &
47    zcube_function_fs2rs,       &
48    fourier_space_op_t,         &
49    fourier_space_op_nullify,   &
50    dfourier_space_op_init,     &
51    dfourier_space_op_apply,    &
52    zfourier_space_op_init,     &
53    zfourier_space_op_apply,    &
54    fourier_space_op_end
55
56  type fourier_space_op_t
57    private
58    FLOAT, pointer :: dop(:, :, :)
59    CMPLX, pointer :: zop(:, :, :)
60    logical :: in_device_memory
61    type(accel_mem_t) :: op_buffer
62    logical :: real_op
63
64    !Parameters used to generate this kernel
65    FLOAT, public :: qq(1:MAX_DIM)
66    FLOAT, public :: singularity
67    FLOAT, public :: mu !< Range separation for the exchange operator
68  end type fourier_space_op_t
69
70contains
71
72  ! ---------------------------------------------------------
73  subroutine fourier_space_op_nullify(this)
74    type(fourier_space_op_t), intent(inout) :: this
75
76    PUSH_SUB(fourier_space_op_end)
77
78    nullify(this%dop)
79    nullify(this%zop)
80    this%in_device_memory = .false.
81
82    !We just set a very large q to guaranty that the kernel is always
83    this%qq = CNST(1e5)
84    this%singularity = M_ZERO
85    this%mu = M_ZERO
86
87    POP_SUB(fourier_space_op_end)
88
89  end subroutine fourier_space_op_nullify
90
91  ! ---------------------------------------------------------
92
93  !> Allocates locally the Fourier space grid, if PFFT library is not used.
94  !! Otherwise, it assigns the PFFT Fourier space grid to the cube Fourier space grid,
95  !! via pointer.
96  subroutine cube_function_alloc_fs(cube, cf, force_alloc)
97    type(cube_t), target,  intent(in)    :: cube
98    type(cube_function_t), intent(inout) :: cf
99    logical, optional,     intent(in)    :: force_alloc
100
101    integer :: n1, n2, n3
102    logical :: allocated
103
104    PUSH_SUB(cube_function_alloc_fs)
105
106    ASSERT(.not. associated(cf%fs))
107    ASSERT(associated(cube%fft))
108
109    cf%forced_alloc = optional_default(force_alloc, .false.)
110
111    n1 = max(1, cube%fs_n(1))
112    n2 = max(1, cube%fs_n(2))
113    n3 = max(1, cube%fs_n(3))
114
115    allocated = .false.
116
117    select case(cube%fft%library)
118    case(FFTLIB_PFFT)
119      if(.not. cf%forced_alloc) then
120        allocated = .true.
121        ASSERT(associated(cube%fft))
122        if(any(cube%fs_n(1:3) == 0)) then
123          cf%fs => cube%fft%fs_data(1:1,1:1,1:1)
124        else
125          cf%fs => cube%fft%fs_data(1:n3,1:n1,1:n2)
126        end if
127      else ! force allocate transposed with PFFT
128        allocated = .true.
129        SAFE_ALLOCATE(cf%fs(1:n3, 1:n1, 1:n2))
130      end if
131    case(FFTLIB_ACCEL)
132      if(cf%in_device_memory) then
133        allocated = .true.
134        call accel_create_buffer(cf%fourier_space_buffer, ACCEL_MEM_READ_WRITE, TYPE_CMPLX, product(cube%fs_n(1:3)))
135      end if
136
137    case(FFTLIB_FFTW)
138      if(.not. cf%forced_alloc) then
139        allocated = .true.
140        ASSERT(associated(cube%fft))
141        cf%fs => cube%fft%fs_data(1:cube%fs_n(1), 1:cube%fs_n(2), 1:cube%fs_n(3))
142      end if
143    end select
144
145    if(.not. allocated) then
146      SAFE_ALLOCATE(cf%fs(1:cube%fs_n(1), 1:cube%fs_n(2), 1:cube%fs_n(3)))
147    end if
148
149    POP_SUB(cube_function_alloc_fs)
150  end subroutine cube_function_alloc_fs
151
152
153  ! ---------------------------------------------------------
154  !> Deallocates the Fourier space grid
155  subroutine cube_function_free_fs(cube, cf)
156    type(cube_t),          intent(in)    :: cube
157    type(cube_function_t), intent(inout) :: cf
158
159    logical :: deallocated
160
161    PUSH_SUB(cube_function_free_fs)
162
163    ASSERT(associated(cube%fft))
164
165    deallocated = .false.
166
167    select case(cube%fft%library)
168    case(FFTLIB_PFFT)
169      if(.not. cf%forced_alloc) then
170        deallocated = .true.
171        nullify(cf%fs)
172      end if
173    case(FFTLIB_ACCEL)
174      if(cf%in_device_memory) then
175        deallocated = .true.
176        call accel_release_buffer(cf%fourier_space_buffer)
177      end if
178    case(FFTLIB_FFTW)
179      if(.not. cf%forced_alloc) then
180        deallocated = .true.
181        nullify(cf%fs)
182      end if
183    end select
184
185    if(.not. deallocated) then
186      ASSERT(associated(cf%fs))
187      SAFE_DEALLOCATE_P(cf%fs)
188    end if
189
190    POP_SUB(cube_function_free_fs)
191  end subroutine cube_function_free_fs
192
193
194  ! ---------------------------------------------------------
195  subroutine fourier_space_op_end(this)
196    type(fourier_space_op_t), intent(inout) :: this
197
198    PUSH_SUB(fourier_space_op_end)
199
200    if(this%in_device_memory) then
201      call accel_release_buffer(this%op_buffer)
202      this%in_device_memory = .false.
203    end if
204    SAFE_DEALLOCATE_P(this%dop)
205    SAFE_DEALLOCATE_P(this%zop)
206
207
208    POP_SUB(fourier_space_op_end)
209  end subroutine fourier_space_op_end
210
211#include "undef.F90"
212#include "real.F90"
213#include "fourier_space_inc.F90"
214
215#include "undef.F90"
216#include "complex.F90"
217#include "fourier_space_inc.F90"
218
219end module fourier_space_oct_m
220
221
222!! Local Variables:
223!! mode: f90
224!! coding: utf-8
225!! End:
226