Epetra Package Browser (Single Doxygen Collection) Development
Epetra_BasicRowMatrix.cpp
Go to the documentation of this file.
00001 
00002 //@HEADER
00003 // ************************************************************************
00004 // 
00005 //               Epetra: Linear Algebra Services Package 
00006 //                 Copyright 2011 Sandia Corporation
00007 // 
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00039 // 
00040 // ************************************************************************
00041 //@HEADER
00042 
00043 #include "Epetra_ConfigDefs.h"
00044 #include "Epetra_BasicRowMatrix.h"
00045 #include "Epetra_Map.h"
00046 #include "Epetra_Import.h"
00047 #include "Epetra_Export.h"
00048 #include "Epetra_Vector.h"
00049 #include "Epetra_MultiVector.h"
00050 #include "Epetra_Comm.h"
00051 #include "Epetra_Util.h"
00052 #include "Epetra_IntSerialDenseVector.h"
00053 #include "Epetra_RowMatrix.h"
00054 #include "Epetra_CrsMatrix.h"
00055 
00056 
00057 //==============================================================================
00058 Epetra_BasicRowMatrix::Epetra_BasicRowMatrix(const Epetra_Comm & comm) 
00059   : Comm_(comm.Clone()),
00060 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
00061     OperatorDomainMap_(Epetra_Map(0,0,comm)),
00062     OperatorRangeMap_(Epetra_Map(0,0,comm)),
00063     RowMatrixRowMap_(Epetra_Map(0,0,comm)),
00064     RowMatrixColMap_(Epetra_Map(0,0,comm)),
00065 #endif
00066     NumMyNonzeros_(0),
00067     NumGlobalNonzeros_(0),
00068     MaxNumEntries_(0),
00069     NormInf_(0.0),
00070     NormOne_(0.0),
00071     UseTranspose_(false),
00072     HasNormInf_(true),
00073     LowerTriangular_(true),
00074     UpperTriangular_(true),
00075     HaveStructureConstants_(false),
00076     HaveNumericConstants_(false),
00077     HaveMaps_(false),
00078     ImportVector_(0),
00079     ExportVector_(0),
00080     Importer_(0),
00081     Exporter_(0)
00082 {
00083   SetLabel("Epetra::BasicRowMatrix");
00084 }
00085 
00086 //==============================================================================
00087 Epetra_BasicRowMatrix::~Epetra_BasicRowMatrix(){
00088 
00089   if (ImportVector_!=0) delete ImportVector_;
00090   ImportVector_=0;
00091   if (ExportVector_!=0) delete ExportVector_;
00092   ExportVector_=0;
00093   if (Importer_!=0) delete Importer_;
00094   Importer_=0;
00095   if (Exporter_!=0) delete Exporter_;
00096   Exporter_=0;
00097   delete Comm_;
00098 }
00099 
00100 //==============================================================================
00101 void Epetra_BasicRowMatrix::SetMaps(const Epetra_Map & RowMap, const Epetra_Map & ColMap) {
00102 
00103   SetMaps(RowMap, ColMap, RowMap, RowMap);
00104 }
00105 
00106 //==============================================================================
00107 void Epetra_BasicRowMatrix::SetMaps(const Epetra_Map & RowMap, const Epetra_Map & ColMap, 
00108       const Epetra_Map & DomainMap, const Epetra_Map & RangeMap) {
00109 
00110   RowMatrixRowMap_ = RowMap;
00111   RowMatrixColMap_ = ColMap;
00112   OperatorDomainMap_ = DomainMap;
00113   OperatorRangeMap_ = RangeMap;
00114   HaveMaps_ = true;
00115   HaveStructureConstants_ = false;
00116   HaveNumericConstants_ = false;
00117 
00118   if (!OperatorDomainMap().UniqueGIDs()) throw RowMatrixRowMap().ReportError("At least one GID is repeated in domain map. Domain and range maps must have unique GIDs", -1);
00119   if (!OperatorRangeMap().UniqueGIDs()) throw RowMatrixRowMap().ReportError("At least one GID is repeated in range map. Domain and range maps must have unique GIDs", -2);
00120   SetImportExport();
00121 }
00122 
00123 //==============================================================================
00124 void Epetra_BasicRowMatrix::SetImportExport() {
00125 
00126   // Check if non-trivial import/export operators
00127   if (!(RowMatrixRowMap().SameAs(OperatorRangeMap()))) 
00128     Exporter_ = new Epetra_Export(RowMatrixRowMap(), OperatorRangeMap());
00129   
00130   if (!(RowMatrixColMap().SameAs(OperatorDomainMap())))
00131     Importer_ = new Epetra_Import(RowMatrixColMap(), OperatorDomainMap());
00132 
00133   NumMyRows_ = RowMatrixRowMap().NumMyPoints();
00134   NumMyCols_ = RowMatrixColMap().NumMyPoints();
00135 }
00136 
00137 //==============================================================================
00138 void Epetra_BasicRowMatrix::ComputeStructureConstants() const {
00139   MaxNumEntries_ = 0;
00140   NumMyNonzeros_ = 0;
00141   NumGlobalNonzeros_ = 0;
00142   int NumEntries = 0;
00143   for (int i=0; i<NumMyRows_; i++) {
00144     NumMyRowEntries(i, NumEntries);
00145     NumMyNonzeros_ += NumEntries;
00146     if (NumEntries>MaxNumEntries_) MaxNumEntries_ = NumEntries;
00147   }
00148 
00149   long long tmp_NumMyNonzeros = NumMyNonzeros_;
00150   RowMatrixRowMap().Comm().SumAll(&tmp_NumMyNonzeros, &NumGlobalNonzeros_, 1);
00151   HaveStructureConstants_ = true;
00152 }
00153 //=============================================================================
00154 void Epetra_BasicRowMatrix::ComputeNumericConstants() const {
00155   Epetra_SerialDenseVector Values(MaxNumEntries());
00156   Epetra_IntSerialDenseVector Indices(MaxNumEntries());
00157   int NumEntries;
00158   Epetra_Vector x1(RowMatrixRowMap()); // Need temp vector for row sums
00159   Epetra_Vector x2(RowMatrixColMap()); // Need temp vector for column sums
00160   for(int i = 0; i < NumMyRows_; i++) {
00161     ExtractMyRowCopy(i, MaxNumEntries(), NumEntries, Values.Values(), Indices.Values());
00162     for(int j = 0; j < NumEntries; j++) {
00163       x1[i] += std::abs(Values[j]);
00164       x2[Indices[j]] += std::abs(Values[j]);
00165       if (Indices[j]<i) UpperTriangular_ = false;
00166       if (Indices[j]>i) LowerTriangular_ = false;
00167     }
00168   }
00169 
00170   // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
00171   if(Exporter() != 0) {
00172     Epetra_Vector xtmp(OperatorRangeMap()); // Create temporary import vector if needed
00173     xtmp.Export(x1,*Exporter(),Add);
00174     xtmp.MaxValue(&NormInf_); // This is the NormInf
00175   }
00176   else
00177     x1.MaxValue(&NormInf_); // Find max
00178 
00179   // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
00180   if(Importer() != 0) {
00181     Epetra_Vector xtmp(OperatorDomainMap()); // Create temporary import vector if needed
00182     xtmp.Export(x2,*Importer(),Add);
00183     xtmp.MaxValue(&NormOne_); // This is the NormOne
00184   }
00185   else
00186     x2.MaxValue(&NormOne_); // Find max
00187 
00188   UpdateFlops(2*NumGlobalNonzeros64());
00189   HaveNumericConstants_ = true;
00190 }
00191 //=============================================================================
00192 int Epetra_BasicRowMatrix::ExtractDiagonalCopy(Epetra_Vector & Diagonal) const {
00193 
00194   if(!RowMatrixRowMap().SameAs(Diagonal.Map())) 
00195     EPETRA_CHK_ERR(-2); // Maps must be the same
00196 
00197   // Crude implementation in terms of ExtractMyRowCopy
00198 
00199   Epetra_SerialDenseVector Values(MaxNumEntries());
00200   Epetra_IntSerialDenseVector Indices(MaxNumEntries());
00201   int NumEntries;
00202 
00203   for(int i = 0; i < NumMyRows_; i++) {
00204     EPETRA_CHK_ERR(ExtractMyRowCopy(i, MaxNumEntries(), NumEntries, Values.Values(), Indices.Values()));
00205     long long ii = RowMatrixRowMap().GID64(i);
00206     
00207     Diagonal[i] = 0.0;
00208     for(int j = 0; j < NumEntries; j++) {
00209       if(ii == RowMatrixColMap().GID64(Indices[j])) {
00210   Diagonal[i] = Values[j];
00211   break;
00212       }
00213     }
00214   }
00215   return(0);
00216 }
00217 //=============================================================================
00218 int Epetra_BasicRowMatrix::InvRowSums(Epetra_Vector & x) const {
00219   int ierr = 0;
00220   int i, j;
00221   Epetra_SerialDenseVector Values(MaxNumEntries());
00222   Epetra_IntSerialDenseVector Indices(MaxNumEntries());
00223   int NumEntries;
00224   x.PutScalar(0.0); // Make sure we sum into a vector of zeros.
00225   double * xp = (double*)x.Values();
00226   if (OperatorRangeMap().SameAs(x.Map()) && Exporter() != 0) {
00227     Epetra_Vector x_tmp(RowMatrixRowMap());
00228     x_tmp.PutScalar(0.0);
00229     double * x_tmp_p = (double*)x_tmp.Values();
00230     for (i=0; i < NumMyRows_; i++) {
00231       EPETRA_CHK_ERR(ExtractMyRowCopy(i, MaxNumEntries(), NumEntries, Values.Values(), Indices.Values()));
00232       for (j=0; j < NumEntries; j++)  x_tmp_p[i] += std::abs(Values[j]);
00233     }
00234     EPETRA_CHK_ERR(x.Export(x_tmp, *Exporter(), Add)); //Export partial row sums to x.
00235     int myLength = x.MyLength();
00236     for (i=0; i<myLength; i++) { 
00237       if (xp[i]<Epetra_MinDouble) {
00238         if (xp[i]==0.0) ierr = 1; // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
00239         else if (ierr!=1) ierr = 2;
00240         xp[i] = Epetra_MaxDouble;
00241       }
00242       else
00243         xp[i] = 1.0/xp[i];
00244     }
00245   }
00246   else if (RowMatrixRowMap().SameAs(x.Map())) {
00247     for (i=0; i < NumMyRows_; i++) {
00248       EPETRA_CHK_ERR(ExtractMyRowCopy(i, MaxNumEntries(), NumEntries, Values.Values(), Indices.Values()));
00249       double scale = 0.0;
00250       for (j=0; j < NumEntries; j++) scale += std::abs(Values[j]);
00251       if (scale<Epetra_MinDouble) {
00252         if (scale==0.0) ierr = 1; // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
00253         else if (ierr!=1) ierr = 2;
00254         xp[i] = Epetra_MaxDouble;
00255       }
00256       else
00257         xp[i] = 1.0/scale;
00258     }
00259   }
00260   else { // x.Map different than both RowMatrixRowMap() and OperatorRangeMap()
00261     EPETRA_CHK_ERR(-2); // The map of x must be the RowMap or RangeMap of A.
00262   }
00263   EPETRA_CHK_ERR(ierr);  
00264   UpdateFlops(NumGlobalNonzeros64());
00265   return(0);
00266 }
00267 //=============================================================================
00268 int Epetra_BasicRowMatrix::LeftScale(const Epetra_Vector & x) {
00269   double *curValue;
00270   int curRowIndex, curColIndex;
00271   if(OperatorRangeMap().SameAs(x.Map()) && Exporter() != 0) {
00272     Epetra_Vector xtmp(RowMatrixRowMap());
00273     xtmp.Import(x,*Exporter(),Insert);
00274     for (int i=0; i<NumMyNonzeros_; i++) {
00275       EPETRA_CHK_ERR(ExtractMyEntryView(i, curValue, curRowIndex, curColIndex));
00276       *curValue *= xtmp[curRowIndex];
00277     }
00278   }
00279   else if (RowMatrixRowMap().SameAs(x.Map()))
00280     for (int i=0; i<NumMyNonzeros_; i++) {
00281       EPETRA_CHK_ERR(ExtractMyEntryView(i, curValue, curRowIndex, curColIndex));
00282       *curValue *= x[curRowIndex];
00283     }
00284   else {
00285     EPETRA_CHK_ERR(-2); // The Map of x must be the RowMap or RangeMap of A.
00286   }
00287   HaveNumericConstants_ = false;
00288   UpdateFlops(NumGlobalNonzeros64());
00289   return(0);
00290 }
00291 //=============================================================================
00292 int Epetra_BasicRowMatrix::InvColSums(Epetra_Vector & x) const {
00293   int ierr = 0;
00294   int i, j;
00295   Epetra_SerialDenseVector Values(MaxNumEntries());
00296   Epetra_IntSerialDenseVector Indices(MaxNumEntries());
00297   int NumEntries;
00298   int MapNumMyElements = x.Map().NumMyElements();
00299   x.PutScalar(0.0); // Make sure we sum into a vector of zeros.
00300   double* xp = (double*)x.Values();
00301   if(OperatorDomainMap().SameAs(x.Map()) && Importer() != 0) {
00302     Epetra_Vector x_tmp(RowMatrixColMap());
00303     x_tmp.PutScalar(0.0);
00304     double * x_tmp_p = (double*)x_tmp.Values();
00305     for(i = 0; i < NumMyRows_; i++) {
00306       EPETRA_CHK_ERR(ExtractMyRowCopy(i, MaxNumEntries(), NumEntries, Values.Values(), Indices.Values()));
00307       for(j = 0; j < NumEntries; j++) 
00308         x_tmp_p[Indices[j]] += std::abs(Values[j]);
00309     }
00310     EPETRA_CHK_ERR(x.Export(x_tmp, *Importer(), Add)); // Fill x with partial column sums
00311   }
00312   else if(RowMatrixColMap().SameAs(x.Map())) {
00313     for(i = 0; i < NumMyRows_; i++) {
00314       EPETRA_CHK_ERR(ExtractMyRowCopy(i, MaxNumEntries(), NumEntries, Values.Values(), Indices.Values()));
00315       for(j = 0; j < NumEntries; j++) 
00316         xp[Indices[j]] += std::abs(Values[j]);
00317     }
00318   }
00319   else { //x.Map different than both RowMatrixColMap() and OperatorDomainMap()
00320     EPETRA_CHK_ERR(-2); // x must have the same distribution as the domain of A
00321   }
00322 
00323   // Invert values, don't allow them to get too large
00324   for(i = 0; i < MapNumMyElements; i++) {
00325     double scale = xp[i];
00326     if(scale < Epetra_MinDouble) {
00327       if(scale == 0.0) 
00328   ierr = 1; // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
00329       else if(ierr != 1) 
00330   ierr = 2;
00331       xp[i] = Epetra_MaxDouble;
00332     }
00333     else
00334       xp[i] = 1.0 / scale;
00335   }
00336 
00337   EPETRA_CHK_ERR(ierr);
00338   UpdateFlops(NumGlobalNonzeros64());
00339   return(0);
00340 }
00341 //=============================================================================
00342 int Epetra_BasicRowMatrix::RightScale(const Epetra_Vector & x) {
00343   double *curValue;
00344   int curRowIndex, curColIndex;
00345   if(OperatorDomainMap().SameAs(x.Map()) && Importer() != 0) {
00346     Epetra_Vector xtmp(RowMatrixColMap());
00347     xtmp.Import(x,*Importer(),Insert);
00348     for (int i=0; i<NumMyNonzeros_; i++) {
00349       EPETRA_CHK_ERR(ExtractMyEntryView(i, curValue, curRowIndex, curColIndex));
00350       *curValue *= xtmp[curColIndex];
00351     }
00352   }
00353   else if (RowMatrixColMap().SameAs(x.Map()))
00354     for (int i=0; i<NumMyNonzeros_; i++) {
00355       EPETRA_CHK_ERR(ExtractMyEntryView(i, curValue, curRowIndex, curColIndex));
00356       *curValue *= x[curColIndex];
00357     }
00358   else {
00359     EPETRA_CHK_ERR(-2); // The Map of x must be the RowMap or RangeMap of A.
00360   }
00361   HaveNumericConstants_ = false;
00362   UpdateFlops(NumGlobalNonzeros64());
00363   return(0);
00364 }
00365 //=============================================================================
00366 int Epetra_BasicRowMatrix::Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
00367   //
00368   // This function forms the product Y = A * Y or Y = A' * X
00369   //
00370 
00371   Epetra_SerialDenseVector Values(MaxNumEntries());
00372   Epetra_IntSerialDenseVector Indices(MaxNumEntries());
00373   int NumEntries;
00374 
00375   int NumVectors = X.NumVectors();
00376   if (NumVectors!=Y.NumVectors()) {
00377     EPETRA_CHK_ERR(-1); // Need same number of vectors in each MV
00378   }
00379 
00380   UpdateImportVector(NumVectors); // Make sure Import and Export Vectors are compatible
00381   UpdateExportVector(NumVectors);
00382 
00383   double ** Xp = (double**) X.Pointers();
00384   double ** Yp = (double**) Y.Pointers();
00385 
00386   if (!TransA) {
00387 
00388     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
00389     if (Importer()!=0) {
00390       EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer(), Insert));
00391       Xp = (double**)ImportVector_->Pointers();
00392     }
00393 
00394     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
00395     if (Exporter()!=0) {
00396       Yp = (double**)ExportVector_->Pointers();
00397     }
00398 
00399     // Do actual computation
00400     for(int i = 0; i < NumMyRows_; i++) {
00401       EPETRA_CHK_ERR(ExtractMyRowCopy(i, MaxNumEntries(), NumEntries, Values.Values(), Indices.Values()));
00402       for (int k=0; k<NumVectors; k++) {
00403   double sum = 0.0;
00404   for(int j = 0; j < NumEntries; j++)
00405     sum += Values[j]*Xp[k][Indices[j]];
00406   Yp[k][i] = sum;
00407       }
00408     }
00409     
00410     if (Exporter()!=0) {
00411       Y.PutScalar(0.0);  // Make sure target is zero
00412       Y.Export(*ExportVector_, *Exporter(), Add); // Fill Y with Values from export vector
00413     }
00414     // Handle case of rangemap being a local replicated map
00415     if (!OperatorRangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
00416   }
00417   else { // Transpose operation
00418     
00419 
00420     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
00421 
00422     if (Exporter()!=0) {
00423       EPETRA_CHK_ERR(ExportVector_->Import(X, *Exporter(), Insert));
00424       Xp = (double**)ExportVector_->Pointers();
00425     }
00426 
00427     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
00428     if (Importer()!=0) {
00429       Yp = (double**)ImportVector_->Pointers();
00430       ImportVector_->PutScalar(0.0);  // Make sure target is zero
00431     }
00432     else Y.PutScalar(0.0); // Make sure target is zero
00433 
00434     // Do actual computation
00435     for(int i = 0; i < NumMyRows_; i++) {
00436       EPETRA_CHK_ERR(ExtractMyRowCopy(i, MaxNumEntries(), NumEntries, Values.Values(), Indices.Values()));
00437       for (int k=0; k<NumVectors; k++) {
00438   double xtmp = Xp[k][i];
00439   for(int j = 0; j < NumEntries; j++)
00440     Yp[k][Indices[j]] += Values[j]*xtmp;
00441       }
00442     }
00443     
00444     if (Importer()!=0) {
00445       Y.PutScalar(0.0);  // Make sure target is zero
00446       EPETRA_CHK_ERR(Y.Export(*ImportVector_, *Importer(), Add)); // Fill Y with Values from export vector
00447     }
00448     // Handle case of rangemap being a local replicated map
00449     if (!OperatorDomainMap().DistributedGlobal() && Comm().NumProc()>1)  EPETRA_CHK_ERR(Y.Reduce());
00450   }
00451 
00452   UpdateFlops(2*NumVectors*NumGlobalNonzeros64());
00453   return(0);
00454 }
00455 //=======================================================================================================
00456 void Epetra_BasicRowMatrix::UpdateImportVector(int NumVectors) const {
00457   if(Importer() != 0) {
00458     if(ImportVector_ != 0) {
00459       if(ImportVector_->NumVectors() != NumVectors) {
00460      delete ImportVector_;
00461      ImportVector_= 0;
00462       }
00463     }
00464     if(ImportVector_ == 0)
00465       ImportVector_ = new Epetra_MultiVector(Importer_->TargetMap(),NumVectors); // Create import vector if needed
00466   }
00467   return;
00468 }
00469 //=======================================================================================================
00470 void Epetra_BasicRowMatrix::UpdateExportVector(int NumVectors) const {
00471   if(Exporter() != 0) {
00472     if(ExportVector_ != 0) {
00473       if(ExportVector_->NumVectors() != NumVectors) {
00474      delete ExportVector_;
00475      ExportVector_= 0;
00476       }
00477     }
00478     if(ExportVector_ == 0)
00479       ExportVector_ = new Epetra_MultiVector(Exporter_->SourceMap(),NumVectors); // Create Export vector if needed
00480   }
00481   return;
00482 }
00483 
00484 //=======================================================================================================
00485 void Epetra_BasicRowMatrix::Print(std::ostream& os) const {
00486 
00487   int MyPID = RowMatrixRowMap().Comm().MyPID();
00488   int NumProc = RowMatrixRowMap().Comm().NumProc();
00489 
00490   for (int iproc=0; iproc < NumProc; iproc++) {
00491     if (MyPID==iproc) {
00492       if (MyPID==0) {
00493     os <<    "\nNumber of Global Rows         = "; os << NumGlobalRows64();    os << std::endl;
00494     os <<    "Number of Global Cols         = "; os << NumGlobalCols64();    os << std::endl;
00495     os <<    "Number of Global Diagonals    = "; os << NumGlobalDiagonals64(); os << std::endl;
00496   os <<    "Number of Global Nonzeros     = "; os << NumGlobalNonzeros_; os << std::endl;
00497       }
00498       
00499       os <<  "\nNumber of My Rows               = "; os << NumMyRows_; os << std::endl;
00500       os <<    "Number of My Cols               = "; os << NumMyCols_; os << std::endl;
00501       os <<    "Number of My Diagonals          = "; os << NumMyDiagonals(); os << std::endl;
00502       os <<    "Number of My Nonzeros           = "; os << NumMyNonzeros_; os << std::endl;
00503       os <<    "My Maximum Num Entries          = "; os << MaxNumEntries_; os << std::endl; os << std::endl;
00504       os << std::flush;
00505       
00506     }
00507     // Do a few global ops to give I/O a chance to complete
00508     Comm().Barrier();
00509     Comm().Barrier();
00510     Comm().Barrier();
00511   }
00512   
00513   for (int iproc=0; iproc < NumProc; iproc++) {
00514     if (MyPID==iproc) {
00515       if (MyPID==0) {
00516   os.width(8);
00517   os <<  "   Processor ";
00518   os.width(10);
00519   os <<  "   Row Index ";
00520   os.width(10);
00521   os <<  "   Col Index ";
00522   os.width(20);
00523   os <<  "   Value     ";
00524   os << std::endl;
00525       }
00526       Epetra_SerialDenseVector Values(MaxNumEntries());
00527       Epetra_IntSerialDenseVector Indices(MaxNumEntries());
00528       int NumEntries;
00529       
00530       for(int i = 0; i < NumMyRows_; i++) {
00531   ExtractMyRowCopy(i, MaxNumEntries(), NumEntries, Values.Values(), Indices.Values());
00532   long long Row = RowMatrixRowMap().GID64(i); // Get global row number
00533   
00534   for (int j = 0; j < NumEntries ; j++) {   
00535     long long Index = RowMatrixColMap().GID64(Indices[j]);
00536     os.width(8);
00537     os <<  MyPID ; os << "    ";  
00538     os.width(10);
00539     os <<  Row ; os << "    ";  
00540     os.width(10);
00541     os <<  Index; os << "    ";
00542     os.width(20);
00543     os <<  Values[j]; os << "    ";
00544     os << std::endl;
00545   }
00546       }
00547     
00548       os << std::flush;
00549       
00550     }
00551     // Do a few global ops to give I/O a chance to complete
00552     Comm().Barrier();
00553     Comm().Barrier();
00554     Comm().Barrier();
00555   }
00556   
00557   return;
00558 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines