Epetra_CrsSingletonFilter.cpp

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

Generated on Wed May 12 21:41:05 2010 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.4.7