Ifpack_Utils.cpp

Go to the documentation of this file.
00001 #include "Ifpack_ConfigDefs.h"
00002 #include "Ifpack_Preconditioner.h"
00003 #include "Ifpack_Utils.h"
00004 #include "Epetra_Comm.h"
00005 #include "Epetra_CrsMatrix.h"
00006 #include "Epetra_CrsGraph.h"
00007 #include "Epetra_Map.h"
00008 #include "Epetra_BlockMap.h"
00009 #include "Epetra_Import.h"
00010 #include "Epetra_MultiVector.h"
00011 #include "Epetra_Vector.h"
00012 
00013 void Ifpack_PrintLine()
00014 {
00015   cout << "================================================================================" << endl;
00016 }
00017 
00018 //============================================================================
00019 void Ifpack_BreakForDebugger(Epetra_Comm& Comm)
00020 {
00021   char hostname[80];
00022   char buf[80];
00023   if (Comm.MyPID()  == 0) cout << "Host and Process Ids for tasks" << endl;
00024   for (int i = 0; i <Comm.NumProc() ; i++) {
00025     if (i == Comm.MyPID() ) {
00026 #if defined(TFLOP) || defined(JANUS_STLPORT)
00027       sprintf(buf, "Host: %s   PID: %d", "janus", getpid());
00028 #elif defined(__INTEL_COMPILER) && defined(_WIN32)
00029       sprintf(buf,"Intel compiler, unknown hostname and PID!");
00030 #else
00031       gethostname(hostname, sizeof(hostname));
00032       sprintf(buf, "Host: %s\tComm.MyPID(): %d\tPID: %d",
00033               hostname, Comm.MyPID(), getpid());
00034 #endif
00035       printf("%s\n",buf);
00036       fflush(stdout);
00037 #if !( defined(__INTEL_COMPILER) && defined(_WIN32) )
00038       sleep(1);
00039 #endif
00040     }
00041   }
00042   if(Comm.MyPID() == 0) {
00043     printf("\n");
00044     printf("** Pausing to attach debugger...\n");
00045     printf("** You may now attach debugger to the processes listed above.\n");
00046     printf( "**\n");
00047     printf( "** Enter a character to continue > "); fflush(stdout);
00048     char go;
00049     scanf("%c",&go);
00050   }
00051 
00052   Comm.Barrier();
00053 
00054 }
00055 
00056 //============================================================================
00057 Epetra_CrsMatrix* Ifpack_CreateOverlappingCrsMatrix(const Epetra_RowMatrix* Matrix,
00058                                                     const int OverlappingLevel)
00059 {
00060 
00061   if (OverlappingLevel == 0) 
00062     return(0); // All done
00063   if (Matrix->Comm().NumProc() == 1) 
00064     return(0); // All done
00065 
00066   Epetra_CrsMatrix* OverlappingMatrix;
00067   OverlappingMatrix = 0;
00068   Epetra_Map* OverlappingMap;
00069   OverlappingMap = (Epetra_Map*)&(Matrix->RowMatrixRowMap());
00070 
00071   const Epetra_RowMatrix* OldMatrix;
00072   const Epetra_Map* DomainMap = &(Matrix->OperatorDomainMap());
00073   const Epetra_Map* RangeMap = &(Matrix->OperatorRangeMap());
00074 
00075   for (int level = 1; level <= OverlappingLevel ; ++level) {
00076 
00077     if (OverlappingMatrix)
00078       OldMatrix = OverlappingMatrix;
00079     else
00080       OldMatrix = Matrix;
00081 
00082     Epetra_Import* OverlappingImporter;
00083     OverlappingImporter = (Epetra_Import*)OldMatrix->RowMatrixImporter();
00084     int NumMyElements = OverlappingImporter->TargetMap().NumMyElements();
00085     int* MyGlobalElements = OverlappingImporter->TargetMap().MyGlobalElements();
00086 
00087     // need to build an Epetra_Map in this way because Epetra_CrsMatrix
00088     // requires Epetra_Map and not Epetra_BlockMap
00089 
00090     OverlappingMap = new Epetra_Map(-1,NumMyElements,MyGlobalElements,
00091                                     0, Matrix->Comm());
00092 
00093     if (level < OverlappingLevel)
00094       OverlappingMatrix = new Epetra_CrsMatrix(Copy, *OverlappingMap, 0);
00095     else
00096       // On last iteration, we want to filter out all columns except 
00097       // those that correspond
00098       // to rows in the graph.  This assures that our matrix is square
00099       OverlappingMatrix = new Epetra_CrsMatrix(Copy, *OverlappingMap, 
00100                                                *OverlappingMap, 0);
00101 
00102     OverlappingMatrix->Import(*OldMatrix, *OverlappingImporter, Insert);
00103     if (level < OverlappingLevel) {
00104       OverlappingMatrix->FillComplete(*DomainMap, *RangeMap);
00105     }
00106     else {
00107       OverlappingMatrix->FillComplete(*DomainMap, *RangeMap);
00108     }
00109 
00110     delete OverlappingMap;
00111 
00112     if (level > 1) {
00113       delete OldMatrix;
00114     }
00115     OverlappingMatrix->FillComplete();
00116 
00117   }
00118 
00119   return(OverlappingMatrix);
00120 }
00121 
00122 //============================================================================
00123 Epetra_CrsGraph* Ifpack_CreateOverlappingCrsMatrix(const Epetra_CrsGraph* Graph,
00124                                                    const int OverlappingLevel)
00125 {
00126 
00127   if (OverlappingLevel == 0) 
00128     return(0); // All done
00129   if (Graph->Comm().NumProc() == 1) 
00130     return(0); // All done
00131 
00132   Epetra_CrsGraph* OverlappingGraph;
00133   Epetra_BlockMap* OverlappingMap;
00134   OverlappingGraph = const_cast<Epetra_CrsGraph*>(Graph);
00135   OverlappingMap = const_cast<Epetra_BlockMap*>(&(Graph->RowMap()));
00136 
00137   Epetra_CrsGraph* OldGraph;
00138   Epetra_BlockMap* OldMap;
00139   const Epetra_BlockMap* DomainMap = &(Graph->DomainMap());
00140   const Epetra_BlockMap* RangeMap = &(Graph->RangeMap());
00141 
00142   for (int level = 1; level <= OverlappingLevel ; ++level) {
00143 
00144     OldGraph = OverlappingGraph;
00145     OldMap = OverlappingMap;
00146 
00147     Epetra_Import* OverlappingImporter;
00148     OverlappingImporter = const_cast<Epetra_Import*>(OldGraph->Importer());
00149     OverlappingMap = new Epetra_BlockMap(OverlappingImporter->TargetMap());
00150 
00151     if (level < OverlappingLevel)
00152       OverlappingGraph = new Epetra_CrsGraph(Copy, *OverlappingMap, 0);
00153     else
00154       // On last iteration, we want to filter out all columns except 
00155       // those that correspond
00156       // to rows in the graph.  This assures that our matrix is square
00157       OverlappingGraph = new Epetra_CrsGraph(Copy, *OverlappingMap, 
00158                                           *OverlappingMap, 0);
00159 
00160     OverlappingGraph->Import(*OldGraph, *OverlappingImporter, Insert);
00161     if (level < OverlappingLevel) 
00162       OverlappingGraph->FillComplete(*DomainMap, *RangeMap);
00163     else {
00164       // Copy last OverlapImporter because we will use it later
00165       OverlappingImporter = new Epetra_Import(*OverlappingMap, *DomainMap);
00166       OverlappingGraph->FillComplete(*DomainMap, *RangeMap);
00167     }
00168 
00169     if (level > 1) {
00170       delete OldGraph;
00171       delete OldMap;
00172     }
00173 
00174     delete OverlappingMap;
00175     OverlappingGraph->FillComplete();
00176 
00177   }
00178 
00179   return(OverlappingGraph);
00180 }
00181 
00182 //============================================================================
00183 string Ifpack_toString(const int& x)
00184 {
00185   char s[100];
00186   sprintf(s, "%d", x);
00187   return string(s);
00188 }
00189 
00190 //============================================================================
00191 string Ifpack_toString(const double& x)
00192 {
00193   char s[100];
00194   sprintf(s, "%g", x);
00195   return string(s);
00196 }
00197 
00198 //============================================================================
00199 int Ifpack_PrintResidual(char* Label, const Epetra_RowMatrix& A,
00200                          const Epetra_MultiVector& X, const Epetra_MultiVector&Y)
00201 {
00202   if (X.Comm().MyPID() == 0) {
00203     cout << "***** " << Label << endl;
00204   }
00205   Ifpack_PrintResidual(0,A,X,Y);
00206 
00207   return(0);
00208 }
00209 
00210 //============================================================================
00211 int Ifpack_PrintResidual(const int iter, const Epetra_RowMatrix& A,
00212                          const Epetra_MultiVector& X, const Epetra_MultiVector&Y)
00213 {
00214   Epetra_MultiVector RHS(X);
00215   std::vector<double> Norm2;
00216   Norm2.resize(X.NumVectors());
00217 
00218   IFPACK_CHK_ERR(A.Multiply(false,X,RHS));
00219   RHS.Update(1.0, Y, -1.0);
00220 
00221   RHS.Norm2(&Norm2[0]);
00222 
00223   if (X.Comm().MyPID() == 0) {
00224     cout << "***** iter: " << iter << ":  ||Ax - b||_2 = " 
00225          << Norm2[0] << endl;
00226   }
00227 
00228   return(0);
00229 }
00230 
00231 //============================================================================
00232 // very simple debugging function that prints on screen the structure
00233 // of the input matrix.
00234 void Ifpack_PrintSparsity_Simple(const Epetra_RowMatrix& A)
00235 {
00236   int MaxEntries = A.MaxNumEntries();
00237   vector<int> Indices(MaxEntries);
00238   vector<double> Values(MaxEntries);
00239   vector<bool> FullRow(A.NumMyRows());
00240 
00241   cout << "+-";
00242   for (int j = 0 ; j < A.NumMyRows() ; ++j)
00243     cout << '-';
00244   cout << "-+" << endl;
00245 
00246   for (int i = 0 ; i < A.NumMyRows() ; ++i) {
00247 
00248     int Length;
00249     A.ExtractMyRowCopy(i,MaxEntries,Length,
00250                        &Values[0], &Indices[0]);
00251 
00252     for (int j = 0 ; j < A.NumMyRows() ; ++j)
00253       FullRow[j] = false;
00254 
00255     for (int j = 0 ; j < Length ; ++j) {
00256       FullRow[Indices[j]] = true;
00257     }
00258 
00259     cout << "| ";
00260     for (int j = 0 ; j < A.NumMyRows() ; ++j) {
00261       if (FullRow[j])
00262         cout << '*';
00263       else
00264         cout << ' ';
00265     }
00266     cout << " |" << endl;
00267   }
00268 
00269   cout << "+-";
00270   for (int j = 0 ; j < A.NumMyRows() ; ++j)
00271     cout << '-';
00272   cout << "-+" << endl << endl;
00273 
00274 }
00275 
00276 //============================================================================
00277 
00278 double Ifpack_FrobeniusNorm(const Epetra_RowMatrix& A)
00279 {
00280   double MyNorm = 0.0, GlobalNorm;
00281 
00282   vector<int> colInd(A.MaxNumEntries());
00283   vector<double> colVal(A.MaxNumEntries());
00284 
00285   for (int i = 0 ; i < A.NumMyRows() ; ++i) {
00286 
00287     int Nnz;
00288     IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
00289                                       &colVal[0],&colInd[0]));
00290 
00291     for (int j = 0 ; j < Nnz ; ++j) {
00292       MyNorm += colVal[j] * colVal[j];
00293     }
00294   }
00295 
00296   A.Comm().SumAll(&MyNorm,&GlobalNorm,1);
00297 
00298   return(sqrt(GlobalNorm));
00299 }
00300 
00301 static void print()
00302 {
00303   printf("\n");
00304 }
00305 
00306 #include <iomanip>
00307 template<class T>
00308 static void print(char* str, T val)
00309 {
00310   cout.width(30); cout.setf(ios::left);
00311   cout << str;
00312   cout << " = " << val << endl;
00313 }
00314 
00315 template<class T>
00316 static void print(char* str, T val, double percentage)
00317 {
00318   cout.width(30); cout.setf(ios::left);
00319   cout << str;
00320   cout << " = ";
00321   cout.width(20); cout.setf(ios::left);
00322   cout << val;
00323   cout << " ( " << percentage << " %)" << endl;
00324 }
00325 template<class T>
00326 static void print(char* str, T one, T two, T three, bool equal = true)
00327 {
00328   cout.width(30); cout.setf(ios::left);
00329   cout << str;
00330   if (equal) 
00331     cout << " = ";
00332   else
00333     cout << "   ";
00334   cout.width(15); cout.setf(ios::left);
00335   cout << one;
00336   cout.width(15); cout.setf(ios::left);
00337   cout << two;
00338   cout.width(15); cout.setf(ios::left);
00339   cout << three;
00340   cout << endl;
00341 }
00342 
00343 //============================================================================
00344 #include "limits.h"
00345 #include "float.h"
00346 #include "Epetra_FECrsMatrix.h"
00347 
00348 int Ifpack_Analyze(const Epetra_RowMatrix& A, const bool Cheap,
00349                    const int NumPDEEqns)
00350 {
00351 
00352   int NumMyRows = A.NumMyRows();
00353   int NumGlobalRows = A.NumGlobalRows();
00354   int NumGlobalCols = A.NumGlobalCols();
00355   int MyBandwidth = 0, GlobalBandwidth;
00356   int MyLowerNonzeros = 0, MyUpperNonzeros = 0;
00357   int GlobalLowerNonzeros, GlobalUpperNonzeros;
00358   int MyDiagonallyDominant = 0, GlobalDiagonallyDominant;
00359   int MyWeaklyDiagonallyDominant = 0, GlobalWeaklyDiagonallyDominant;
00360   double MyMin, MyAvg, MyMax;
00361   double GlobalMin, GlobalAvg, GlobalMax;
00362   int GlobalStorage;
00363 
00364   bool verbose = (A.Comm().MyPID() == 0);
00365 
00366   GlobalStorage = sizeof(int*) * NumGlobalRows + 
00367     sizeof(int) * A.NumGlobalNonzeros() + 
00368     sizeof(double) * A.NumGlobalNonzeros();
00369 
00370   if (verbose) {
00371     print();
00372     Ifpack_PrintLine();
00373     print<const char*>("Label", A.Label());
00374     print<int>("Global rows", NumGlobalRows);
00375     print<int>("Global columns", NumGlobalCols);
00376     print<int>("Stored nonzeros", A.NumGlobalNonzeros());
00377     print<int>("Nonzeros / row", A.NumGlobalNonzeros() / NumGlobalRows);
00378     print<double>("Estimated storage (Mbytes)", 1.0e-6 * GlobalStorage);
00379   }
00380 
00381   int NumMyActualNonzeros = 0, NumGlobalActualNonzeros;
00382   int NumMyEmptyRows = 0, NumGlobalEmptyRows;
00383   int NumMyDirichletRows = 0, NumGlobalDirichletRows;
00384 
00385   vector<int> colInd(A.MaxNumEntries());
00386   vector<double> colVal(A.MaxNumEntries());
00387 
00388   Epetra_Vector Diag(A.RowMatrixRowMap());
00389   Epetra_Vector RowSum(A.RowMatrixRowMap());
00390   Diag.PutScalar(0.0);
00391   RowSum.PutScalar(0.0);
00392 
00393   for (int i = 0 ; i < NumMyRows ; ++i) {
00394 
00395     int GRID = A.RowMatrixRowMap().GID(i);
00396     int Nnz;
00397     IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
00398                                       &colVal[0],&colInd[0]));
00399 
00400     if (Nnz == 0)
00401       NumMyEmptyRows++;
00402 
00403     if (Nnz == 1)
00404       NumMyDirichletRows++;
00405 
00406     for (int j = 0 ; j < Nnz ; ++j) {
00407 
00408       double v = colVal[j];
00409       if (v < 0) v = -v;
00410       if (colVal[j] != 0.0)
00411         NumMyActualNonzeros++;
00412 
00413       int GCID = A.RowMatrixColMap().GID(colInd[j]);
00414 
00415       if (GCID != GRID)
00416         RowSum[i] += v;
00417       else
00418         Diag[i] = v;
00419 
00420       if (GCID < GRID) 
00421         MyLowerNonzeros++;
00422       else if (GCID > GRID) 
00423         MyUpperNonzeros++;
00424       int b = GCID - GRID;
00425       if (b < 0) b = -b;
00426       if (b > MyBandwidth)
00427         MyBandwidth = b;
00428     }
00429 
00430     if (Diag[i] > RowSum[i])
00431       MyDiagonallyDominant++;
00432 
00433     if (Diag[i] >= RowSum[i])
00434       MyWeaklyDiagonallyDominant++;
00435 
00436     RowSum[i] += Diag[i];
00437   }
00438 
00439   // ======================== //
00440   // summing up global values //
00441   // ======================== //
00442  
00443   A.Comm().SumAll(&MyDiagonallyDominant,&GlobalDiagonallyDominant,1);
00444   A.Comm().SumAll(&MyWeaklyDiagonallyDominant,&GlobalWeaklyDiagonallyDominant,1);
00445   A.Comm().SumAll(&NumMyActualNonzeros, &NumGlobalActualNonzeros, 1);
00446   A.Comm().SumAll(&NumMyEmptyRows, &NumGlobalEmptyRows, 1);
00447   A.Comm().SumAll(&NumMyDirichletRows, &NumGlobalDirichletRows, 1);
00448   A.Comm().SumAll(&MyBandwidth, &GlobalBandwidth, 1);
00449   A.Comm().SumAll(&MyLowerNonzeros, &GlobalLowerNonzeros, 1);
00450   A.Comm().SumAll(&MyUpperNonzeros, &GlobalUpperNonzeros, 1);
00451   A.Comm().SumAll(&MyDiagonallyDominant, &GlobalDiagonallyDominant, 1);
00452   A.Comm().SumAll(&MyWeaklyDiagonallyDominant, &GlobalWeaklyDiagonallyDominant, 1);
00453  
00454   double NormOne = A.NormOne();
00455   double NormInf = A.NormInf();
00456   double NormF   = Ifpack_FrobeniusNorm(A);
00457 
00458   if (verbose) {
00459     print();
00460     print<int>("Actual nonzeros", NumGlobalActualNonzeros);
00461     print<int>("Nonzeros in strict lower part", GlobalLowerNonzeros);
00462     print<int>("Nonzeros in strict upper part", GlobalUpperNonzeros);
00463     print();
00464     print<int>("Empty rows", NumGlobalEmptyRows,
00465                100.0 * NumGlobalEmptyRows / NumGlobalRows);
00466     print<int>("Dirichlet rows", NumGlobalDirichletRows,
00467                100.0 * NumGlobalDirichletRows / NumGlobalRows);
00468     print<int>("Diagonally dominant rows", GlobalDiagonallyDominant,
00469                100.0 * GlobalDiagonallyDominant / NumGlobalRows);
00470     print<int>("Weakly diag. dominant rows", 
00471                GlobalWeaklyDiagonallyDominant,
00472                100.0 * GlobalWeaklyDiagonallyDominant / NumGlobalRows);
00473     print();
00474     print<int>("Maximum bandwidth", GlobalBandwidth);
00475 
00476     print();
00477     print("", "one-norm", "inf-norm", "Frobenius", false);
00478     print("", "========", "========", "=========", false);
00479     print();
00480 
00481     print<double>("A", NormOne, NormInf, NormF);
00482   }
00483 
00484   if (Cheap == false) {
00485 
00486     // create A + A^T and A - A^T
00487 
00488     Epetra_FECrsMatrix AplusAT(Copy, A.RowMatrixRowMap(), 0);
00489     Epetra_FECrsMatrix AminusAT(Copy, A.RowMatrixRowMap(), 0);
00490 
00491     for (int i = 0 ; i < NumMyRows ; ++i) {
00492 
00493       int GRID = A.RowMatrixRowMap().GID(i);
00494       assert (GRID != -1);
00495 
00496       int Nnz;
00497       IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
00498                                         &colVal[0],&colInd[0]));
00499 
00500       for (int j = 0 ; j < Nnz ; ++j) {
00501 
00502         int GCID         = A.RowMatrixColMap().GID(colInd[j]);
00503         assert (GCID != -1);
00504 
00505         double plus_val  = colVal[j];
00506         double minus_val = -colVal[j];
00507 
00508         if (AplusAT.SumIntoGlobalValues(1,&GRID,1,&GCID,&plus_val) != 0) {
00509           IFPACK_CHK_ERR(AplusAT.InsertGlobalValues(1,&GRID,1,&GCID,&plus_val));
00510         }
00511 
00512         if (AplusAT.SumIntoGlobalValues(1,&GCID,1,&GRID,&plus_val) != 0) {
00513           IFPACK_CHK_ERR(AplusAT.InsertGlobalValues(1,&GCID,1,&GRID,&plus_val));
00514         }
00515 
00516         if (AminusAT.SumIntoGlobalValues(1,&GRID,1,&GCID,&plus_val) != 0) {
00517           IFPACK_CHK_ERR(AminusAT.InsertGlobalValues(1,&GRID,1,&GCID,&plus_val));
00518         }
00519 
00520         if (AminusAT.SumIntoGlobalValues(1,&GCID,1,&GRID,&minus_val) != 0) {
00521           IFPACK_CHK_ERR(AminusAT.InsertGlobalValues(1,&GCID,1,&GRID,&minus_val));
00522         }
00523 
00524       }
00525     }
00526 
00527     AplusAT.FillComplete();
00528     AminusAT.FillComplete();
00529 
00530     AplusAT.Scale(0.5);
00531     AminusAT.Scale(0.5);
00532 
00533     NormOne = AplusAT.NormOne();
00534     NormInf = AplusAT.NormInf();
00535     NormF   = Ifpack_FrobeniusNorm(AplusAT);
00536 
00537     if (verbose) {
00538       print<double>("A + A^T", NormOne, NormInf, NormF);
00539     }
00540 
00541     NormOne = AminusAT.NormOne();
00542     NormInf = AminusAT.NormInf();
00543     NormF   = Ifpack_FrobeniusNorm(AminusAT);
00544 
00545     if (verbose) {
00546       print<double>("A - A^T", NormOne, NormInf, NormF);
00547     }
00548   }
00549 
00550   if (verbose) {
00551     print();
00552     print<char*>("", "min", "avg", "max", false);
00553     print<char*>("", "===", "===", "===", false);
00554   }
00555 
00556   MyMax = -DBL_MAX;
00557   MyMin = DBL_MAX;
00558   MyAvg = 0.0;
00559 
00560   for (int i = 0 ; i < NumMyRows ; ++i) {
00561 
00562     int Nnz;
00563     IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
00564                                       &colVal[0],&colInd[0]));
00565 
00566     for (int j = 0 ; j < Nnz ; ++j) {
00567       MyAvg += colVal[j];
00568       if (colVal[j] > MyMax) MyMax = colVal[j];
00569       if (colVal[j] < MyMin) MyMin = colVal[j];
00570     }
00571   }
00572 
00573   A.Comm().MaxAll(&MyMax, &GlobalMax, 1);
00574   A.Comm().MinAll(&MyMin, &GlobalMin, 1);
00575   A.Comm().SumAll(&MyAvg, &GlobalAvg, 1);
00576   GlobalAvg /= A.NumGlobalNonzeros();
00577 
00578   if (verbose) {
00579     print();
00580     print<double>(" A(i,j)", GlobalMin, GlobalAvg, GlobalMax);
00581   }
00582 
00583   MyMax = 0.0;
00584   MyMin = DBL_MAX;
00585   MyAvg = 0.0;
00586 
00587   for (int i = 0 ; i < NumMyRows ; ++i) {
00588 
00589     int Nnz;
00590     IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
00591                                       &colVal[0],&colInd[0]));
00592 
00593     for (int j = 0 ; j < Nnz ; ++j) {
00594       double v = colVal[j];
00595       if (v < 0) v = -v;
00596       MyAvg += v;
00597       if (colVal[j] > MyMax) MyMax = v;
00598       if (colVal[j] < MyMin) MyMin = v;
00599     }
00600   }
00601 
00602   A.Comm().MaxAll(&MyMax, &GlobalMax, 1);
00603   A.Comm().MinAll(&MyMin, &GlobalMin, 1);
00604   A.Comm().SumAll(&MyAvg, &GlobalAvg, 1);
00605   GlobalAvg /= A.NumGlobalNonzeros();
00606 
00607   if (verbose) {
00608     print<double>("|A(i,j)|", GlobalMin, GlobalAvg, GlobalMax);
00609   }
00610 
00611   // ================= //
00612   // diagonal elements //
00613   // ================= //
00614 
00615   Diag.MinValue(&GlobalMin);
00616   Diag.MaxValue(&GlobalMax);
00617   Diag.MeanValue(&GlobalAvg);
00618 
00619   if (verbose) {
00620     print();
00621     print<double>(" A(k,k)", GlobalMin, GlobalAvg, GlobalMax);
00622   }
00623 
00624   Diag.Abs(Diag);
00625   Diag.MinValue(&GlobalMin);
00626   Diag.MaxValue(&GlobalMax);
00627   Diag.MeanValue(&GlobalAvg);
00628   if (verbose) {
00629     print<double>("|A(k,k)|", GlobalMin, GlobalAvg, GlobalMax);
00630   }
00631   
00632   // ============================================== //
00633   // cycle over all equations for diagonal elements //
00634   // ============================================== //
00635 
00636   if (NumPDEEqns > 1 ) {
00637 
00638     if (verbose) print();
00639 
00640     for (int ie = 0 ; ie < NumPDEEqns ; ie++) {
00641 
00642       MyMin = DBL_MAX;
00643       MyMax = -DBL_MAX;
00644       MyAvg = 0.0;
00645 
00646       for (int i = ie ; i < Diag.MyLength() ; i += NumPDEEqns) {
00647         double d = Diag[i];
00648         MyAvg += d;
00649         if (d < MyMin)
00650           MyMin = d;
00651         if (d > MyMax)
00652           MyMax = d;
00653       }
00654       A.Comm().MinAll(&MyMin, &GlobalMin, 1);
00655       A.Comm().MaxAll(&MyMax, &GlobalMax, 1);
00656       A.Comm().SumAll(&MyAvg, &GlobalAvg, 1);
00657       // does not really work fine if the number of global
00658       // elements is not a multiple of NumPDEEqns
00659       GlobalAvg /= (Diag.GlobalLength() / NumPDEEqns);
00660 
00661       if (verbose) {
00662         char str[80];
00663         sprintf(str, " A(k,k), eq %d", ie);
00664         print<double>(str, GlobalMin, GlobalAvg, GlobalMax);
00665       }
00666     }
00667   }
00668 
00669   // ======== //
00670   // row sums //
00671   // ======== //
00672   
00673   RowSum.MinValue(&GlobalMin);
00674   RowSum.MaxValue(&GlobalMax);
00675   RowSum.MeanValue(&GlobalAvg);
00676 
00677   if (verbose) {
00678     print();
00679     print<double>(" sum_j A(k,j)", GlobalMin, GlobalAvg, GlobalMax);
00680   }
00681 
00682   // ===================================== //
00683   // cycle over all equations for row sums //
00684   // ===================================== //
00685 
00686   if (NumPDEEqns > 1 ) {
00687 
00688     if (verbose) print();
00689 
00690     for (int ie = 0 ; ie < NumPDEEqns ; ie++) {
00691 
00692       MyMin = DBL_MAX;
00693       MyMax = -DBL_MAX;
00694       MyAvg = 0.0;
00695 
00696       for (int i = ie ; i < Diag.MyLength() ; i += NumPDEEqns) {
00697         double d = RowSum[i];
00698         MyAvg += d;
00699         if (d < MyMin)
00700           MyMin = d;
00701         if (d > MyMax)
00702           MyMax = d;
00703       }
00704       A.Comm().MinAll(&MyMin, &GlobalMin, 1);
00705       A.Comm().MaxAll(&MyMax, &GlobalMax, 1);
00706       A.Comm().SumAll(&MyAvg, &GlobalAvg, 1);
00707       // does not really work fine if the number of global
00708       // elements is not a multiple of NumPDEEqns
00709       GlobalAvg /= (Diag.GlobalLength() / NumPDEEqns);
00710 
00711       if (verbose) {
00712         char str[80];
00713         sprintf(str, " sum_j A(k,j), eq %d", ie);
00714         print<double>(str, GlobalMin, GlobalAvg, GlobalMax);
00715       }
00716     }
00717   }
00718 
00719   if (verbose)
00720     Ifpack_PrintLine();
00721 
00722   return(0);
00723 }
00724 
00725 int Ifpack_AnalyzeVectorElements(const Epetra_Vector& Diagonal,
00726                                  const bool abs, const int steps)
00727 {
00728 
00729   bool verbose = (Diagonal.Comm().MyPID() == 0);
00730   double min_val =  DBL_MAX;
00731   double max_val = -DBL_MAX;
00732 
00733   for (int i = 0 ; i < Diagonal.MyLength() ; ++i) {
00734     double v = Diagonal[i];
00735     if (abs)
00736       if (v < 0) v = -v;
00737     if (v > max_val)
00738       max_val = v;
00739     if (v < min_val)
00740       min_val = v;
00741   }
00742 
00743   if (verbose) {
00744     cout << endl;
00745     Ifpack_PrintLine();
00746     cout << "Vector label = " << Diagonal.Label() << endl;
00747     cout << endl;
00748   }
00749 
00750   double delta = (max_val - min_val) / steps;
00751   for (int k = 0 ; k < steps ; ++k) {
00752 
00753     double below = delta * k + min_val;
00754     double above = below + delta;
00755     int MyBelow = 0, GlobalBelow;
00756 
00757     for (int i = 0 ; i < Diagonal.MyLength() ; ++i) {
00758       double v = Diagonal[i];
00759       if (v < 0) v = -v;
00760       if (v >= below && v < above) MyBelow++;
00761     }
00762 
00763     Diagonal.Comm().SumAll(&MyBelow, &GlobalBelow, 1);
00764 
00765     if (verbose) {
00766       printf("Elements in [%+7e, %+7e) = %10d ( = %5.2f %%)\n",
00767              below, above, GlobalBelow,
00768              100.0 * GlobalBelow / Diagonal.GlobalLength());
00769     }
00770   }
00771   
00772   if (verbose) {
00773     Ifpack_PrintLine();
00774     cout << endl;
00775   }
00776 
00777   return(0);
00778 }
00779 
00780 // ====================================================================== 
00781 
00782 int Ifpack_AnalyzeMatrixElements(const Epetra_RowMatrix& A,
00783                                  const bool abs, const int steps)
00784 {
00785 
00786   bool verbose = (A.Comm().MyPID() == 0);
00787   double min_val =  DBL_MAX;
00788   double max_val = -DBL_MAX;
00789 
00790   vector<int>    colInd(A.MaxNumEntries());
00791   vector<double> colVal(A.MaxNumEntries());
00792   
00793   for (int i = 0 ; i < A.NumMyRows() ; ++i) {
00794 
00795     int Nnz;
00796     IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
00797                                       &colVal[0],&colInd[0]));
00798 
00799     for (int j = 0 ; j < Nnz ; ++j) {
00800       double v = colVal[j];
00801       if (abs)
00802         if (v < 0) v = -v;
00803       if (v < min_val)
00804         min_val = v;
00805       if (v > max_val)
00806         max_val = v;
00807     }
00808   }
00809 
00810   if (verbose) {
00811     cout << endl;
00812     Ifpack_PrintLine();
00813     cout << "Label of matrix = " << A.Label() << endl;
00814     cout << endl;
00815   }
00816 
00817   double delta = (max_val - min_val) / steps;
00818   for (int k = 0 ; k < steps ; ++k) {
00819 
00820     double below = delta * k + min_val;
00821     double above = below + delta;
00822     int MyBelow = 0, GlobalBelow;
00823 
00824     for (int i = 0 ; i < A.NumMyRows() ; ++i) {
00825 
00826       int Nnz;
00827       IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
00828                                         &colVal[0],&colInd[0]));
00829 
00830       for (int j = 0 ; j < Nnz ; ++j) {
00831         double v = colVal[j];
00832         if (abs)
00833           if (v < 0) v = -v;
00834         if (v >= below && v < above) MyBelow++;
00835       }
00836 
00837     }
00838     A.Comm().SumAll(&MyBelow, &GlobalBelow, 1);
00839     if (verbose) {
00840       printf("Elements in [%+7e, %+7e) = %10d ( = %5.2f %%)\n",
00841              below, above, GlobalBelow,
00842              100.0 * GlobalBelow / A.NumGlobalNonzeros());
00843     }
00844   }
00845 
00846   if (verbose) {
00847     Ifpack_PrintLine();
00848     cout << endl;
00849   }
00850 
00851   return(0);
00852 }
00853 
00854 // ====================================================================== 
00855 int Ifpack_PrintSparsity(const Epetra_RowMatrix& A, const char* InputFileName, 
00856                          const int NumPDEEqns)
00857 {
00858 
00859   int m,nc,nr,maxdim,ltit;
00860   double lrmrgn,botmrgn,xtit,ytit,ytitof,fnstit,siz = 0.0;
00861   double xl,xr, yb,yt, scfct,u2dot,frlw,delt,paperx;
00862   bool square = false;
00863   /*change square to .true. if you prefer a square frame around
00864     a rectangular matrix */
00865   double conv = 2.54;
00866   char munt = 'E'; /* put 'E' for centimeters, 'U' for inches */
00867   int ptitle = 0; /* position of the title, 0 under the drawing,
00868                      else above */
00869   FILE* fp = NULL;
00870   int NumMyRows;
00871   //int NumMyCols;
00872   int NumGlobalRows;
00873   int NumGlobalCols;
00874   int MyPID;
00875   int NumProc;
00876   char FileName[1024];
00877   char title[1024];
00878   
00879   const Epetra_Comm& Comm = A.Comm();
00880 
00881   /* --------------------- execution begins ---------------------- */
00882 
00883   if (strlen(A.Label()) != 0)
00884     strcpy(title, A.Label());
00885   else
00886     sprintf(title, "%s", "matrix");
00887 
00888   if (InputFileName == 0)
00889     sprintf(FileName, "%s.ps", title);
00890   else
00891     strcpy(FileName, InputFileName);
00892 
00893   MyPID = Comm.MyPID();
00894   NumProc = Comm.NumProc();
00895 
00896   NumMyRows = A.NumMyRows();
00897   //NumMyCols = A.NumMyCols();
00898 
00899   NumGlobalRows = A.NumGlobalRows();
00900   NumGlobalCols = A.NumGlobalCols();
00901 
00902   if (NumGlobalRows != NumGlobalCols)
00903     IFPACK_CHK_ERR(-1); // never tested
00904 
00905   /* to be changed for rect matrices */
00906   maxdim = (NumGlobalRows>NumGlobalCols)?NumGlobalRows:NumGlobalCols;
00907   maxdim /= NumPDEEqns;
00908 
00909   m = 1 + maxdim;
00910   nr = NumGlobalRows / NumPDEEqns + 1;
00911   nc = NumGlobalCols / NumPDEEqns + 1;
00912 
00913   if (munt == 'E') {
00914     u2dot = 72.0/conv;
00915     paperx = 21.0;
00916     siz = 10.0;
00917   }
00918   else {
00919     u2dot = 72.0;
00920     paperx = 8.5*conv;
00921     siz = siz*conv;
00922   }
00923 
00924   /* left and right margins (drawing is centered) */
00925 
00926   lrmrgn = (paperx-siz)/2.0;
00927 
00928   /* bottom margin : 2 cm */
00929 
00930   botmrgn = 2.0;
00931   /* c scaling factor */
00932   scfct = siz*u2dot/m;
00933   /* matrix frame line witdh */
00934   frlw = 0.25;
00935   /* font size for title (cm) */
00936   fnstit = 0.5;
00937   if (title) ltit = strlen(title);
00938   else       ltit = 0;
00939   /* position of title : centered horizontally */
00940   /*                     at 1.0 cm vertically over the drawing */
00941   ytitof = 1.0;
00942   xtit = paperx/2.0;
00943   ytit = botmrgn+siz*nr/m + ytitof;
00944   /* almost exact bounding box */
00945   xl = lrmrgn*u2dot - scfct*frlw/2;
00946   xr = (lrmrgn+siz)*u2dot + scfct*frlw/2;
00947   yb = botmrgn*u2dot - scfct*frlw/2;
00948   yt = (botmrgn+siz*nr/m)*u2dot + scfct*frlw/2;
00949   if (ltit == 0) {
00950     yt = yt + (ytitof+fnstit*0.70)*u2dot;
00951   } 
00952   /* add some room to bounding box */
00953   delt = 10.0;
00954   xl = xl-delt;
00955   xr = xr+delt;
00956   yb = yb-delt;
00957   yt = yt+delt;
00958 
00959   /* correction for title under the drawing */
00960   if ((ptitle == 0) && (ltit == 0)) {
00961     ytit = botmrgn + fnstit*0.3;
00962     botmrgn = botmrgn + ytitof + fnstit*0.7;
00963   }
00964 
00965   /* begin of output */
00966 
00967   if (MyPID == 0) {
00968 
00969     fp = fopen(FileName,"w");
00970 
00971     fprintf(fp,"%%%%!PS-Adobe-2.0\n");
00972     fprintf(fp,"%%%%Creator: IFPACK\n");
00973     fprintf(fp,"%%%%BoundingBox: %f %f %f %f\n",
00974             xl,yb,xr,yt);
00975     fprintf(fp,"%%%%EndComments\n");
00976     fprintf(fp,"/cm {72 mul 2.54 div} def\n");
00977     fprintf(fp,"/mc {72 div 2.54 mul} def\n");
00978     fprintf(fp,"/pnum { 72 div 2.54 mul 20 string ");
00979     fprintf(fp,"cvs print ( ) print} def\n");
00980     fprintf(fp,"/Cshow {dup stringwidth pop -2 div 0 rmoveto show} def\n");
00981 
00982     /* we leave margins etc. in cm so it is easy to modify them if
00983        needed by editing the output file */
00984     fprintf(fp,"gsave\n");
00985     if (ltit != 0) {
00986       fprintf(fp,"/Helvetica findfont %e cm scalefont setfont\n",
00987               fnstit);
00988       fprintf(fp,"%f cm %f cm moveto\n",
00989               xtit,ytit);
00990       fprintf(fp,"(%s) Cshow\n", title);
00991       fprintf(fp,"%f cm %f cm translate\n",
00992               lrmrgn,botmrgn);
00993     }
00994     fprintf(fp,"%f cm %d div dup scale \n",
00995             siz,m);
00996     /* draw a frame around the matrix */
00997 
00998     fprintf(fp,"%f setlinewidth\n",
00999             frlw);
01000     fprintf(fp,"newpath\n");
01001     fprintf(fp,"0 0 moveto ");
01002     if (square) {
01003       printf("------------------- %d\n",m);
01004       fprintf(fp,"%d %d lineto\n",
01005               m,0);
01006       fprintf(fp,"%d %d lineto\n",
01007               m, m);
01008       fprintf(fp,"%d %d lineto\n",
01009               0, m);
01010     } 
01011     else {
01012       fprintf(fp,"%d %d lineto\n",
01013               nc, 0);
01014       fprintf(fp,"%d %d lineto\n",
01015               nc, nr);
01016       fprintf(fp,"%d %d lineto\n",
01017               0, nr);
01018     }
01019     fprintf(fp,"closepath stroke\n");
01020 
01021     /* plotting loop */
01022 
01023     fprintf(fp,"1 1 translate\n");
01024     fprintf(fp,"0.8 setlinewidth\n");
01025     fprintf(fp,"/p {moveto 0 -.40 rmoveto \n");
01026     fprintf(fp,"           0  .80 rlineto stroke} def\n");
01027 
01028     fclose(fp);
01029   }
01030   
01031   int MaxEntries = A.MaxNumEntries();
01032   vector<int> Indices(MaxEntries);
01033   vector<double> Values(MaxEntries);
01034 
01035   for (int pid = 0 ; pid < NumProc ; ++pid) {
01036 
01037     if (pid == MyPID) {
01038 
01039       fp = fopen(FileName,"a");
01040       if( fp == NULL ) {
01041         fprintf(stderr,"ERROR\n");
01042         exit(EXIT_FAILURE);
01043       }
01044 
01045       for (int i = 0 ; i < NumMyRows ; ++i) {
01046 
01047         if (i % NumPDEEqns) continue;
01048 
01049         int Nnz;
01050         A.ExtractMyRowCopy(i,MaxEntries,Nnz,&Values[0],&Indices[0]);
01051 
01052         int grow = A.RowMatrixRowMap().GID(i);
01053 
01054         for (int j = 0 ; j < Nnz ; ++j) {
01055           int col = Indices[j];
01056           if (col % NumPDEEqns == 0) {
01057             int gcol = A.RowMatrixColMap().GID(Indices[j]);
01058             grow /= NumPDEEqns;
01059             gcol /= NumPDEEqns;
01060             fprintf(fp,"%d %d p\n",
01061                     gcol, NumGlobalRows - grow - 1); 
01062           }
01063         }
01064       }
01065 
01066       fprintf(fp,"%%end of data for this process\n");
01067 
01068       if( pid == NumProc - 1 )
01069         fprintf(fp,"showpage\n");
01070 
01071       fclose(fp);
01072     }
01073     Comm.Barrier();
01074   }
01075 
01076   return(0);
01077 }

Generated on Thu Sep 18 12:37:22 2008 for Ifpack Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1