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 (2001) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include <EpetraExt_StaticCondensation_LinearProblem.h>
00030 
00031 #include <Epetra_Export.h>
00032 #include <Epetra_LinearProblem.h>
00033 #include <Epetra_CrsGraph.h>
00034 #include <Epetra_CrsMatrix.h>
00035 #include <Epetra_MultiVector.h>
00036 #include <Epetra_Vector.h>
00037 #include <Epetra_IntVector.h>
00038 #include <Epetra_Map.h>
00039 #include <Epetra_Comm.h>
00040 
00041 #include <vector>
00042 #include <map>
00043 #include <set>
00044 
00045 namespace EpetraExt {
00046 
00047 LinearProblem_StaticCondensation::
00048 ~LinearProblem_StaticCondensation()
00049 {
00050   if( Exporter_ ) delete Exporter_;
00051 
00052   if( NewProblem_ ) delete NewProblem_;
00053   if( NewRHS_ ) delete NewRHS_;
00054   if( NewLHS_ ) delete NewLHS_;
00055   if( NewMatrix_ ) delete NewMatrix_;
00056   if( NewGraph_ ) delete NewGraph_;
00057   if( NewRowMap_ ) delete NewRowMap_;
00058   if( NewColMap_ ) delete NewColMap_;
00059 
00060   if( ULHS_ ) delete ULHS_;
00061   if( RLHS_ ) delete RLHS_;
00062   if( LLHS_ ) delete LLHS_;
00063   if( URHS_ ) delete URHS_;
00064   if( RRHS_ ) delete RRHS_;
00065   if( LRHS_ ) delete LRHS_;
00066 
00067   if( UUMatrix_ ) delete UUMatrix_;
00068   if( URMatrix_ ) delete URMatrix_;
00069   if( ULMatrix_ ) delete ULMatrix_;
00070   if( RRMatrix_ ) delete RRMatrix_;
00071   if( RLMatrix_ ) delete RLMatrix_;
00072   if( LLMatrix_ ) delete LLMatrix_;
00073 
00074   if( UUGraph_ ) delete UUGraph_;
00075   if( URGraph_ ) delete URGraph_;
00076   if( ULGraph_ ) delete ULGraph_;
00077   if( RRGraph_ ) delete RRGraph_;
00078   if( RLGraph_ ) delete RLGraph_;
00079   if( LLGraph_ ) delete LLGraph_;
00080 
00081   if( UExporter_ ) delete UExporter_;
00082   if( RExporter_ ) delete RExporter_;
00083   if( LExporter_ ) delete LExporter_;
00084 
00085   if( UMap_ ) delete UMap_;
00086   if( RMap_ ) delete RMap_;
00087   if( LMap_ ) delete LMap_;
00088 }
00089 
00090 LinearProblem_StaticCondensation::NewTypeRef
00091 LinearProblem_StaticCondensation::
00092 operator()( OriginalTypeRef orig )
00093 {
00094   origObj_ = &orig;
00095 
00096   int ierr = 0;
00097 
00098   OldMatrix_ = dynamic_cast<Epetra_CrsMatrix*>( orig.GetMatrix() );
00099   OldGraph_  = &OldMatrix_->Graph();
00100   OldRHS_    = orig.GetRHS();
00101   OldLHS_    = orig.GetLHS();
00102   OldRowMap_ = &OldMatrix_->RowMap();
00103 
00104   const Epetra_Comm & CommObj = OldRowMap_->Comm();
00105 
00106   if( !OldMatrix_ ) ierr = -2;
00107   if( !OldRHS_ )    ierr = -3;
00108   if( !OldLHS_ )    ierr = -4;
00109 
00110   if( OldRowMap_->DistributedGlobal() ) ierr = -5;
00111   if( degree_ != 1 ) ierr = -6;
00112 
00113   int NRows = OldGraph_->NumMyRows();
00114   int IndexBase = OldRowMap_->IndexBase();
00115 
00116   vector<int> ColNZCnt( NRows );
00117   vector<int> CS_RowIndices( NRows );
00118   map<int,int> RS_Map;
00119   map<int,int> CS_Map;
00120 
00121   int NumIndices;
00122   int * Indices;
00123   for( int i = 0; i < NRows; ++i )
00124   {
00125     ierr = OldGraph_->ExtractMyRowView( i, NumIndices, Indices );
00126 
00127     for( int j = 0; j < NumIndices; ++j )
00128     {
00129       ++ColNZCnt[ Indices[j] ];
00130       CS_RowIndices[ Indices[j] ] = i;
00131     }
00132 
00133     if( NumIndices == 1 ) RS_Map[i] = Indices[0];
00134   }
00135 
00136   if( verbose_ )
00137   {
00138     cout << "-------------------------\n";
00139     cout << "Row Singletons\n";
00140     for( map<int,int>::iterator itM = RS_Map.begin(); itM != RS_Map.end(); ++itM )
00141       cout << (*itM).first << "\t" << (*itM).second << endl;
00142     cout << "Col Counts\n";
00143     for( int i = 0; i < NRows; ++i )
00144       cout << i << "\t" << ColNZCnt[i] << "\t" << CS_RowIndices[i] << endl;
00145     cout << "-------------------------\n";
00146   }
00147 
00148   set<int> RS_Set;
00149   set<int> CS_Set;
00150 
00151   for( int i = 0; i < NRows; ++i )
00152     if( ColNZCnt[i] == 1 )
00153     {
00154       int RowIndex = CS_RowIndices[i];
00155       if( RS_Map.count(i) && RS_Map[i] == RowIndex )
00156       {
00157         CS_Set.insert(i);
00158         RS_Set.insert( RowIndex );
00159       }
00160     }
00161 
00162   if( verbose_ )
00163   {
00164     cout << "-------------------------\n";
00165     cout << "Singletons: " << CS_Set.size() << endl;
00166     set<int>::iterator itRS = RS_Set.begin();
00167     set<int>::iterator itCS = CS_Set.begin();
00168     set<int>::iterator endRS = RS_Set.end();
00169     set<int>::iterator endCS = CS_Set.end();
00170     for( ; itRS != endRS; ++itRS, ++itCS )
00171       cout << *itRS << "\t" << *itCS << endl;
00172     cout << "-------------------------\n";
00173   }
00174 
00175   //build new maps
00176   int NSingletons = CS_Set.size();
00177   int NReducedRows = NRows - 2* NSingletons; 
00178   vector<int> ReducedIndices( NReducedRows );
00179   vector<int> CSIndices( NSingletons );
00180   vector<int> RSIndices( NSingletons );
00181   int Reduced_Loc = 0;
00182   int RS_Loc = 0;
00183   int CS_Loc = 0;
00184   for( int i = 0; i < NRows; ++i )
00185   {
00186     int GlobalIndex = OldRowMap_->GID(i);
00187     if     ( RS_Set.count(i) ) RSIndices[RS_Loc++] = GlobalIndex;
00188     else if( CS_Set.count(i) ) CSIndices[CS_Loc++] = GlobalIndex;
00189     else                       ReducedIndices[Reduced_Loc++] = GlobalIndex;
00190   }
00191 
00192   vector<int> RowPermutedIndices( NRows );
00193   copy( RSIndices.begin(), RSIndices.end(), RowPermutedIndices.begin() );
00194   copy( ReducedIndices.begin(), ReducedIndices.end(), RowPermutedIndices.begin() + NSingletons );
00195   copy( CSIndices.begin(), CSIndices.end(), RowPermutedIndices.begin() + NReducedRows + NSingletons );
00196 
00197   vector<int> ColPermutedIndices( NRows );
00198   copy( CSIndices.begin(), CSIndices.end(), ColPermutedIndices.begin() );
00199   copy( ReducedIndices.begin(), ReducedIndices.end(), ColPermutedIndices.begin() + NSingletons );
00200   copy( RSIndices.begin(), RSIndices.end(), ColPermutedIndices.begin() + NReducedRows + NSingletons );
00201 
00202   UMap_ = new Epetra_Map( NSingletons, NSingletons, &RSIndices[0], IndexBase, CommObj );
00203   RMap_ = new Epetra_Map( NReducedRows, NReducedRows, &ReducedIndices[0], IndexBase, CommObj );
00204   LMap_ = new Epetra_Map( NSingletons, NSingletons, &CSIndices[0], IndexBase, CommObj );
00205 
00206   NewRowMap_ = new Epetra_Map( NRows, NRows, &RowPermutedIndices[0], IndexBase, CommObj );
00207   NewColMap_ = new Epetra_Map( NRows, NRows, &ColPermutedIndices[0], IndexBase, CommObj );
00208 
00209   //Construct Permuted System
00210   Exporter_ = new Epetra_Export( *OldRowMap_, *NewRowMap_ );
00211 
00212   NewRHS_ = new Epetra_MultiVector( *NewRowMap_, OldRHS_->NumVectors() );
00213   NewRHS_->Export( *OldRHS_, *Exporter_, Insert );
00214 
00215   NewLHS_ = new Epetra_MultiVector( *NewRowMap_, OldLHS_->NumVectors() );
00216   NewLHS_->Export( *OldLHS_, *Exporter_, Insert );
00217 
00218   NewGraph_ = new Epetra_CrsGraph( Copy, *NewRowMap_, *NewColMap_, 0 );
00219   NewGraph_->Export( *OldGraph_, *Exporter_, Insert );
00220   NewGraph_->FillComplete();
00221   cout << "Permuted Graph:\n" << *NewGraph_;
00222 
00223   NewMatrix_ = new Epetra_CrsMatrix( Copy, *NewGraph_ );
00224   NewMatrix_->Export( *OldMatrix_, *Exporter_, Insert );
00225   NewMatrix_->FillComplete();
00226   cout << "Permuted Matrix:\n" << *NewMatrix_;
00227 
00228   UExporter_ = new Epetra_Export( *OldRowMap_, *UMap_ );
00229   RExporter_ = new Epetra_Export( *OldRowMap_, *RMap_ );
00230   LExporter_ = new Epetra_Export( *OldRowMap_, *LMap_ );
00231 cout << "UExporter:\n" << *UExporter_;
00232 cout << "RExporter:\n" << *RExporter_;
00233 cout << "LExporter:\n" << *LExporter_;
00234 
00235   ULHS_ = new Epetra_MultiVector( *LMap_, OldLHS_->NumVectors() );
00236   ULHS_->Export( *OldLHS_, *LExporter_, Insert );
00237   cout << "ULHS:\n" << *ULHS_;
00238 
00239   RLHS_ = new Epetra_MultiVector( *RMap_, OldLHS_->NumVectors() );
00240   RLHS_->Export( *OldLHS_, *RExporter_, Insert );
00241   cout << "RLHS:\n" << *RLHS_;
00242 
00243   LLHS_ = new Epetra_MultiVector( *UMap_, OldLHS_->NumVectors() );
00244   LLHS_->Export( *OldLHS_, *UExporter_, Insert );
00245   cout << "LLHS:\n" << *LLHS_;
00246 
00247   URHS_ = new Epetra_MultiVector( *UMap_, OldRHS_->NumVectors() );
00248   URHS_->Export( *OldRHS_, *UExporter_, Insert );
00249   cout << "URHS:\n" << *URHS_;
00250 
00251   RRHS_ = new Epetra_MultiVector( *RMap_, OldRHS_->NumVectors() );
00252   RRHS_->Export( *OldRHS_, *RExporter_, Insert );
00253   cout << "RRHS:\n" << *RRHS_;
00254 
00255   LRHS_ = new Epetra_MultiVector( *LMap_, OldRHS_->NumVectors() );
00256   LRHS_->Export( *OldRHS_, *LExporter_, Insert );
00257   cout << "LRHS:\n" << *LRHS_;
00258 
00259   UUGraph_ = new Epetra_CrsGraph( Copy, *UMap_, *LMap_, 0 );
00260   UUGraph_->Export( *OldGraph_, *UExporter_, Insert );
00261   UUGraph_->FillComplete( LMap_, UMap_ );
00262   cout << "UUGraph:\n" << *UUGraph_;
00263 
00264   UUMatrix_ = new Epetra_CrsMatrix( Copy, *UUGraph_ );
00265   UUMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
00266   UUMatrix_->FillComplete();
00267   cout << "UUMatrix:\n" << *UUMatrix_;
00268 
00269   URGraph_ = new Epetra_CrsGraph( Copy, *UMap_, *RMap_, 0 );
00270   URGraph_->Export( *OldGraph_, *UExporter_, Insert );
00271   URGraph_->FillComplete( RMap_, UMap_ );
00272   cout << "URGraph:\n" << *URGraph_;
00273 
00274   URMatrix_ = new Epetra_CrsMatrix( Copy, *URGraph_ );
00275   URMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
00276   URMatrix_->FillComplete();
00277   cout << "URMatrix:\n" << *URMatrix_;
00278 
00279   ULGraph_ = new Epetra_CrsGraph( Copy, *UMap_, *UMap_, 0 );
00280   ULGraph_->Export( *OldGraph_, *UExporter_, Insert );
00281   ULGraph_->FillComplete();
00282   cout << "ULGraph:\n" << *ULGraph_;
00283 
00284   ULMatrix_ = new Epetra_CrsMatrix( Copy, *ULGraph_ );
00285   ULMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
00286   ULMatrix_->FillComplete();
00287   cout << "ULMatrix:\n" << *ULMatrix_;
00288 
00289   RRGraph_ = new Epetra_CrsGraph( Copy, *RMap_, *RMap_, 0 );
00290   RRGraph_->Export( *OldGraph_, *RExporter_, Insert );
00291   RRGraph_->FillComplete();
00292   cout << "RRGraph:\n" << *RRGraph_;
00293 
00294   RRMatrix_ = new Epetra_CrsMatrix( Copy, *RRGraph_ );
00295   RRMatrix_->Export( *OldMatrix_, *RExporter_, Insert );
00296   RRMatrix_->FillComplete();
00297   cout << "RRMatrix:\n" << *RRMatrix_;
00298 
00299   RLGraph_ = new Epetra_CrsGraph( Copy, *RMap_, *UMap_, 0 );
00300   RLGraph_->Export( *OldGraph_, *RExporter_, Insert );
00301   RLGraph_->FillComplete( UMap_, RMap_ );
00302   cout << "RLGraph:\n" << *RLGraph_;
00303 
00304   RLMatrix_ = new Epetra_CrsMatrix( Copy, *RLGraph_ );
00305   RLMatrix_->Export( *OldMatrix_, *RExporter_, Insert );
00306   RLMatrix_->FillComplete();
00307   cout << "RLMatrix:\n" << *RLMatrix_;
00308 
00309   LLGraph_ = new Epetra_CrsGraph( Copy, *LMap_, *UMap_, 0 );
00310   LLGraph_->Export( *OldGraph_, *LExporter_, Insert );
00311   LLGraph_->FillComplete( UMap_, LMap_ );
00312   cout << "LLGraph:\n" << *LLGraph_;
00313 
00314   LLMatrix_ = new Epetra_CrsMatrix( Copy, *LLGraph_ );
00315   LLMatrix_->Export( *OldMatrix_, *LExporter_, Insert );
00316   LLMatrix_->FillComplete();
00317   cout << "LLMatrix:\n" << *LLMatrix_;
00318 
00319   if( verbose_ )
00320   {
00321     cout << "Full System Characteristics:" << endl
00322          << "----------------------------" << endl
00323          << " Dimension                   = " << NRows << endl
00324          << " NNZs                        = " << OldMatrix_->NumGlobalNonzeros() << endl
00325          << " Max Num Row Entries         = " << OldMatrix_->GlobalMaxNumEntries() << endl << endl
00326          << "Reduced System Characteristics:" << endl
00327          << " Dimension                   = " << NReducedRows << endl
00328          << " NNZs                        = " << RRMatrix_->NumGlobalNonzeros() << endl
00329          << " Max Num Row Entries         = " << RRMatrix_->GlobalMaxNumEntries() << endl
00330          << "Singleton Info:" << endl
00331          << " Num Singleton                 = " << NSingletons << endl
00332          << "Ratios:" << endl
00333          << " % Reduction in Dimension    = " << 100.0*(NRows-NReducedRows)/NRows << endl
00334          << " % Reduction in NNZs         = " << (OldMatrix_->NumGlobalNonzeros()-RRMatrix_->NumGlobalNonzeros())/100.0 << endl
00335          << "-------------------------------" << endl;
00336   }
00337 
00338   NewProblem_ = new Epetra_LinearProblem( RRMatrix_, RLHS_, RRHS_ );
00339 
00340   newObj_ = NewProblem_;
00341 
00342   cout << "done with SC\n";
00343 
00344   return *NewProblem_;
00345 }
00346 
00347 bool
00348 LinearProblem_StaticCondensation::
00349 fwd()
00350 {
00351   if( verbose_ ) cout << "LP_SC::fwd()\n";
00352   if( verbose_ ) cout << "LP_SC::fwd() : Exporting to New LHS\n";
00353   ULHS_->Export( *OldLHS_, *LExporter_, Insert );
00354   RLHS_->Export( *OldLHS_, *RExporter_, Insert );
00355   LLHS_->Export( *OldLHS_, *UExporter_, Insert );
00356 
00357   if( verbose_ ) cout << "LP_SC::fwd() : Exporting to New RHS\n";
00358   URHS_->Export( *OldRHS_, *UExporter_, Insert );
00359   RRHS_->Export( *OldRHS_, *RExporter_, Insert );
00360   LRHS_->Export( *OldRHS_, *LExporter_, Insert );
00361 
00362   UUMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
00363   URMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
00364   ULMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
00365   RRMatrix_->Export( *OldMatrix_, *RExporter_, Insert );
00366   RLMatrix_->Export( *OldMatrix_, *RExporter_, Insert );
00367   LLMatrix_->Export( *OldMatrix_, *LExporter_, Insert );
00368 
00369   Epetra_Vector LLDiag( *LMap_ );
00370   LLMatrix_->ExtractDiagonalCopy( LLDiag );
00371   Epetra_Vector LLRecipDiag( *LMap_ );
00372   LLRecipDiag.Reciprocal( LLDiag );
00373 
00374   if( verbose_ ) cout << "LP_SC::fwd() : Forming LLHS\n";
00375   LLDiag.Multiply( 1.0, LLRecipDiag, *LRHS_, 0.0 );
00376   int LSize = LMap_->NumMyElements();
00377   for( int i = 0; i < LSize; ++i ) (*LLHS_)[0][i] = LLDiag[i];
00378 
00379   if( verbose_ ) cout << "LP_SC::fwd() : Updating RRHS\n";
00380   Epetra_Vector RUpdate( *RMap_ );
00381   RLMatrix_->Multiply( false, *LLHS_, RUpdate );
00382   RRHS_->Update( -1.0, RUpdate, 1.0 );
00383 
00384   if( verbose_ ) cout << "LP_SC::fwd() : Updating URHS\n";
00385   Epetra_Vector UUpdate( *UMap_ );
00386   ULMatrix_->Multiply( false, *LLHS_, UUpdate );
00387   URHS_->Update( -1.0, UUpdate, 1.0 );
00388   
00389   return true;
00390 }
00391 
00392 bool
00393 LinearProblem_StaticCondensation::
00394 rvs()
00395 {
00396   if( verbose_ ) cout << "LP_SC::rvs()\n";
00397   if( verbose_ ) cout << "LP_SC::rvs() : Updating URHS\n";
00398   Epetra_Vector UUpdate( *UMap_ );
00399   URMatrix_->Multiply( false, *RLHS_, UUpdate );
00400   URHS_->Update( -1.0, UUpdate, 1.0 );
00401 
00402   Epetra_Vector UUDiag( *UMap_ );
00403   UUMatrix_->ExtractDiagonalCopy( UUDiag );
00404   Epetra_Vector UURecipDiag( *UMap_ );
00405   UURecipDiag.Reciprocal( UUDiag );
00406 
00407   if( verbose_ ) cout << "LP_SC::rvs() : Forming ULHS\n";
00408   UUDiag.Multiply( 1.0, UURecipDiag, *URHS_, 0.0 );
00409   int USize = UMap_->NumMyElements();
00410   for( int i = 0; i < USize; ++i ) (*ULHS_)[0][i] = UUDiag[i];
00411 
00412   if( verbose_ ) cout << "LP_SC::rvs() : Importing to Old LHS\n";
00413   OldLHS_->Import( *ULHS_, *LExporter_, Insert );
00414   OldLHS_->Import( *RLHS_, *RExporter_, Insert );
00415   OldLHS_->Import( *LLHS_, *UExporter_, Insert );
00416 
00417   return true;
00418 }
00419 
00420 } // namespace EpetraExt
00421 

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