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