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