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

Generated on Thu Sep 18 12:37:57 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1