Line data Source code
1 : //# tFFT2D.cc: this tests FFT2D
2 : //# Copyright (C) 2016
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 :
29 : #include <casacore/casa/aips.h>
30 : #include <casacore/casa/Exceptions/Error.h>
31 : #include <casacore/casa/BasicSL/String.h>
32 : #include <casacore/casa/Logging/LogIO.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/images/Images/PagedImage.h>
35 : #include <casacore/images/Images/ImageConcat.h>
36 : #include <casacore/casa/namespace.h>
37 : #include <casacore/scimath/Mathematics/FFTPack.h>
38 : #include <casacore/images/Images/TempImage.h>
39 : #include <casacore/images/Images/PagedImage.h>
40 : #include <casacore/coordinates/Coordinates/CoordinateUtil.h>
41 :
42 : #include <synthesis/Utilities/FFT2D.h>
43 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
44 : #ifdef _OPENMP
45 : #include <omp.h>
46 : #endif
47 : using namespace casa;
48 1 : int main(int argc, char **argv)
49 : {
50 :
51 : try{
52 :
53 :
54 1 : Matrix<Double> xform(2,2);
55 1 : xform = 0.0;
56 1 : xform.diagonal() = 1.0;
57 2 : DirectionCoordinate dc(MDirection::J2000, Projection::SIN, Quantity(20.0,"deg"), Quantity(20.0, "deg"),
58 2 : Quantity(0.5, "arcsec"), Quantity(0.5,"arcsec"),
59 0 : xform, 50.0, 50.0, 999.0,
60 4 : 999.0);
61 1 : Vector<Int> whichStokes(1, Stokes::I);
62 1 : StokesCoordinate stc(whichStokes);
63 1 : SpectralCoordinate spc(MFrequency::LSRK, 1.5e9, 1e6, 0.0 , 1.420405752E9);
64 1 : CoordinateSystem cs;
65 1 : cs.addCoordinate(dc); cs.addCoordinate(stc); cs.addCoordinate(spc);
66 1 : Int x=10000; Int y=12000;
67 2 : PagedImage<Complex> im(IPosition(4,x,y,1,1), cs, "gulu0.image");
68 1 : im.set(0.0);
69 1 : im.putAt(Complex(3.0, 1.5),IPosition(4,x/100*20,y/100*80,0,0));
70 1 : im.putAt (Complex(4.0, 2.0),IPosition(4,x/100*60,y/100*70,0,0));
71 1 : im.putAt(Complex(1.0, 0.5),IPosition(4,x/100*10,y/100*10,0,0));
72 1 : im.putAt(Complex(2.0, 1.0),IPosition(4,x/100*75,y/100*25,0,0));
73 2 : PagedImage<Complex> im3(IPosition(4,x,y,1,1), cs, "gulu_lat.image");
74 1 : im3.copyData(im);
75 1 : Double wtime0=0.0;
76 1 : Double wtime1=0.0;
77 1 : Double wtime2=0.0;
78 :
79 1 : Int numthreads=-1;
80 1 : if(argc >1 ){
81 0 : numthreads=atoi(argv[1]);
82 : #ifdef _OPENMP
83 0 : omp_set_num_threads(numthreads);
84 : #endif
85 : }
86 1 : FFT2D ftw(true);
87 : {
88 1 : Array<Complex> arr0;
89 1 : Array<Complex> arr1;
90 : //////Lets do FFT via FFTPack
91 : {
92 1 : FFT2D ft(false);
93 1 : im.get(arr0, true);
94 : #ifdef _OPENMP
95 1 : wtime0=omp_get_wtime();
96 : #endif
97 : Bool del;
98 1 : Complex *scr= arr0.getStorage(del);
99 : //ft.fftShift(scr, Long(x), Long(y), true);
100 1 : ft.doFFT(scr, Long(x), Long(y), true);
101 : #ifdef _OPENMP
102 1 : cerr << "FFTPack forward " << x << " by " << y << " complex takes " << omp_get_wtime()-wtime0 << endl;
103 1 : arr0.putStorage(scr, del);
104 1 : arr0=Complex(x*y)*arr0;
105 1 : scr= arr0.getStorage(del);
106 1 : wtime0=omp_get_wtime();
107 : #endif
108 :
109 1 : ft.doFFT(scr, Long(x), Long(y), false);
110 : #ifdef _OPENMP
111 1 : cerr << "FFTPack backward " << x << " by " << y << " complex takes " << omp_get_wtime()-wtime0 << endl;
112 : #endif
113 1 : arr0.putStorage(scr, del);
114 1 : }
115 : //////Lets do FFT via FFTW
116 :
117 : {
118 :
119 1 : im.get(arr1, true);
120 : #ifdef _OPENMP
121 1 : wtime0=omp_get_wtime();
122 : #endif
123 :
124 : Bool del;
125 1 : Complex *scr= arr1.getStorage(del);
126 1 : ftw.doFFT(scr, Long(x), Long(y), true);
127 :
128 : #ifdef _OPENMP
129 1 : cerr << "FFTW forward " << x << " by " << y << " complex takes " << omp_get_wtime()-wtime0 << endl;
130 1 : wtime0=omp_get_wtime();
131 : #endif
132 :
133 1 : ftw.doFFT(scr, Long(x), Long(y), false);
134 :
135 : #ifdef _OPENMP
136 1 : cerr << "FFTW backward " << x << " by " << y << " complex takes " << omp_get_wtime()-wtime0 << endl;
137 : #endif
138 :
139 :
140 1 : arr1.putStorage(scr, del);
141 : }
142 :
143 : //////Lets do FFT via FFTW 1 thread
144 : {
145 :
146 :
147 1 : im.get(arr1, true);
148 :
149 : #ifdef _OPENMP
150 1 : numthreads=omp_get_max_threads();
151 1 : omp_set_num_threads(1);
152 1 : wtime0=omp_get_wtime();
153 : #endif
154 :
155 1 : FFT2D ft(true);
156 : Bool del;
157 1 : Complex *scr= arr1.getStorage(del);
158 1 : ft.doFFT(scr, Long(x), Long(y), true);
159 :
160 : #ifdef _OPENMP
161 1 : cerr << "FFTW 1-thread forward " << x << " by " << y << " complex takes " << omp_get_wtime()-wtime0 << endl;
162 1 : omp_set_num_threads(numthreads);
163 1 : wtime0=omp_get_wtime();
164 : #endif
165 :
166 1 : ft.doFFT(scr, Long(x), Long(y), false);
167 :
168 : #ifdef _OPENMP
169 1 : cerr << "FFTW 1-thread backward " << x << " by " << y << " complex takes " << omp_get_wtime()-wtime0 << endl;
170 1 : arr1.putStorage(scr, del);
171 : #endif
172 1 : }
173 : //Bool del;
174 : //Complex *scr0=arr0.getStorage(del);
175 : //Complex *scr1=arr1.getStorage(del);
176 1 : cerr << std::setprecision(9) << "max bet FFTW/FFTPack " << max(arr0) << " " << max(arr1) << endl;
177 1 : }
178 : #ifdef _OPENMP
179 1 : wtime0=omp_get_wtime();
180 : #endif
181 1 : Array<Complex> arr;
182 1 : Bool isRef=im.get(arr, true);
183 : #ifdef _OPENMP
184 1 : wtime1=omp_get_wtime();
185 : #endif
186 : Bool del;
187 1 : Complex *scr= arr.getStorage(del);
188 1 : cerr << "isRef " << isRef << " storage is copy " << del << endl;
189 : //FFT2D smp(true);
190 1 : ftw.c2cFFT(scr, Long(x), Long(y));
191 : #ifdef _OPENMP
192 1 : wtime2=omp_get_wtime();
193 : #endif
194 1 : if(!isRef)
195 1 : im.put(arr.reform(IPosition(4, x, y,1,1)));
196 1 : Double totTimeFFT2D=0.0;
197 : #ifdef _OPENMP
198 1 : totTimeFFT2D= -wtime0 + omp_get_wtime();
199 1 : cerr << "getting array " << wtime1-wtime0 << " fft " << wtime2-wtime1 << " put " << omp_get_wtime()-wtime2 << endl;
200 1 : wtime0=omp_get_wtime();
201 : #endif
202 1 : ftw.c2cFFT(scr, Long(x), Long(y), false);
203 : #ifdef _OPENMP
204 1 : wtime1=omp_get_wtime();
205 : #endif
206 2 : PagedImage<Complex> im2(IPosition(4,x,y,1,1), cs, "gulu_out.image");
207 1 : im2.put(arr.reform(IPosition(4, x, y,1,1)));
208 : #ifdef _OPENMP
209 1 : cerr << "putting array " << omp_get_wtime()-wtime1 << " fft " << wtime1 -wtime0 << endl;
210 : #endif
211 : {
212 :
213 : #ifdef _OPENMP
214 1 : wtime1=omp_get_wtime();
215 : #endif
216 1 : LatticeFFT::cfft2d(im3, true);
217 : }
218 : #ifdef _OPENMP
219 1 : cerr << "LatticeFFT::cfft2d " << omp_get_wtime()-wtime1 << " as compared to disk array based FFT2D " << totTimeFFT2D << endl;
220 : #endif
221 1 : Array<Complex> arr2;
222 1 : im.get(arr, true);
223 1 : im3.get(arr2, true);
224 1 : cerr << "max diff lattFFT and FFT2D " << max(arr-arr2) << endl;
225 :
226 :
227 1 : }catch( AipsError e ){
228 0 : cout << "Exception ocurred." << endl;
229 0 : cout << e.getMesg() << endl;
230 0 : }
231 1 : cout << "OK"<< endl;
232 1 : return 0;
233 : };
|
