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