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_;
111 OldMatrix_ =
dynamic_cast<Epetra_CrsMatrix*
>( orig.GetMatrix() );
112 OldGraph_ = &OldMatrix_->Graph();
113 OldRHS_ = orig.GetRHS();
114 OldLHS_ = orig.GetLHS();
115 OldRowMap_ = &OldMatrix_->RowMap();
117 const Epetra_Comm & CommObj = OldRowMap_->Comm();
119 if( !OldMatrix_ ) ierr = -2;
120 if( !OldRHS_ ) ierr = -3;
121 if( !OldLHS_ ) ierr = -4;
123 if( OldRowMap_->DistributedGlobal() ) ierr = -5;
124 if( degree_ != 1 ) ierr = -6;
126 int NRows = OldGraph_->NumMyRows();
127 int IndexBase = OldRowMap_->IndexBase();
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 )
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 );
223 Exporter_ =
new Epetra_Export( *OldRowMap_, *NewRowMap_ );
225 NewRHS_ =
new Epetra_MultiVector( *NewRowMap_, OldRHS_->NumVectors() );
226 NewRHS_->Export( *OldRHS_, *Exporter_, Insert );
228 NewLHS_ =
new Epetra_MultiVector( *NewRowMap_, OldLHS_->NumVectors() );
229 NewLHS_->Export( *OldLHS_, *Exporter_, Insert );
231 NewGraph_ =
new Epetra_CrsGraph( Copy, *NewRowMap_, *NewColMap_, 0 );
232 NewGraph_->Export( *OldGraph_, *Exporter_, Insert );
233 NewGraph_->FillComplete();
234 cout <<
"Permuted Graph:\n" << *NewGraph_;
236 NewMatrix_ =
new Epetra_CrsMatrix( Copy, *NewGraph_ );
237 NewMatrix_->Export( *OldMatrix_, *Exporter_, Insert );
238 NewMatrix_->FillComplete();
239 cout <<
"Permuted Matrix:\n" << *NewMatrix_;
241 UExporter_ =
new Epetra_Export( *OldRowMap_, *UMap_ );
242 RExporter_ =
new Epetra_Export( *OldRowMap_, *RMap_ );
243 LExporter_ =
new Epetra_Export( *OldRowMap_, *LMap_ );
244 cout <<
"UExporter:\n" << *UExporter_;
245 cout <<
"RExporter:\n" << *RExporter_;
246 cout <<
"LExporter:\n" << *LExporter_;
248 ULHS_ =
new Epetra_MultiVector( *LMap_, OldLHS_->NumVectors() );
249 ULHS_->Export( *OldLHS_, *LExporter_, Insert );
250 cout <<
"ULHS:\n" << *ULHS_;
252 RLHS_ =
new Epetra_MultiVector( *RMap_, OldLHS_->NumVectors() );
253 RLHS_->Export( *OldLHS_, *RExporter_, Insert );
254 cout <<
"RLHS:\n" << *RLHS_;
256 LLHS_ =
new Epetra_MultiVector( *UMap_, OldLHS_->NumVectors() );
257 LLHS_->Export( *OldLHS_, *UExporter_, Insert );
258 cout <<
"LLHS:\n" << *LLHS_;
260 URHS_ =
new Epetra_MultiVector( *UMap_, OldRHS_->NumVectors() );
261 URHS_->Export( *OldRHS_, *UExporter_, Insert );
262 cout <<
"URHS:\n" << *URHS_;
264 RRHS_ =
new Epetra_MultiVector( *RMap_, OldRHS_->NumVectors() );
265 RRHS_->Export( *OldRHS_, *RExporter_, Insert );
266 cout <<
"RRHS:\n" << *RRHS_;
268 LRHS_ =
new Epetra_MultiVector( *LMap_, OldRHS_->NumVectors() );
269 LRHS_->Export( *OldRHS_, *LExporter_, Insert );
270 cout <<
"LRHS:\n" << *LRHS_;
272 UUGraph_ =
new Epetra_CrsGraph( Copy, *UMap_, *LMap_, 0 );
273 UUGraph_->Export( *OldGraph_, *UExporter_, Insert );
274 UUGraph_->FillComplete( LMap_, UMap_ );
275 cout <<
"UUGraph:\n" << *UUGraph_;
277 UUMatrix_ =
new Epetra_CrsMatrix( Copy, *UUGraph_ );
278 UUMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
279 UUMatrix_->FillComplete();
280 cout <<
"UUMatrix:\n" << *UUMatrix_;
282 URGraph_ =
new Epetra_CrsGraph( Copy, *UMap_, *RMap_, 0 );
283 URGraph_->Export( *OldGraph_, *UExporter_, Insert );
284 URGraph_->FillComplete( RMap_, UMap_ );
285 cout <<
"URGraph:\n" << *URGraph_;
287 URMatrix_ =
new Epetra_CrsMatrix( Copy, *URGraph_ );
288 URMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
289 URMatrix_->FillComplete();
290 cout <<
"URMatrix:\n" << *URMatrix_;
292 ULGraph_ =
new Epetra_CrsGraph( Copy, *UMap_, *UMap_, 0 );
293 ULGraph_->Export( *OldGraph_, *UExporter_, Insert );
294 ULGraph_->FillComplete();
295 cout <<
"ULGraph:\n" << *ULGraph_;
297 ULMatrix_ =
new Epetra_CrsMatrix( Copy, *ULGraph_ );
298 ULMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
299 ULMatrix_->FillComplete();
300 cout <<
"ULMatrix:\n" << *ULMatrix_;
302 RRGraph_ =
new Epetra_CrsGraph( Copy, *RMap_, *RMap_, 0 );
303 RRGraph_->Export( *OldGraph_, *RExporter_, Insert );
304 RRGraph_->FillComplete();
305 cout <<
"RRGraph:\n" << *RRGraph_;
307 RRMatrix_ =
new Epetra_CrsMatrix( Copy, *RRGraph_ );
308 RRMatrix_->Export( *OldMatrix_, *RExporter_, Insert );
309 RRMatrix_->FillComplete();
310 cout <<
"RRMatrix:\n" << *RRMatrix_;
312 RLGraph_ =
new Epetra_CrsGraph( Copy, *RMap_, *UMap_, 0 );
313 RLGraph_->Export( *OldGraph_, *RExporter_, Insert );
314 RLGraph_->FillComplete( UMap_, RMap_ );
315 cout <<
"RLGraph:\n" << *RLGraph_;
317 RLMatrix_ =
new Epetra_CrsMatrix( Copy, *RLGraph_ );
318 RLMatrix_->Export( *OldMatrix_, *RExporter_, Insert );
319 RLMatrix_->FillComplete();
320 cout <<
"RLMatrix:\n" << *RLMatrix_;
322 LLGraph_ =
new Epetra_CrsGraph( Copy, *LMap_, *UMap_, 0 );
323 LLGraph_->Export( *OldGraph_, *LExporter_, Insert );
324 LLGraph_->FillComplete( UMap_, LMap_ );
325 cout <<
"LLGraph:\n" << *LLGraph_;
327 LLMatrix_ =
new Epetra_CrsMatrix( Copy, *LLGraph_ );
328 LLMatrix_->Export( *OldMatrix_, *LExporter_, Insert );
329 LLMatrix_->FillComplete();
330 cout <<
"LLMatrix:\n" << *LLMatrix_;
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;
351 NewProblem_ =
new Epetra_LinearProblem( RRMatrix_, RLHS_, RRHS_ );
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";
366 ULHS_->Export( *OldLHS_, *LExporter_, Insert );
367 RLHS_->Export( *OldLHS_, *RExporter_, Insert );
368 LLHS_->Export( *OldLHS_, *UExporter_, Insert );
370 if( verbose_ ) cout <<
"LP_SC::fwd() : Exporting to New RHS\n";
371 URHS_->Export( *OldRHS_, *UExporter_, Insert );
372 RRHS_->Export( *OldRHS_, *RExporter_, Insert );
373 LRHS_->Export( *OldRHS_, *LExporter_, Insert );
375 UUMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
376 URMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
377 ULMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
378 RRMatrix_->Export( *OldMatrix_, *RExporter_, Insert );
379 RLMatrix_->Export( *OldMatrix_, *RExporter_, Insert );
380 LLMatrix_->Export( *OldMatrix_, *LExporter_, Insert );
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_ );
394 RLMatrix_->Multiply(
false, *LLHS_, RUpdate );
395 RRHS_->Update( -1.0, RUpdate, 1.0 );
397 if( verbose_ ) cout <<
"LP_SC::fwd() : Updating URHS\n";
398 Epetra_Vector UUpdate( *UMap_ );
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";
411 Epetra_Vector UUpdate( *UMap_ );
412 URMatrix_->Multiply(
false, *RLHS_, UUpdate );
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";
426 OldLHS_->Import( *ULHS_, *LExporter_, Insert );
427 OldLHS_->Import( *RLHS_, *RExporter_, Insert );
428 OldLHS_->Import( *LLHS_, *UExporter_, Insert );
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...
NewTypeRef operator()(OriginalTypeRef orig)
~LinearProblem_StaticCondensation()
bool rvs()
Reverse transfer of data from new object created in the operator() method call to the orig object inp...