EpetraExt Development
EpetraExt_CrsSingletonFilter_LinearProblem.cpp
Go to the documentation of this file.
00001 
00002 //@HEADER
00003 // ***********************************************************************
00004 //
00005 //     EpetraExt: Epetra Extended - 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_Map.h"
00045 #include "Epetra_Util.h"
00046 #include "Epetra_Export.h"
00047 #include "Epetra_Import.h"
00048 #include "Epetra_MultiVector.h"
00049 #include "Epetra_Vector.h"
00050 #include "Epetra_GIDTypeVector.h"
00051 #include "Epetra_Comm.h"
00052 #include "Epetra_LinearProblem.h"
00053 #include "Epetra_MapColoring.h"
00054 #include "EpetraExt_CrsSingletonFilter_LinearProblem.h"
00055 
00056 namespace EpetraExt {
00057 
00058 //==============================================================================
00059 LinearProblem_CrsSingletonFilter::
00060 LinearProblem_CrsSingletonFilter( bool verbose )
00061 : verbose_(verbose)
00062 {
00063   InitializeDefaults();
00064 }
00065 //==============================================================================
00066 LinearProblem_CrsSingletonFilter::
00067 ~LinearProblem_CrsSingletonFilter()
00068 {
00069   if (ReducedLHS_!=0) delete ReducedLHS_;
00070   if (ReducedRHS_!=0) delete ReducedRHS_;
00071   if (ReducedMatrixDomainMap_!=ReducedMatrixColMap_) delete ReducedMatrixDomainMap_;
00072   if (OrigReducedMatrixDomainMap_!=ReducedMatrixColMap_ &&
00073       OrigReducedMatrixDomainMap_!=0) delete OrigReducedMatrixDomainMap_;
00074   if (ReducedMatrixRangeMap_!=ReducedMatrixRowMap_) delete ReducedMatrixRangeMap_;
00075   if (ReducedMatrixRowMap_!=0) delete ReducedMatrixRowMap_;
00076   if (ReducedMatrixColMap_!=0) delete ReducedMatrixColMap_;
00077   if (Full2ReducedRHSImporter_!=0) delete Full2ReducedRHSImporter_;
00078   if (Full2ReducedLHSImporter_!=0) delete Full2ReducedLHSImporter_;
00079   if (RedistributeDomainExporter_!=0) delete RedistributeDomainExporter_;
00080   if (RowMapColors_!=0) delete RowMapColors_;
00081   if (ColMapColors_!=0) delete ColMapColors_;
00082 
00083   if (ColSingletonRowLIDs_ != 0) delete [] ColSingletonRowLIDs_;
00084   if (ColSingletonColLIDs_ != 0) delete [] ColSingletonColLIDs_;
00085   if (ColSingletonPivotLIDs_ != 0) delete [] ColSingletonPivotLIDs_;
00086   if (ColSingletonPivots_ != 0) delete [] ColSingletonPivots_;
00087   if (tempExportX_ != 0) delete tempExportX_;
00088   if (Indices_int_ != 0) delete [] Indices_int_;
00089 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00090   if (Indices_LL_ != 0) delete [] Indices_LL_;
00091 #endif
00092   if (tempX_ != 0) delete tempX_;
00093   if (tempB_ != 0) delete tempB_;
00094 
00095 }
00096 
00097 //==============================================================================
00098 LinearProblem_CrsSingletonFilter::NewTypeRef
00099 LinearProblem_CrsSingletonFilter::
00100 operator()( LinearProblem_CrsSingletonFilter::OriginalTypeRef orig )
00101 {
00102   analyze( orig );
00103   return construct();
00104 }
00105 
00106 //==============================================================================
00107 bool
00108 LinearProblem_CrsSingletonFilter::
00109 analyze( LinearProblem_CrsSingletonFilter::OriginalTypeRef orig )
00110 {
00111   origObj_ = &orig;
00112 
00113   FullMatrix_ = orig.GetMatrix();
00114 
00115   int flag = Analyze( FullMatrix_ );
00116   assert( flag >= 0 );
00117 
00118   if ( verbose_ && FullMatrix_->Comm().MyPID()==0 ) {
00119     std::cout << "\nAnalyzed Singleton Problem:\n";
00120     std::cout << "---------------------------\n";
00121   }
00122   if ( SingletonsDetected() ) {
00123     if ( verbose_ && FullMatrix_->Comm().MyPID()==0 ) {
00124       std::cout << "Singletons Detected!" << std::endl;;
00125       std::cout << "Num Singletons:      " << NumSingletons() << std::endl;
00126     }
00127   }
00128   else {
00129     if ( verbose_ && FullMatrix_->Comm().MyPID()==0 ) 
00130         std::cout << "No Singletons Detected!" << std::endl;
00131   }
00132 /*
00133     std::cout << "List of Row Singletons:\n";
00134     int * slist = RowMapColors_->ColorLIDList(1);
00135     for( int i = 0; i < RowMapColors_->NumElementsWithColor(1); ++i )
00136       std::cout << slist[i] << " ";
00137     std::cout << "\n";
00138     std::cout << "List of Col Singletons:\n";
00139     slist = ColMapColors_->ColorLIDList(1);
00140     for( int i = 0; i < ColMapColors_->NumElementsWithColor(1); ++i )
00141       std::cout << slist[i] << " ";
00142     std::cout << "\n";
00143 */
00144   if ( verbose_ && FullMatrix_->Comm().MyPID()==0 )
00145     std::cout << "---------------------------\n\n";
00146  
00147   return true;
00148 }
00149 
00150 //==============================================================================
00151 LinearProblem_CrsSingletonFilter::NewTypeRef
00152 LinearProblem_CrsSingletonFilter::
00153 construct()
00154 {
00155   if( !origObj_ ) abort();
00156 
00157   int flag = ConstructReducedProblem( origObj_ );
00158   assert( flag >= 0 );
00159 
00160   newObj_ = ReducedProblem();
00161 
00162   if( verbose_ && SingletonsDetected() && FullMatrix_->Comm().MyPID()==0 )
00163   {
00164     std::cout << "\nConstructedSingleton Problem:\n";
00165     std::cout << "---------------------------\n";
00166     std::cout << "RatioOfDimensions:   " << RatioOfDimensions() << std::endl;
00167     std::cout << "RatioOfNonzeros:     " << RatioOfNonzeros() << std::endl;
00168     std::cout << "---------------------------\n\n";
00169   }
00170 
00171   return *newObj_;
00172 }
00173 
00174 //==============================================================================
00175 bool LinearProblem_CrsSingletonFilter::fwd()
00176 {
00177   int ierr = UpdateReducedProblem( FullProblem_ );
00178   if( ierr ) std::cout << "EDT_LinearProblem_CrsSingletonFilter::UpdateReducedProblem FAILED!\n";
00179 
00180   return (ierr==0);
00181 }
00182 
00183 //==============================================================================
00184 bool LinearProblem_CrsSingletonFilter::rvs()
00185 {
00186   int ierr = ComputeFullSolution();
00187   if( ierr ) std::cout << "EDT_LinearProblem_CrsSingletonFilter::ComputeFullSolution FAILED!\n";
00188 
00189   return (ierr==0);
00190 }
00191 
00192 //==============================================================================
00193 void LinearProblem_CrsSingletonFilter::InitializeDefaults() { 
00194 
00195 // Initialize all attributes that have trivial default values
00196 
00197   FullProblem_ = 0;
00198   FullMatrix_ = 0;
00199   ReducedRHS_ = 0;
00200   ReducedLHS_ = 0;
00201   ReducedMatrixRowMap_ = 0;
00202   ReducedMatrixColMap_ = 0;
00203   ReducedMatrixDomainMap_ = 0;
00204   ReducedMatrixRangeMap_ = 0;
00205   OrigReducedMatrixDomainMap_ = 0;
00206   Full2ReducedRHSImporter_ = 0;
00207   Full2ReducedLHSImporter_ = 0;
00208   RedistributeDomainExporter_ = 0;
00209 
00210   ColSingletonRowLIDs_ = 0;
00211   ColSingletonColLIDs_ = 0;
00212   ColSingletonPivotLIDs_ = 0;
00213   ColSingletonPivots_ = 0;
00214 
00215   AbsoluteThreshold_ = 0;
00216   RelativeThreshold_ = 0;
00217 
00218   NumMyRowSingletons_ = -1;
00219   NumMyColSingletons_ = -1;
00220   NumGlobalRowSingletons_ = -1;
00221   NumGlobalColSingletons_ = -1;
00222   RatioOfDimensions_ = -1.0;
00223   RatioOfNonzeros_ = -1.0;
00224 
00225   HaveReducedProblem_ = false;
00226   UserDefinedEliminateMaps_ = false;
00227   AnalysisDone_ = false;
00228   SymmetricElimination_ = true;
00229 
00230   tempExportX_ = 0;
00231   tempX_ = 0;
00232   tempB_ = 0;
00233 
00234   Indices_int_ = 0;
00235 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00236   Indices_LL_ = 0;
00237 #endif
00238 
00239   RowMapColors_ = 0;
00240   ColMapColors_ = 0;
00241 
00242   FullMatrixIsCrsMatrix_ = false;
00243   MaxNumMyEntries_ = 0;
00244   return;
00245 }
00246 //==============================================================================
00247 int LinearProblem_CrsSingletonFilter::Analyze(Epetra_RowMatrix * FullMatrix) {
00248 
00249   int i, j, jj;
00250 
00251   FullMatrix_ = FullMatrix; 
00252 
00253   if (AnalysisDone_) EPETRA_CHK_ERR(-1); // Analysis already done once.  Cannot do it again
00254   if (FullMatrix==0) EPETRA_CHK_ERR(-2); // Input matrix pointer is zero
00255   if (FullMatrix->NumGlobalRows64()==0) EPETRA_CHK_ERR(-3); // Full matrix has zero dimension.
00256   if (FullMatrix->NumGlobalNonzeros64()==0) EPETRA_CHK_ERR(-4); // Full matrix has no nonzero terms.
00257 
00258   // First check for columns with single entries and find columns with singleton rows
00259   Epetra_IntVector ColProfiles(FullMatrixColMap()); ColProfiles.PutValue(0);
00260   Epetra_IntVector ColHasRowWithSingleton(FullMatrixColMap()); ColHasRowWithSingleton.PutValue(0);
00261 
00262   // RowIDs[j] will contain the local row ID associated with the jth column, 
00263   // if the jth col has a single entry
00264   Epetra_IntVector RowIDs(FullMatrixColMap()); RowIDs.PutValue(-1);
00265 
00266   // Define MapColoring objects
00267   RowMapColors_ = new Epetra_MapColoring(FullMatrixRowMap());  // Initial colors are all 0
00268   ColMapColors_ = new Epetra_MapColoring(FullMatrixColMap());
00269   Epetra_MapColoring & RowMapColors = *RowMapColors_;
00270   Epetra_MapColoring & ColMapColors = *ColMapColors_;
00271 
00272   int NumMyRows = FullMatrix->NumMyRows();
00273   int NumMyCols = FullMatrix->NumMyCols();
00274 
00275   // Set up for accessing full matrix.  Will do so row-by-row.
00276   EPETRA_CHK_ERR(InitFullMatrixAccess()); 
00277 
00278   // Scan matrix for singleton rows, build up column profiles
00279   int NumIndices;
00280   int * Indices;
00281   NumMyRowSingletons_ = 0;
00282   for (i=0; i<NumMyRows; i++) {
00283     // Get ith row
00284     EPETRA_CHK_ERR(GetRow(i, NumIndices, Indices));
00285     for (j=0; j<NumIndices; j++) {
00286       int ColumnIndex = Indices[j];
00287       ColProfiles[ColumnIndex]++; // Increment column count
00288       // Record local row ID for current column
00289       // will use to identify row to eliminate if column is a singleton
00290       RowIDs[ColumnIndex] = i;
00291     }
00292     // If row has single entry, color it and associated column with color=1
00293     if (NumIndices==1) {
00294       int j = Indices[0];
00295       ColHasRowWithSingleton[j]++;
00296       RowMapColors[i] = 1;
00297       ColMapColors[j] = 1;
00298       NumMyRowSingletons_++;
00299     }
00300   }
00301 
00302   // 1) The vector ColProfiles has column nonzero counts for each processor's contribution
00303   // Combine these to get total column profile information and then redistribute to processors
00304   // so each can determine if it is the owner of the row associated with the singleton column
00305   // 2) The vector ColHasRowWithSingleton[i] contain count of singleton rows  that are associated with 
00306   // the ith column on this processor.  Must tell other processors that they should also eliminate 
00307   // these columns.
00308 
00309   // Make a copy of ColProfiles for later use when detecting columns that disappear locally
00310 
00311   Epetra_IntVector NewColProfiles(ColProfiles);
00312 
00313   // If RowMatrixImporter is non-trivial, we need to perform a gather/scatter to accumulate results
00314 
00315   if (FullMatrix->RowMatrixImporter()!=0) {
00316     Epetra_IntVector tmpVec(FullMatrixDomainMap()); // Use for gather/scatter of column vectors
00317     EPETRA_CHK_ERR(tmpVec.PutValue(0));
00318     EPETRA_CHK_ERR(tmpVec.Export(ColProfiles, *FullMatrix->RowMatrixImporter(), Add));
00319     EPETRA_CHK_ERR(ColProfiles.Import(tmpVec, *FullMatrix->RowMatrixImporter(), Insert));
00320     
00321     EPETRA_CHK_ERR(tmpVec.PutValue(0));
00322     EPETRA_CHK_ERR(tmpVec.Export(ColHasRowWithSingleton, *FullMatrix->RowMatrixImporter(), Add));
00323     EPETRA_CHK_ERR(ColHasRowWithSingleton.Import(tmpVec, *FullMatrix->RowMatrixImporter(), Insert));
00324   }
00325   // ColProfiles now contains the nonzero column entry count for all columns that have
00326   // an entry on this processor.
00327   // ColHasRowWithSingleton now contains a count of singleton rows associated with the corresponding
00328   // local column.  Next we check to make sure no column is associated with more than one singleton row.
00329 
00330   if (ColHasRowWithSingleton.MaxValue()>1) {
00331     EPETRA_CHK_ERR(-2); // At least one col is associated with two singleton rows, can't handle it.
00332   }
00333 
00334   Epetra_IntVector RowHasColWithSingleton(FullMatrix->RowMatrixRowMap()); // Use to check for errors
00335   RowHasColWithSingleton.PutValue(0);
00336  
00337   NumMyColSingletons_ = 0;
00338   // Count singleton columns (that were not already counted as singleton rows)
00339   for (j=0; j<NumMyCols; j++) {
00340     int i = RowIDs[j];
00341     // Check if column is a singleton
00342     if (ColProfiles[j]==1 ) {
00343       // Check to make sure RowID is not invalid
00344 //      assert(i!=-1);
00345       // Check to see if this column already eliminated by the row check above
00346       if (RowMapColors[i]!=1) {
00347   RowHasColWithSingleton[i]++; // Increment col singleton counter for ith row
00348   RowMapColors[i] = 2; // Use 2 for now, to distinguish between row eliminated directly or via column singletons
00349   ColMapColors[j] = 1;
00350   NumMyColSingletons_++;
00351   // If we delete a row, we need to keep track of associated column entries that were also deleted 
00352   // in case all entries in a column are eventually deleted, in which case the column should
00353   // also be deleted.
00354   EPETRA_CHK_ERR(GetRow(i, NumIndices, Indices));
00355   for (jj=0; jj<NumIndices; jj++) NewColProfiles[Indices[jj]]--;
00356   
00357       }
00358     }
00359     // Check if some other processor eliminated this column    
00360     else if (ColHasRowWithSingleton[j]==1 && RowMapColors[i]!=1) { 
00361   ColMapColors[j] = 1;
00362     }
00363   }
00364   if (RowHasColWithSingleton.MaxValue()>1) {
00365     EPETRA_CHK_ERR(-3); // At lease one row is associated with two singleton cols, can't handle it.
00366   }
00367 
00368  // Generate arrays that keep track of column singleton row, col and pivot info needed for post-solve phase
00369   EPETRA_CHK_ERR(CreatePostSolveArrays(RowIDs, RowMapColors, ColProfiles, NewColProfiles,
00370                ColHasRowWithSingleton));
00371 
00372   for (i=0; i<NumMyRows; i++) if (RowMapColors[i]==2) RowMapColors[i] = 1; // Convert all eliminated rows to same color
00373 
00374   FullMatrix->RowMatrixRowMap().Comm().SumAll(&NumMyRowSingletons_, &NumGlobalRowSingletons_, 1);
00375   FullMatrix->RowMatrixRowMap().Comm().SumAll(&NumMyColSingletons_, &NumGlobalColSingletons_, 1);
00376   AnalysisDone_ = true;
00377   return(0);
00378 }
00379 //==============================================================================
00380 template<typename int_type>
00381 int LinearProblem_CrsSingletonFilter::TConstructReducedProblem(Epetra_LinearProblem * Problem) {
00382   int i, j;
00383   if (HaveReducedProblem_) EPETRA_CHK_ERR(-1); // Setup already done once.  Cannot do it again
00384   if (Problem==0) EPETRA_CHK_ERR(-2); // Null problem pointer
00385 
00386   FullProblem_ = Problem;
00387   FullMatrix_ = dynamic_cast<Epetra_RowMatrix *>(Problem->GetMatrix());
00388   if (FullMatrix_==0) EPETRA_CHK_ERR(-3); // Need a RowMatrix
00389   if (Problem->GetRHS()==0) EPETRA_CHK_ERR(-4); // Need a RHS
00390   if (Problem->GetLHS()==0) EPETRA_CHK_ERR(-5); // Need a LHS
00391   // Generate reduced row and column maps
00392 
00393   if ( SingletonsDetected() ) {
00394 
00395     Epetra_MapColoring & RowMapColors = *RowMapColors_;
00396     Epetra_MapColoring & ColMapColors = *ColMapColors_;
00397 
00398     ReducedMatrixRowMap_ = RowMapColors.GenerateMap(0);
00399     ReducedMatrixColMap_ = ColMapColors.GenerateMap(0);
00400 
00401     // Create domain and range map colorings by exporting map coloring of column and row maps
00402 
00403     if (FullMatrix()->RowMatrixImporter()!=0) {
00404       Epetra_MapColoring DomainMapColors(FullMatrixDomainMap());
00405       EPETRA_CHK_ERR(DomainMapColors.Export(*ColMapColors_, *FullMatrix()->RowMatrixImporter(), AbsMax));
00406       OrigReducedMatrixDomainMap_ = DomainMapColors.GenerateMap(0);
00407     }
00408     else
00409       OrigReducedMatrixDomainMap_ = ReducedMatrixColMap_;
00410 
00411     if (FullMatrixIsCrsMatrix_) {
00412       if (FullCrsMatrix()->Exporter()!=0) { // Non-trivial exporter
00413         Epetra_MapColoring RangeMapColors(FullMatrixRangeMap());
00414         EPETRA_CHK_ERR(RangeMapColors.Export(*RowMapColors_, *FullCrsMatrix()->Exporter(), 
00415              AbsMax));
00416         ReducedMatrixRangeMap_ = RangeMapColors.GenerateMap(0);
00417       }
00418       else
00419         ReducedMatrixRangeMap_ = ReducedMatrixRowMap_;
00420     }
00421     else
00422       ReducedMatrixRangeMap_ = ReducedMatrixRowMap_;
00423 
00424     // Check to see if the reduced system domain and range maps are the same.
00425     // If not, we need to remap entries of the LHS multivector so that they are distributed
00426     // conformally with the rows of the reduced matrix and the RHS multivector
00427     SymmetricElimination_ = ReducedMatrixRangeMap_->SameAs(*OrigReducedMatrixDomainMap_);
00428     if (!SymmetricElimination_) 
00429       ConstructRedistributeExporter(OrigReducedMatrixDomainMap_, ReducedMatrixRangeMap_, 
00430               RedistributeDomainExporter_, ReducedMatrixDomainMap_);
00431     else {
00432       ReducedMatrixDomainMap_ = OrigReducedMatrixDomainMap_;
00433       OrigReducedMatrixDomainMap_ = 0;
00434       RedistributeDomainExporter_ = 0;
00435     }
00436   
00437     // Create pointer to Full RHS, LHS
00438     Epetra_MultiVector * FullRHS = FullProblem()->GetRHS();
00439     Epetra_MultiVector * FullLHS = FullProblem()->GetLHS();
00440     int NumVectors = FullLHS->NumVectors();
00441   
00442     // Create importers
00443     Full2ReducedLHSImporter_ = new Epetra_Import(*ReducedMatrixDomainMap(), FullMatrixDomainMap());
00444     Full2ReducedRHSImporter_ = new Epetra_Import(*ReducedMatrixRowMap(), FullRHS->Map());
00445 
00446     // Construct Reduced Matrix
00447     ReducedMatrix_ = Teuchos::rcp( new Epetra_CrsMatrix(Copy, *ReducedMatrixRowMap(), *ReducedMatrixColMap(), 0) );
00448 
00449     // Create storage for temporary X values due to explicit elimination of rows
00450     tempExportX_ = new Epetra_MultiVector(FullMatrixColMap(), NumVectors);
00451 
00452     int NumEntries;
00453     double * Values;
00454     int NumMyRows = FullMatrix()->NumMyRows();
00455     int ColSingletonCounter = 0;
00456     for (i=0; i<NumMyRows; i++) {
00457       int_type curGRID = (int_type) FullMatrixRowMap().GID64(i);
00458       if (ReducedMatrixRowMap()->MyGID(curGRID)) { // Check if this row should go into reduced matrix
00459         int_type * Indices;
00460         EPETRA_CHK_ERR(GetRowGCIDs(i, NumEntries, Values, Indices)); // Get current row (Indices are global)
00461       
00462         int ierr = ReducedMatrix()->InsertGlobalValues(curGRID, NumEntries, 
00463                  Values, Indices); // Insert into reduce matrix
00464         // Positive errors will occur because we are submitting col entries that are not part of
00465         // reduced system.  However, because we specified a column map to the ReducedMatrix constructor
00466         // these extra column entries will be ignored and we will be politely reminded by a positive
00467         // error code
00468         if (ierr<0) EPETRA_CHK_ERR(ierr); 
00469       }
00470       else {
00471         int * Indices;
00472         EPETRA_CHK_ERR(GetRow(i, NumEntries, Values, Indices)); // Get current row
00473         if (NumEntries==1) {
00474       double pivot = Values[0];
00475       if (pivot==0.0) EPETRA_CHK_ERR(-1); // Encountered zero row, unable to continue
00476     int indX = Indices[0];
00477     for (j=0; j<NumVectors; j++)
00478       (*tempExportX_)[j][indX] = (*FullRHS)[j][i]/pivot;
00479         }
00480         // Otherwise, this is a singleton column and we will scan for the pivot element needed 
00481         // for post-solve equations
00482         else {
00483     int targetCol = ColSingletonColLIDs_[ColSingletonCounter];
00484     for (j=0; j<NumEntries; j++) {
00485         if (Indices[j]==targetCol) {
00486         double pivot = Values[j];
00487         if (pivot==0.0) EPETRA_CHK_ERR(-2); // Encountered zero column, unable to continue
00488         ColSingletonPivotLIDs_[ColSingletonCounter] = j; // Save for later use
00489         ColSingletonPivots_[ColSingletonCounter] = pivot;
00490         ColSingletonCounter++;
00491         break;
00492       }
00493     }
00494         }
00495       }
00496     }
00497 
00498     // Now convert to local indexing.  We have constructed things so that the domain and range of the
00499     // matrix will have the same map.  If the reduced matrix domain and range maps were not the same, the
00500     // differences were addressed in the ConstructRedistributeExporter() method
00501     EPETRA_CHK_ERR(ReducedMatrix()->FillComplete(*ReducedMatrixDomainMap(), *ReducedMatrixRangeMap()));
00502 
00503     // 1) The vector ColProfiles has column nonzero counts for each processor's contribution
00504     // Construct Reduced LHS (Puts any initial guess values into reduced system)
00505 
00506     ReducedLHS_ = new Epetra_MultiVector(*ReducedMatrixDomainMap(), NumVectors);
00507     EPETRA_CHK_ERR(ReducedLHS_->Import(*FullLHS, *Full2ReducedLHSImporter_, Insert));
00508     FullLHS->PutScalar(0.0); // zero out Full LHS since we will inject values as we get them
00509 
00510     // Construct Reduced RHS
00511 
00512     // First compute influence of already-known values of X on RHS
00513     tempX_ = new Epetra_MultiVector(FullMatrixDomainMap(), NumVectors);
00514     tempB_ = new Epetra_MultiVector(FullRHS->Map(), NumVectors);
00515   
00516     //Inject known X values into tempX for purpose of computing tempB = FullMatrix*tempX
00517     // Also inject into full X since we already know the solution
00518 
00519     if (FullMatrix()->RowMatrixImporter()!=0) {
00520       EPETRA_CHK_ERR(tempX_->Export(*tempExportX_, *FullMatrix()->RowMatrixImporter(), Add));
00521       EPETRA_CHK_ERR(FullLHS->Export(*tempExportX_, *FullMatrix()->RowMatrixImporter(), Add));
00522     }
00523     else {
00524       tempX_->Update(1.0, *tempExportX_, 0.0);
00525       FullLHS->Update(1.0, *tempExportX_, 0.0);
00526     }
00527 
00528     EPETRA_CHK_ERR(FullMatrix()->Multiply(false, *tempX_, *tempB_));
00529 
00530     EPETRA_CHK_ERR(tempB_->Update(1.0, *FullRHS, -1.0)); // tempB now has influence of already-known X values
00531 
00532     ReducedRHS_ = new Epetra_MultiVector(*ReducedMatrixRowMap(), FullRHS->NumVectors());
00533     EPETRA_CHK_ERR(ReducedRHS_->Import(*tempB_, *Full2ReducedRHSImporter_, Insert));
00534 
00535     // Finally construct Reduced Linear Problem
00536     ReducedProblem_ = Teuchos::rcp( new Epetra_LinearProblem(ReducedMatrix_.get(), ReducedLHS_, ReducedRHS_) );
00537   }
00538   else {
00539 
00540     // There are no singletons, so don't bother building a reduced problem. 
00541     ReducedProblem_ = Teuchos::rcp( Problem, false );
00542     ReducedMatrix_ = Teuchos::rcp( dynamic_cast<Epetra_CrsMatrix *>(Problem->GetMatrix()), false );
00543   }
00544 
00545   double fn = (double) FullMatrix()->NumGlobalRows64();
00546   double fnnz = (double) FullMatrix()->NumGlobalNonzeros64();
00547   double rn = (double) ReducedMatrix()->NumGlobalRows64();
00548   double rnnz = (double) ReducedMatrix()->NumGlobalNonzeros64();
00549 
00550   RatioOfDimensions_ = rn/fn;
00551   RatioOfNonzeros_ = rnnz/fnnz;
00552   HaveReducedProblem_ = true;
00553   
00554   return(0);
00555 }
00556 
00557 int LinearProblem_CrsSingletonFilter::ConstructReducedProblem(Epetra_LinearProblem * Problem) {
00558 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00559   if(Problem->GetMatrix()->RowMatrixRowMap().GlobalIndicesInt()) {
00560     return TConstructReducedProblem<int>(Problem);
00561   }
00562   else
00563 #endif
00564 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00565   if(Problem->GetMatrix()->RowMatrixRowMap().GlobalIndicesLongLong()) {
00566     return TConstructReducedProblem<long long>(Problem);
00567   }
00568   else
00569 #endif
00570     throw "LinearProblem_CrsSingletonFilter::ConstructReducedProblem: ERROR, GlobalIndices type unknown.";
00571 }
00572 
00573 //==============================================================================
00574 template<typename int_type>
00575 int LinearProblem_CrsSingletonFilter::TUpdateReducedProblem(Epetra_LinearProblem * Problem) {
00576 
00577   int i, j;
00578 
00579   if (Problem==0) EPETRA_CHK_ERR(-1); // Null problem pointer
00580 
00581   FullProblem_ = Problem;
00582   FullMatrix_ = dynamic_cast<Epetra_RowMatrix *>(Problem->GetMatrix());
00583   if (FullMatrix_==0) EPETRA_CHK_ERR(-2); // Need a RowMatrix
00584   if (Problem->GetRHS()==0) EPETRA_CHK_ERR(-3); // Need a RHS
00585   if (Problem->GetLHS()==0) EPETRA_CHK_ERR(-4); // Need a LHS
00586   if (!HaveReducedProblem_) EPETRA_CHK_ERR(-5); // Must have set up reduced problem
00587 
00588   if ( SingletonsDetected() ) {
00589 
00590     // Create pointer to Full RHS, LHS
00591     Epetra_MultiVector * FullRHS = FullProblem()->GetRHS();
00592     Epetra_MultiVector * FullLHS = FullProblem()->GetLHS();
00593 
00594     int NumVectors = FullLHS->NumVectors();
00595     tempExportX_->PutScalar(0.0);
00596 
00597     int NumEntries;
00598     double * Values;
00599     int NumMyRows = FullMatrix()->NumMyRows();
00600     int ColSingletonCounter = 0;
00601     for (i=0; i<NumMyRows; i++) {
00602       int_type curGRID = (int_type) FullMatrixRowMap().GID64(i);
00603       if (ReducedMatrixRowMap()->MyGID(curGRID)) { // Check if this row should go into reduced matrix
00604         int_type * Indices;
00605         EPETRA_CHK_ERR(GetRowGCIDs(i, NumEntries, Values, Indices)); // Get current row (indices global)
00606         int ierr = ReducedMatrix()->ReplaceGlobalValues(curGRID, NumEntries, 
00607                   Values, Indices);
00608         // Positive errors will occur because we are submitting col entries that are not part of
00609         // reduced system.  However, because we specified a column map to the ReducedMatrix constructor
00610         // these extra column entries will be ignored and we will be politely reminded by a positive
00611         // error code
00612         if (ierr<0) EPETRA_CHK_ERR(ierr); 
00613       }
00614       // Otherwise if singleton row we explicitly eliminate this row and solve for corresponding X value
00615       else {
00616         int * Indices;
00617         EPETRA_CHK_ERR(GetRow(i, NumEntries, Values, Indices)); // Get current row
00618         if (NumEntries==1) {
00619       double pivot = Values[0];
00620     if (pivot==0.0) EPETRA_CHK_ERR(-1); // Encountered zero row, unable to continue
00621     int indX = Indices[0];
00622     for (j=0; j<NumVectors; j++)
00623       (*tempExportX_)[j][indX] = (*FullRHS)[j][i]/pivot;
00624         }
00625         // Otherwise, this is a singleton column and we will scan for the pivot element needed 
00626         // for post-solve equations
00627         else {
00628     j = ColSingletonPivotLIDs_[ColSingletonCounter];
00629     double pivot = Values[j];
00630     if (pivot==0.0) EPETRA_CHK_ERR(-2); // Encountered zero column, unable to continue
00631     ColSingletonPivots_[ColSingletonCounter] = pivot;
00632     ColSingletonCounter++;
00633         }
00634       }
00635     }
00636 
00637     assert(ColSingletonCounter==NumMyColSingletons_); // Sanity test
00638 
00639     // Update Reduced LHS (Puts any initial guess values into reduced system)
00640 
00641     ReducedLHS_->PutScalar(0.0); // zero out Reduced LHS
00642     EPETRA_CHK_ERR(ReducedLHS_->Import(*FullLHS, *Full2ReducedLHSImporter_, Insert));
00643     FullLHS->PutScalar(0.0); // zero out Full LHS since we will inject values as we get them
00644 
00645     // Construct Reduced RHS
00646 
00647     // Zero out temp space
00648     tempX_->PutScalar(0.0);
00649     tempB_->PutScalar(0.0);
00650   
00651     //Inject known X values into tempX for purpose of computing tempB = FullMatrix*tempX
00652     // Also inject into full X since we already know the solution
00653 
00654     if (FullMatrix()->RowMatrixImporter()!=0) {
00655       EPETRA_CHK_ERR(tempX_->Export(*tempExportX_, *FullMatrix()->RowMatrixImporter(), Add));
00656       EPETRA_CHK_ERR(FullLHS->Export(*tempExportX_, *FullMatrix()->RowMatrixImporter(), Add));
00657     }
00658     else {
00659       tempX_->Update(1.0, *tempExportX_, 0.0);
00660       FullLHS->Update(1.0, *tempExportX_, 0.0);
00661     }
00662 
00663     EPETRA_CHK_ERR(FullMatrix()->Multiply(false, *tempX_, *tempB_));
00664 
00665     EPETRA_CHK_ERR(tempB_->Update(1.0, *FullRHS, -1.0)); // tempB now has influence of already-known X values
00666 
00667     ReducedRHS_->PutScalar(0.0);
00668     EPETRA_CHK_ERR(ReducedRHS_->Import(*tempB_, *Full2ReducedRHSImporter_, Insert));
00669   }  
00670   else {
00671 
00672     // There are no singletons, so don't bother building a reduced problem. 
00673     ReducedProblem_ = Teuchos::rcp( Problem, false );
00674     ReducedMatrix_ = Teuchos::rcp( dynamic_cast<Epetra_CrsMatrix *>(Problem->GetMatrix()), false );
00675   }
00676 
00677   return(0);
00678 }
00679 
00680 int LinearProblem_CrsSingletonFilter::UpdateReducedProblem(Epetra_LinearProblem * Problem) {
00681 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00682   if(Problem->GetMatrix()->RowMatrixRowMap().GlobalIndicesInt()) {
00683     return TUpdateReducedProblem<int>(Problem);
00684   }
00685   else
00686 #endif
00687 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00688   if(Problem->GetMatrix()->RowMatrixRowMap().GlobalIndicesLongLong()) {
00689     return TUpdateReducedProblem<long long>(Problem);
00690   }
00691   else
00692 #endif
00693     throw "LinearProblem_CrsSingletonFilter::UpdateReducedProblem: ERROR, GlobalIndices type unknown.";
00694 }
00695 //==============================================================================
00696 template<typename int_type>
00697 int LinearProblem_CrsSingletonFilter::TConstructRedistributeExporter(Epetra_Map * SourceMap, Epetra_Map * TargetMap,
00698                    Epetra_Export * & RedistributeExporter,
00699                    Epetra_Map * & RedistributeMap) {
00700 
00701   int_type IndexBase = (int_type) SourceMap->IndexBase64();
00702   if (IndexBase!=(int_type) TargetMap->IndexBase64()) EPETRA_CHK_ERR(-1);
00703 
00704   const Epetra_Comm & Comm = TargetMap->Comm();
00705 
00706   int TargetNumMyElements = TargetMap->NumMyElements();
00707   int SourceNumMyElements = SourceMap->NumMyElements();
00708 
00709   // ContiguousTargetMap has same number of elements per PE as TargetMap, but uses contigious indexing 
00710   Epetra_Map ContiguousTargetMap((int_type) -1, TargetNumMyElements, IndexBase,Comm);
00711 
00712   // Same for ContiguousSourceMap
00713   Epetra_Map ContiguousSourceMap((int_type) -1, SourceNumMyElements, IndexBase, Comm);
00714 
00715   assert(ContiguousSourceMap.NumGlobalElements64()==ContiguousTargetMap.NumGlobalElements64());
00716 
00717   // Now create a vector that contains the global indices of the Source Epetra_MultiVector
00718   int_type* SourceMapMyGlobalElements = 0;
00719   SourceMap->MyGlobalElementsPtr(SourceMapMyGlobalElements);
00720   typename Epetra_GIDTypeVector<int_type>::impl SourceIndices(View, ContiguousSourceMap, SourceMapMyGlobalElements);
00721 
00722   // Create an exporter to send the SourceMap global IDs to the target distribution
00723   Epetra_Export Exporter(ContiguousSourceMap, ContiguousTargetMap);
00724   
00725   // Create a vector to catch the global IDs in the target distribution
00726   typename Epetra_GIDTypeVector<int_type>::impl TargetIndices(ContiguousTargetMap);
00727   TargetIndices.Export(SourceIndices, Exporter, Insert);
00728 
00729   // Create a new map that describes how the Source MultiVector should be laid out so that it has
00730   // the same number of elements on each processor as the TargetMap
00731   RedistributeMap = new Epetra_Map((int_type) -1, TargetNumMyElements, TargetIndices.Values(), IndexBase, Comm);
00732 
00733   // This exporter will finally redistribute the Source MultiVector to the same layout as the TargetMap
00734   RedistributeExporter = new Epetra_Export(*SourceMap, *RedistributeMap);
00735   return(0);
00736 }
00737 
00738 int LinearProblem_CrsSingletonFilter::ConstructRedistributeExporter(Epetra_Map * SourceMap, Epetra_Map * TargetMap,
00739                    Epetra_Export * & RedistributeExporter,
00740                    Epetra_Map * & RedistributeMap) {
00741 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00742   if(SourceMap->GlobalIndicesInt() && TargetMap->GlobalIndicesInt()) {
00743     return TConstructRedistributeExporter<int>(SourceMap, TargetMap, RedistributeExporter, RedistributeMap);
00744   }
00745   else
00746 #endif
00747 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00748   if(SourceMap->GlobalIndicesLongLong() && TargetMap->GlobalIndicesLongLong()) {
00749     return TConstructRedistributeExporter<long long>(SourceMap, TargetMap, RedistributeExporter, RedistributeMap);
00750   }
00751   else
00752 #endif
00753     throw "LinearProblem_CrsSingletonFilter::ConstructRedistributeExporter: ERROR, GlobalIndices type unknown.";
00754 }
00755 //==============================================================================
00756 int LinearProblem_CrsSingletonFilter::ComputeFullSolution() {
00757 
00758   if ( SingletonsDetected() ) {
00759     int jj, k;
00760 
00761     Epetra_MultiVector * FullLHS = FullProblem()->GetLHS(); 
00762     Epetra_MultiVector * FullRHS = FullProblem()->GetRHS(); 
00763 
00764     tempX_->PutScalar(0.0); tempExportX_->PutScalar(0.0);
00765     // Inject values that the user computed for the reduced problem into the full solution vector
00766     EPETRA_CHK_ERR(tempX_->Export(*ReducedLHS_, *Full2ReducedLHSImporter_, Add));
00767 
00768     FullLHS->Update(1.0, *tempX_, 1.0);
00769 
00770     // Next we will use our full solution vector which is populated with pre-filter solution
00771     // values and reduced system solution values to compute the sum of the row contributions
00772     // that must be subtracted to get the post-filter solution values
00773 
00774     EPETRA_CHK_ERR(FullMatrix()->Multiply(false, *FullLHS, *tempB_));
00775 
00776     // Finally we loop through the local rows that were associated with column singletons and compute the
00777     // solution for these equations.
00778 
00779     int NumVectors = tempB_->NumVectors();
00780     for (k=0; k<NumMyColSingletons_; k++) {
00781       int i = ColSingletonRowLIDs_[k];
00782       int j = ColSingletonColLIDs_[k];
00783       double pivot = ColSingletonPivots_[k];
00784       for (jj=0; jj<NumVectors; jj++)
00785         (*tempExportX_)[jj][j]= ((*FullRHS)[jj][i] - (*tempB_)[jj][i])/pivot;
00786     }
00787 
00788     // Finally, insert values from post-solve step and we are done!!!!
00789   
00790     if (FullMatrix()->RowMatrixImporter()!=0) {
00791       EPETRA_CHK_ERR(tempX_->Export(*tempExportX_, *FullMatrix()->RowMatrixImporter(), Add));
00792     }
00793     else {
00794       tempX_->Update(1.0, *tempExportX_, 0.0);
00795     }
00796 
00797     FullLHS->Update(1.0, *tempX_, 1.0);
00798   }
00799      
00800   return(0);
00801 }
00802 //==============================================================================
00803 int LinearProblem_CrsSingletonFilter::InitFullMatrixAccess() {
00804 
00805   MaxNumMyEntries_ = FullMatrix()->MaxNumEntries(); 
00806 
00807   // Cast to CrsMatrix, if possible.  Can save some work.
00808   FullCrsMatrix_ = dynamic_cast<Epetra_CrsMatrix *>(FullMatrix());
00809   FullMatrixIsCrsMatrix_ = (FullCrsMatrix_!=0); // Pointer is non-zero if cast worked
00810   Indices_int_ = new int[MaxNumMyEntries_];
00811 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00812   if(FullMatrix()->RowMatrixRowMap().GlobalIndicesLongLong())
00813     Indices_LL_ = new long long[MaxNumMyEntries_];
00814 #endif
00815   Values_.Size(MaxNumMyEntries_);
00816 
00817   return(0);
00818 }
00819 //==============================================================================
00820 int LinearProblem_CrsSingletonFilter::GetRow(int Row, int & NumIndices, int * & Indices) {
00821 
00822   if (FullMatrixIsCrsMatrix_) { // View of current row
00823     EPETRA_CHK_ERR(FullCrsMatrix()->Graph().ExtractMyRowView(Row, NumIndices, Indices)); 
00824   }
00825   else { // Copy of current row (we must get the values, but we ignore them)
00826     EPETRA_CHK_ERR(FullMatrix()->ExtractMyRowCopy(Row, MaxNumMyEntries_, NumIndices, 
00827               Values_.Values(), Indices_int_));
00828     Indices = Indices_int_;
00829   } 
00830   return(0);
00831 }
00832 //==============================================================================
00833 int LinearProblem_CrsSingletonFilter::GetRow(int Row, int & NumIndices, 
00834               double * & Values, int * & Indices) {
00835 
00836   if (FullMatrixIsCrsMatrix_) { // View of current row
00837     EPETRA_CHK_ERR(FullCrsMatrix_->ExtractMyRowView(Row, NumIndices, Values, Indices)); 
00838   }
00839   else { // Copy of current row (we must get the values, but we ignore them)
00840     EPETRA_CHK_ERR(FullMatrix()->ExtractMyRowCopy(Row, MaxNumMyEntries_, NumIndices, 
00841               Values_.Values(), Indices_int_));
00842     Values = Values_.Values();
00843     Indices = Indices_int_;
00844   } 
00845   return(0);
00846 }
00847 //==============================================================================
00848 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00849 int LinearProblem_CrsSingletonFilter::GetRowGCIDs(int Row, int & NumIndices, 
00850              double * & Values, int * & GlobalIndices) {
00851 
00852     EPETRA_CHK_ERR(FullMatrix()->ExtractMyRowCopy(Row, MaxNumMyEntries_, NumIndices, 
00853               Values_.Values(), Indices_int_));
00854     for (int j=0; j<NumIndices; j++) Indices_int_[j] = FullMatrixColMap().GID(Indices_int_[j]);
00855     Values = Values_.Values();
00856     GlobalIndices = Indices_int_;
00857   return(0);
00858 }
00859 #endif
00860 
00861 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00862 int LinearProblem_CrsSingletonFilter::GetRowGCIDs(int Row, int & NumIndices, 
00863              double * & Values, long long * & GlobalIndices) {
00864     EPETRA_CHK_ERR(FullMatrix()->ExtractMyRowCopy(Row, MaxNumMyEntries_, NumIndices, 
00865               Values_.Values(), Indices_int_));
00866     for (int j=0; j<NumIndices; j++) Indices_LL_[j] = FullMatrixColMap().GID64(Indices_int_[j]);
00867     Values = Values_.Values();
00868     GlobalIndices = Indices_LL_;
00869   return(0);
00870 }
00871 #endif
00872 //==============================================================================
00873 int LinearProblem_CrsSingletonFilter::CreatePostSolveArrays(const Epetra_IntVector & RowIDs,
00874                  const Epetra_MapColoring & RowMapColors,
00875                  const Epetra_IntVector & ColProfiles,
00876                  const Epetra_IntVector & NewColProfiles,
00877                  const Epetra_IntVector & ColHasRowWithSingleton) {
00878 
00879   int j;
00880 
00881   if (NumMyColSingletons_==0) return(0); // Nothing to do
00882 
00883   Epetra_MapColoring & ColMapColors = *ColMapColors_;
00884 
00885   int NumMyCols = FullMatrix()->NumMyCols();
00886 
00887   // We will need these arrays for the post-solve phase
00888   ColSingletonRowLIDs_ = new int[NumMyColSingletons_];
00889   ColSingletonColLIDs_ = new int[NumMyColSingletons_];
00890   ColSingletonPivotLIDs_ = new int[NumMyColSingletons_];
00891   ColSingletonPivots_ = new double[NumMyColSingletons_];
00892   
00893   // Register singleton columns (that were not already counted as singleton rows)
00894   // Check to see if any columns disappeared because all associated rows were eliminated
00895   int NumMyColSingletonstmp = 0;
00896   for (j=0; j<NumMyCols; j++) {
00897     int i = RowIDs[j];
00898     if ( ColProfiles[j]==1 && RowMapColors[i]!=1 ) {
00899       ColSingletonRowLIDs_[NumMyColSingletonstmp] = i;
00900       ColSingletonColLIDs_[NumMyColSingletonstmp] = j;
00901       NumMyColSingletonstmp++;
00902     }
00903     // Also check for columns that were eliminated implicitly by 
00904     // having all associated row eliminated
00905     else if (NewColProfiles[j]==0 && ColHasRowWithSingleton[j]!=1 && RowMapColors[i]==0) {
00906     ColMapColors[j] = 1;
00907     }
00908   }
00909 
00910   assert(NumMyColSingletonstmp==NumMyColSingletons_); //Sanity check
00911   Epetra_Util sorter;
00912   sorter.Sort(true, NumMyColSingletons_, ColSingletonRowLIDs_, 0, 0, 1, &ColSingletonColLIDs_);
00913     
00914   return(0);
00915 }
00916 
00917 } //namespace EpetraExt
00918 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines