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

Generated on Tue Oct 20 12:45:29 2009 for EpetraExt by doxygen 1.4.7