source: LMDZ4/branches/LMDZ4_AR5/libf/filtrez/mod_fft_mkl.F90

Last change on this file was 1601, checked in by Ehouarn Millour, 13 years ago

Updates and upgrades of filter:

  • add minor corrections (r1591 of trunk) in filtreg_mod.F90: some arrays were oversized and the computation of the indexes from which the filter is applied could go wrong in extreme cases.
  • adapt filtreg_p.F (r1597 of trunk) so that BLAS routine DGEMM is only used if 'BLAS' preprocessing flag is set.
  • update FFT filter routines to match current 'trunk' versions (mostly cosmetic changes, exept for the implementation of use of FFTW in mod_fft_fftw.F90).
  • update "arch-PW6_VARGAS.fcm" to enable use of FFTW

NB: implemented FFTs assume working precision to be double precision ; trying to use them when working precision is single precision will clearly end in tragedy.
EM
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.0 KB
Line 
1MODULE mod_fft_mkl
2#ifdef FFT_MKL
3
4  USE MKL_DFTI
5 
6  REAL,SAVE                :: scale_factor
7  INTEGER,SAVE             :: vsize
8  INTEGER,PARAMETER        :: inc=1
9 
10!  TYPE FFT_HANDLE
11!    TYPE(DFTI_DESCRIPTOR), POINTER :: Pt
12!    LOGICAL :: IsAllocated
13!  END TYPE FFT_HANDLE
14 
15!  TYPE(FFT_HANDLE),SAVE,ALLOCATABLE :: Forward_Handle(:)
16!  TYPE(FFT_HANDLE),SAVE,ALLOCATABLE :: Backward_Handle(:)
17!!$OMP THREADPRIVATE(Forward_Handle,Backward_Handle) 
18 
19CONTAINS
20 
21  SUBROUTINE Init_fft(iim,nb_vect_max)
22    IMPLICIT NONE
23    INTEGER :: iim
24    INTEGER :: nb_vect_max
25    REAL    :: rtmp=1.
26    COMPLEX :: ctmp
27    INTEGER :: itmp=1
28    INTEGER :: isign=0
29    INTEGER :: ierr
30   
31    vsize=iim
32    scale_factor=1./SQRT(1.*vsize)
33!    ALLOCATE(Forward_Handle(nb_vect_max))
34!    ALLOCATE(Backward_Handle(nb_vect_max))
35   
36!    Forward_Handle(:)%IsAllocated=.FALSE.
37!    Backward_Handle(:)%IsAllocated=.FALSE.
38   
39!    ALLOCATE(Table_forward(2*vsize+64))
40!    ALLOCATE(Table_backward(2*vsize+64))
41!   
42!    CALL DZFFTM(isign,vsize,itmp,scale_factor,rtmp,vsize+inc,ctmp,vsize/2+1,table_forward,rtmp,ierr)
43!   
44!    CALL ZDFFTM(isign,vsize,itmp,scale_factor,ctmp,vsize/2+1,rtmp,vsize+inc,table_backward,rtmp,ierr)
45
46!    ierr = DftiCreateDescriptor( FFT_Handle, DFTI_DOUBLE, DFTI_REAL, 1, vsize )
47!    ierr = DftiSetValue(FFT_Handle,DFTI_NUMBER_OF_TRANSFORMS,nb_vect)
48!    ierr = DftiSetValue(FFT_Handle,DFTI_FORWARD_SCALE,scale_factor)
49!    ierr = DftiSetValue(FFT_Handle,DFTI_BACKWARD_SCALE,scale_factor)
50!    ierr = DftiSetValue(FFT_Handle,DFTI_PLACEMENT,DFTI_NOT_INPLACE)
51!    ierr = DftiSetValue(Desc_Handle, DFTI_INPUT_DISTANCE, vsize+inc)
52!    ierr = DftiSetValue(Desc_Handle, DFTI_OUTPUT_DISTANCE, vsize)
53!    ierr = DftiCommitDescriptor( FFT_HANDLE )
54
55  END SUBROUTINE Init_fft
56 
57 
58  SUBROUTINE fft_forward(vect,TF_vect,nb_vect)
59    IMPLICIT NONE
60    INTEGER,INTENT(IN)  :: nb_vect
61    REAL,INTENT(IN)     :: vect((vsize+inc)*nb_vect)
62    COMPLEX,INTENT(OUT) :: TF_vect((vsize/2+1)*nb_vect)
63    REAL                :: work(4*vsize*nb_vect)
64    INTEGER             :: ierr
65    INTEGER, PARAMETER :: isign=-1
66    REAL               :: vect_out((vsize+inc)*nb_vect)
67    TYPE(DFTI_DESCRIPTOR), POINTER :: FFT_Handle
68   
69!    IF ( .NOT. Forward_handle(nb_vect)%IsAllocated) THEN
70      ierr = DftiCreateDescriptor( FFT_Handle, DFTI_DOUBLE, DFTI_REAL, 1, vsize )
71      ierr = DftiSetValue(FFT_Handle,DFTI_NUMBER_OF_TRANSFORMS,nb_vect)
72      ierr = DftiSetValue(FFT_Handle,DFTI_FORWARD_SCALE,scale_factor)
73      ierr = DftiSetValue(FFT_Handle,DFTI_BACKWARD_SCALE,scale_factor)
74      ierr = DftiSetValue(FFT_Handle,DFTI_PLACEMENT,DFTI_NOT_INPLACE)
75      ierr = DftiSetValue(FFT_Handle, DFTI_INPUT_DISTANCE, vsize+inc)
76      ierr = DftiSetValue(FFT_Handle, DFTI_OUTPUT_DISTANCE, (vsize/2+1)*2)
77      ierr = DftiCommitDescriptor( FFT_Handle )
78!      Forward_handle(nb_vect)%IsAllocated=.TRUE.
79!    ENDIF
80   
81    ierr = DftiComputeForward( FFT_Handle, vect, TF_vect )
82
83    ierr = DftiFreeDescriptor( FFT_Handle )
84
85!    ierr = DftiCreateDescriptor( FFT_Handle, DFTI_DOUBLE, DFTI_REAL, 1, vsize )
86!    ierr = DftiSetValue(FFT_Handle,DFTI_NUMBER_OF_TRANSFORMS,nb_vect)
87!    ierr = DftiSetValue(FFT_Handle,DFTI_FORWARD_SCALE,scale_factor)
88!    ierr = DftiSetValue(FFT_Handle,DFTI_BACKWARD_SCALE,scale_factor)
89!    ierr = DftiSetValue(FFT_Handle,DFTI_PLACEMENT,DFTI_NOT_INPLACE)
90!    ierr = DftiSetValue(FFT_HANDLE, DFTI_INPUT_DISTANCE, vsize/2+1)
91!    ierr = DftiSetValue(FFT_HANDLE, DFTI_OUTPUT_DISTANCE, vsize+inc)
92!    ierr = DftiCommitDescriptor( FFT_HANDLE )
93!    ierr = DftiComputeBackward( FFT_HANDLE, TF_vect, vect_out )
94!    ierr = DftiFreeDescriptor( FFT_HANDLE )
95
96!    CALL DZFFTM(isign,vsize,nb_vect,scale_factor,vect,vsize+inc,TF_vect,vsize/2+1,table_forward,work,ierr)
97 
98  END SUBROUTINE fft_forward
99 
100  SUBROUTINE fft_backward(TF_vect,vect,nb_vect)
101    IMPLICIT NONE
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    REAL                :: work(4*vsize*nb_vect)
106    INTEGER             :: ierr
107    INTEGER, PARAMETER :: isign=1
108    TYPE(DFTI_DESCRIPTOR),POINTER :: FFT_Handle
109!    CALL ZDFFTM(isign,vsize,nb_vect,scale_factor,TF_vect,vsize/2+1,vect,vsize+inc,table_backward,work,ierr)
110!    IF ( .NOT. Backward_handle(nb_vect)%IsAllocated) THEN
111      ierr = DftiCreateDescriptor( FFT_Handle, DFTI_DOUBLE, DFTI_REAL, 1, vsize )
112      ierr = DftiSetValue(FFT_Handle,DFTI_NUMBER_OF_TRANSFORMS,nb_vect)
113      ierr = DftiSetValue(FFT_Handle,DFTI_FORWARD_SCALE,scale_factor)
114      ierr = DftiSetValue(FFT_Handle,DFTI_BACKWARD_SCALE,scale_factor)
115      ierr = DftiSetValue(FFT_Handle,DFTI_PLACEMENT,DFTI_NOT_INPLACE)
116      ierr = DftiSetValue(FFT_Handle, DFTI_INPUT_DISTANCE,  (vsize/2+1)*2)
117      ierr = DftiSetValue(FFT_Handle, DFTI_OUTPUT_DISTANCE, vsize+inc)
118      ierr = DftiCommitDescriptor( FFT_Handle )
119!      Backward_handle(nb_vect)%IsAllocated=.TRUE.
120!    ENDIF
121    ierr = DftiComputeBackward( FFT_Handle, TF_vect, vect )
122    ierr = DftiFreeDescriptor( FFT_Handle)
123 
124  END SUBROUTINE fft_backward
125
126#endif
127 
128END MODULE mod_fft_mkl
Note: See TracBrowser for help on using the repository browser.