44#include <Epetra_Export.h>
45#include <Epetra_LinearProblem.h>
46#include <Epetra_CrsGraph.h>
47#include <Epetra_CrsMatrix.h>
48#include <Epetra_MultiVector.h>
49#include <Epetra_Vector.h>
50#include <Epetra_IntVector.h>
51#include <Epetra_Map.h>
52#include <Epetra_Comm.h>
63 if( Exporter_ )
delete Exporter_;
65 if( NewProblem_ )
delete NewProblem_;
66 if( NewRHS_ )
delete NewRHS_;
67 if( NewLHS_ )
delete NewLHS_;
68 if( NewMatrix_ )
delete NewMatrix_;
69 if( NewGraph_ )
delete NewGraph_;
70 if( NewRowMap_ )
delete NewRowMap_;
71 if( NewColMap_ )
delete NewColMap_;
73 if( ULHS_ )
delete ULHS_;
74 if( RLHS_ )
delete RLHS_;
75 if( LLHS_ )
delete LLHS_;
76 if( URHS_ )
delete URHS_;
77 if( RRHS_ )
delete RRHS_;
78 if( LRHS_ )
delete LRHS_;
80 if( UUMatrix_ )
delete UUMatrix_;
81 if( URMatrix_ )
delete URMatrix_;
82 if( ULMatrix_ )
delete ULMatrix_;
83 if( RRMatrix_ )
delete RRMatrix_;
84 if( RLMatrix_ )
delete RLMatrix_;
85 if( LLMatrix_ )
delete LLMatrix_;
87 if( UUGraph_ )
delete UUGraph_;
88 if( URGraph_ )
delete URGraph_;
89 if( ULGraph_ )
delete ULGraph_;
90 if( RRGraph_ )
delete RRGraph_;
91 if( RLGraph_ )
delete RLGraph_;
92 if( LLGraph_ )
delete LLGraph_;
94 if( UExporter_ )
delete UExporter_;
95 if( RExporter_ )
delete RExporter_;
96 if( LExporter_ )
delete LExporter_;
98 if( UMap_ )
delete UMap_;
99 if( RMap_ )
delete RMap_;
100 if( LMap_ )
delete LMap_;
112 OldGraph_ = &OldMatrix_->
Graph();
113 OldRHS_ = orig.GetRHS();
114 OldLHS_ = orig.GetLHS();
115 OldRowMap_ = &OldMatrix_->
RowMap();
119 if( !OldMatrix_ ) ierr = -2;
120 if( !OldRHS_ ) ierr = -3;
121 if( !OldLHS_ ) ierr = -4;
124 if( degree_ != 1 ) ierr = -6;
129 vector<int> ColNZCnt( NRows );
130 vector<int> CS_RowIndices( NRows );
136 for(
int i = 0; i < NRows; ++i )
140 for(
int j = 0; j < NumIndices; ++j )
142 ++ColNZCnt[ Indices[j] ];
143 CS_RowIndices[ Indices[j] ] = i;
146 if( NumIndices == 1 ) RS_Map[i] = Indices[0];
151 cout <<
"-------------------------\n";
152 cout <<
"Row Singletons\n";
153 for( map<int,int>::iterator itM = RS_Map.begin(); itM != RS_Map.end(); ++itM )
154 cout << (*itM).first <<
"\t" << (*itM).second << endl;
155 cout <<
"Col Counts\n";
156 for(
int i = 0; i < NRows; ++i )
157 cout << i <<
"\t" << ColNZCnt[i] <<
"\t" << CS_RowIndices[i] << endl;
158 cout <<
"-------------------------\n";
164 for(
int i = 0; i < NRows; ++i )
165 if( ColNZCnt[i] == 1 )
167 int RowIndex = CS_RowIndices[i];
168 if( RS_Map.count(i) && RS_Map[i] == RowIndex )
171 RS_Set.insert( RowIndex );
177 cout <<
"-------------------------\n";
178 cout <<
"Singletons: " << CS_Set.size() << endl;
179 set<int>::iterator itRS = RS_Set.begin();
180 set<int>::iterator itCS = CS_Set.begin();
181 set<int>::iterator endRS = RS_Set.end();
182 set<int>::iterator endCS = CS_Set.end();
183 for( ; itRS != endRS; ++itRS, ++itCS )
184 cout << *itRS <<
"\t" << *itCS << endl;
185 cout <<
"-------------------------\n";
189 int NSingletons = CS_Set.size();
190 int NReducedRows = NRows - 2* NSingletons;
191 vector<int> ReducedIndices( NReducedRows );
192 vector<int> CSIndices( NSingletons );
193 vector<int> RSIndices( NSingletons );
197 for(
int i = 0; i < NRows; ++i )
199 int GlobalIndex = OldRowMap_->
GID(i);
200 if ( RS_Set.count(i) ) RSIndices[RS_Loc++] = GlobalIndex;
201 else if( CS_Set.count(i) ) CSIndices[CS_Loc++] = GlobalIndex;
202 else ReducedIndices[Reduced_Loc++] = GlobalIndex;
205 vector<int> RowPermutedIndices( NRows );
206 copy( RSIndices.begin(), RSIndices.end(), RowPermutedIndices.begin() );
207 copy( ReducedIndices.begin(), ReducedIndices.end(), RowPermutedIndices.begin() + NSingletons );
208 copy( CSIndices.begin(), CSIndices.end(), RowPermutedIndices.begin() + NReducedRows + NSingletons );
210 vector<int> ColPermutedIndices( NRows );
211 copy( CSIndices.begin(), CSIndices.end(), ColPermutedIndices.begin() );
212 copy( ReducedIndices.begin(), ReducedIndices.end(), ColPermutedIndices.begin() + NSingletons );
213 copy( RSIndices.begin(), RSIndices.end(), ColPermutedIndices.begin() + NReducedRows + NSingletons );
215 UMap_ =
new Epetra_Map( NSingletons, NSingletons, &RSIndices[0], IndexBase, CommObj );
216 RMap_ =
new Epetra_Map( NReducedRows, NReducedRows, &ReducedIndices[0], IndexBase, CommObj );
217 LMap_ =
new Epetra_Map( NSingletons, NSingletons, &CSIndices[0], IndexBase, CommObj );
219 NewRowMap_ =
new Epetra_Map( NRows, NRows, &RowPermutedIndices[0], IndexBase, CommObj );
220 NewColMap_ =
new Epetra_Map( NRows, NRows, &ColPermutedIndices[0], IndexBase, CommObj );
234 cout <<
"Permuted Graph:\n" << *NewGraph_;
239 cout <<
"Permuted Matrix:\n" << *NewMatrix_;
244cout <<
"UExporter:\n" << *UExporter_;
245cout <<
"RExporter:\n" << *RExporter_;
246cout <<
"LExporter:\n" << *LExporter_;
250 cout <<
"ULHS:\n" << *ULHS_;
254 cout <<
"RLHS:\n" << *RLHS_;
258 cout <<
"LLHS:\n" << *LLHS_;
262 cout <<
"URHS:\n" << *URHS_;
266 cout <<
"RRHS:\n" << *RRHS_;
270 cout <<
"LRHS:\n" << *LRHS_;
275 cout <<
"UUGraph:\n" << *UUGraph_;
280 cout <<
"UUMatrix:\n" << *UUMatrix_;
285 cout <<
"URGraph:\n" << *URGraph_;
290 cout <<
"URMatrix:\n" << *URMatrix_;
295 cout <<
"ULGraph:\n" << *ULGraph_;
300 cout <<
"ULMatrix:\n" << *ULMatrix_;
305 cout <<
"RRGraph:\n" << *RRGraph_;
310 cout <<
"RRMatrix:\n" << *RRMatrix_;
315 cout <<
"RLGraph:\n" << *RLGraph_;
320 cout <<
"RLMatrix:\n" << *RLMatrix_;
325 cout <<
"LLGraph:\n" << *LLGraph_;
330 cout <<
"LLMatrix:\n" << *LLMatrix_;
334 cout <<
"Full System Characteristics:" << endl
335 <<
"----------------------------" << endl
336 <<
" Dimension = " << NRows << endl
339 <<
"Reduced System Characteristics:" << endl
340 <<
" Dimension = " << NReducedRows << endl
341 <<
" NNZs = " << RRMatrix_->NumGlobalNonzeros() << endl
342 <<
" Max Num Row Entries = " << RRMatrix_->GlobalMaxNumEntries() << endl
343 <<
"Singleton Info:" << endl
344 <<
" Num Singleton = " << NSingletons << endl
346 <<
" % Reduction in Dimension = " << 100.0*(NRows-NReducedRows)/NRows << endl
347 <<
" % Reduction in NNZs = " << (OldMatrix_->
NumGlobalNonzeros()-RRMatrix_->NumGlobalNonzeros())/100.0 << endl
348 <<
"-------------------------------" << endl;
355 cout <<
"done with SC\n";
364 if( verbose_ ) cout <<
"LP_SC::fwd()\n";
365 if( verbose_ ) cout <<
"LP_SC::fwd() : Exporting to New LHS\n";
370 if( verbose_ ) cout <<
"LP_SC::fwd() : Exporting to New RHS\n";
387 if( verbose_ ) cout <<
"LP_SC::fwd() : Forming LLHS\n";
388 LLDiag.
Multiply( 1.0, LLRecipDiag, *LRHS_, 0.0 );
390 for(
int i = 0; i < LSize; ++i ) (*LLHS_)[0][i] = LLDiag[i];
392 if( verbose_ ) cout <<
"LP_SC::fwd() : Updating RRHS\n";
394 RLMatrix_->
Multiply(
false, *LLHS_, RUpdate );
395 RRHS_->
Update( -1.0, RUpdate, 1.0 );
397 if( verbose_ ) cout <<
"LP_SC::fwd() : Updating URHS\n";
399 ULMatrix_->
Multiply(
false, *LLHS_, UUpdate );
400 URHS_->
Update( -1.0, UUpdate, 1.0 );
409 if( verbose_ ) cout <<
"LP_SC::rvs()\n";
410 if( verbose_ ) cout <<
"LP_SC::rvs() : Updating URHS\n";
412 URMatrix_->
Multiply(
false, *RLHS_, UUpdate );
413 URHS_->
Update( -1.0, UUpdate, 1.0 );
420 if( verbose_ ) cout <<
"LP_SC::rvs() : Forming ULHS\n";
421 UUDiag.
Multiply( 1.0, UURecipDiag, *URHS_, 0.0 );
423 for(
int i = 0; i < USize; ++i ) (*ULHS_)[0][i] = UUDiag[i];
425 if( verbose_ ) cout <<
"LP_SC::rvs() : Importing to Old LHS\n";
bool rvs()
Reverse transfer of data from new object created in the operator() method call to the orig object inp...
bool fwd()
Forward transfer of data from orig object input in the operator() method call to the new object creat...
NewTypeRef operator()(OriginalTypeRef orig)
Analysis of transform operation on original object and construction of new object.
~LinearProblem_StaticCondensation()
bool DistributedGlobal() const
const Epetra_Comm & Comm() const
int NumMyElements() const
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
const Epetra_Map & RowMap() const
int GlobalMaxNumEntries() const
int FillComplete(bool OptimizeDataStorage=true)
int NumGlobalNonzeros() const
const Epetra_CrsGraph & Graph() const
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int Reciprocal(const Epetra_MultiVector &A)
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.