| 1 | ! |
|---|
| 2 | ! $Id: mod_fft_fftw.F90 1403 2010-07-01 09:02:53Z fairhead $ |
|---|
| 3 | ! |
|---|
| 4 | |
|---|
| 5 | MODULE mod_fft_fftw |
|---|
| 6 | |
|---|
| 7 | #ifdef FFT_FFTW |
|---|
| 8 | |
|---|
| 9 | REAL, SAVE :: scale_factor |
|---|
| 10 | INTEGER, SAVE :: vsize |
|---|
| 11 | INTEGER, PARAMETER :: inc=1 |
|---|
| 12 | |
|---|
| 13 | INTEGER*8, ALLOCATABLE, DIMENSION(:), SAVE :: plan_forward |
|---|
| 14 | INTEGER*8, ALLOCATABLE, DIMENSION(:), SAVE :: plan_backward |
|---|
| 15 | |
|---|
| 16 | CONTAINS |
|---|
| 17 | |
|---|
| 18 | SUBROUTINE Init_fft(iim,nvectmax) |
|---|
| 19 | IMPLICIT NONE |
|---|
| 20 | #include <fftw3.f> |
|---|
| 21 | INTEGER :: iim |
|---|
| 22 | INTEGER :: nvectmax |
|---|
| 23 | |
|---|
| 24 | INTEGER :: itmp |
|---|
| 25 | |
|---|
| 26 | INTEGER :: rank |
|---|
| 27 | INTEGER :: howmany |
|---|
| 28 | INTEGER :: istride, idist |
|---|
| 29 | INTEGER :: ostride, odist |
|---|
| 30 | INTEGER, DIMENSION(1) :: n_array, inembed, onembed |
|---|
| 31 | |
|---|
| 32 | REAL, DIMENSION(iim+1,nvectmax) :: dbidon |
|---|
| 33 | COMPLEX, DIMENSION(iim/2+1,nvectmax) :: cbidon |
|---|
| 34 | |
|---|
| 35 | vsize = iim |
|---|
| 36 | scale_factor = 1./SQRT(1.*vsize) |
|---|
| 37 | |
|---|
| 38 | dbidon = 0 |
|---|
| 39 | cbidon = 0 |
|---|
| 40 | |
|---|
| 41 | ALLOCATE(plan_forward(nvectmax)) |
|---|
| 42 | ALLOCATE(plan_backward(nvectmax)) |
|---|
| 43 | |
|---|
| 44 | WRITE(*,*)"!---------------------!" |
|---|
| 45 | WRITE(*,*)"! !" |
|---|
| 46 | WRITE(*,*)"! INITIALISATION FFTW !" |
|---|
| 47 | WRITE(*,*)"! !" |
|---|
| 48 | WRITE(*,*)"!---------------------!" |
|---|
| 49 | |
|---|
| 50 | ! On initialise tous les plans |
|---|
| 51 | DO itmp = 1, nvectmax |
|---|
| 52 | rank = 1 |
|---|
| 53 | n_array(1) = iim |
|---|
| 54 | howmany = itmp |
|---|
| 55 | inembed(1) = iim + 1 ; onembed(1) = iim/2 + 1 |
|---|
| 56 | istride = 1 ; ostride = 1 |
|---|
| 57 | idist = iim + 1 ; odist = iim/2 + 1 |
|---|
| 58 | |
|---|
| 59 | CALL dfftw_plan_many_dft_r2c(plan_forward(itmp), rank, n_array, howmany, & |
|---|
| 60 | & dbidon, inembed, istride, idist, & |
|---|
| 61 | & cbidon, onembed, ostride, odist, & |
|---|
| 62 | & FFTW_ESTIMATE) |
|---|
| 63 | |
|---|
| 64 | rank = 1 |
|---|
| 65 | n_array(1) = iim |
|---|
| 66 | howmany = itmp |
|---|
| 67 | inembed(1) = iim/2 + 1 ; onembed(1) = iim + 1 |
|---|
| 68 | istride = 1 ; ostride = 1 |
|---|
| 69 | idist = iim/2 + 1 ; odist = iim + 1 |
|---|
| 70 | CALL dfftw_plan_many_dft_c2r(plan_backward(itmp), rank, n_array, howmany, & |
|---|
| 71 | & cbidon, inembed, istride, idist, & |
|---|
| 72 | & dbidon, onembed, ostride, odist, & |
|---|
| 73 | & FFTW_ESTIMATE) |
|---|
| 74 | |
|---|
| 75 | ENDDO |
|---|
| 76 | |
|---|
| 77 | WRITE(*,*)"!-------------------------!" |
|---|
| 78 | WRITE(*,*)"! !" |
|---|
| 79 | WRITE(*,*)"! FIN INITIALISATION FFTW !" |
|---|
| 80 | WRITE(*,*)"! !" |
|---|
| 81 | WRITE(*,*)"!-------------------------!" |
|---|
| 82 | |
|---|
| 83 | END SUBROUTINE Init_fft |
|---|
| 84 | |
|---|
| 85 | |
|---|
| 86 | SUBROUTINE fft_forward(vect,TF_vect,nb_vect) |
|---|
| 87 | IMPLICIT NONE |
|---|
| 88 | #include <fftw3.f> |
|---|
| 89 | INTEGER,INTENT(IN) :: nb_vect |
|---|
| 90 | REAL,INTENT(IN) :: vect(vsize+inc,nb_vect) |
|---|
| 91 | COMPLEX,INTENT(OUT) :: TF_vect(vsize/2+1,nb_vect) |
|---|
| 92 | |
|---|
| 93 | CALL dfftw_execute_dft_r2c(plan_forward(nb_vect),vect,TF_vect) |
|---|
| 94 | |
|---|
| 95 | TF_vect = scale_factor * TF_vect |
|---|
| 96 | |
|---|
| 97 | END SUBROUTINE fft_forward |
|---|
| 98 | |
|---|
| 99 | SUBROUTINE fft_backward(TF_vect,vect,nb_vect) |
|---|
| 100 | IMPLICIT NONE |
|---|
| 101 | #include <fftw3.f> |
|---|
| 102 | INTEGER,INTENT(IN) :: nb_vect |
|---|
| 103 | REAL,INTENT(OUT) :: vect(vsize+inc,nb_vect) |
|---|
| 104 | COMPLEX,INTENT(IN ) :: TF_vect(vsize/2+1,nb_vect) |
|---|
| 105 | |
|---|
| 106 | CALL dfftw_execute_dft_c2r(plan_backward(nb_vect),TF_vect,vect) |
|---|
| 107 | |
|---|
| 108 | vect = scale_factor * vect |
|---|
| 109 | |
|---|
| 110 | END SUBROUTINE fft_backward |
|---|
| 111 | |
|---|
| 112 | #endif |
|---|
| 113 | |
|---|
| 114 | END MODULE mod_fft_fftw |
|---|
