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> 111 OldMatrix_ =
dynamic_cast<Epetra_CrsMatrix*
>( orig.GetMatrix() );
117 const Epetra_Comm & CommObj =
OldRowMap_->Comm();
123 if(
OldRowMap_->DistributedGlobal() ) ierr = -5;
129 vector<int> ColNZCnt( NRows );
130 vector<int> CS_RowIndices( NRows );
136 for(
int i = 0; i < NRows; ++i )
138 ierr =
OldGraph_->ExtractMyRowView( i, NumIndices, Indices );
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 )
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_;
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_;
334 cout <<
"Full System Characteristics:" << endl
335 <<
"----------------------------" << endl
336 <<
" Dimension = " << NRows << endl
337 <<
" NNZs = " <<
OldMatrix_->NumGlobalNonzeros() << endl
338 <<
" Max Num Row Entries = " <<
OldMatrix_->GlobalMaxNumEntries() << endl << 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";
382 Epetra_Vector LLDiag( *
LMap_ );
383 LLMatrix_->ExtractDiagonalCopy( LLDiag );
384 Epetra_Vector LLRecipDiag( *
LMap_ );
385 LLRecipDiag.Reciprocal( LLDiag );
387 if(
verbose_ ) cout <<
"LP_SC::fwd() : Forming LLHS\n";
388 LLDiag.Multiply( 1.0, LLRecipDiag, *
LRHS_, 0.0 );
389 int LSize =
LMap_->NumMyElements();
390 for(
int i = 0; i < LSize; ++i ) (*
LLHS_)[0][i] = LLDiag[i];
392 if(
verbose_ ) cout <<
"LP_SC::fwd() : Updating RRHS\n";
393 Epetra_Vector RUpdate( *
RMap_ );
395 RRHS_->Update( -1.0, RUpdate, 1.0 );
397 if(
verbose_ ) cout <<
"LP_SC::fwd() : Updating URHS\n";
398 Epetra_Vector UUpdate( *
UMap_ );
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";
411 Epetra_Vector UUpdate( *
UMap_ );
413 URHS_->Update( -1.0, UUpdate, 1.0 );
415 Epetra_Vector UUDiag( *
UMap_ );
416 UUMatrix_->ExtractDiagonalCopy( UUDiag );
417 Epetra_Vector UURecipDiag( *
UMap_ );
418 UURecipDiag.Reciprocal( UUDiag );
420 if(
verbose_ ) cout <<
"LP_SC::rvs() : Forming ULHS\n";
421 UUDiag.Multiply( 1.0, UURecipDiag, *
URHS_, 0.0 );
422 int USize =
UMap_->NumMyElements();
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";
Epetra_MultiVector * OldLHS_
Epetra_MultiVector * ULHS_
Epetra_CrsGraph * URGraph_
Epetra_CrsMatrix * UUMatrix_
Epetra_CrsGraph * ULGraph_
Epetra_CrsMatrix * OldMatrix_
Epetra_CrsGraph * RRGraph_
Epetra_CrsMatrix * RRMatrix_
const Epetra_Map * OldRowMap_
Epetra_MultiVector * URHS_
Epetra_Export * RExporter_
const Epetra_CrsGraph * OldGraph_
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
bool fwd()
Forward transfer of data from orig object input in the operator() method call to the new object creat...
Epetra_CrsGraph * LLGraph_
Epetra_CrsMatrix * URMatrix_
Epetra_MultiVector * NewLHS_
NewTypeRef operator()(OriginalTypeRef orig)
Epetra_CrsGraph * UUGraph_
Epetra_MultiVector * NewRHS_
Epetra_MultiVector * LRHS_
Epetra_MultiVector * RRHS_
Epetra_LinearProblem * NewProblem_
Epetra_CrsGraph * RLGraph_
Epetra_CrsGraph * NewGraph_
Epetra_Export * UExporter_
Epetra_MultiVector * LLHS_
~LinearProblem_StaticCondensation()
Epetra_CrsMatrix * NewMatrix_
Epetra_CrsMatrix * RLMatrix_
Epetra_MultiVector * RLHS_
Epetra_CrsMatrix * ULMatrix_
Epetra_Export * Exporter_
Epetra_CrsMatrix * LLMatrix_
Epetra_Export * LExporter_
bool rvs()
Reverse transfer of data from new object created in the operator() method call to the orig object inp...
Epetra_MultiVector * OldRHS_