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