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