FFTW and Subroutines, copying arrays
in [Fortran]
|
From: Ralph Kube on 4 Nov 2009 12:06 Richard Maine wrote: > Ideal, since you are talking about a run-time error, would be a small, > self-contained, runable example. Lacking that, it at least needs more > information than this. > Okay, here we go. My problem was that I got different results from the routines csfield1 and csfield2. I found the error by now, I was overwriting a variable along the way in the real program. Thank you everybody for your debugging help, I am going home now. Cheers, Ralph This is a stripped down version of the code that was bugging me. subroutine csfield1(My, Nx, Ahat, A, planbw) integer My, Nx double complex Ahat(0:My/2+1, 0:Nx+1) double precision A(0:My+1, 0:Nx+1) integer*8 planbw double complex carr(0:My/2+1, 0:Nx+1) double precision rarr(0:My+1, 0:Nx+1) call dfftw_execute_dft_c2r(planbw, Ahat, A) end subroutine subroutine csfield2(My, Nx, Ahat, A, planbw) integer My, Nx double complex Ahat(0:My/2+1, 0:Nx+1) double precision A(0:My+1, 0:Nx+1) integer*8 planbw double complex carr(0:My/2+1, 0:Nx+1) double precision rarr(0:My+1, 0:Nx+1) carr = Ahat call dfftw_execute_dft_c2r(planbw, carr, rarr) Ahat = rarr end subroutine program fftwtest integer, parameter :: Nx = 8 integer, parameter :: My = 8 integer, parameter :: tlevs = 3 real(kind=8), allocatable :: theta(:,:) double complex, allocatable :: thetahat(:,:,:) integer*8 planfw, planbw integer i,j real(kind=8) pi, Lx, Ly allocate(theta(0:My+1, 0:Nx+1)) allocate(thetahat(0:My/2+1, 0:Nx+1, 0:tlevs-1)) print*, 'planning fftw' call dfftw_plan_dft_r2c_2d(planfw, My+2, Nx+2, theta, thetahat(0,0,0), FFTW_UNALIGNED + FFTW_ESTIMATE) call dfftw_plan_dft_c2r_2d(planbw, My+2, Nx+2, thetahat(0,0,0), theta, FFTW_UNALIGNED + FFTW_ESTIMATE) pi = acos(-1.0) Lx = real(Nx, 8) Ly = real(My, 8) do i = 0, My+1 do j = 0, Nx+1 theta(i,j) = sin(2*pi*j/Ly)*sin(2*pi*i/Lx) end do end do call dfftw_execute_dft_r2c(planfw, theta, thetahat(0,0,0)) print*, 'calling subroutine' print*, theta call csfield1(My, Nx, thetahat(0,0,0), theta, planbw) print*, '-------------------------' print*, theta/(real((My+1)*(Nx+1), 8)) call csfield2(My, Nx, thetahat(0,0,0), theta, planbw) print*, '-------------------------' print*, theta/(real((My+1)*(Nx+1), 8)) call dfftw_destroy_plan(planfw) call dfftw_destroy_plan(planbw) end program
From: Steven G. Johnson on 4 Nov 2009 14:52 On Nov 4, 9:44 am, Ralph Kube <rku...(a)post.uit.no> wrote: > dpb wrote: > > In the code as shown here, sure... > > > call dfftw_execute_dft_r2c(planfw, akw, thrhshat) > > This is the call to thefftwsubroutine that gives me a headache. When > I pass an array defined in the same subroutine to it, it works okay. > But when I pass an array, that has been passed to the subroutine where > I call dfftw_execute_dft_r2c, it segfaults. > > Any ideas what causes this behaviour and how I can avoid declaring > an extra array? One thing you might look at is the fine manual. It has a section explicitly discussing the dfftw_execute* functions in Fortran: http://fftw.org/doc/FFTW-Execution-in-Fortran.html In particular, it says: If you pass different input/output arrays compared to those used when creating the plan, you must abide by all the restrictions of the new-array execute functions (see New-array Execute Functions). The most difficult of these, in Fortran, is the requirement that the new arrays have the same alignment as the original arrays, because there seems to be no way in Fortran to obtain guaranteed-aligned arrays (analogous to fftw_malloc in C). You can, of course, use the FFTW_UNALIGNED flag when creating the plan, in which case the plan does not depend on the alignment, but this may sacrifice substantial performance on architectures (like x86) with SIMD instructions (see SIMD alignment and fftw_malloc). Regards, Steven G. Johnson
From: Steven G. Johnson on 4 Nov 2009 15:16 Ah, nevermind, I see by your posted code that you are indeed using FFTW_UNALIGNED; I'm glad you found your bug elsewhere. Steven
From: robin on 5 Nov 2009 10:47 "Ralph Kube" <rku000(a)post.uit.no> wrote in message news:hcs40n$2np0$1(a)inn.uit.no... | dpb wrote: | > In the code as shown here, sure... | > | > call dfftw_execute_dft_r2c(planfw, akw, thrhshat) | | This is the call to the fftw subroutine that gives me a headache. When | I pass an array defined in the same subroutine to it, it works okay. | But when I pass an array, that has been passed to the subroutine where | I call dfftw_execute_dft_r2c, it segfaults. | | Any ideas what causes this behaviour and how I can avoid declaring | an extra array? Do you have an interface for it? Do you have subscript bounds checking enabled?
From: robin on 5 Nov 2009 10:47
"Ralph Kube" <rku000(a)post.uit.no> wrote in message news:hcromn$2h4h$1(a)inn.uit.no... | Hello everybody, | I am using FFTW with fortran and wonder about the way I have to pass | arrays to subroutines. In every subroutine I use FFTW in, I have | to locally declare the array fftw stores the array in. If I use an array | passed to the subroutine as the target for the fftw routine, my | program crashes. | Here is an example: This isn't an example. It's not even a complete subroutine. What is resol.h? Anyway, why not just define theta etc as double precision theta (0: , 0: ) ? | subroutine thetaeq(theta, strmf, thrhshat, planfw) | include 'resol.h' | double precision theta(0:My+1, 0:Nx+1) | double precision strmf(0:My+1, 0:Nx+1) | double complex thrhshat(0:My/2+1, 0:Nx+1) | integer*8 planfw | | double precision akw(0:My+1, 0:Nx+1) | double complex akwhat(0:My/2+1, 0:Nx+1) | | call arakw(theta, strmf, akw) | call dfftw_execute_dft_r2c(planfw, akw, akwhat) | thrhshat = akwhat | end subroutine | | The routine dfftw_execute_dft_r2c stores the result in akwhat. Is there | any way I can store the result in thrhshat so that I can ommit copying | the array? |