EpetraExt Development
Loading...
Searching...
No Matches
EpetraExt_StaticCondensation_LinearProblem.cpp
Go to the documentation of this file.
1//@HEADER
2// ***********************************************************************
3//
4// EpetraExt: Epetra Extended - Linear Algebra Services Package
5// Copyright (2011) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41
43
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>
53
54#include <vector>
55#include <map>
56#include <set>
57
58namespace EpetraExt {
59
62{
63 if( Exporter_ ) delete Exporter_;
64
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_;
72
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_;
79
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_;
86
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_;
93
94 if( UExporter_ ) delete UExporter_;
95 if( RExporter_ ) delete RExporter_;
96 if( LExporter_ ) delete LExporter_;
97
98 if( UMap_ ) delete UMap_;
99 if( RMap_ ) delete RMap_;
100 if( LMap_ ) delete LMap_;
101}
102
106{
107 origObj_ = &orig;
108
109 int ierr = 0;
110
111 OldMatrix_ = dynamic_cast<Epetra_CrsMatrix*>( orig.GetMatrix() );
112 OldGraph_ = &OldMatrix_->Graph();
113 OldRHS_ = orig.GetRHS();
114 OldLHS_ = orig.GetLHS();
115 OldRowMap_ = &OldMatrix_->RowMap();
116
117 const Epetra_Comm & CommObj = OldRowMap_->Comm();
118
119 if( !OldMatrix_ ) ierr = -2;
120 if( !OldRHS_ ) ierr = -3;
121 if( !OldLHS_ ) ierr = -4;
122
123 if( OldRowMap_->DistributedGlobal() ) ierr = -5;
124 if( degree_ != 1 ) ierr = -6;
125
126 int NRows = OldGraph_->NumMyRows();
127 int IndexBase = OldRowMap_->IndexBase();
128
129 vector<int> ColNZCnt( NRows );
130 vector<int> CS_RowIndices( NRows );
131 map<int,int> RS_Map;
132 map<int,int> CS_Map;
133
134 int NumIndices;
135 int * Indices;
136 for( int i = 0; i < NRows; ++i )
137 {
138 ierr = OldGraph_->ExtractMyRowView( i, NumIndices, Indices );
139
140 for( int j = 0; j < NumIndices; ++j )
141 {
142 ++ColNZCnt[ Indices[j] ];
143 CS_RowIndices[ Indices[j] ] = i;
144 }
145
146 if( NumIndices == 1 ) RS_Map[i] = Indices[0];
147 }
148
149 if( verbose_ )
150 {
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";
159 }
160
161 set<int> RS_Set;
162 set<int> CS_Set;
163
164 for( int i = 0; i < NRows; ++i )
165 if( ColNZCnt[i] == 1 )
166 {
167 int RowIndex = CS_RowIndices[i];
168 if( RS_Map.count(i) && RS_Map[i] == RowIndex )
169 {
170 CS_Set.insert(i);
171 RS_Set.insert( RowIndex );
172 }
173 }
174
175 if( verbose_ )
176 {
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";
186 }
187
188 //build new maps
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 );
194 int Reduced_Loc = 0;
195 int RS_Loc = 0;
196 int CS_Loc = 0;
197 for( int i = 0; i < NRows; ++i )
198 {
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;
203 }
204
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 );
209
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 );
214
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 );
218
219 NewRowMap_ = new Epetra_Map( NRows, NRows, &RowPermutedIndices[0], IndexBase, CommObj );
220 NewColMap_ = new Epetra_Map( NRows, NRows, &ColPermutedIndices[0], IndexBase, CommObj );
221
222 //Construct Permuted System
223 Exporter_ = new Epetra_Export( *OldRowMap_, *NewRowMap_ );
224
225 NewRHS_ = new Epetra_MultiVector( *NewRowMap_, OldRHS_->NumVectors() );
226 NewRHS_->Export( *OldRHS_, *Exporter_, Insert );
227
228 NewLHS_ = new Epetra_MultiVector( *NewRowMap_, OldLHS_->NumVectors() );
229 NewLHS_->Export( *OldLHS_, *Exporter_, Insert );
230
231 NewGraph_ = new Epetra_CrsGraph( Copy, *NewRowMap_, *NewColMap_, 0 );
232 NewGraph_->Export( *OldGraph_, *Exporter_, Insert );
233 NewGraph_->FillComplete();
234 cout << "Permuted Graph:\n" << *NewGraph_;
235
236 NewMatrix_ = new Epetra_CrsMatrix( Copy, *NewGraph_ );
237 NewMatrix_->Export( *OldMatrix_, *Exporter_, Insert );
238 NewMatrix_->FillComplete();
239 cout << "Permuted Matrix:\n" << *NewMatrix_;
240
241 UExporter_ = new Epetra_Export( *OldRowMap_, *UMap_ );
242 RExporter_ = new Epetra_Export( *OldRowMap_, *RMap_ );
243 LExporter_ = new Epetra_Export( *OldRowMap_, *LMap_ );
244cout << "UExporter:\n" << *UExporter_;
245cout << "RExporter:\n" << *RExporter_;
246cout << "LExporter:\n" << *LExporter_;
247
248 ULHS_ = new Epetra_MultiVector( *LMap_, OldLHS_->NumVectors() );
249 ULHS_->Export( *OldLHS_, *LExporter_, Insert );
250 cout << "ULHS:\n" << *ULHS_;
251
252 RLHS_ = new Epetra_MultiVector( *RMap_, OldLHS_->NumVectors() );
253 RLHS_->Export( *OldLHS_, *RExporter_, Insert );
254 cout << "RLHS:\n" << *RLHS_;
255
256 LLHS_ = new Epetra_MultiVector( *UMap_, OldLHS_->NumVectors() );
257 LLHS_->Export( *OldLHS_, *UExporter_, Insert );
258 cout << "LLHS:\n" << *LLHS_;
259
260 URHS_ = new Epetra_MultiVector( *UMap_, OldRHS_->NumVectors() );
261 URHS_->Export( *OldRHS_, *UExporter_, Insert );
262 cout << "URHS:\n" << *URHS_;
263
264 RRHS_ = new Epetra_MultiVector( *RMap_, OldRHS_->NumVectors() );
265 RRHS_->Export( *OldRHS_, *RExporter_, Insert );
266 cout << "RRHS:\n" << *RRHS_;
267
268 LRHS_ = new Epetra_MultiVector( *LMap_, OldRHS_->NumVectors() );
269 LRHS_->Export( *OldRHS_, *LExporter_, Insert );
270 cout << "LRHS:\n" << *LRHS_;
271
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_;
276
277 UUMatrix_ = new Epetra_CrsMatrix( Copy, *UUGraph_ );
278 UUMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
279 UUMatrix_->FillComplete();
280 cout << "UUMatrix:\n" << *UUMatrix_;
281
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_;
286
287 URMatrix_ = new Epetra_CrsMatrix( Copy, *URGraph_ );
288 URMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
289 URMatrix_->FillComplete();
290 cout << "URMatrix:\n" << *URMatrix_;
291
292 ULGraph_ = new Epetra_CrsGraph( Copy, *UMap_, *UMap_, 0 );
293 ULGraph_->Export( *OldGraph_, *UExporter_, Insert );
294 ULGraph_->FillComplete();
295 cout << "ULGraph:\n" << *ULGraph_;
296
297 ULMatrix_ = new Epetra_CrsMatrix( Copy, *ULGraph_ );
298 ULMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
299 ULMatrix_->FillComplete();
300 cout << "ULMatrix:\n" << *ULMatrix_;
301
302 RRGraph_ = new Epetra_CrsGraph( Copy, *RMap_, *RMap_, 0 );
303 RRGraph_->Export( *OldGraph_, *RExporter_, Insert );
304 RRGraph_->FillComplete();
305 cout << "RRGraph:\n" << *RRGraph_;
306
307 RRMatrix_ = new Epetra_CrsMatrix( Copy, *RRGraph_ );
308 RRMatrix_->Export( *OldMatrix_, *RExporter_, Insert );
309 RRMatrix_->FillComplete();
310 cout << "RRMatrix:\n" << *RRMatrix_;
311
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_;
316
317 RLMatrix_ = new Epetra_CrsMatrix( Copy, *RLGraph_ );
318 RLMatrix_->Export( *OldMatrix_, *RExporter_, Insert );
319 RLMatrix_->FillComplete();
320 cout << "RLMatrix:\n" << *RLMatrix_;
321
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_;
326
327 LLMatrix_ = new Epetra_CrsMatrix( Copy, *LLGraph_ );
328 LLMatrix_->Export( *OldMatrix_, *LExporter_, Insert );
329 LLMatrix_->FillComplete();
330 cout << "LLMatrix:\n" << *LLMatrix_;
331
332 if( verbose_ )
333 {
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
345 << "Ratios:" << endl
346 << " % Reduction in Dimension = " << 100.0*(NRows-NReducedRows)/NRows << endl
347 << " % Reduction in NNZs = " << (OldMatrix_->NumGlobalNonzeros()-RRMatrix_->NumGlobalNonzeros())/100.0 << endl
348 << "-------------------------------" << endl;
349 }
350
351 NewProblem_ = new Epetra_LinearProblem( RRMatrix_, RLHS_, RRHS_ );
352
353 newObj_ = NewProblem_;
354
355 cout << "done with SC\n";
356
357 return *NewProblem_;
358}
359
360bool
362fwd()
363{
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 );
369
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 );
374
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 );
381
382 Epetra_Vector LLDiag( *LMap_ );
383 LLMatrix_->ExtractDiagonalCopy( LLDiag );
384 Epetra_Vector LLRecipDiag( *LMap_ );
385 LLRecipDiag.Reciprocal( LLDiag );
386
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];
391
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 );
396
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 );
401
402 return true;
403}
404
405bool
407rvs()
408{
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 );
414
415 Epetra_Vector UUDiag( *UMap_ );
416 UUMatrix_->ExtractDiagonalCopy( UUDiag );
417 Epetra_Vector UURecipDiag( *UMap_ );
418 UURecipDiag.Reciprocal( UUDiag );
419
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];
424
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 );
429
430 return true;
431}
432
433} // namespace EpetraExt
434
Insert
Copy
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.
bool DistributedGlobal() const
int GID(int LID) const
int IndexBase() const
const Epetra_Comm & Comm() const
int NumMyElements() const
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
int NumMyRows() 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 NumVectors() const
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.