1 /*
2 Copyright (c) 2009-2014, Jack Poulson
3 All rights reserved.
4
5 This file is part of Elemental and is under the BSD 2-Clause License,
6 which can be found in the LICENSE file in the root directory, or at
7 http://opensource.org/licenses/BSD-2-Clause
8 */
9 // NOTE: It is possible to simply include "elemental.hpp" instead
10 #include "elemental-lite.hpp"
11 #include ELEM_ENTRYWISEMAP_INC
12 #include ELEM_FROBENIUSNORM_INC
13 #include ELEM_PSEUDOSPECTRUM_INC
14
15 #include ELEM_BULLSHEAD_INC
16 #include ELEM_FOXLI_INC
17 #include ELEM_GRCAR_INC
18 #include ELEM_HATANONELSON_INC
19 #include ELEM_HELMHOLTZPML_INC
20 #include ELEM_LOTKIN_INC
21 #include ELEM_TREFETHEN_INC
22 #include ELEM_TRIANGLE_INC
23 #include ELEM_UNIFORM_INC
24 #include ELEM_UNIFORMHELMHOLTZGREENS_INC
25 #include ELEM_WHALE_INC
26 using namespace std;
27 using namespace elem;
28
29 typedef double Real;
30 typedef Complex<Real> C;
31
32 int
main(int argc,char * argv[])33 main( int argc, char* argv[] )
34 {
35 Initialize( argc, argv );
36
37 try
38 {
39 Int r = Input("--gridHeight","process grid height",0);
40 const bool colMajor = Input("--colMajor","column-major ordering?",true);
41 const Int matType =
42 Input("--matType","0:uniform,1:Haar,2:Lotkin,3:Grcar,4:FoxLi,"
43 "5:HelmholtzPML1D,6:HelmholtzPML2D,7:Trefethen,"
44 "8:Bull's head,9:Triangle,10:Whale,"
45 "11:UniformHelmholtzGreen's,12:HatanoNelson",4);
46 const Int n = Input("--size","height of matrix",100);
47 const Int nbAlg = Input("--nbAlg","algorithmic blocksize",96);
48 #ifdef ELEM_HAVE_SCALAPACK
49 // QR algorithm options
50 const Int nbDist = Input("--nbDist","distribution blocksize",32);
51 #else
52 // Spectral Divide and Conquer options
53 const Int cutoff = Input("--cutoff","problem size for QR",256);
54 const Int maxInnerIts = Input("--maxInnerIts","SDC limit",2);
55 const Int maxOuterIts = Input("--maxOuterIts","SDC limit",10);
56 const bool random = Input("--random","Random RRQR in SDC",true);
57 const Real sdcTol = Input("--sdcTol","Rel. tol. for SDC",1e-6);
58 const Real spreadFactor = Input("--spreadFactor","median pert.",1e-6);
59 const Real signTol = Input("--signTol","Sign tolerance for SDC",1e-9);
60 #endif
61 const Real realCenter = Input("--realCenter","real center",0.);
62 const Real imagCenter = Input("--imagCenter","imag center",0.);
63 Real realWidth = Input("--realWidth","x width of image",0.);
64 Real imagWidth = Input("--imagWidth","y width of image",0.);
65 const Real numReal = Input("--numReal","num real chunks",2);
66 const Real numImag = Input("--numImag","num imag chunks",2);
67 const Int realSize = Input("--realSize","number of x samples",100);
68 const Int imagSize = Input("--imagSize","number of y samples",100);
69 const bool arnoldi = Input("--arnoldi","use Arnoldi?",true);
70 const Int basisSize = Input("--basisSize","num basis vectors",10);
71 const Int maxIts = Input("--maxIts","maximum pseudospec iter's",200);
72 const Real psTol = Input("--psTol","tolerance for pseudospectra",1e-6);
73 // Uniform options
74 const Real uniformRealCenter =
75 Input("--uniformRealCenter","real center of uniform dist",0.);
76 const Real uniformImagCenter =
77 Input("--uniformImagCenter","imag center of uniform dist",0.);
78 const Real uniformRadius =
79 Input("--uniformRadius","radius of uniform dist",1.);
80 // Grcar options
81 const Int numBands = Input("--numBands","num bands for Grcar",3);
82 // Fox-Li options
83 const Real omega = Input("--omega","frequency for Fox-Li/Helm",16*M_PI);
84 // Helmholtz-PML options [also uses omega from Fox-Li]
85 const Int mx = Input("--mx","number of x points for HelmholtzPML",30);
86 const Int my = Input("--my","number of y points for HelmholtzPML",30);
87 const Int numPmlPoints = Input("--numPml","num PML points for Helm",5);
88 const double sigma = Input("--sigma","PML amplitude",1.5);
89 const double pmlExp = Input("--pmlExp","PML takeoff exponent",3.);
90 // Uniform Helmholtz Green's options
91 const double lambda = Input("--lambda","wavelength of U.H.Green's",0.1);
92 // Hatano-Nelson options
93 const double gHatano = Input("--gHatano","g in Hatano-Nelson",0.5);
94 const bool periodic = Input("--periodic","periodic HatanoNelson?",true);
95 // Input/Output options
96 const bool progress = Input("--progress","print progress?",true);
97 const bool deflate = Input("--deflate","deflate?",true);
98 const bool display = Input("--display","display matrices?",false);
99 const bool write = Input("--write","write matrices?",false);
100 const bool saveSchur = Input("--saveSchur","save Schur factor?",true);
101 const Int numSaveFreq =
102 Input("--numSaveFreq","numerical save frequency",-1);
103 const Int imgSaveFreq =
104 Input("--imgSaveFreq","image save frequency",-1);
105 const Int imgDispFreq =
106 Input("--imgDispFreq","image display frequency",-1);
107 const std::string numBase =
108 Input("--numBase","numerical save basename",std::string("num"));
109 const std::string imgBase =
110 Input("--imgBase","image save basename",std::string("img"));
111 const Int numFormatInt = Input("--numFormat","numerical format",2);
112 const Int imgFormatInt = Input("--imgFormat","image format",8);
113 const Int colorMapInt = Input("--colorMap","color map",0);
114 const bool itCounts = Input("--itCounts","display iter. counts?",true);
115 ProcessInput();
116 PrintInputReport();
117
118 if( r == 0 )
119 r = Grid::FindFactor( mpi::Size(mpi::COMM_WORLD) );
120 const GridOrder order = ( colMajor ? COLUMN_MAJOR : ROW_MAJOR );
121 const Grid g( mpi::COMM_WORLD, r, order );
122 SetBlocksize( nbAlg );
123 if( numFormatInt < 1 || numFormatInt >= FileFormat_MAX )
124 LogicError("Invalid numerical format integer, should be in [1,",
125 FileFormat_MAX,")");
126 if( imgFormatInt < 1 || imgFormatInt >= FileFormat_MAX )
127 LogicError("Invalid image format integer, should be in [1,",
128 FileFormat_MAX,")");
129
130 const FileFormat numFormat = static_cast<FileFormat>(numFormatInt);
131 const FileFormat imgFormat = static_cast<FileFormat>(imgFormatInt);
132 const ColorMap colorMap = static_cast<ColorMap>(colorMapInt);
133 SetColorMap( colorMap );
134 const C center(realCenter,imagCenter);
135 const C uniformCenter(uniformRealCenter,uniformImagCenter);
136
137 bool isReal = true;
138 std::string matName;
139 DistMatrix<Real> AReal(g);
140 DistMatrix<C> ACpx(g);
141 switch( matType )
142 {
143 case 0: matName="uniform";
144 Uniform( ACpx, n, n, uniformCenter, uniformRadius );
145 isReal = false;
146 break;
147 case 1: matName="Haar";
148 Haar( ACpx, n );
149 isReal = false;
150 break;
151 case 2: matName="Lotkin";
152 Lotkin( AReal, n );
153 isReal = true;
154 break;
155 case 3: matName="Grcar";
156 Grcar( AReal, n, numBands );
157 isReal = true;
158 break;
159 case 4: matName="FoxLi";
160 FoxLi( ACpx, n, omega );
161 isReal = false;
162 break;
163 case 5: matName="HelmholtzPML";
164 HelmholtzPML
165 ( ACpx, n, C(omega), numPmlPoints, sigma, pmlExp );
166 isReal = false;
167 break;
168 case 6: matName="HelmholtzPML2D";
169 HelmholtzPML
170 ( ACpx, mx, my, C(omega), numPmlPoints, sigma, pmlExp );
171 isReal = false;
172 break;
173 case 7: matName="Trefethen";
174 Trefethen( ACpx, n );
175 isReal = false;
176 break;
177 case 8: matName="BullsHead";
178 BullsHead( ACpx, n );
179 isReal = false;
180 break;
181 case 9: matName="Triangle";
182 Triangle( AReal, n );
183 isReal = true;
184 break;
185 case 10: matName="Whale";
186 Whale( ACpx, n );
187 isReal = false;
188 break;
189 case 11: matName="UniformHelmholtzGreens";
190 UniformHelmholtzGreens( ACpx, n, lambda );
191 isReal = false;
192 break;
193 case 12: matName="HatanoNelson";
194 HatanoNelson
195 ( AReal, n, realCenter, uniformRadius, gHatano, periodic );
196 isReal = true;
197 break;
198 default: LogicError("Invalid matrix type");
199 }
200 if( display )
201 {
202 if( isReal )
203 Display( AReal, "A" );
204 else
205 Display( ACpx, "A" );
206 }
207 if( write )
208 {
209 if( isReal )
210 {
211 Write( AReal, "A", numFormat );
212 Write( AReal, "A", imgFormat );
213 }
214 else
215 {
216 Write( ACpx, "A", numFormat );
217 Write( ACpx, "A", imgFormat );
218 }
219 }
220
221 // Begin by computing the Schur decomposition
222 Timer timer;
223 DistMatrix<C,VR,STAR> w(g);
224 mpi::Barrier( mpi::COMM_WORLD );
225 const bool formATR = true;
226 #ifdef ELEM_HAVE_SCALAPACK
227 SetDefaultBlockHeight( nbDist );
228 SetDefaultBlockWidth( nbDist );
229 timer.Start();
230 if( isReal )
231 schur::QR( AReal, w, formATR );
232 else
233 schur::QR( ACpx, w, formATR );
234 mpi::Barrier( mpi::COMM_WORLD );
235 const double qrTime = timer.Stop();
236 if( mpi::WorldRank() == 0 )
237 std::cout << "QR algorithm took " << qrTime << " seconds"
238 << std::endl;
239 #else
240 SdcCtrl<Real> sdcCtrl;
241 sdcCtrl.cutoff = cutoff;
242 sdcCtrl.maxInnerIts = maxInnerIts;
243 sdcCtrl.maxOuterIts = maxOuterIts;
244 sdcCtrl.tol = sdcTol;
245 sdcCtrl.spreadFactor = spreadFactor;
246 sdcCtrl.random = random;
247 sdcCtrl.progress = progress;
248 sdcCtrl.signCtrl.tol = signTol;
249 sdcCtrl.signCtrl.progress = progress;
250 timer.Start();
251 if( isReal )
252 {
253 DistMatrix<Real> XReal(g);
254 schur::SDC( AReal, w, XReal, formATR, sdcCtrl );
255 }
256 else
257 {
258 DistMatrix<C> XCpx(g);
259 schur::SDC( ACpx, w, XCpx, formATR, sdcCtrl );
260 }
261 mpi::Barrier( mpi::COMM_WORLD );
262 const double sdcTime = timer.Stop();
263 if( mpi::WorldRank() == 0 )
264 std::cout << "SDC took " << sdcTime << " seconds" << std::endl;
265 #endif
266 if( saveSchur )
267 {
268 if( mpi::WorldRank() == 0 )
269 {
270 std::cout << "Writing Schur decomposition to file...";
271 std::cout.flush();
272 }
273 timer.Start();
274 if( isReal )
275 {
276 std::ostringstream os;
277 os << matName << "-"
278 << AReal.ColStride() << "x" << AReal.RowStride()
279 << "-" << AReal.DistRank();
280 write::Binary( AReal.LockedMatrix(), os.str() );
281 }
282 else
283 {
284 std::ostringstream os;
285 os << matName << "-"
286 << ACpx.ColStride() << "x" << ACpx.RowStride()
287 << "-" << ACpx.DistRank();
288 write::Binary( ACpx.LockedMatrix(), os.str() );
289 }
290 mpi::Barrier( mpi::COMM_WORLD );
291 const double saveSchurTime = timer.Stop();
292 if( mpi::WorldRank() == 0 )
293 std::cout << "DONE. " << saveSchurTime << " seconds"
294 << std::endl;
295 }
296
297 // Find a window if none is specified
298 if( realWidth == 0. || imagWidth == 0. )
299 {
300 const Real radius = MaxNorm( w );
301 const Real oneNorm = ( isReal ? OneNorm(AReal) : OneNorm(ACpx) );
302 Real width;
303 if( oneNorm == 0. && radius == 0. )
304 {
305 width = 1;
306 if( mpi::WorldRank() == 0 )
307 std::cout << "Setting width to 1 to handle zero matrix"
308 << std::endl;
309 }
310 else if( radius >= 0.2*oneNorm )
311 {
312 width = 2.5*radius;
313 if( mpi::WorldRank() == 0 )
314 std::cout << "Setting width to " << width
315 << " based on the spectral radius, "
316 << radius << std::endl;
317 }
318 else
319 {
320 width = 0.8*oneNorm;
321 if( mpi::WorldRank() == 0 )
322 std::cout << "Setting width to " << width
323 << " based on the one norm, " << oneNorm
324 << std::endl;
325 }
326 realWidth = width;
327 imagWidth = width;
328 }
329
330 PseudospecCtrl<Real> psCtrl;
331 psCtrl.schur = true;
332 psCtrl.maxIts = maxIts;
333 psCtrl.tol = psTol;
334 psCtrl.deflate = deflate;
335 psCtrl.arnoldi = arnoldi;
336 psCtrl.basisSize = basisSize;
337 psCtrl.progress = progress;
338 psCtrl.snapCtrl.imgSaveFreq = imgSaveFreq;
339 psCtrl.snapCtrl.numSaveFreq = numSaveFreq;
340 psCtrl.snapCtrl.imgDispFreq = imgDispFreq;
341 psCtrl.snapCtrl.imgFormat = imgFormat;
342 psCtrl.snapCtrl.numFormat = numFormat;
343 psCtrl.snapCtrl.itCounts = itCounts;
344
345 // Visualize/write the pseudospectrum within each window
346 DistMatrix<Real> invNormMap(g);
347 DistMatrix<Int> itCountMap(g);
348 const Int xBlock = realSize / numReal;
349 const Int yBlock = imagSize / numImag;
350 const Int xLeftover = realSize - (numReal-1)*xBlock;
351 const Int yLeftover = imagSize - (numImag-1)*yBlock;
352 const Real realStep = realWidth/realSize;
353 const Real imagStep = imagWidth/imagSize;
354 const C corner = center - C(realWidth/2,imagWidth/2);
355 for( Int realChunk=0; realChunk<numReal; ++realChunk )
356 {
357 const Int realChunkSize =
358 ( realChunk==numReal-1 ? xLeftover : xBlock );
359 const Real realChunkWidth = realStep*realChunkSize;
360 for( Int imagChunk=0; imagChunk<numImag; ++imagChunk )
361 {
362 std::ostringstream chunkStream;
363 chunkStream << "_" << realChunk << "_" << imagChunk;
364 const std::string chunkTag = chunkStream.str();
365
366 const Int imagChunkSize =
367 ( imagChunk==numImag-1 ? yLeftover : yBlock );
368 const Real imagChunkWidth = imagStep*imagChunkSize;
369
370 const C chunkCorner = corner +
371 C(realStep*realChunk*xBlock,imagStep*imagChunk*yBlock);
372 const C chunkCenter = chunkCorner +
373 0.5*C(realStep*realChunkSize,imagStep*imagChunkSize);
374
375 if( mpi::WorldRank() == 0 )
376 std::cout << "Starting computation for chunk centered at "
377 << chunkCenter << std::endl;
378 mpi::Barrier( mpi::COMM_WORLD );
379 timer.Start();
380 psCtrl.snapCtrl.numBase = matName+"-"+numBase+chunkTag;
381 psCtrl.snapCtrl.imgBase = matName+"-"+imgBase+chunkTag;
382 if( isReal )
383 {
384 itCountMap = QuasiTriangularPseudospectrum
385 ( AReal, invNormMap, chunkCenter,
386 realChunkWidth, imagChunkWidth,
387 realChunkSize, imagChunkSize, psCtrl );
388 }
389 else
390 {
391 itCountMap = TriangularPseudospectrum
392 ( ACpx, invNormMap, chunkCenter,
393 realChunkWidth, imagChunkWidth,
394 realChunkSize, imagChunkSize, psCtrl );
395 }
396 mpi::Barrier( mpi::COMM_WORLD );
397 const double pseudoTime = timer.Stop();
398 const Int numIts = MaxNorm( itCountMap );
399 if( mpi::WorldRank() == 0 )
400 {
401 std::cout << "num seconds=" << pseudoTime << "\n"
402 << "num iterations=" << numIts << std::endl;
403 }
404 }
405 }
406 }
407 catch( exception& e ) { ReportException(e); }
408
409 Finalize();
410 return 0;
411 }
412