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