1!! Copyright (C) 2011 U. De Giovannini 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 20#include "global.h" 21 22module nfft_oct_m 23 use fftw_params_oct_m 24 use global_oct_m 25 use, intrinsic :: iso_c_binding 26 use loct_math_oct_m 27 use messages_oct_m 28 use namespace_oct_m 29 use parser_oct_m 30 use varinfo_oct_m 31 implicit none 32 33 private 34 35 public :: & 36 nfft_t, & 37 nfft_copy_info, & 38 nfft_init, & 39 nfft_end, & 40 nfft_precompute, & 41 nfft_write_info, & 42 nfft_guru_options, & 43 znfft_forward, & 44 znfft_backward, & 45 dnfft_forward, & 46 dnfft_backward 47 48 ! global constants 49 integer, public, parameter :: & 50 nfft_real = 0, & 51 nfft_complex = 1 52 53 !NFFT flags 54 integer, public, parameter :: & 55 NFFT_PRE_PHI_HUT = 0, & 56 NFFT_FG_PSI = 2, & 57 NFFT_PRE_LIN_PSI = 4, & 58 NFFT_PRE_FG_PSI = 8, & 59 NFFT_PRE_PSI = 16, & 60 NFFT_PRE_FULL_PSI = 32, & 61 NFFT_MALLOC_X = 64, & 62 NFFT_MALLOC_F_HAT = 128, & 63 NFFT_MALLOC_F = 256, & 64 NFFT_FFT_OUT_OF_PLACE = 512, & 65 NFFT_FFTW_INIT = 1024 66 67 68 type nfft_t 69 private 70 71 integer :: N(3) !> size of the nfft bandwidths 72 integer :: M(3) !> Number of the nfft nodes 73 integer :: dim !> the dimension 74 integer :: fftN(3) !> size of the fft used 75 FLOAT, public :: norm !> Normalization 76 77 ! Guru options 78 logical, public :: set_defaults = .false. !> the defaults can be overriden 79 logical, public :: guru !> use guru options? 80 integer, public :: precompute !> precompute strategy 81 integer, public :: mm !> Window function cut-off parameter 82 FLOAT, public :: sigma !> Oversampling factor 83 84 type(c_ptr) :: plan !> the NFFT plan 85 86 end type nfft_t 87 88 89contains 90 91 92 93 ! --------------------------------------------------------- 94 ! GURU options 95 subroutine nfft_guru_options(nfft, namespace) 96 type(nfft_t), intent(inout) :: nfft 97 type(namespace_t), intent(in) :: namespace 98 99 PUSH_SUB(nfft_guru_options) 100 101 !%Variable NFFTGuruInterface 102 !%Type logical 103 !%Default false 104 !%Section Mesh::FFTs 105 !%Description 106 !% Perform NFFT with guru interface. This permits the fine tuning of several critical parameters. 107 !%End 108 call parse_variable(namespace, 'NFFTGuruInterface', .false., nfft%guru) 109 110 111 !%Variable NFFTCutoff 112 !%Type integer 113 !%Default 6 114 !%Section Mesh::FFTs 115 !%Description 116 !% Cut-off parameter of the window function. 117 !% See NFFT manual for details. 118 !%End 119 call parse_variable(namespace, 'NFFTCutoff', 6, nfft%mm) 120 121 122 !%Variable NFFTOversampling 123 !%Type float 124 !%Default 2 125 !%Section Mesh::FFTs 126 !%Description 127 !% NFFT oversampling factor (sigma). This will rule the size of the FFT under the hood. 128 !%End 129 call parse_variable(namespace, 'NFFTOversampling', M_TWO, nfft%sigma) 130 131 !%Variable NFFTPrecompute 132 !%Type integer 133 !%Default NFFT_PRE_PSI 134 !%Section Mesh::FFTs 135 !%Description 136 !% NFFT precomputation strategy. 137 !%Option NFFT_PRE_LIN_PSI 4 138 !% This method implements a linear interpolation from a lookup table. 139 !%Option NFFT_PRE_PSI 16 140 !% This method uses a medium amount of memory to store d*(2*m+1)*M real numbers and requires at most 141 !% 2(2m + 1)d extra multiplications for each node. 142 !% This is the default option. 143 !%Option NFFT_PRE_FULL_PSI 32 144 !% Is the fastest method but requires a large amount of memory as it requires to store (2*m+1)^d*M 145 !% real numbers. No extra operations are needed during matrix vector multiplication. 146 !%End 147 call parse_variable(namespace, 'NFFTPrecompute', NFFT_PRE_PSI, nfft%precompute) 148 if(.not.varinfo_valid_option('NFFTPrecompute', nfft%precompute)) call messages_input_error(namespace, 'NFFTPrecompute') 149! call messages_print_var_option(stdout, "NFFTPrecompute", nfft%precompute) 150 151! if(.not.varinfo_valid_option('NFFTPrecompute', nfft%precompute, is_flag=.true.)) then 152! call messages_input_error('NFFTPrecompute') 153! end if 154 155 156 POP_SUB(nfft_guru_options) 157 end subroutine nfft_guru_options 158 159 160 161 ! --------------------------------------------------------- 162 subroutine nfft_init(nfft, nfft_options, N, dim, M, optimize) 163 type(nfft_t), intent(inout) :: nfft 164 type(nfft_t), intent(in) :: nfft_options 165 integer, intent(inout) :: N(3) !> nfft bandwidths 166 integer, intent(inout) :: M(3) !> nfft nodes 167 integer, intent(in) :: dim 168 logical, optional, intent(in) :: optimize 169 170 integer :: ii, my_N(3) 171 logical :: optimize_ 172 integer :: nfft_flags 173 174 175 PUSH_SUB(nfft_init) 176 177 optimize_ = optional_default(optimize, .true.) 178 179 nfft%dim = dim 180 nfft%M(:) = M(:) 181 nfft%N(:) = N(:) 182 183 if(.not. nfft%set_defaults) then 184 nfft%guru = nfft_options%guru 185 nfft%mm = nfft_options%mm 186 nfft%sigma = nfft_options%sigma 187 nfft%precompute = nfft_options%precompute 188 end if 189 190 ! set unused dimensions to 1 191 nfft%M(dim+1:3) = 1 192 193 my_N = 0 194 do ii = 1, dim 195 my_N(ii) = N(ii)*nfft%sigma 196 if(optimize_ .or. (.not. nfft%guru)) call loct_fft_optimize(my_N(ii), 1) ! ask for an odd number 197 end do 198 199 nfft%fftN(1:dim) = my_N(1:dim) 200 201 if(nfft%guru) then ! Guru interface 202 nfft_flags = nfft_PRE_PHI_HUT + nfft_MALLOC_X +nfft_MALLOC_F_HAT +& 203 nfft_MALLOC_F + nfft_FFTW_INIT + nfft_FFT_OUT_OF_PLACE 204 205 nfft_flags = nfft_flags + nfft%precompute 206 207 call oct_nfft_init_guru(nfft%plan, dim, N, M(1)*M(2)*M(3), my_N, nfft%mm, & 208 nfft_flags, FFTW_MEASURE + FFTW_DESTROY_INPUT) 209 210 else ! Default interfaces 211 212 select case(dim) 213 case(3) 214 call oct_nfft_init_3d(nfft%plan, N(1), N(2),N(3), M(1)*M(2)*M(3)) 215 case(2) 216 call oct_nfft_init_2d(nfft%plan, N(1), N(2), M(1)*M(2)) 217 case(1) 218 call oct_nfft_init_1d(nfft%plan,N(1),M(1)) 219 end select 220 221 end if 222 223 224 POP_SUB(nfft_init) 225 end subroutine nfft_init 226 227 ! --------------------------------------------------------- 228 subroutine nfft_write_info(nfft) 229 type(nfft_t), intent(inout) :: nfft 230 231 integer :: idir 232! integer :: mm 233 234 PUSH_SUB(nfft_write_info) 235 236 call messages_write("Info: NFFT parameters") 237 call messages_new_line() 238 call messages_write(" Fourier coefficients N = ") 239 do idir = 1, nfft%dim 240 call messages_write(nfft%N(idir)) 241 if(idir < nfft%dim) call messages_write(" x ") 242 end do 243 call messages_new_line() 244 245 call messages_write(" Spatial nodes M = ") 246 247! mm = nfft%M(1)*nfft%M(2)*nfft%M(3) 248! 249! call messages_write(mm) 250! call messages_new_line() 251 do idir = 1, nfft%dim 252 call messages_write(nfft%M(idir)) 253 if(idir < nfft%dim) call messages_write(" x ") 254 end do 255 call messages_new_line() 256 257 258 call messages_write(" Oversampling factor sigma = ") 259 call messages_write(nfft%sigma) 260 call messages_new_line() 261 262 call messages_write(" FFT grid size n = ") 263 do idir = 1, nfft%dim 264 call messages_write(nfft%fftN(idir)) 265 if(idir < nfft%dim) call messages_write(" x ") 266 end do 267 call messages_new_line() 268 269 call messages_write(" Real space cutoff m = ") 270 call messages_write(nfft%mm) 271 call messages_new_line() 272 273 call messages_write(" Pre-computation strategy = ") 274 select case(nfft%precompute) 275 case(nfft_PRE_LIN_PSI) 276 call messages_write(" NFFT_PRE_LIN_PSI") 277 case(nfft_PRE_PSI) 278 call messages_write(" NFFT_PRE_PSI") 279 case(nfft_PRE_FULL_PSI) 280 call messages_write(" NFFT_PRE_FULL_PSI") 281 end select 282 283 call messages_info() 284 285 286 POP_SUB(nfft_write_info) 287 end subroutine nfft_write_info 288 289 290 ! --------------------------------------------------------- 291 subroutine nfft_end(nfft) 292 type(nfft_t), intent(inout) :: nfft 293 294 PUSH_SUB(nfft_end) 295 296 call oct_nfft_finalize(nfft%plan); 297 298 POP_SUB(nfft_end) 299 end subroutine nfft_end 300 301 ! --------------------------------------------------------- 302 ! This routine is intend to copy the configuration parameters 303 ! rather the whole structure. 304 ! --------------------------------------------------------- 305 subroutine nfft_copy_info(in, out) 306 type(nfft_t), intent(in) :: in 307 type(nfft_t), intent(out) :: out 308 309 310 PUSH_SUB(nfft_copy_info) 311 312 out%N = in%N 313 out%M = in%M 314 out%dim = in%dim 315 out%fftN = in%fftN 316 out%norm = in%norm 317 318 out%guru = in%guru 319 out%precompute = in%precompute 320 out%mm = in%mm 321 out%sigma = in%sigma 322 323 324 POP_SUB(nfft_copy_info) 325 end subroutine nfft_copy_info 326 327 328 !---------------------------------------------------------- 329 ! Precompute the plan according to the position the grid nodes in real space 330 ! x axis is X1, y axis is X2, z axis is X3 331 ! NOTE: We only allow different spacing for each direction x,y,z 332 ! the NFFT interface however is more general 333 ! --------------------------------------------------------- 334 subroutine nfft_precompute(nfft, X1, X2, X3) 335 type(nfft_t), intent(inout) :: nfft 336 FLOAT, intent(in) :: X1(:) 337 FLOAT, optional, intent(in) :: X2(:) 338 FLOAT, optional, intent(in) :: X3(:) 339 340 341 342 FLOAT :: x1_(1:nfft%M(1)), x2_(1:nfft%M(2)), x3_(1:nfft%M(3)) 343 FLOAT :: length, cc, eps, dX1(1:nfft%M(1)-1), dX2(1:nfft%M(2)-1), dX3(1:nfft%M(3)-1) 344 345 integer :: ii 346 347 PUSH_SUB(nfft_precompute) 348 349! eps = 1.000001 ! the sample nodes must be in [0.5,0.5) 350 eps = M_ONE+M_EPSILON ! the sample nodes must be in [0.5,0.5) 351 352 select case(nfft%dim) 353 case(3) 354 length = (maxval(X1)-minval(X1))*eps 355 cc = (minval(X1)+maxval(X1))/M_TWO 356 x1_ =(X1-cc)/length 357 length = (maxval(X2)-minval(X2))*eps 358 cc = (minval(X2)+maxval(X2))/M_TWO 359 x2_ =(X2-cc)/length 360 length = (maxval(X3)-minval(X3))*eps 361 cc = (minval(X3)+maxval(X3))/M_TWO 362 x3_ =(X3-cc)/length 363 call oct_nfft_precompute_one_psi_3d(nfft%plan, nfft%M, x1_, x2_, x3_) 364 365 ! Set the normalization factor 366 do ii = 1, nfft%M(1)-1 367 dX1(ii)= abs(x1_(ii+1)-x1_(ii)) 368 end do 369 do ii = 1, nfft%M(2)-1 370 dX2(ii)= abs(x2_(ii+1)-x2_(ii)) 371 end do 372 do ii = 1, nfft%M(3)-1 373 dX3(ii)= abs(x3_(ii+1)-x3_(ii)) 374 end do 375 nfft%norm = M_ONE/(minval(dX1(:)) * minval(dX2(:)) * minval(dX3(:))) 376 377 case(2) 378 length = (maxval(X1)-minval(X1))*eps 379 cc = (minval(X1)+maxval(X1))/M_TWO 380 x1_ =(X1-cc)/length 381 length = (maxval(X2)-minval(X2))*eps 382 cc = (minval(X2)+maxval(X2))/M_TWO 383 x2_ =(X2-cc)/length 384 call oct_nfft_precompute_one_psi_2d(nfft%plan, nfft%M, x1_, x2_) 385 386 ! Set the normalization factor 387 do ii = 1, nfft%M(1)-1 388 dX1(ii)= abs(x1_(ii+1)-x1_(ii)) 389 end do 390 do ii = 1, nfft%M(2)-1 391 dX2(ii)= abs(x2_(ii+1)-x2_(ii)) 392 end do 393 nfft%norm = M_ONE/(minval(dX1(:)) * minval(dX2(:))) 394 395 396 case(1) 397 length = (maxval(X1)-minval(X1))*eps 398 cc = (minval(X1)+maxval(X1))/M_TWO 399 x1_ =(X1-cc)/length 400 call oct_nfft_precompute_one_psi_1d(nfft%plan,nfft%M(1),x1_) 401 402 ! Set the normalization factor 403 do ii = 1, nfft%M(1)-1 404 dX1(ii)= abs(x1_(ii+1)-x1_(ii)) 405 end do 406 nfft%norm = M_ONE/(minval(dX1(:))) 407 408 end select 409 410 ! check the plan 411 call oct_nfft_check(nfft%plan) 412 413 414 write(message(1), '(a)') "Info: NFFT plan precomputed." 415 call messages_info(1) 416 417 418 POP_SUB(nfft_precompute) 419 end subroutine nfft_precompute 420 421 !-------------------------------------------- 422 subroutine dnfft_forward(nfft, in, out) 423 type(nfft_t), intent(in) :: nfft 424 FLOAT, intent(in) :: in(:,:,:) 425 CMPLX, intent(out) :: out(:,:,:) 426 427 CMPLX, allocatable :: zin(:,:,:) 428 integer:: b(6) 429 430 PUSH_SUB(dnfft_forward) 431 432 b(1) = lbound(in, dim=1) 433 b(2) = ubound(in, dim=1) 434 b(3) = lbound(in, dim=2) 435 b(4) = ubound(in, dim=3) 436 b(5) = lbound(in, dim=3) 437 b(6) = ubound(in, dim=3) 438 439! SAxFE_ALLOCATE(zin(b(1):b(2),b(3):b(4),b(5):b(6))) 440 allocate(zin(b(1):b(2),b(3):b(4),b(5):b(6))) 441 zin = in 442 call znfft_forward(nfft, zin, out) 443 444 deallocate(zin) 445! SAxFE_DEALLOCATE_A(zin) 446 447 POP_SUB(dnfft_forward) 448 end subroutine dnfft_forward 449 450 !-------------------------------------------- 451 subroutine dnfft_backward(nfft, in, out) 452 type(nfft_t), intent(in) :: nfft 453 CMPLX, intent(in) :: in (:,:,:) 454 FLOAT, intent(out) :: out(:,:,:) 455 456 CMPLX, allocatable :: zout(:,:,:) 457 integer:: b(6) 458 459 PUSH_SUB(dnfft_backward) 460 461 b(1) = lbound(out, dim=1) 462 b(2) = ubound(out, dim=1) 463 b(3) = lbound(out, dim=2) 464 b(4) = ubound(out, dim=3) 465 b(5) = lbound(out, dim=3) 466 b(6) = ubound(out, dim=3) 467 468 allocate(zout(b(1):b(2),b(3):b(4),b(5):b(6))) 469 470 call znfft_backward(nfft, in, zout) 471 out = zout 472 deallocate(zout) 473 474 POP_SUB(dnfft_backward) 475 end subroutine dnfft_backward 476 477 !-------------------------------------------- 478 subroutine znfft_forward(nfft, in, out) 479 type(nfft_t), intent(in) :: nfft 480 CMPLX, intent(in) :: in(:,:,:) 481 CMPLX, intent(out) :: out(:,:,:) 482 483 integer :: ix, iy, iz 484 485 PUSH_SUB(znfft_forward) 486 487 do ix = 1, nfft%N(1) 488 do iy = 1, nfft%N(2) 489 do iz = 1, nfft%N(3) 490 call zoct_set_f_hat(nfft%plan, nfft%dim, in(ix,iy,iz), ix, iy, iz) 491 end do 492 end do 493 end do 494 495 call oct_nfft_trafo(nfft%plan) 496 497 do ix = 1, nfft%M(1) 498 do iy = 1, nfft%M(2) 499 do iz = 1, nfft%M(3) 500 call zoct_get_f(nfft%plan, nfft%M, nfft%dim, out(ix,iy,iz), ix, iy, iz) 501 end do 502 end do 503 end do 504 505 POP_SUB(znfft_forward) 506 end subroutine znfft_forward 507 508 ! --------------------------------------------------------- 509 subroutine znfft_backward(nfft, in, out) 510 type(nfft_t), intent(in) :: nfft 511 CMPLX, intent(in) :: in (:,:,:) 512 CMPLX, intent(out) :: out(:,:,:) 513 514 integer :: ix, iy, iz 515 516 PUSH_SUB(znfft_backward) 517 518 do ix = 1, nfft%M(1) 519 do iy = 1, nfft%M(2) 520 do iz = 1, nfft%M(3) 521 call zoct_set_f(nfft%plan, nfft%M, nfft%dim, in(ix,iy,iz), ix, iy, iz) 522 end do 523 end do 524 end do 525 526 call oct_nfft_adjoint(nfft%plan) 527 528 do ix = 1,nfft%N(1) 529 do iy = 1, nfft%N(2) 530 do iz = 1, nfft%N(3) 531 call zoct_get_f_hat(nfft%plan, nfft%dim, out(ix,iy,iz), ix, iy, iz) 532 end do 533 end do 534 end do 535 536 out = out/nfft%norm 537 538 POP_SUB(znfft_backward) 539 end subroutine znfft_backward 540 541end module nfft_oct_m 542 543!! Local Variables: 544!! mode: f90 545!! coding: utf-8 546!! End: 547