EpetraExt Development
EpetraExt_StaticCondensation_LinearProblem.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 
00042 #include <EpetraExt_StaticCondensation_LinearProblem.h>
00043 
00044 #include <Epetra_Export.h>
00045 #include <Epetra_LinearProblem.h>
00046 #include <Epetra_CrsGraph.h>
00047 #include <Epetra_CrsMatrix.h>
00048 #include <Epetra_MultiVector.h>
00049 #include <Epetra_Vector.h>
00050 #include <Epetra_IntVector.h>
00051 #include <Epetra_Map.h>
00052 #include <Epetra_Comm.h>
00053 
00054 #include <vector>
00055 #include <map>
00056 #include <set>
00057 
00058 namespace EpetraExt {
00059 
00060 LinearProblem_StaticCondensation::
00061 ~LinearProblem_StaticCondensation()
00062 {
00063   if( Exporter_ ) delete Exporter_;
00064 
00065   if( NewProblem_ ) delete NewProblem_;
00066   if( NewRHS_ ) delete NewRHS_;
00067   if( NewLHS_ ) delete NewLHS_;
00068   if( NewMatrix_ ) delete NewMatrix_;
00069   if( NewGraph_ ) delete NewGraph_;
00070   if( NewRowMap_ ) delete NewRowMap_;
00071   if( NewColMap_ ) delete NewColMap_;
00072 
00073   if( ULHS_ ) delete ULHS_;
00074   if( RLHS_ ) delete RLHS_;
00075   if( LLHS_ ) delete LLHS_;
00076   if( URHS_ ) delete URHS_;
00077   if( RRHS_ ) delete RRHS_;
00078   if( LRHS_ ) delete LRHS_;
00079 
00080   if( UUMatrix_ ) delete UUMatrix_;
00081   if( URMatrix_ ) delete URMatrix_;
00082   if( ULMatrix_ ) delete ULMatrix_;
00083   if( RRMatrix_ ) delete RRMatrix_;
00084   if( RLMatrix_ ) delete RLMatrix_;
00085   if( LLMatrix_ ) delete LLMatrix_;
00086 
00087   if( UUGraph_ ) delete UUGraph_;
00088   if( URGraph_ ) delete URGraph_;
00089   if( ULGraph_ ) delete ULGraph_;
00090   if( RRGraph_ ) delete RRGraph_;
00091   if( RLGraph_ ) delete RLGraph_;
00092   if( LLGraph_ ) delete LLGraph_;
00093 
00094   if( UExporter_ ) delete UExporter_;
00095   if( RExporter_ ) delete RExporter_;
00096   if( LExporter_ ) delete LExporter_;
00097 
00098   if( UMap_ ) delete UMap_;
00099   if( RMap_ ) delete RMap_;
00100   if( LMap_ ) delete LMap_;
00101 }
00102 
00103 LinearProblem_StaticCondensation::NewTypeRef
00104 LinearProblem_StaticCondensation::
00105 operator()( OriginalTypeRef orig )
00106 {
00107   origObj_ = &orig;
00108 
00109   int ierr = 0;
00110 
00111   OldMatrix_ = dynamic_cast<Epetra_CrsMatrix*>( orig.GetMatrix() );
00112   OldGraph_  = &OldMatrix_->Graph();
00113   OldRHS_    = orig.GetRHS();
00114   OldLHS_    = orig.GetLHS();
00115   OldRowMap_ = &OldMatrix_->RowMap();
00116 
00117   const Epetra_Comm & CommObj = OldRowMap_->Comm();
00118 
00119   if( !OldMatrix_ ) ierr = -2;
00120   if( !OldRHS_ )    ierr = -3;
00121   if( !OldLHS_ )    ierr = -4;
00122 
00123   if( OldRowMap_->DistributedGlobal() ) ierr = -5;
00124   if( degree_ != 1 ) ierr = -6;
00125 
00126   int NRows = OldGraph_->NumMyRows();
00127   int IndexBase = OldRowMap_->IndexBase();
00128 
00129   vector<int> ColNZCnt( NRows );
00130   vector<int> CS_RowIndices( NRows );
00131   map<int,int> RS_Map;
00132   map<int,int> CS_Map;
00133 
00134   int NumIndices;
00135   int * Indices;
00136   for( int i = 0; i < NRows; ++i )
00137   {
00138     ierr = OldGraph_->ExtractMyRowView( i, NumIndices, Indices );
00139 
00140     for( int j = 0; j < NumIndices; ++j )
00141     {
00142       ++ColNZCnt[ Indices[j] ];
00143       CS_RowIndices[ Indices[j] ] = i;
00144     }
00145 
00146     if( NumIndices == 1 ) RS_Map[i] = Indices[0];
00147   }
00148 
00149   if( verbose_ )
00150   {
00151     cout << "-------------------------\n";
00152     cout << "Row Singletons\n";
00153     for( map<int,int>::iterator itM = RS_Map.begin(); itM != RS_Map.end(); ++itM )
00154       cout << (*itM).first << "\t" << (*itM).second << endl;
00155     cout << "Col Counts\n";
00156     for( int i = 0; i < NRows; ++i )
00157       cout << i << "\t" << ColNZCnt[i] << "\t" << CS_RowIndices[i] << endl;
00158     cout << "-------------------------\n";
00159   }
00160 
00161   set<int> RS_Set;
00162   set<int> CS_Set;
00163 
00164   for( int i = 0; i < NRows; ++i )
00165     if( ColNZCnt[i] == 1 )
00166     {
00167       int RowIndex = CS_RowIndices[i];
00168       if( RS_Map.count(i) && RS_Map[i] == RowIndex )
00169       {
00170         CS_Set.insert(i);
00171         RS_Set.insert( RowIndex );
00172       }
00173     }
00174 
00175   if( verbose_ )
00176   {
00177     cout << "-------------------------\n";
00178     cout << "Singletons: " << CS_Set.size() << endl;
00179     set<int>::iterator itRS = RS_Set.begin();
00180     set<int>::iterator itCS = CS_Set.begin();
00181     set<int>::iterator endRS = RS_Set.end();
00182     set<int>::iterator endCS = CS_Set.end();
00183     for( ; itRS != endRS; ++itRS, ++itCS )
00184       cout << *itRS << "\t" << *itCS << endl;
00185     cout << "-------------------------\n";
00186   }
00187 
00188   //build new maps
00189   int NSingletons = CS_Set.size();
00190   int NReducedRows = NRows - 2* NSingletons; 
00191   vector<int> ReducedIndices( NReducedRows );
00192   vector<int> CSIndices( NSingletons );
00193   vector<int> RSIndices( NSingletons );
00194   int Reduced_Loc = 0;
00195   int RS_Loc = 0;
00196   int CS_Loc = 0;
00197   for( int i = 0; i < NRows; ++i )
00198   {
00199     int GlobalIndex = OldRowMap_->GID(i);
00200     if     ( RS_Set.count(i) ) RSIndices[RS_Loc++] = GlobalIndex;
00201     else if( CS_Set.count(i) ) CSIndices[CS_Loc++] = GlobalIndex;
00202     else                       ReducedIndices[Reduced_Loc++] = GlobalIndex;
00203   }
00204 
00205   vector<int> RowPermutedIndices( NRows );
00206   copy( RSIndices.begin(), RSIndices.end(), RowPermutedIndices.begin() );
00207   copy( ReducedIndices.begin(), ReducedIndices.end(), RowPermutedIndices.begin() + NSingletons );
00208   copy( CSIndices.begin(), CSIndices.end(), RowPermutedIndices.begin() + NReducedRows + NSingletons );
00209 
00210   vector<int> ColPermutedIndices( NRows );
00211   copy( CSIndices.begin(), CSIndices.end(), ColPermutedIndices.begin() );
00212   copy( ReducedIndices.begin(), ReducedIndices.end(), ColPermutedIndices.begin() + NSingletons );
00213   copy( RSIndices.begin(), RSIndices.end(), ColPermutedIndices.begin() + NReducedRows + NSingletons );
00214 
00215   UMap_ = new Epetra_Map( NSingletons, NSingletons, &RSIndices[0], IndexBase, CommObj );
00216   RMap_ = new Epetra_Map( NReducedRows, NReducedRows, &ReducedIndices[0], IndexBase, CommObj );
00217   LMap_ = new Epetra_Map( NSingletons, NSingletons, &CSIndices[0], IndexBase, CommObj );
00218 
00219   NewRowMap_ = new Epetra_Map( NRows, NRows, &RowPermutedIndices[0], IndexBase, CommObj );
00220   NewColMap_ = new Epetra_Map( NRows, NRows, &ColPermutedIndices[0], IndexBase, CommObj );
00221 
00222   //Construct Permuted System
00223   Exporter_ = new Epetra_Export( *OldRowMap_, *NewRowMap_ );
00224 
00225   NewRHS_ = new Epetra_MultiVector( *NewRowMap_, OldRHS_->NumVectors() );
00226   NewRHS_->Export( *OldRHS_, *Exporter_, Insert );
00227 
00228   NewLHS_ = new Epetra_MultiVector( *NewRowMap_, OldLHS_->NumVectors() );
00229   NewLHS_->Export( *OldLHS_, *Exporter_, Insert );
00230 
00231   NewGraph_ = new Epetra_CrsGraph( Copy, *NewRowMap_, *NewColMap_, 0 );
00232   NewGraph_->Export( *OldGraph_, *Exporter_, Insert );
00233   NewGraph_->FillComplete();
00234   cout << "Permuted Graph:\n" << *NewGraph_;
00235 
00236   NewMatrix_ = new Epetra_CrsMatrix( Copy, *NewGraph_ );
00237   NewMatrix_->Export( *OldMatrix_, *Exporter_, Insert );
00238   NewMatrix_->FillComplete();
00239   cout << "Permuted Matrix:\n" << *NewMatrix_;
00240 
00241   UExporter_ = new Epetra_Export( *OldRowMap_, *UMap_ );
00242   RExporter_ = new Epetra_Export( *OldRowMap_, *RMap_ );
00243   LExporter_ = new Epetra_Export( *OldRowMap_, *LMap_ );
00244 cout << "UExporter:\n" << *UExporter_;
00245 cout << "RExporter:\n" << *RExporter_;
00246 cout << "LExporter:\n" << *LExporter_;
00247 
00248   ULHS_ = new Epetra_MultiVector( *LMap_, OldLHS_->NumVectors() );
00249   ULHS_->Export( *OldLHS_, *LExporter_, Insert );
00250   cout << "ULHS:\n" << *ULHS_;
00251 
00252   RLHS_ = new Epetra_MultiVector( *RMap_, OldLHS_->NumVectors() );
00253   RLHS_->Export( *OldLHS_, *RExporter_, Insert );
00254   cout << "RLHS:\n" << *RLHS_;
00255 
00256   LLHS_ = new Epetra_MultiVector( *UMap_, OldLHS_->NumVectors() );
00257   LLHS_->Export( *OldLHS_, *UExporter_, Insert );
00258   cout << "LLHS:\n" << *LLHS_;
00259 
00260   URHS_ = new Epetra_MultiVector( *UMap_, OldRHS_->NumVectors() );
00261   URHS_->Export( *OldRHS_, *UExporter_, Insert );
00262   cout << "URHS:\n" << *URHS_;
00263 
00264   RRHS_ = new Epetra_MultiVector( *RMap_, OldRHS_->NumVectors() );
00265   RRHS_->Export( *OldRHS_, *RExporter_, Insert );
00266   cout << "RRHS:\n" << *RRHS_;
00267 
00268   LRHS_ = new Epetra_MultiVector( *LMap_, OldRHS_->NumVectors() );
00269   LRHS_->Export( *OldRHS_, *LExporter_, Insert );
00270   cout << "LRHS:\n" << *LRHS_;
00271 
00272   UUGraph_ = new Epetra_CrsGraph( Copy, *UMap_, *LMap_, 0 );
00273   UUGraph_->Export( *OldGraph_, *UExporter_, Insert );
00274   UUGraph_->FillComplete( LMap_, UMap_ );
00275   cout << "UUGraph:\n" << *UUGraph_;
00276 
00277   UUMatrix_ = new Epetra_CrsMatrix( Copy, *UUGraph_ );
00278   UUMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
00279   UUMatrix_->FillComplete();
00280   cout << "UUMatrix:\n" << *UUMatrix_;
00281 
00282   URGraph_ = new Epetra_CrsGraph( Copy, *UMap_, *RMap_, 0 );
00283   URGraph_->Export( *OldGraph_, *UExporter_, Insert );
00284   URGraph_->FillComplete( RMap_, UMap_ );
00285   cout << "URGraph:\n" << *URGraph_;
00286 
00287   URMatrix_ = new Epetra_CrsMatrix( Copy, *URGraph_ );
00288   URMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
00289   URMatrix_->FillComplete();
00290   cout << "URMatrix:\n" << *URMatrix_;
00291 
00292   ULGraph_ = new Epetra_CrsGraph( Copy, *UMap_, *UMap_, 0 );
00293   ULGraph_->Export( *OldGraph_, *UExporter_, Insert );
00294   ULGraph_->FillComplete();
00295   cout << "ULGraph:\n" << *ULGraph_;
00296 
00297   ULMatrix_ = new Epetra_CrsMatrix( Copy, *ULGraph_ );
00298   ULMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
00299   ULMatrix_->FillComplete();
00300   cout << "ULMatrix:\n" << *ULMatrix_;
00301 
00302   RRGraph_ = new Epetra_CrsGraph( Copy, *RMap_, *RMap_, 0 );
00303   RRGraph_->Export( *OldGraph_, *RExporter_, Insert );
00304   RRGraph_->FillComplete();
00305   cout << "RRGraph:\n" << *RRGraph_;
00306 
00307   RRMatrix_ = new Epetra_CrsMatrix( Copy, *RRGraph_ );
00308   RRMatrix_->Export( *OldMatrix_, *RExporter_, Insert );
00309   RRMatrix_->FillComplete();
00310   cout << "RRMatrix:\n" << *RRMatrix_;
00311 
00312   RLGraph_ = new Epetra_CrsGraph( Copy, *RMap_, *UMap_, 0 );
00313   RLGraph_->Export( *OldGraph_, *RExporter_, Insert );
00314   RLGraph_->FillComplete( UMap_, RMap_ );
00315   cout << "RLGraph:\n" << *RLGraph_;
00316 
00317   RLMatrix_ = new Epetra_CrsMatrix( Copy, *RLGraph_ );
00318   RLMatrix_->Export( *OldMatrix_, *RExporter_, Insert );
00319   RLMatrix_->FillComplete();
00320   cout << "RLMatrix:\n" << *RLMatrix_;
00321 
00322   LLGraph_ = new Epetra_CrsGraph( Copy, *LMap_, *UMap_, 0 );
00323   LLGraph_->Export( *OldGraph_, *LExporter_, Insert );
00324   LLGraph_->FillComplete( UMap_, LMap_ );
00325   cout << "LLGraph:\n" << *LLGraph_;
00326 
00327   LLMatrix_ = new Epetra_CrsMatrix( Copy, *LLGraph_ );
00328   LLMatrix_->Export( *OldMatrix_, *LExporter_, Insert );
00329   LLMatrix_->FillComplete();
00330   cout << "LLMatrix:\n" << *LLMatrix_;
00331 
00332   if( verbose_ )
00333   {
00334     cout << "Full System Characteristics:" << endl
00335          << "----------------------------" << endl
00336          << " Dimension                   = " << NRows << endl
00337          << " NNZs                        = " << OldMatrix_->NumGlobalNonzeros() << endl
00338          << " Max Num Row Entries         = " << OldMatrix_->GlobalMaxNumEntries() << endl << endl
00339          << "Reduced System Characteristics:" << endl
00340          << " Dimension                   = " << NReducedRows << endl
00341          << " NNZs                        = " << RRMatrix_->NumGlobalNonzeros() << endl
00342          << " Max Num Row Entries         = " << RRMatrix_->GlobalMaxNumEntries() << endl
00343          << "Singleton Info:" << endl
00344          << " Num Singleton                 = " << NSingletons << endl
00345          << "Ratios:" << endl
00346          << " % Reduction in Dimension    = " << 100.0*(NRows-NReducedRows)/NRows << endl
00347          << " % Reduction in NNZs         = " << (OldMatrix_->NumGlobalNonzeros()-RRMatrix_->NumGlobalNonzeros())/100.0 << endl
00348          << "-------------------------------" << endl;
00349   }
00350 
00351   NewProblem_ = new Epetra_LinearProblem( RRMatrix_, RLHS_, RRHS_ );
00352 
00353   newObj_ = NewProblem_;
00354 
00355   cout << "done with SC\n";
00356 
00357   return *NewProblem_;
00358 }
00359 
00360 bool
00361 LinearProblem_StaticCondensation::
00362 fwd()
00363 {
00364   if( verbose_ ) cout << "LP_SC::fwd()\n";
00365   if( verbose_ ) cout << "LP_SC::fwd() : Exporting to New LHS\n";
00366   ULHS_->Export( *OldLHS_, *LExporter_, Insert );
00367   RLHS_->Export( *OldLHS_, *RExporter_, Insert );
00368   LLHS_->Export( *OldLHS_, *UExporter_, Insert );
00369 
00370   if( verbose_ ) cout << "LP_SC::fwd() : Exporting to New RHS\n";
00371   URHS_->Export( *OldRHS_, *UExporter_, Insert );
00372   RRHS_->Export( *OldRHS_, *RExporter_, Insert );
00373   LRHS_->Export( *OldRHS_, *LExporter_, Insert );
00374 
00375   UUMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
00376   URMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
00377   ULMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
00378   RRMatrix_->Export( *OldMatrix_, *RExporter_, Insert );
00379   RLMatrix_->Export( *OldMatrix_, *RExporter_, Insert );
00380   LLMatrix_->Export( *OldMatrix_, *LExporter_, Insert );
00381 
00382   Epetra_Vector LLDiag( *LMap_ );
00383   LLMatrix_->ExtractDiagonalCopy( LLDiag );
00384   Epetra_Vector LLRecipDiag( *LMap_ );
00385   LLRecipDiag.Reciprocal( LLDiag );
00386 
00387   if( verbose_ ) cout << "LP_SC::fwd() : Forming LLHS\n";
00388   LLDiag.Multiply( 1.0, LLRecipDiag, *LRHS_, 0.0 );
00389   int LSize = LMap_->NumMyElements();
00390   for( int i = 0; i < LSize; ++i ) (*LLHS_)[0][i] = LLDiag[i];
00391 
00392   if( verbose_ ) cout << "LP_SC::fwd() : Updating RRHS\n";
00393   Epetra_Vector RUpdate( *RMap_ );
00394   RLMatrix_->Multiply( false, *LLHS_, RUpdate );
00395   RRHS_->Update( -1.0, RUpdate, 1.0 );
00396 
00397   if( verbose_ ) cout << "LP_SC::fwd() : Updating URHS\n";
00398   Epetra_Vector UUpdate( *UMap_ );
00399   ULMatrix_->Multiply( false, *LLHS_, UUpdate );
00400   URHS_->Update( -1.0, UUpdate, 1.0 );
00401   
00402   return true;
00403 }
00404 
00405 bool
00406 LinearProblem_StaticCondensation::
00407 rvs()
00408 {
00409   if( verbose_ ) cout << "LP_SC::rvs()\n";
00410   if( verbose_ ) cout << "LP_SC::rvs() : Updating URHS\n";
00411   Epetra_Vector UUpdate( *UMap_ );
00412   URMatrix_->Multiply( false, *RLHS_, UUpdate );
00413   URHS_->Update( -1.0, UUpdate, 1.0 );
00414 
00415   Epetra_Vector UUDiag( *UMap_ );
00416   UUMatrix_->ExtractDiagonalCopy( UUDiag );
00417   Epetra_Vector UURecipDiag( *UMap_ );
00418   UURecipDiag.Reciprocal( UUDiag );
00419 
00420   if( verbose_ ) cout << "LP_SC::rvs() : Forming ULHS\n";
00421   UUDiag.Multiply( 1.0, UURecipDiag, *URHS_, 0.0 );
00422   int USize = UMap_->NumMyElements();
00423   for( int i = 0; i < USize; ++i ) (*ULHS_)[0][i] = UUDiag[i];
00424 
00425   if( verbose_ ) cout << "LP_SC::rvs() : Importing to Old LHS\n";
00426   OldLHS_->Import( *ULHS_, *LExporter_, Insert );
00427   OldLHS_->Import( *RLHS_, *RExporter_, Insert );
00428   OldLHS_->Import( *LLHS_, *UExporter_, Insert );
00429 
00430   return true;
00431 }
00432 
00433 } // namespace EpetraExt
00434 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines