EpetraExt_CrsSingletonFilter_LinearProblem.cpp

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

Generated on Tue Jul 13 09:23:05 2010 for EpetraExt by  doxygen 1.4.7