Index: /LMDZ4/branches/LMDZ4V5.0-dev/arch/arch-PW6_VARGAS.fcm =================================================================== --- /LMDZ4/branches/LMDZ4V5.0-dev/arch/arch-PW6_VARGAS.fcm (revision 1382) +++ /LMDZ4/branches/LMDZ4V5.0-dev/arch/arch-PW6_VARGAS.fcm (revision 1383) @@ -3,6 +3,6 @@ %AR ar %MAKE gmake -%FPP_FLAGS -P -%FPP_DEF NC_DOUBLE BLAS SGEMV=DGEMV SGEMM=DGEMM +%FPP_FLAGS -P -I/usr/local/pub/FFTW/3.2/include +%FPP_DEF NC_DOUBLE BLAS SGEMV=DGEMV SGEMM=DGEMM FFT_FFTW %BASE_FFLAGS -qautodbl=dbl4 -qxlf90=autodealloc -qmaxmem=-1 -qzerosize %PROD_FFLAGS -O5 Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/filtrez/mod_fft_fftw.F90 =================================================================== --- /LMDZ4/branches/LMDZ4V5.0-dev/libf/filtrez/mod_fft_fftw.F90 (revision 1382) +++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/filtrez/mod_fft_fftw.F90 (revision 1383) @@ -1,39 +1,84 @@ +! +! $Id$ +! + MODULE mod_fft_fftw #ifdef FFT_FFTW - REAL,SAVE,ALLOCATABLE :: Table_forward(:) - REAL,SAVE,ALLOCATABLE :: Table_backward(:) - REAL,SAVE :: scale_factor - INTEGER,SAVE :: vsize - INTEGER,PARAMETER :: inc=1 + REAL, SAVE :: scale_factor + INTEGER, SAVE :: vsize + INTEGER, PARAMETER :: inc=1 - INTEGER,SAVE :: plan_forward - INTEGER,SAVE :: plan_backward + INTEGER*8, ALLOCATABLE, DIMENSION(:), SAVE :: plan_forward + INTEGER*8, ALLOCATABLE, DIMENSION(:), SAVE :: plan_backward CONTAINS - SUBROUTINE Init_fft(iim) + SUBROUTINE Init_fft(iim,nvectmax) IMPLICIT NONE -#include +#include INTEGER :: iim - REAL :: rtmp=1. - COMPLEX*16 :: ctmp - INTEGER :: itmp=1 - INTEGER :: isign=0 - INTEGER :: ierr + INTEGER :: nvectmax + + INTEGER :: itmp + + INTEGER :: rank + INTEGER :: howmany + INTEGER :: istride, idist + INTEGER :: ostride, odist + INTEGER, DIMENSION(1) :: n_array, inembed, onembed + + REAL, DIMENSION(iim+1,nvectmax) :: dbidon + COMPLEX, DIMENSION(iim/2+1,nvectmax) :: cbidon + + vsize = iim + scale_factor = 1./SQRT(1.*vsize) + + dbidon = 0 + cbidon = 0 + + ALLOCATE(plan_forward(nvectmax)) + ALLOCATE(plan_backward(nvectmax)) - vsize=iim - scale_factor=1./SQRT(1.*vsize) - ALLOCATE(Table_forward(2*vsize+64)) - ALLOCATE(Table_backward(2*vsize+64)) + WRITE(*,*)"!---------------------!" + WRITE(*,*)"! !" + WRITE(*,*)"! INITIALISATION FFTW !" + WRITE(*,*)"! !" + WRITE(*,*)"!---------------------!" -! CALL DZFFTM(isign,vsize,itmp,scale_factor,rtmp,vsize+inc,ctmp,vsize/2+1,table_forward,rtmp,ierr) - -! CALL ZDFFTM(isign,vsize,itmp,scale_factor,ctmp,vsize/2+1,rtmp,vsize+inc,table_backward,rtmp,ierr) +! On initialise tous les plans + DO itmp = 1, nvectmax + rank = 1 + n_array(1) = iim + howmany = itmp + inembed(1) = iim + 1 ; onembed(1) = iim/2 + 1 + istride = 1 ; ostride = 1 + idist = iim + 1 ; odist = iim/2 + 1 - CALL rfftw_f77_create_plan(plan_forward,iim,FFTW_REAL_TO_COMPLEX,FFTW_ESTIMATE) - CALL rfftw_f77_create_plan(plan_backward,iim,FFTW_COMPLEX_TO_REAL,FFTW_ESTIMATE) - + CALL dfftw_plan_many_dft_r2c(plan_forward(itmp), rank, n_array, howmany, & + & dbidon, inembed, istride, idist, & + & cbidon, onembed, ostride, odist, & + & FFTW_ESTIMATE) + + rank = 1 + n_array(1) = iim + howmany = itmp + inembed(1) = iim/2 + 1 ; onembed(1) = iim + 1 + istride = 1 ; ostride = 1 + idist = iim/2 + 1 ; odist = iim + 1 + CALL dfftw_plan_many_dft_c2r(plan_backward(itmp), rank, n_array, howmany, & + & cbidon, inembed, istride, idist, & + & dbidon, onembed, ostride, odist, & + & FFTW_ESTIMATE) + + ENDDO + + WRITE(*,*)"!-------------------------!" + WRITE(*,*)"! !" + WRITE(*,*)"! FIN INITIALISATION FFTW !" + WRITE(*,*)"! !" + WRITE(*,*)"!-------------------------!" + END SUBROUTINE Init_fft @@ -41,29 +86,25 @@ SUBROUTINE fft_forward(vect,TF_vect,nb_vect) IMPLICIT NONE -#include - INTEGER,INTENT(IN) :: nb_vect - REAL,INTENT(IN) :: vect(vsize+inc,nb_vect) +#include + INTEGER,INTENT(IN) :: nb_vect + REAL,INTENT(IN) :: vect(vsize+inc,nb_vect) COMPLEX*16,INTENT(OUT) :: TF_vect(vsize/2+1,nb_vect) - REAL :: work(4*vsize*nb_vect) - INTEGER :: ierr - INTEGER, PARAMETER :: isign=-1 - -! CALL DZFFTM(isign,vsize,nb_vect,scale_factor,vect,vsize+inc,TF_vect,vsize/2+1,table_forward,work,ierr) - CALL rfftwnd_f77_real_to_complex(plan_forward,nb_vect,vect, 1, vsize+inc , TF_vect, 1, vsize/2+1); - + + CALL dfftw_execute_dft_r2c(plan_forward(nb_vect),vect,TF_vect) + + TF_vect = scale_factor * TF_vect + END SUBROUTINE fft_forward SUBROUTINE fft_backward(TF_vect,vect,nb_vect) IMPLICIT NONE -#include - INTEGER,INTENT(IN) :: nb_vect - REAL,INTENT(OUT) :: vect(vsize+inc,nb_vect) +#include + INTEGER,INTENT(IN) :: nb_vect + REAL,INTENT(OUT) :: vect(vsize+inc,nb_vect) COMPLEX*16,INTENT(IN ) :: TF_vect(vsize/2+1,nb_vect) - REAL :: work(4*vsize*nb_vect) - INTEGER :: ierr - INTEGER, PARAMETER :: isign=1 - -! CALL ZDFFTM(isign,vsize,nb_vect,scale_factor,TF_vect,vsize/2+1,vect,vsize+inc,table_backward,work,ierr) - CALL rfftwnd_f77_complex_to_real(plan_forward,nb_vect,TF_vect, 1, vsize/2+1 , vect, 1, vsize+inc); + + CALL dfftw_execute_dft_c2r(plan_backward(nb_vect),TF_vect,vect) + + vect = scale_factor * vect END SUBROUTINE fft_backward @@ -72,3 +113,2 @@ END MODULE mod_fft_fftw - Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/filtrez/mod_filtre_fft.F90 =================================================================== --- /LMDZ4/branches/LMDZ4V5.0-dev/libf/filtrez/mod_filtre_fft.F90 (revision 1382) +++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/filtrez/mod_filtre_fft.F90 (revision 1383) @@ -1,2 +1,6 @@ +! +! $Id$ +! + MODULE mod_filtre_fft @@ -23,5 +27,6 @@ INTEGER :: index_vp(iim) INTEGER :: i,j - + INTEGER :: l,ll_nb + index_vp(1)=1 DO i=1,iim/2 @@ -98,7 +103,13 @@ ENDDO - +#ifdef FFT_FFTW + + WRITE (*,*)"COTH jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv" + WRITE (*,*)jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv + WRITE (*,*)MAX(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1 + CALL Init_FFT(iim,(llm+1)*(MAX(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1)) +#else CALL Init_FFT(iim,(jjm+1)*(llm+1)) - +#endif END SUBROUTINE Init_filtre_fft @@ -118,11 +129,8 @@ REAL :: vect(iim+inc,jj_end-jj_begin+1,nbniv) -! REAL :: vect_test(iim+inc,jj_end-jj_begin+1,nbniv) - COMPLEX*16 :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv) -! COMPLEX*16 :: TF_vect_test(iim/2+1,jj_end-jj_begin+1,nbniv) + COMPLEX :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv) INTEGER :: nb_vect INTEGER :: i,j,l INTEGER :: ll_nb -! REAL :: vect_tmp(iim+inc,jj_end-jj_begin+1,nbniv) ll_nb=0 @@ -140,15 +148,5 @@ nb_vect=(jj_end-jj_begin+1)*ll_nb -! vect_tmp=vect - CALL FFT_forward(vect,TF_vect,nb_vect) - -! CALL FFT_forward(vect,TF_vect_test,nb_vect) -! PRINT *,"XXXXXXXXXXXXX Filtre_u_FFT xxxxxxxxxxxx" -! DO j=1,jj_end-jj_begin+1 -! DO i=1,iim/2+1 -! PRINT *,"====",i,j,"----->",TF_vect_test(i,j,1) -! ENDDO -! ENDDO DO l=1,ll_nb @@ -159,15 +157,8 @@ ENDDO ENDDO - + CALL FFT_backward(TF_vect,vect,nb_vect) -! CALL FFT_backward(TF_vect_test,vect_test,nb_vect) - -! PRINT *,"XXXXXXXXXXXXX Filtre_u_FFT xxxxxxxxxxxx" -! DO j=1,jj_end-jj_begin+1 -! DO i=1,iim -! PRINT *,"====",i,j,"----->",vect_test(i,j,1) -! ENDDO -! ENDDO - + + ll_nb=0 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) @@ -199,5 +190,5 @@ REAL :: vect(iim+inc,jj_end-jj_begin+1,nbniv) - COMPLEX*16 :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv) + COMPLEX :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv) INTEGER :: nb_vect INTEGER :: i,j,l @@ -260,6 +251,6 @@ REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv) - REAL :: vect(iim+inc,jj_end-jj_begin+1,nbniv) - COMPLEX*16 :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv) + REAL :: vect(iim+inc,jj_end-jj_begin+1,nbniv) + COMPLEX :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv) INTEGER :: nb_vect INTEGER :: i,j,l @@ -305,37 +296,4 @@ END SUBROUTINE Filtre_inv_fft - - -! SUBROUTINE get_ll_index(nbniv,ll_index,ll_nb) -! IMPLICIT NONE -! INTEGER,INTENT(IN) :: nbniv -! INTEGER,INTENT(OUT) :: ll_index(nbniv) -! INTEGER,INTENT(OUT) :: ll_nb -! -! INTEGER :: l,ll_begin, ll_end -! INTEGER :: omp_rank,omp_size -! INTEGER :: OMP_GET_NUM_THREADS -! INTEGER :: omp_chunk -! EXTERNAL OMP_GET_NUM_THREADS -! INTEGER :: OMP_GET_THREAD_NUM -! EXTERNAL OMP_GET_THREAD_NUM -! -! -! omp_size=OMP_GET_NUM_THREADS() -! omp_rank=OMP_GET_THREAD_NUM() -! omp_chunk=nbniv/omp_size+min(1,MOD(nbniv,omp_size)) -! -! ll_begin=omp_rank*OMP_CHUNK+1 -! ll_nb=0 -! DO WHILE (ll_begin<=nbniv) -! ll_end=min(ll_begin+OMP_CHUNK-1,nbniv) -! DO l=ll_begin,ll_end -! ll_nb=ll_nb+1 -! ll_index(ll_nb)=l -! ENDDO -! ll_begin=ll_begin+omp_size*OMP_CHUNK -! ENDDO -! -! END SUBROUTINE get_ll_index END MODULE mod_filtre_fft