IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
Ifpack_Utils.cpp
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack: Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2002) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #include "Ifpack_ConfigDefs.h"
44 #include "Ifpack_Preconditioner.h"
45 #include "Ifpack_Utils.h"
46 #include "Epetra_Comm.h"
47 #include "Epetra_CrsMatrix.h"
48 #include "Epetra_CrsGraph.h"
49 #include "Epetra_Map.h"
50 #include "Epetra_BlockMap.h"
51 #include "Epetra_Import.h"
52 #include "Epetra_MultiVector.h"
53 #include "Epetra_Vector.h"
54 
56 {
57  std::cout << "================================================================================" << std::endl;
58 }
59 
60 //============================================================================
62 {
63  using std::cout;
64  using std::endl;
65 
66  char hostname[80];
67  char buf[160];
68  if (Comm.MyPID() == 0) cout << "Host and Process Ids for tasks" << endl;
69  for (int i = 0; i <Comm.NumProc() ; i++) {
70  if (i == Comm.MyPID() ) {
71 #if defined(TFLOP) || defined(JANUS_STLPORT)
72  sprintf(buf, "Host: %s PID: %d", "janus", getpid());
73 #elif defined(_WIN32)
74  sprintf(buf,"Windows compiler, unknown hostname and PID!");
75 #else
76  gethostname(hostname, sizeof(hostname));
77  sprintf(buf, "Host: %s\tComm.MyPID(): %d\tPID: %d",
78  hostname, Comm.MyPID(), getpid());
79 #endif
80  printf("%s\n",buf);
81  fflush(stdout);
82 #if !( defined(_WIN32) )
83  sleep(1);
84 #endif
85  }
86  }
87  if(Comm.MyPID() == 0) {
88  printf("\n");
89  printf("** Pausing to attach debugger...\n");
90  printf("** You may now attach debugger to the processes listed above.\n");
91  printf( "**\n");
92  printf( "** Enter a character to continue > "); fflush(stdout);
93  char go;
94  TEUCHOS_ASSERT(scanf("%c",&go) != EOF);
95  }
96 
97  Comm.Barrier();
98 
99 }
100 
101 //============================================================================
103  const int OverlappingLevel)
104 {
105 
106  if (OverlappingLevel == 0)
107  return(0); // All done
108  if (Matrix->Comm().NumProc() == 1)
109  return(0); // All done
110 
111  Epetra_CrsMatrix* OverlappingMatrix;
112  OverlappingMatrix = 0;
113  Epetra_Map* OverlappingMap;
114  OverlappingMap = (Epetra_Map*)&(Matrix->RowMatrixRowMap());
115 
116  const Epetra_RowMatrix* OldMatrix;
117  const Epetra_Map* DomainMap = &(Matrix->OperatorDomainMap());
118  const Epetra_Map* RangeMap = &(Matrix->OperatorRangeMap());
119 
120  for (int level = 1; level <= OverlappingLevel ; ++level) {
121 
122  if (OverlappingMatrix)
123  OldMatrix = OverlappingMatrix;
124  else
125  OldMatrix = Matrix;
126 
127  Epetra_Import* OverlappingImporter;
128  OverlappingImporter = (Epetra_Import*)OldMatrix->RowMatrixImporter();
129  int NumMyElements = OverlappingImporter->TargetMap().NumMyElements();
130 
131  // need to build an Epetra_Map in this way because Epetra_CrsMatrix
132  // requires Epetra_Map and not Epetra_BlockMap
133 
134 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
135  if(OverlappingImporter->TargetMap().GlobalIndicesInt()) {
136  int* MyGlobalElements = OverlappingImporter->TargetMap().MyGlobalElements();
137  OverlappingMap = new Epetra_Map(-1,NumMyElements,MyGlobalElements,
138  0, Matrix->Comm());
139  }
140  else
141 #endif
142 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
143  if(OverlappingImporter->TargetMap().GlobalIndicesLongLong()) {
144  long long* MyGlobalElements = OverlappingImporter->TargetMap().MyGlobalElements64();
145  OverlappingMap = new Epetra_Map((long long) -1,NumMyElements,MyGlobalElements,
146  0, Matrix->Comm());
147  }
148  else
149 #endif
150  throw "Ifpack_CreateOverlappingCrsMatrix: GlobalIndices type unknown";
151 
152  if (level < OverlappingLevel)
153  OverlappingMatrix = new Epetra_CrsMatrix(Copy, *OverlappingMap, 0);
154  else
155  // On last iteration, we want to filter out all columns except
156  // those that correspond
157  // to rows in the graph. This assures that our matrix is square
158  OverlappingMatrix = new Epetra_CrsMatrix(Copy, *OverlappingMap,
159  *OverlappingMap, 0);
160 
161  OverlappingMatrix->Import(*OldMatrix, *OverlappingImporter, Insert);
162  if (level < OverlappingLevel) {
163  OverlappingMatrix->FillComplete(*DomainMap, *RangeMap);
164  }
165  else {
166  OverlappingMatrix->FillComplete(*DomainMap, *RangeMap);
167  }
168 
169  delete OverlappingMap;
170 
171  if (level > 1) {
172  delete OldMatrix;
173  }
174  OverlappingMatrix->FillComplete();
175 
176  }
177 
178  return(OverlappingMatrix);
179 }
180 
181 //============================================================================
183  const int OverlappingLevel)
184 {
185 
186  if (OverlappingLevel == 0)
187  return(0); // All done
188  if (Graph->Comm().NumProc() == 1)
189  return(0); // All done
190 
191  Epetra_CrsGraph* OverlappingGraph;
192  Epetra_BlockMap* OverlappingMap;
193  OverlappingGraph = const_cast<Epetra_CrsGraph*>(Graph);
194  OverlappingMap = const_cast<Epetra_BlockMap*>(&(Graph->RowMap()));
195 
196  Epetra_CrsGraph* OldGraph;
197  Epetra_BlockMap* OldMap;
198  const Epetra_BlockMap* DomainMap = &(Graph->DomainMap());
199  const Epetra_BlockMap* RangeMap = &(Graph->RangeMap());
200 
201  for (int level = 1; level <= OverlappingLevel ; ++level) {
202 
203  OldGraph = OverlappingGraph;
204  OldMap = OverlappingMap;
205 
206  Epetra_Import* OverlappingImporter;
207  OverlappingImporter = const_cast<Epetra_Import*>(OldGraph->Importer());
208  OverlappingMap = new Epetra_BlockMap(OverlappingImporter->TargetMap());
209 
210  if (level < OverlappingLevel)
211  OverlappingGraph = new Epetra_CrsGraph(Copy, *OverlappingMap, 0);
212  else
213  // On last iteration, we want to filter out all columns except
214  // those that correspond
215  // to rows in the graph. This assures that our matrix is square
216  OverlappingGraph = new Epetra_CrsGraph(Copy, *OverlappingMap,
217  *OverlappingMap, 0);
218 
219  OverlappingGraph->Import(*OldGraph, *OverlappingImporter, Insert);
220  if (level < OverlappingLevel)
221  OverlappingGraph->FillComplete(*DomainMap, *RangeMap);
222  else {
223  // Copy last OverlapImporter because we will use it later
224  OverlappingImporter = new Epetra_Import(*OverlappingMap, *DomainMap);
225  OverlappingGraph->FillComplete(*DomainMap, *RangeMap);
226  }
227 
228  if (level > 1) {
229  delete OldGraph;
230  delete OldMap;
231  }
232 
233  delete OverlappingMap;
234  OverlappingGraph->FillComplete();
235 
236  }
237 
238  return(OverlappingGraph);
239 }
240 
241 //============================================================================
242 std::string Ifpack_toString(const int& x)
243 {
244  char s[100];
245  sprintf(s, "%d", x);
246  return std::string(s);
247 }
248 
249 //============================================================================
250 std::string Ifpack_toString(const double& x)
251 {
252  char s[100];
253  sprintf(s, "%g", x);
254  return std::string(s);
255 }
256 
257 //============================================================================
258 int Ifpack_PrintResidual(char* Label, const Epetra_RowMatrix& A,
259  const Epetra_MultiVector& X, const Epetra_MultiVector&Y)
260 {
261  using std::cout;
262  using std::endl;
263 
264  if (X.Comm().MyPID() == 0) {
265  cout << "***** " << Label << endl;
266  }
267  Ifpack_PrintResidual(0,A,X,Y);
268 
269  return(0);
270 }
271 
272 //============================================================================
273 int Ifpack_PrintResidual(const int iter, const Epetra_RowMatrix& A,
274  const Epetra_MultiVector& X, const Epetra_MultiVector&Y)
275 {
276  using std::cout;
277  using std::endl;
278 
279  Epetra_MultiVector RHS(X);
280  std::vector<double> Norm2;
281  Norm2.resize(X.NumVectors());
282 
283  IFPACK_CHK_ERR(A.Multiply(false,X,RHS));
284  RHS.Update(1.0, Y, -1.0);
285 
286  RHS.Norm2(&Norm2[0]);
287 
288  if (X.Comm().MyPID() == 0) {
289  cout << "***** iter: " << iter << ": ||Ax - b||_2 = "
290  << Norm2[0] << endl;
291  }
292 
293  return(0);
294 }
295 
296 //============================================================================
297 // very simple debugging function that prints on screen the structure
298 // of the input matrix.
299 void Ifpack_PrintSparsity_Simple(const Epetra_RowMatrix& A)
300 {
301  using std::cout;
302  using std::endl;
303 
304  int MaxEntries = A.MaxNumEntries();
305  std::vector<int> Indices(MaxEntries);
306  std::vector<double> Values(MaxEntries);
307  std::vector<bool> FullRow(A.NumMyRows());
308 
309  cout << "+-";
310  for (int j = 0 ; j < A.NumMyRows() ; ++j)
311  cout << '-';
312  cout << "-+" << endl;
313 
314  for (int i = 0 ; i < A.NumMyRows() ; ++i) {
315 
316  int Length;
317  A.ExtractMyRowCopy(i,MaxEntries,Length,
318  &Values[0], &Indices[0]);
319 
320  for (int j = 0 ; j < A.NumMyRows() ; ++j)
321  FullRow[j] = false;
322 
323  for (int j = 0 ; j < Length ; ++j) {
324  FullRow[Indices[j]] = true;
325  }
326 
327  cout << "| ";
328  for (int j = 0 ; j < A.NumMyRows() ; ++j) {
329  if (FullRow[j])
330  cout << '*';
331  else
332  cout << ' ';
333  }
334  cout << " |" << endl;
335  }
336 
337  cout << "+-";
338  for (int j = 0 ; j < A.NumMyRows() ; ++j)
339  cout << '-';
340  cout << "-+" << endl << endl;
341 
342 }
343 
344 //============================================================================
345 
346 double Ifpack_FrobeniusNorm(const Epetra_RowMatrix& A)
347 {
348  double MyNorm = 0.0, GlobalNorm;
349 
350  std::vector<int> colInd(A.MaxNumEntries());
351  std::vector<double> colVal(A.MaxNumEntries());
352 
353  for (int i = 0 ; i < A.NumMyRows() ; ++i) {
354 
355  int Nnz;
356  IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
357  &colVal[0],&colInd[0]));
358 
359  for (int j = 0 ; j < Nnz ; ++j) {
360  MyNorm += colVal[j] * colVal[j];
361  }
362  }
363 
364  A.Comm().SumAll(&MyNorm,&GlobalNorm,1);
365 
366  return(sqrt(GlobalNorm));
367 }
368 
369 static void print()
370 {
371  printf("\n");
372 }
373 
374 #include <iomanip>
375 template<class T>
376 static void print(const char str[], T val)
377 {
378  std::cout.width(30); std::cout.setf(std::ios::left);
379  std::cout << str;
380  std::cout << " = " << val << std::endl;
381 }
382 
383 template<class T>
384 static void print(const char str[], T val, double percentage)
385 {
386  std::cout.width(30); std::cout.setf(std::ios::left);
387  std::cout << str;
388  std::cout << " = ";
389  std::cout.width(20); std::cout.setf(std::ios::left);
390  std::cout << val;
391  std::cout << " ( " << percentage << " %)" << std::endl;
392 }
393 template<class T>
394 static void print(const char str[], T one, T two, T three, bool equal = true)
395 {
396  using std::endl;
397 
398  std::cout.width(30); std::cout.setf(std::ios::left);
399  std::cout << str;
400  if (equal)
401  std::cout << " = ";
402  else
403  std::cout << " ";
404  std::cout.width(15); std::cout.setf(std::ios::left);
405  std::cout << one;
406  std::cout.width(15); std::cout.setf(std::ios::left);
407  std::cout << two;
408  std::cout.width(15); std::cout.setf(std::ios::left);
409  std::cout << three;
410  std::cout << endl;
411 }
412 
413 //============================================================================
414 #include "limits.h"
415 #include "float.h"
416 #include "Epetra_FECrsMatrix.h"
417 
418 int Ifpack_Analyze(const Epetra_RowMatrix& A, const bool Cheap,
419  const int NumPDEEqns)
420 {
421 
422  int NumMyRows = A.NumMyRows();
423  long long NumGlobalRows = A.NumGlobalRows64();
424  long long NumGlobalCols = A.NumGlobalCols64();
425  long long MyBandwidth = 0, GlobalBandwidth;
426  long long MyLowerNonzeros = 0, MyUpperNonzeros = 0;
427  long long GlobalLowerNonzeros, GlobalUpperNonzeros;
428  long long MyDiagonallyDominant = 0, GlobalDiagonallyDominant;
429  long long MyWeaklyDiagonallyDominant = 0, GlobalWeaklyDiagonallyDominant;
430  double MyMin, MyAvg, MyMax;
431  double GlobalMin, GlobalAvg, GlobalMax;
432  long long GlobalStorage;
433 
434  bool verbose = (A.Comm().MyPID() == 0);
435 
436  GlobalStorage = sizeof(int*) * NumGlobalRows +
437  sizeof(int) * A.NumGlobalNonzeros64() +
438  sizeof(double) * A.NumGlobalNonzeros64();
439 
440  if (verbose) {
441  print();
443  print<const char*>("Label", A.Label());
444  print<long long>("Global rows", NumGlobalRows);
445  print<long long>("Global columns", NumGlobalCols);
446  print<long long>("Stored nonzeros", A.NumGlobalNonzeros64());
447  print<long long>("Nonzeros / row", A.NumGlobalNonzeros64() / NumGlobalRows);
448  print<double>("Estimated storage (Mbytes)", 1.0e-6 * GlobalStorage);
449  }
450 
451  long long NumMyActualNonzeros = 0, NumGlobalActualNonzeros;
452  long long NumMyEmptyRows = 0, NumGlobalEmptyRows;
453  long long NumMyDirichletRows = 0, NumGlobalDirichletRows;
454 
455  std::vector<int> colInd(A.MaxNumEntries());
456  std::vector<double> colVal(A.MaxNumEntries());
457 
458  Epetra_Vector Diag(A.RowMatrixRowMap());
459  Epetra_Vector RowSum(A.RowMatrixRowMap());
460  Diag.PutScalar(0.0);
461  RowSum.PutScalar(0.0);
462 
463  for (int i = 0 ; i < NumMyRows ; ++i) {
464 
465  long long GRID = A.RowMatrixRowMap().GID64(i);
466  int Nnz;
467  IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
468  &colVal[0],&colInd[0]));
469 
470  if (Nnz == 0)
471  NumMyEmptyRows++;
472 
473  if (Nnz == 1)
474  NumMyDirichletRows++;
475 
476  for (int j = 0 ; j < Nnz ; ++j) {
477 
478  double v = colVal[j];
479  if (v < 0) v = -v;
480  if (colVal[j] != 0.0)
481  NumMyActualNonzeros++;
482 
483  long long GCID = A.RowMatrixColMap().GID64(colInd[j]);
484 
485  if (GCID != GRID)
486  RowSum[i] += v;
487  else
488  Diag[i] = v;
489 
490  if (GCID < GRID)
491  MyLowerNonzeros++;
492  else if (GCID > GRID)
493  MyUpperNonzeros++;
494  long long b = GCID - GRID;
495  if (b < 0) b = -b;
496  if (b > MyBandwidth)
497  MyBandwidth = b;
498  }
499 
500  if (Diag[i] > RowSum[i])
501  MyDiagonallyDominant++;
502 
503  if (Diag[i] >= RowSum[i])
504  MyWeaklyDiagonallyDominant++;
505 
506  RowSum[i] += Diag[i];
507  }
508 
509  // ======================== //
510  // summing up global values //
511  // ======================== //
512 
513  A.Comm().SumAll(&MyDiagonallyDominant,&GlobalDiagonallyDominant,1);
514  A.Comm().SumAll(&MyWeaklyDiagonallyDominant,&GlobalWeaklyDiagonallyDominant,1);
515  A.Comm().SumAll(&NumMyActualNonzeros, &NumGlobalActualNonzeros, 1);
516  A.Comm().SumAll(&NumMyEmptyRows, &NumGlobalEmptyRows, 1);
517  A.Comm().SumAll(&NumMyDirichletRows, &NumGlobalDirichletRows, 1);
518  A.Comm().SumAll(&MyBandwidth, &GlobalBandwidth, 1);
519  A.Comm().SumAll(&MyLowerNonzeros, &GlobalLowerNonzeros, 1);
520  A.Comm().SumAll(&MyUpperNonzeros, &GlobalUpperNonzeros, 1);
521  A.Comm().SumAll(&MyDiagonallyDominant, &GlobalDiagonallyDominant, 1);
522  A.Comm().SumAll(&MyWeaklyDiagonallyDominant, &GlobalWeaklyDiagonallyDominant, 1);
523 
524  double NormOne = A.NormOne();
525  double NormInf = A.NormInf();
526  double NormF = Ifpack_FrobeniusNorm(A);
527 
528  if (verbose) {
529  print();
530  print<long long>("Actual nonzeros", NumGlobalActualNonzeros);
531  print<long long>("Nonzeros in strict lower part", GlobalLowerNonzeros);
532  print<long long>("Nonzeros in strict upper part", GlobalUpperNonzeros);
533  print();
534  print<long long>("Empty rows", NumGlobalEmptyRows,
535  100.0 * NumGlobalEmptyRows / NumGlobalRows);
536  print<long long>("Dirichlet rows", NumGlobalDirichletRows,
537  100.0 * NumGlobalDirichletRows / NumGlobalRows);
538  print<long long>("Diagonally dominant rows", GlobalDiagonallyDominant,
539  100.0 * GlobalDiagonallyDominant / NumGlobalRows);
540  print<long long>("Weakly diag. dominant rows",
541  GlobalWeaklyDiagonallyDominant,
542  100.0 * GlobalWeaklyDiagonallyDominant / NumGlobalRows);
543  print();
544  print<long long>("Maximum bandwidth", GlobalBandwidth);
545 
546  print();
547  print("", "one-norm", "inf-norm", "Frobenius", false);
548  print("", "========", "========", "=========", false);
549  print();
550 
551  print<double>("A", NormOne, NormInf, NormF);
552  }
553 
554  if (Cheap == false) {
555 
556  // create A + A^T and A - A^T
557 
558  Epetra_FECrsMatrix AplusAT(Copy, A.RowMatrixRowMap(), 0);
559  Epetra_FECrsMatrix AminusAT(Copy, A.RowMatrixRowMap(), 0);
560 
561 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
563  for (int i = 0 ; i < NumMyRows ; ++i) {
564 
565  int GRID = A.RowMatrixRowMap().GID(i);
566  assert (GRID != -1);
567 
568  int Nnz;
569  IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
570  &colVal[0],&colInd[0]));
571 
572  for (int j = 0 ; j < Nnz ; ++j) {
573 
574  int GCID = A.RowMatrixColMap().GID(colInd[j]);
575  assert (GCID != -1);
576 
577  double plus_val = colVal[j];
578  double minus_val = -colVal[j];
579 
580  if (AplusAT.SumIntoGlobalValues(1,&GRID,1,&GCID,&plus_val) != 0) {
581  IFPACK_CHK_ERR(AplusAT.InsertGlobalValues(1,&GRID,1,&GCID,&plus_val));
582  }
583 
584  if (AplusAT.SumIntoGlobalValues(1,&GCID,1,&GRID,&plus_val) != 0) {
585  IFPACK_CHK_ERR(AplusAT.InsertGlobalValues(1,&GCID,1,&GRID,&plus_val));
586  }
587 
588  if (AminusAT.SumIntoGlobalValues(1,&GRID,1,&GCID,&plus_val) != 0) {
589  IFPACK_CHK_ERR(AminusAT.InsertGlobalValues(1,&GRID,1,&GCID,&plus_val));
590  }
591 
592  if (AminusAT.SumIntoGlobalValues(1,&GCID,1,&GRID,&minus_val) != 0) {
593  IFPACK_CHK_ERR(AminusAT.InsertGlobalValues(1,&GCID,1,&GRID,&minus_val));
594  }
595 
596  }
597  }
598  }
599  else
600 #endif
601 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
603  for (int i = 0 ; i < NumMyRows ; ++i) {
604 
605  long long GRID = A.RowMatrixRowMap().GID64(i);
606  assert (GRID != -1);
607 
608  int Nnz;
609  IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
610  &colVal[0],&colInd[0]));
611 
612  for (int j = 0 ; j < Nnz ; ++j) {
613 
614  long long GCID = A.RowMatrixColMap().GID64(colInd[j]);
615  assert (GCID != -1);
616 
617  double plus_val = colVal[j];
618  double minus_val = -colVal[j];
619 
620  if (AplusAT.SumIntoGlobalValues(1,&GRID,1,&GCID,&plus_val) != 0) {
621  IFPACK_CHK_ERR(AplusAT.InsertGlobalValues(1,&GRID,1,&GCID,&plus_val));
622  }
623 
624  if (AplusAT.SumIntoGlobalValues(1,&GCID,1,&GRID,&plus_val) != 0) {
625  IFPACK_CHK_ERR(AplusAT.InsertGlobalValues(1,&GCID,1,&GRID,&plus_val));
626  }
627 
628  if (AminusAT.SumIntoGlobalValues(1,&GRID,1,&GCID,&plus_val) != 0) {
629  IFPACK_CHK_ERR(AminusAT.InsertGlobalValues(1,&GRID,1,&GCID,&plus_val));
630  }
631 
632  if (AminusAT.SumIntoGlobalValues(1,&GCID,1,&GRID,&minus_val) != 0) {
633  IFPACK_CHK_ERR(AminusAT.InsertGlobalValues(1,&GCID,1,&GRID,&minus_val));
634  }
635 
636  }
637  }
638  }
639  else
640 #endif
641  throw "Ifpack_Analyze: GlobalIndices type unknown";
642 
643  AplusAT.FillComplete();
644  AminusAT.FillComplete();
645 
646  AplusAT.Scale(0.5);
647  AminusAT.Scale(0.5);
648 
649  NormOne = AplusAT.NormOne();
650  NormInf = AplusAT.NormInf();
651  NormF = Ifpack_FrobeniusNorm(AplusAT);
652 
653  if (verbose) {
654  print<double>("A + A^T", NormOne, NormInf, NormF);
655  }
656 
657  NormOne = AminusAT.NormOne();
658  NormInf = AminusAT.NormInf();
659  NormF = Ifpack_FrobeniusNorm(AminusAT);
660 
661  if (verbose) {
662  print<double>("A - A^T", NormOne, NormInf, NormF);
663  }
664  }
665 
666  if (verbose) {
667  print();
668  print<const char*>("", "min", "avg", "max", false);
669  print<const char*>("", "===", "===", "===", false);
670  }
671 
672  MyMax = -DBL_MAX;
673  MyMin = DBL_MAX;
674  MyAvg = 0.0;
675 
676  for (int i = 0 ; i < NumMyRows ; ++i) {
677 
678  int Nnz;
679  IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
680  &colVal[0],&colInd[0]));
681 
682  for (int j = 0 ; j < Nnz ; ++j) {
683  MyAvg += colVal[j];
684  if (colVal[j] > MyMax) MyMax = colVal[j];
685  if (colVal[j] < MyMin) MyMin = colVal[j];
686  }
687  }
688 
689  A.Comm().MaxAll(&MyMax, &GlobalMax, 1);
690  A.Comm().MinAll(&MyMin, &GlobalMin, 1);
691  A.Comm().SumAll(&MyAvg, &GlobalAvg, 1);
692  GlobalAvg /= A.NumGlobalNonzeros64();
693 
694  if (verbose) {
695  print();
696  print<double>(" A(i,j)", GlobalMin, GlobalAvg, GlobalMax);
697  }
698 
699  MyMax = 0.0;
700  MyMin = DBL_MAX;
701  MyAvg = 0.0;
702 
703  for (int i = 0 ; i < NumMyRows ; ++i) {
704 
705  int Nnz;
706  IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
707  &colVal[0],&colInd[0]));
708 
709  for (int j = 0 ; j < Nnz ; ++j) {
710  double v = colVal[j];
711  if (v < 0) v = -v;
712  MyAvg += v;
713  if (colVal[j] > MyMax) MyMax = v;
714  if (colVal[j] < MyMin) MyMin = v;
715  }
716  }
717 
718  A.Comm().MaxAll(&MyMax, &GlobalMax, 1);
719  A.Comm().MinAll(&MyMin, &GlobalMin, 1);
720  A.Comm().SumAll(&MyAvg, &GlobalAvg, 1);
721  GlobalAvg /= A.NumGlobalNonzeros64();
722 
723  if (verbose) {
724  print<double>("|A(i,j)|", GlobalMin, GlobalAvg, GlobalMax);
725  }
726 
727  // ================= //
728  // diagonal elements //
729  // ================= //
730 
731  Diag.MinValue(&GlobalMin);
732  Diag.MaxValue(&GlobalMax);
733  Diag.MeanValue(&GlobalAvg);
734 
735  if (verbose) {
736  print();
737  print<double>(" A(k,k)", GlobalMin, GlobalAvg, GlobalMax);
738  }
739 
740  Diag.Abs(Diag);
741  Diag.MinValue(&GlobalMin);
742  Diag.MaxValue(&GlobalMax);
743  Diag.MeanValue(&GlobalAvg);
744  if (verbose) {
745  print<double>("|A(k,k)|", GlobalMin, GlobalAvg, GlobalMax);
746  }
747 
748  // ============================================== //
749  // cycle over all equations for diagonal elements //
750  // ============================================== //
751 
752  if (NumPDEEqns > 1 ) {
753 
754  if (verbose) print();
755 
756  for (int ie = 0 ; ie < NumPDEEqns ; ie++) {
757 
758  MyMin = DBL_MAX;
759  MyMax = -DBL_MAX;
760  MyAvg = 0.0;
761 
762  for (int i = ie ; i < Diag.MyLength() ; i += NumPDEEqns) {
763  double d = Diag[i];
764  MyAvg += d;
765  if (d < MyMin)
766  MyMin = d;
767  if (d > MyMax)
768  MyMax = d;
769  }
770  A.Comm().MinAll(&MyMin, &GlobalMin, 1);
771  A.Comm().MaxAll(&MyMax, &GlobalMax, 1);
772  A.Comm().SumAll(&MyAvg, &GlobalAvg, 1);
773  // does not really work fine if the number of global
774  // elements is not a multiple of NumPDEEqns
775  GlobalAvg /= (Diag.GlobalLength64() / NumPDEEqns);
776 
777  if (verbose) {
778  char str[80];
779  sprintf(str, " A(k,k), eq %d", ie);
780  print<double>(str, GlobalMin, GlobalAvg, GlobalMax);
781  }
782  }
783  }
784 
785  // ======== //
786  // row sums //
787  // ======== //
788 
789  RowSum.MinValue(&GlobalMin);
790  RowSum.MaxValue(&GlobalMax);
791  RowSum.MeanValue(&GlobalAvg);
792 
793  if (verbose) {
794  print();
795  print<double>(" sum_j A(k,j)", GlobalMin, GlobalAvg, GlobalMax);
796  }
797 
798  // ===================================== //
799  // cycle over all equations for row sums //
800  // ===================================== //
801 
802  if (NumPDEEqns > 1 ) {
803 
804  if (verbose) print();
805 
806  for (int ie = 0 ; ie < NumPDEEqns ; ie++) {
807 
808  MyMin = DBL_MAX;
809  MyMax = -DBL_MAX;
810  MyAvg = 0.0;
811 
812  for (int i = ie ; i < Diag.MyLength() ; i += NumPDEEqns) {
813  double d = RowSum[i];
814  MyAvg += d;
815  if (d < MyMin)
816  MyMin = d;
817  if (d > MyMax)
818  MyMax = d;
819  }
820  A.Comm().MinAll(&MyMin, &GlobalMin, 1);
821  A.Comm().MaxAll(&MyMax, &GlobalMax, 1);
822  A.Comm().SumAll(&MyAvg, &GlobalAvg, 1);
823  // does not really work fine if the number of global
824  // elements is not a multiple of NumPDEEqns
825  GlobalAvg /= (Diag.GlobalLength64() / NumPDEEqns);
826 
827  if (verbose) {
828  char str[80];
829  sprintf(str, " sum_j A(k,j), eq %d", ie);
830  print<double>(str, GlobalMin, GlobalAvg, GlobalMax);
831  }
832  }
833  }
834 
835  if (verbose)
837 
838  return(0);
839 }
840 
842  const bool abs, const int steps)
843 {
844  using std::cout;
845  using std::endl;
846 
847  bool verbose = (Diagonal.Comm().MyPID() == 0);
848  double min_val = DBL_MAX;
849  double max_val = -DBL_MAX;
850 
851  for (int i = 0 ; i < Diagonal.MyLength() ; ++i) {
852  double v = Diagonal[i];
853  if (abs)
854  if (v < 0) v = -v;
855  if (v > max_val)
856  max_val = v;
857  if (v < min_val)
858  min_val = v;
859  }
860 
861  if (verbose) {
862  cout << endl;
864  cout << "Vector label = " << Diagonal.Label() << endl;
865  cout << endl;
866  }
867 
868  double delta = (max_val - min_val) / steps;
869  for (int k = 0 ; k < steps ; ++k) {
870 
871  double below = delta * k + min_val;
872  double above = below + delta;
873  int MyBelow = 0, GlobalBelow;
874 
875  for (int i = 0 ; i < Diagonal.MyLength() ; ++i) {
876  double v = Diagonal[i];
877  if (v < 0) v = -v;
878  if (v >= below && v < above) MyBelow++;
879  }
880 
881  Diagonal.Comm().SumAll(&MyBelow, &GlobalBelow, 1);
882 
883  if (verbose) {
884  printf("Elements in [%+7e, %+7e) = %10d ( = %5.2f %%)\n",
885  below, above, GlobalBelow,
886  100.0 * GlobalBelow / Diagonal.GlobalLength64());
887  }
888  }
889 
890  if (verbose) {
892  cout << endl;
893  }
894 
895  return(0);
896 }
897 
898 // ======================================================================
899 
901  const bool abs, const int steps)
902 {
903  using std::cout;
904  using std::endl;
905 
906  bool verbose = (A.Comm().MyPID() == 0);
907  double min_val = DBL_MAX;
908  double max_val = -DBL_MAX;
909 
910  std::vector<int> colInd(A.MaxNumEntries());
911  std::vector<double> colVal(A.MaxNumEntries());
912 
913  for (int i = 0 ; i < A.NumMyRows() ; ++i) {
914 
915  int Nnz;
916  IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
917  &colVal[0],&colInd[0]));
918 
919  for (int j = 0 ; j < Nnz ; ++j) {
920  double v = colVal[j];
921  if (abs)
922  if (v < 0) v = -v;
923  if (v < min_val)
924  min_val = v;
925  if (v > max_val)
926  max_val = v;
927  }
928  }
929 
930  if (verbose) {
931  cout << endl;
933  cout << "Label of matrix = " << A.Label() << endl;
934  cout << endl;
935  }
936 
937  double delta = (max_val - min_val) / steps;
938  for (int k = 0 ; k < steps ; ++k) {
939 
940  double below = delta * k + min_val;
941  double above = below + delta;
942  int MyBelow = 0, GlobalBelow;
943 
944  for (int i = 0 ; i < A.NumMyRows() ; ++i) {
945 
946  int Nnz;
947  IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
948  &colVal[0],&colInd[0]));
949 
950  for (int j = 0 ; j < Nnz ; ++j) {
951  double v = colVal[j];
952  if (abs)
953  if (v < 0) v = -v;
954  if (v >= below && v < above) MyBelow++;
955  }
956 
957  }
958  A.Comm().SumAll(&MyBelow, &GlobalBelow, 1);
959  if (verbose) {
960  printf("Elements in [%+7e, %+7e) = %10d ( = %5.2f %%)\n",
961  below, above, GlobalBelow,
962  100.0 * GlobalBelow / A.NumGlobalNonzeros64());
963  }
964  }
965 
966  if (verbose) {
968  cout << endl;
969  }
970 
971  return(0);
972 }
973 
974 // ======================================================================
975 int Ifpack_PrintSparsity(const Epetra_RowMatrix& A, const char* InputFileName,
976  const int NumPDEEqns)
977 {
978 
979  int ltit;
980  long long m,nc,nr,maxdim;
981  double lrmrgn,botmrgn,xtit,ytit,ytitof,fnstit,siz = 0.0;
982  double xl,xr, yb,yt, scfct,u2dot,frlw,delt,paperx;
983  bool square = false;
984  /*change square to .true. if you prefer a square frame around
985  a rectangular matrix */
986  double conv = 2.54;
987  char munt = 'E'; /* put 'E' for centimeters, 'U' for inches */
988  int ptitle = 0; /* position of the title, 0 under the drawing,
989  else above */
990  FILE* fp = NULL;
991  int NumMyRows;
992  //int NumMyCols;
993  long long NumGlobalRows;
994  long long NumGlobalCols;
995  int MyPID;
996  int NumProc;
997  char FileName[1024];
998  char title[1020];
999 
1000  const Epetra_Comm& Comm = A.Comm();
1001 
1002  /* --------------------- execution begins ---------------------- */
1003 
1004  if (strlen(A.Label()) != 0)
1005  strcpy(title, A.Label());
1006  else
1007  sprintf(title, "%s", "matrix");
1008 
1009  if (InputFileName == 0)
1010  sprintf(FileName, "%s.ps", title);
1011  else
1012  strcpy(FileName, InputFileName);
1013 
1014  MyPID = Comm.MyPID();
1015  NumProc = Comm.NumProc();
1016 
1017  NumMyRows = A.NumMyRows();
1018  //NumMyCols = A.NumMyCols();
1019 
1020  NumGlobalRows = A.NumGlobalRows64();
1021  NumGlobalCols = A.NumGlobalCols64();
1022 
1023  if (NumGlobalRows != NumGlobalCols)
1024  IFPACK_CHK_ERR(-1); // never tested
1025 
1026  /* to be changed for rect matrices */
1027  maxdim = (NumGlobalRows>NumGlobalCols)?NumGlobalRows:NumGlobalCols;
1028  maxdim /= NumPDEEqns;
1029 
1030  m = 1 + maxdim;
1031  nr = NumGlobalRows / NumPDEEqns + 1;
1032  nc = NumGlobalCols / NumPDEEqns + 1;
1033 
1034  if (munt == 'E') {
1035  u2dot = 72.0/conv;
1036  paperx = 21.0;
1037  siz = 10.0;
1038  }
1039  else {
1040  u2dot = 72.0;
1041  paperx = 8.5*conv;
1042  siz = siz*conv;
1043  }
1044 
1045  /* left and right margins (drawing is centered) */
1046 
1047  lrmrgn = (paperx-siz)/2.0;
1048 
1049  /* bottom margin : 2 cm */
1050 
1051  botmrgn = 2.0;
1052  /* c scaling factor */
1053  scfct = siz*u2dot/m;
1054  /* matrix frame line witdh */
1055  frlw = 0.25;
1056  /* font size for title (cm) */
1057  fnstit = 0.5;
1058  /* mfh 23 Jan 2013: title is always nonnull, since it's an array of
1059  fixed nonzero length. The 'if' test thus results in a compiler
1060  warning. */
1061  /*if (title) ltit = strlen(title);*/
1062  /*else ltit = 0;*/
1063  ltit = strlen(title);
1064 
1065  /* position of title : centered horizontally */
1066  /* at 1.0 cm vertically over the drawing */
1067  ytitof = 1.0;
1068  xtit = paperx/2.0;
1069  ytit = botmrgn+siz*nr/m + ytitof;
1070  /* almost exact bounding box */
1071  xl = lrmrgn*u2dot - scfct*frlw/2;
1072  xr = (lrmrgn+siz)*u2dot + scfct*frlw/2;
1073  yb = botmrgn*u2dot - scfct*frlw/2;
1074  yt = (botmrgn+siz*nr/m)*u2dot + scfct*frlw/2;
1075  if (ltit == 0) {
1076  yt = yt + (ytitof+fnstit*0.70)*u2dot;
1077  }
1078  /* add some room to bounding box */
1079  delt = 10.0;
1080  xl = xl-delt;
1081  xr = xr+delt;
1082  yb = yb-delt;
1083  yt = yt+delt;
1084 
1085  /* correction for title under the drawing */
1086  if ((ptitle == 0) && (ltit == 0)) {
1087  ytit = botmrgn + fnstit*0.3;
1088  botmrgn = botmrgn + ytitof + fnstit*0.7;
1089  }
1090 
1091  /* begin of output */
1092 
1093  if (MyPID == 0) {
1094 
1095  fp = fopen(FileName,"w");
1096 
1097  fprintf(fp,"%s","%%!PS-Adobe-2.0\n");
1098  fprintf(fp,"%s","%%Creator: IFPACK\n");
1099  fprintf(fp,"%%%%BoundingBox: %f %f %f %f\n",
1100  xl,yb,xr,yt);
1101  fprintf(fp,"%s","%%EndComments\n");
1102  fprintf(fp,"%s","/cm {72 mul 2.54 div} def\n");
1103  fprintf(fp,"%s","/mc {72 div 2.54 mul} def\n");
1104  fprintf(fp,"%s","/pnum { 72 div 2.54 mul 20 std::string ");
1105  fprintf(fp,"%s","cvs print ( ) print} def\n");
1106  fprintf(fp,"%s","/Cshow {dup stringwidth pop -2 div 0 rmoveto show} def\n");
1107 
1108  /* we leave margins etc. in cm so it is easy to modify them if
1109  needed by editing the output file */
1110  fprintf(fp,"%s","gsave\n");
1111  if (ltit != 0) {
1112  fprintf(fp,"/Helvetica findfont %e cm scalefont setfont\n",
1113  fnstit);
1114  fprintf(fp,"%f cm %f cm moveto\n",
1115  xtit,ytit);
1116  fprintf(fp,"(%s) Cshow\n", title);
1117  fprintf(fp,"%f cm %f cm translate\n",
1118  lrmrgn,botmrgn);
1119  }
1120  fprintf(fp,"%f cm %d div dup scale \n",
1121  siz, (int) m);
1122  /* draw a frame around the matrix */
1123 
1124  fprintf(fp,"%f setlinewidth\n",
1125  frlw);
1126  fprintf(fp,"%s","newpath\n");
1127  fprintf(fp,"%s","0 0 moveto ");
1128  if (square) {
1129  printf("------------------- %d\n", (int) m);
1130  fprintf(fp,"%d %d lineto\n",
1131  (int) m, 0);
1132  fprintf(fp,"%d %d lineto\n",
1133  (int) m, (int) m);
1134  fprintf(fp,"%d %d lineto\n",
1135  0, (int) m);
1136  }
1137  else {
1138  fprintf(fp,"%d %d lineto\n",
1139  (int) nc, 0);
1140  fprintf(fp,"%d %d lineto\n",
1141  (int) nc, (int) nr);
1142  fprintf(fp,"%d %d lineto\n",
1143  0, (int) nr);
1144  }
1145  fprintf(fp,"%s","closepath stroke\n");
1146 
1147  /* plotting loop */
1148 
1149  fprintf(fp,"%s","1 1 translate\n");
1150  fprintf(fp,"%s","0.8 setlinewidth\n");
1151  fprintf(fp,"%s","/p {moveto 0 -.40 rmoveto \n");
1152  fprintf(fp,"%s"," 0 .80 rlineto stroke} def\n");
1153 
1154  fclose(fp);
1155  }
1156 
1157  int MaxEntries = A.MaxNumEntries();
1158  std::vector<int> Indices(MaxEntries);
1159  std::vector<double> Values(MaxEntries);
1160 
1161  for (int pid = 0 ; pid < NumProc ; ++pid) {
1162 
1163  if (pid == MyPID) {
1164 
1165  fp = fopen(FileName,"a");
1166  TEUCHOS_ASSERT(fp != NULL);
1167 
1168  for (int i = 0 ; i < NumMyRows ; ++i) {
1169 
1170  if (i % NumPDEEqns) continue;
1171 
1172  int Nnz;
1173  A.ExtractMyRowCopy(i,MaxEntries,Nnz,&Values[0],&Indices[0]);
1174 
1175  long long grow = A.RowMatrixRowMap().GID64(i);
1176 
1177  for (int j = 0 ; j < Nnz ; ++j) {
1178  int col = Indices[j];
1179  if (col % NumPDEEqns == 0) {
1180  long long gcol = A.RowMatrixColMap().GID64(Indices[j]);
1181  grow /= NumPDEEqns;
1182  gcol /= NumPDEEqns;
1183  fprintf(fp,"%lld %lld p\n",
1184  gcol, NumGlobalRows - grow - 1);
1185  }
1186  }
1187  }
1188 
1189  fprintf(fp,"%s","%end of data for this process\n");
1190 
1191  if( pid == NumProc - 1 )
1192  fprintf(fp,"%s","showpage\n");
1193 
1194  fclose(fp);
1195  }
1196  Comm.Barrier();
1197  }
1198 
1199  return(0);
1200 }
const Epetra_BlockMap & RangeMap() const
virtual const Epetra_Map & RowMatrixRowMap() const =0
const Epetra_Comm & Comm() const
virtual double NormOne() const =0
bool GlobalIndicesLongLong() const
int MyGlobalElements(int *MyGlobalElementList) const
const Epetra_Import * Importer() const
int Ifpack_Analyze(const Epetra_RowMatrix &A, const bool Cheap=false, const int NumPDEEqns=1)
Analyzes the basic properties of the input matrix A; see Usage of Ifpack_Analyze()..
virtual const Epetra_Map & OperatorDomainMap() const =0
virtual const char * Label() const =0
bool GlobalIndicesInt() const
const Epetra_BlockMap & DomainMap() const
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
virtual void Barrier() const =0
virtual int MyPID() const =0
int FillComplete(bool OptimizeDataStorage=true)
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
int NumMyElements() const
virtual int MaxNumEntries() const =0
virtual const Epetra_Map & OperatorRangeMap() const =0
virtual const Epetra_Comm & Comm() const =0
virtual const Epetra_Import * RowMatrixImporter() const =0
const Epetra_BlockMap & TargetMap() const
int GID(int LID) const
virtual double NormInf() const =0
virtual int NumMyRows() const =0
const Epetra_BlockMap & RowMap() const
int Ifpack_AnalyzeVectorElements(const Epetra_Vector &Diagonal, const bool abs=false, const int steps=10)
Analyzes the distribution of values of the input vector Diagonal.
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
virtual int NumProc() const =0
int Ifpack_AnalyzeMatrixElements(const Epetra_RowMatrix &A, const bool abs=false, const int steps=10)
Analyzes the distribution of values of the input matrix A.
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
virtual const Epetra_Map & RowMatrixColMap() const =0
void Ifpack_PrintLine()
Prints a line of `=&#39; on cout.
Epetra_CrsMatrix * Ifpack_CreateOverlappingCrsMatrix(const Epetra_RowMatrix *Matrix, const int OverlappingLevel)
Creates an overlapping Epetra_CrsMatrix. Returns 0 if OverlappingLevel is 0.
void Ifpack_BreakForDebugger(Epetra_Comm &Comm)
Stops the execution of code, so that a debugger can be attached.
int Ifpack_PrintResidual(char *Label, const Epetra_RowMatrix &A, const Epetra_MultiVector &X, const Epetra_MultiVector &Y)
Prints on cout the true residual.
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.