Changeset 1383


Ignore:
Timestamp:
May 6, 2010, 4:18:50 PM (15 years ago)
Author:
yann meurdesoif
Message:

ajout filtre FFT - FFTW pour vargas (entre autre)

YM + EM (pour le support moral)

Location:
LMDZ4/branches/LMDZ4V5.0-dev
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • LMDZ4/branches/LMDZ4V5.0-dev/arch/arch-PW6_VARGAS.fcm

    r1357 r1383  
    33%AR                  ar
    44%MAKE                gmake
    5 %FPP_FLAGS           -P
    6 %FPP_DEF             NC_DOUBLE BLAS SGEMV=DGEMV SGEMM=DGEMM
     5%FPP_FLAGS           -P -I/usr/local/pub/FFTW/3.2/include
     6%FPP_DEF             NC_DOUBLE BLAS SGEMV=DGEMV SGEMM=DGEMM FFT_FFTW
    77%BASE_FFLAGS         -qautodbl=dbl4 -qxlf90=autodealloc -qmaxmem=-1 -qzerosize
    88%PROD_FFLAGS         -O5
  • LMDZ4/branches/LMDZ4V5.0-dev/libf/filtrez/mod_fft_fftw.F90

    r986 r1383  
     1!
     2! $Id$
     3!
     4
    15MODULE mod_fft_fftw
    26
    37#ifdef FFT_FFTW
    48
    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
    1012 
    11   INTEGER,SAVE            :: plan_forward
    12   INTEGER,SAVE            :: plan_backward
     13  INTEGER*8, ALLOCATABLE, DIMENSION(:), SAVE :: plan_forward
     14  INTEGER*8, ALLOCATABLE, DIMENSION(:), SAVE :: plan_backward
    1315 
    1416CONTAINS
    1517 
    16   SUBROUTINE Init_fft(iim)
     18  SUBROUTINE Init_fft(iim,nvectmax)
    1719  IMPLICIT NONE
    18 #include <rfftw.h>
     20#include <fftw3.f>
    1921    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))
    2543   
    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(*,*)"!---------------------!"
    3049   
    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
    3458
    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
    3883  END SUBROUTINE Init_fft
    3984 
     
    4186  SUBROUTINE fft_forward(vect,TF_vect,nb_vect)
    4287    IMPLICIT NONE
    43 #include <rfftw.h>
    44     INTEGER,INTENT(IN)  :: nb_vect
    45     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)
    4691    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
    5497  END SUBROUTINE fft_forward
    5598 
    5699  SUBROUTINE fft_backward(TF_vect,vect,nb_vect)
    57100    IMPLICIT NONE
    58 #include <rfftw.h>
    59     INTEGER,INTENT(IN)  :: nb_vect
    60     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)
    61104    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
    68109
    69110  END SUBROUTINE fft_backward
     
    72113 
    73114END MODULE mod_fft_fftw
    74 
  • LMDZ4/branches/LMDZ4V5.0-dev/libf/filtrez/mod_filtre_fft.F90

    r1279 r1383  
     1!
     2! $Id$
     3!
     4
    15MODULE mod_filtre_fft
    26
     
    2327    INTEGER            :: index_vp(iim)
    2428    INTEGER            :: i,j
    25    
     29    INTEGER            :: l,ll_nb
     30
    2631    index_vp(1)=1
    2732    DO i=1,iim/2
     
    98103    ENDDO
    99104   
    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   
    101112    CALL Init_FFT(iim,(jjm+1)*(llm+1))
    102        
     113#endif       
    103114   
    104115  END SUBROUTINE Init_filtre_fft
     
    118129
    119130    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)
    123132    INTEGER            :: nb_vect
    124133    INTEGER :: i,j,l
    125134    INTEGER :: ll_nb
    126 !    REAL               :: vect_tmp(iim+inc,jj_end-jj_begin+1,nbniv)
    127135   
    128136    ll_nb=0
     
    140148    nb_vect=(jj_end-jj_begin+1)*ll_nb
    141149
    142 !    vect_tmp=vect
    143 
    144150    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+1
    149 !      DO i=1,iim/2+1
    150 !         PRINT *,"====",i,j,"----->",TF_vect_test(i,j,1)
    151 !       ENDDO
    152 !      ENDDO
    153151
    154152    DO l=1,ll_nb
     
    159157      ENDDO
    160158    ENDDO
    161        
     159 
    162160    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     
    172163    ll_nb=0
    173164!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     
    199190
    200191    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)
    202193    INTEGER            :: nb_vect
    203194    INTEGER :: i,j,l
     
    260251    REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
    261252
    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)
    264255    INTEGER            :: nb_vect
    265256    INTEGER :: i,j,l
     
    305296
    306297  END SUBROUTINE Filtre_inv_fft 
    307  
    308  
    309 !  SUBROUTINE get_ll_index(nbniv,ll_index,ll_nb)
    310 !  IMPLICIT NONE
    311 !    INTEGER,INTENT(IN)  :: nbniv
    312 !    INTEGER,INTENT(OUT) :: ll_index(nbniv)
    313 !    INTEGER,INTENT(OUT) :: ll_nb
    314 !
    315 !    INTEGER :: l,ll_begin, ll_end
    316 !   INTEGER :: omp_rank,omp_size
    317 !   INTEGER :: OMP_GET_NUM_THREADS
    318 !   INTEGER :: omp_chunk
    319 !   EXTERNAL OMP_GET_NUM_THREADS
    320 !   INTEGER :: OMP_GET_THREAD_NUM
    321 !   EXTERNAL OMP_GET_THREAD_NUM
    322 !
    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+1
    329 !   ll_nb=0
    330 !   DO WHILE (ll_begin<=nbniv)
    331 !     ll_end=min(ll_begin+OMP_CHUNK-1,nbniv)
    332 !     DO l=ll_begin,ll_end
    333 !       ll_nb=ll_nb+1
    334 !       ll_index(ll_nb)=l
    335 !     ENDDO
    336 !     ll_begin=ll_begin+omp_size*OMP_CHUNK
    337 !   ENDDO
    338 
    339 !  END SUBROUTINE get_ll_index
    340298   
    341299END MODULE mod_filtre_fft
Note: See TracChangeset for help on using the changeset viewer.