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_GRCAR_INC
17 #include ELEM_HATANONELSON_INC
18 #include ELEM_FOXLI_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 const Real realWidth = Input("--realWidth","x width of image",0.);
64 const Real imagWidth = Input("--imagWidth","y width of image",0.);
65 const Int realSize = Input("--realSize","number of x samples",100);
66 const Int imagSize = Input("--imagSize","number of y samples",100);
67 const bool schur = Input("--schur","Schur decomposition?",true);
68 const bool forceComplexSchur =
69 Input("--forceComplexSchur",
70 "switch to complex arithmetic for QR alg.",false);
71 const bool forceComplexPs =
72 Input("--forceComplexPs",
73 "switch to complex arithmetic for PS iter's",true);
74 const bool arnoldi = Input("--arnoldi","use Arnoldi?",true);
75 const Int basisSize = Input("--basisSize","num Arnoldi vectors",10);
76 const Int maxIts = Input("--maxIts","maximum pseudospec iter's",200);
77 const Real psTol = Input("--psTol","tolerance for pseudospectra",1e-6);
78 // Uniform options
79 const Real uniformRealCenter =
80 Input("--uniformRealCenter","real center of uniform dist",0.);
81 const Real uniformImagCenter =
82 Input("--uniformImagCenter","imag center of uniform dist",0.);
83 const Real uniformRadius =
84 Input("--uniformRadius","radius of uniform dist",1.);
85 // Grcar options
86 const Int numBands = Input("--numBands","num bands for Grcar",3);
87 // Fox-Li options
88 const Real omega = Input("--omega","frequency for Fox-Li/Helm",16*M_PI);
89 // Helmholtz-PML options [also uses Fox-Li omega]
90 const Int mx = Input("--mx","number of x points for HelmholtzPML",30);
91 const Int my = Input("--my","number of y points for HelmholtzPML",30);
92 const Int numPmlPoints = Input("--numPml","num PML points for Helm",5);
93 const double sigma = Input("--sigma","PML amplitude",1.5);
94 const double pmlExp = Input("--pmlExp","PML takeoff exponent",3.);
95 // Uniform Helmholtz Green's options
96 const double lambda = Input("--lambda","wavelength of U.H.Green's",0.1);
97 // Hatano-Nelson options [also uses uniform real center]
98 const double gHatano = Input("--gHatano","g in Hatano-Nelson",0.5);
99 const bool periodic = Input("--periodic","periodic HatanoNelson?",true);
100 // Input/Output options
101 const bool progress = Input("--progress","print progress?",true);
102 const bool deflate = Input("--deflate","deflate?",true);
103 const bool display = Input("--display","display matrices?",false);
104 const bool write = Input("--write","write matrices?",false);
105 const Int numSaveFreq =
106 Input("--numSaveFreq","numerical save frequency",-1);
107 const Int imgSaveFreq =
108 Input("--imgSaveFreq","image save frequency",-1);
109 const Int imgDispFreq =
110 Input("--imgDispFreq","image display frequency",-1);
111 const std::string numBase =
112 Input("--numBase","numerical save basename",std::string("num"));
113 const std::string imgBase =
114 Input("--imgBase","image save basename",std::string("img"));
115 const Int numFormatInt = Input("--numFormat","numerical format",2);
116 const Int imgFormatInt = Input("--imgFormat","image format",8);
117 const Int colorMapInt = Input("--colorMap","color map",0);
118 const bool itCounts = Input("--itCounts","display iter. counts?",true);
119 ProcessInput();
120 PrintInputReport();
121
122 if( r == 0 )
123 r = Grid::FindFactor( mpi::Size(mpi::COMM_WORLD) );
124 const GridOrder order = ( colMajor ? COLUMN_MAJOR : ROW_MAJOR );
125 const Grid g( mpi::COMM_WORLD, r, order );
126 SetBlocksize( nbAlg );
127 #ifdef ELEM_HAVE_SCALAPACK
128 SetDefaultBlockHeight( nbDist );
129 SetDefaultBlockWidth( nbDist );
130 #endif
131 if( numFormatInt < 1 || numFormatInt >= FileFormat_MAX )
132 LogicError("Invalid numerical format integer, should be in [1,",
133 FileFormat_MAX,")");
134 if( imgFormatInt < 1 || imgFormatInt >= FileFormat_MAX )
135 LogicError("Invalid image format integer, should be in [1,",
136 FileFormat_MAX,")");
137
138 const FileFormat numFormat = static_cast<FileFormat>(numFormatInt);
139 const FileFormat imgFormat = static_cast<FileFormat>(imgFormatInt);
140 const ColorMap colorMap = static_cast<ColorMap>(colorMapInt);
141 SetColorMap( colorMap );
142 const C center(realCenter,imagCenter);
143 const C uniformCenter(uniformRealCenter,uniformImagCenter);
144
145 bool isReal = true;
146 std::string matName;
147 DistMatrix<Real> AReal(g);
148 DistMatrix<C> ACpx(g);
149 switch( matType )
150 {
151 case 0: matName="uniform";
152 Uniform( ACpx, n, n, uniformCenter, uniformRadius );
153 isReal = false;
154 break;
155 case 1: matName="Haar";
156 Haar( ACpx, n );
157 isReal = false;
158 break;
159 case 2: matName="Lotkin";
160 Lotkin( AReal, n );
161 isReal = true;
162 break;
163 case 3: matName="Grcar";
164 Grcar( AReal, n, numBands );
165 isReal = true;
166 break;
167 case 4: matName="FoxLi";
168 FoxLi( ACpx, n, omega );
169 isReal = false;
170 break;
171 case 5: matName="HelmholtzPML";
172 HelmholtzPML
173 ( ACpx, n, C(omega), numPmlPoints, sigma, pmlExp );
174 isReal = false;
175 break;
176 case 6: matName="HelmholtzPML2D";
177 HelmholtzPML
178 ( ACpx, mx, my, C(omega), numPmlPoints, sigma, pmlExp );
179 isReal = false;
180 break;
181 case 7: matName="Trefethen";
182 Trefethen( ACpx, n );
183 isReal = false;
184 break;
185 case 8: matName="BullsHead";
186 BullsHead( ACpx, n );
187 isReal = false;
188 break;
189 case 9: matName="Triangle";
190 Triangle( AReal, n );
191 isReal = true;
192 break;
193 case 10: matName="Whale";
194 Whale( ACpx, n );
195 isReal = false;
196 break;
197 case 11: matName="UniformHelmholtzGreens";
198 UniformHelmholtzGreens( ACpx, n, lambda );
199 isReal = false;
200 break;
201 case 12: matName="HatanoNelson";
202 HatanoNelson
203 ( AReal, n, realCenter, uniformRadius, gHatano, periodic );
204 isReal = true;
205 break;
206 default: LogicError("Invalid matrix type");
207 }
208 if( display )
209 {
210 if( isReal )
211 Display( AReal, "A" );
212 else
213 Display( ACpx, "A" );
214 }
215 if( write )
216 {
217 if( isReal )
218 {
219 Write( AReal, "A", numFormat );
220 Write( AReal, "A", imgFormat );
221 }
222 else
223 {
224 Write( ACpx, "A", numFormat );
225 Write( ACpx, "A", imgFormat );
226 }
227 }
228
229 PseudospecCtrl<Real> psCtrl;
230 psCtrl.schur = schur;
231 psCtrl.forceComplexSchur = forceComplexSchur;
232 psCtrl.forceComplexPs = forceComplexPs;
233 psCtrl.maxIts = maxIts;
234 psCtrl.tol = psTol;
235 psCtrl.deflate = deflate;
236 psCtrl.arnoldi = arnoldi;
237 psCtrl.basisSize = basisSize;
238 psCtrl.progress = progress;
239 #ifndef ELEM_HAVE_SCALAPACK
240 psCtrl.sdcCtrl.cutoff = cutoff;
241 psCtrl.sdcCtrl.maxInnerIts = maxInnerIts;
242 psCtrl.sdcCtrl.maxOuterIts = maxOuterIts;
243 psCtrl.sdcCtrl.tol = sdcTol;
244 psCtrl.sdcCtrl.spreadFactor = spreadFactor;
245 psCtrl.sdcCtrl.random = random;
246 psCtrl.sdcCtrl.progress = progress;
247 psCtrl.sdcCtrl.signCtrl.tol = signTol;
248 psCtrl.sdcCtrl.signCtrl.progress = progress;
249 #endif
250 psCtrl.snapCtrl.imgSaveFreq = imgSaveFreq;
251 psCtrl.snapCtrl.numSaveFreq = numSaveFreq;
252 psCtrl.snapCtrl.imgDispFreq = imgDispFreq;
253 psCtrl.snapCtrl.imgFormat = imgFormat;
254 psCtrl.snapCtrl.numFormat = numFormat;
255 psCtrl.snapCtrl.imgBase = matName+"-"+imgBase;
256 psCtrl.snapCtrl.numBase = matName+"-"+numBase;
257 psCtrl.snapCtrl.itCounts = itCounts;
258
259 // Visualize the pseudospectrum by evaluating ||inv(A-sigma I)||_2
260 // for a grid of complex sigma's.
261 DistMatrix<Real> invNormMap(g);
262 DistMatrix<Int> itCountMap(g);
263 if( realWidth != 0. && imagWidth != 0. )
264 {
265 if( isReal )
266 itCountMap = Pseudospectrum
267 ( AReal, invNormMap, center, realWidth, imagWidth,
268 realSize, imagSize, psCtrl );
269 else
270 itCountMap = Pseudospectrum
271 ( ACpx, invNormMap, center, realWidth, imagWidth,
272 realSize, imagSize, psCtrl );
273 }
274 else
275 {
276 if( isReal )
277 itCountMap = Pseudospectrum
278 ( AReal, invNormMap, center, realSize, imagSize, psCtrl );
279 else
280 itCountMap = Pseudospectrum
281 ( ACpx, invNormMap, center, realSize, imagSize, psCtrl );
282 }
283 const Int numIts = MaxNorm( itCountMap );
284 if( mpi::WorldRank() == 0 )
285 std::cout << "num iterations=" << numIts << std::endl;
286 }
287 catch( exception& e ) { ReportException(e); }
288
289 Finalize();
290 return 0;
291 }
292