Changeset 1383 for LMDZ4/branches/LMDZ4V5.0-dev/libf
- Timestamp:
- May 6, 2010, 4:18:50 PM (15 years ago)
- Location:
- LMDZ4/branches/LMDZ4V5.0-dev/libf/filtrez
- Files:
-
- 2 edited
-
mod_fft_fftw.F90 (modified) (3 diffs)
-
mod_filtre_fft.F90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ4/branches/LMDZ4V5.0-dev/libf/filtrez/mod_fft_fftw.F90
r986 r1383 1 ! 2 ! $Id$ 3 ! 4 1 5 MODULE mod_fft_fftw 2 6 3 7 #ifdef FFT_FFTW 4 8 5 REAL,SAVE,ALLOCATABLE :: Table_forward(:) 6 REAL,SAVE,ALLOCATABLE :: Table_backward(:) 7 REAL,SAVE :: scale_factor 8 INTEGER,SAVE :: vsize 9 INTEGER,PARAMETER :: inc=1 9 REAL, SAVE :: scale_factor 10 INTEGER, SAVE :: vsize 11 INTEGER, PARAMETER :: inc=1 10 12 11 INTEGER ,SAVE:: plan_forward12 INTEGER ,SAVE:: plan_backward13 INTEGER*8, ALLOCATABLE, DIMENSION(:), SAVE :: plan_forward 14 INTEGER*8, ALLOCATABLE, DIMENSION(:), SAVE :: plan_backward 13 15 14 16 CONTAINS 15 17 16 SUBROUTINE Init_fft(iim )18 SUBROUTINE Init_fft(iim,nvectmax) 17 19 IMPLICIT NONE 18 #include < rfftw.h>20 #include <fftw3.f> 19 21 INTEGER :: iim 20 REAL :: rtmp=1. 21 COMPLEX*16 :: ctmp 22 INTEGER :: itmp=1 23 INTEGER :: isign=0 24 INTEGER :: ierr 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)) 25 43 26 vsize=iim 27 scale_factor=1./SQRT(1.*vsize) 28 ALLOCATE(Table_forward(2*vsize+64)) 29 ALLOCATE(Table_backward(2*vsize+64)) 44 WRITE(*,*)"!---------------------!" 45 WRITE(*,*)"! !" 46 WRITE(*,*)"! INITIALISATION FFTW !" 47 WRITE(*,*)"! !" 48 WRITE(*,*)"!---------------------!" 30 49 31 ! CALL DZFFTM(isign,vsize,itmp,scale_factor,rtmp,vsize+inc,ctmp,vsize/2+1,table_forward,rtmp,ierr) 32 33 ! CALL ZDFFTM(isign,vsize,itmp,scale_factor,ctmp,vsize/2+1,rtmp,vsize+inc,table_backward,rtmp,ierr) 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 34 58 35 CALL rfftw_f77_create_plan(plan_forward,iim,FFTW_REAL_TO_COMPLEX,FFTW_ESTIMATE) 36 CALL rfftw_f77_create_plan(plan_backward,iim,FFTW_COMPLEX_TO_REAL,FFTW_ESTIMATE) 37 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 38 83 END SUBROUTINE Init_fft 39 84 … … 41 86 SUBROUTINE fft_forward(vect,TF_vect,nb_vect) 42 87 IMPLICIT NONE 43 #include < rfftw.h>44 INTEGER,INTENT(IN) :: nb_vect45 REAL,INTENT(IN) :: vect(vsize+inc,nb_vect)88 #include <fftw3.f> 89 INTEGER,INTENT(IN) :: nb_vect 90 REAL,INTENT(IN) :: vect(vsize+inc,nb_vect) 46 91 COMPLEX*16,INTENT(OUT) :: TF_vect(vsize/2+1,nb_vect) 47 REAL :: work(4*vsize*nb_vect) 48 INTEGER :: ierr 49 INTEGER, PARAMETER :: isign=-1 50 51 ! CALL DZFFTM(isign,vsize,nb_vect,scale_factor,vect,vsize+inc,TF_vect,vsize/2+1,table_forward,work,ierr) 52 CALL rfftwnd_f77_real_to_complex(plan_forward,nb_vect,vect, 1, vsize+inc , TF_vect, 1, vsize/2+1); 53 92 93 CALL dfftw_execute_dft_r2c(plan_forward(nb_vect),vect,TF_vect) 94 95 TF_vect = scale_factor * TF_vect 96 54 97 END SUBROUTINE fft_forward 55 98 56 99 SUBROUTINE fft_backward(TF_vect,vect,nb_vect) 57 100 IMPLICIT NONE 58 #include < rfftw.h>59 INTEGER,INTENT(IN) :: nb_vect60 REAL,INTENT(OUT) :: vect(vsize+inc,nb_vect)101 #include <fftw3.f> 102 INTEGER,INTENT(IN) :: nb_vect 103 REAL,INTENT(OUT) :: vect(vsize+inc,nb_vect) 61 104 COMPLEX*16,INTENT(IN ) :: TF_vect(vsize/2+1,nb_vect) 62 REAL :: work(4*vsize*nb_vect) 63 INTEGER :: ierr 64 INTEGER, PARAMETER :: isign=1 65 66 ! CALL ZDFFTM(isign,vsize,nb_vect,scale_factor,TF_vect,vsize/2+1,vect,vsize+inc,table_backward,work,ierr) 67 CALL rfftwnd_f77_complex_to_real(plan_forward,nb_vect,TF_vect, 1, vsize/2+1 , vect, 1, vsize+inc); 105 106 CALL dfftw_execute_dft_c2r(plan_backward(nb_vect),TF_vect,vect) 107 108 vect = scale_factor * vect 68 109 69 110 END SUBROUTINE fft_backward … … 72 113 73 114 END MODULE mod_fft_fftw 74 -
LMDZ4/branches/LMDZ4V5.0-dev/libf/filtrez/mod_filtre_fft.F90
r1279 r1383 1 ! 2 ! $Id$ 3 ! 4 1 5 MODULE mod_filtre_fft 2 6 … … 23 27 INTEGER :: index_vp(iim) 24 28 INTEGER :: i,j 25 29 INTEGER :: l,ll_nb 30 26 31 index_vp(1)=1 27 32 DO i=1,iim/2 … … 98 103 ENDDO 99 104 100 105 #ifdef FFT_FFTW 106 107 WRITE (*,*)"COTH jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv" 108 WRITE (*,*)jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv 109 WRITE (*,*)MAX(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1 110 CALL Init_FFT(iim,(llm+1)*(MAX(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1)) 111 #else 101 112 CALL Init_FFT(iim,(jjm+1)*(llm+1)) 102 113 #endif 103 114 104 115 END SUBROUTINE Init_filtre_fft … … 118 129 119 130 REAL :: vect(iim+inc,jj_end-jj_begin+1,nbniv) 120 ! REAL :: vect_test(iim+inc,jj_end-jj_begin+1,nbniv) 121 COMPLEX*16 :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv) 122 ! COMPLEX*16 :: TF_vect_test(iim/2+1,jj_end-jj_begin+1,nbniv) 131 COMPLEX :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv) 123 132 INTEGER :: nb_vect 124 133 INTEGER :: i,j,l 125 134 INTEGER :: ll_nb 126 ! REAL :: vect_tmp(iim+inc,jj_end-jj_begin+1,nbniv)127 135 128 136 ll_nb=0 … … 140 148 nb_vect=(jj_end-jj_begin+1)*ll_nb 141 149 142 ! vect_tmp=vect143 144 150 CALL FFT_forward(vect,TF_vect,nb_vect) 145 146 ! CALL FFT_forward(vect,TF_vect_test,nb_vect)147 ! PRINT *,"XXXXXXXXXXXXX Filtre_u_FFT xxxxxxxxxxxx"148 ! DO j=1,jj_end-jj_begin+1149 ! DO i=1,iim/2+1150 ! PRINT *,"====",i,j,"----->",TF_vect_test(i,j,1)151 ! ENDDO152 ! ENDDO153 151 154 152 DO l=1,ll_nb … … 159 157 ENDDO 160 158 ENDDO 161 159 162 160 CALL FFT_backward(TF_vect,vect,nb_vect) 163 ! CALL FFT_backward(TF_vect_test,vect_test,nb_vect) 164 165 ! PRINT *,"XXXXXXXXXXXXX Filtre_u_FFT xxxxxxxxxxxx" 166 ! DO j=1,jj_end-jj_begin+1 167 ! DO i=1,iim 168 ! PRINT *,"====",i,j,"----->",vect_test(i,j,1) 169 ! ENDDO 170 ! ENDDO 171 161 162 172 163 ll_nb=0 173 164 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) … … 199 190 200 191 REAL :: vect(iim+inc,jj_end-jj_begin+1,nbniv) 201 COMPLEX *16:: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)192 COMPLEX :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv) 202 193 INTEGER :: nb_vect 203 194 INTEGER :: i,j,l … … 260 251 REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv) 261 252 262 REAL :: vect(iim+inc,jj_end-jj_begin+1,nbniv)263 COMPLEX *16:: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)253 REAL :: vect(iim+inc,jj_end-jj_begin+1,nbniv) 254 COMPLEX :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv) 264 255 INTEGER :: nb_vect 265 256 INTEGER :: i,j,l … … 305 296 306 297 END SUBROUTINE Filtre_inv_fft 307 308 309 ! SUBROUTINE get_ll_index(nbniv,ll_index,ll_nb)310 ! IMPLICIT NONE311 ! INTEGER,INTENT(IN) :: nbniv312 ! INTEGER,INTENT(OUT) :: ll_index(nbniv)313 ! INTEGER,INTENT(OUT) :: ll_nb314 !315 ! INTEGER :: l,ll_begin, ll_end316 ! INTEGER :: omp_rank,omp_size317 ! INTEGER :: OMP_GET_NUM_THREADS318 ! INTEGER :: omp_chunk319 ! EXTERNAL OMP_GET_NUM_THREADS320 ! INTEGER :: OMP_GET_THREAD_NUM321 ! EXTERNAL OMP_GET_THREAD_NUM322 !323 !324 ! omp_size=OMP_GET_NUM_THREADS()325 ! omp_rank=OMP_GET_THREAD_NUM()326 ! omp_chunk=nbniv/omp_size+min(1,MOD(nbniv,omp_size))327 !328 ! ll_begin=omp_rank*OMP_CHUNK+1329 ! ll_nb=0330 ! DO WHILE (ll_begin<=nbniv)331 ! ll_end=min(ll_begin+OMP_CHUNK-1,nbniv)332 ! DO l=ll_begin,ll_end333 ! ll_nb=ll_nb+1334 ! ll_index(ll_nb)=l335 ! ENDDO336 ! ll_begin=ll_begin+omp_size*OMP_CHUNK337 ! ENDDO338 !339 ! END SUBROUTINE get_ll_index340 298 341 299 END MODULE mod_filtre_fft
Note: See TracChangeset
for help on using the changeset viewer.
