EpetraExt Package Browser (Single Doxygen Collection) 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>
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() );
113 OldRHS_ = orig.GetRHS();
114 OldLHS_ = orig.GetLHS();
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
224
227
230
234 cout << "Permuted Graph:\n" << *NewGraph_;
235
239 cout << "Permuted Matrix:\n" << *NewMatrix_;
240
244cout << "UExporter:\n" << *UExporter_;
245cout << "RExporter:\n" << *RExporter_;
246cout << "LExporter:\n" << *LExporter_;
247
250 cout << "ULHS:\n" << *ULHS_;
251
254 cout << "RLHS:\n" << *RLHS_;
255
258 cout << "LLHS:\n" << *LLHS_;
259
262 cout << "URHS:\n" << *URHS_;
263
266 cout << "RRHS:\n" << *RRHS_;
267
270 cout << "LRHS:\n" << *LRHS_;
271
272 UUGraph_ = new Epetra_CrsGraph( Copy, *UMap_, *LMap_, 0 );
275 cout << "UUGraph:\n" << *UUGraph_;
276
280 cout << "UUMatrix:\n" << *UUMatrix_;
281
282 URGraph_ = new Epetra_CrsGraph( Copy, *UMap_, *RMap_, 0 );
285 cout << "URGraph:\n" << *URGraph_;
286
290 cout << "URMatrix:\n" << *URMatrix_;
291
292 ULGraph_ = new Epetra_CrsGraph( Copy, *UMap_, *UMap_, 0 );
295 cout << "ULGraph:\n" << *ULGraph_;
296
300 cout << "ULMatrix:\n" << *ULMatrix_;
301
302 RRGraph_ = new Epetra_CrsGraph( Copy, *RMap_, *RMap_, 0 );
305 cout << "RRGraph:\n" << *RRGraph_;
306
310 cout << "RRMatrix:\n" << *RRMatrix_;
311
312 RLGraph_ = new Epetra_CrsGraph( Copy, *RMap_, *UMap_, 0 );
315 cout << "RLGraph:\n" << *RLGraph_;
316
320 cout << "RLMatrix:\n" << *RLMatrix_;
321
322 LLGraph_ = new Epetra_CrsGraph( Copy, *LMap_, *UMap_, 0 );
325 cout << "LLGraph:\n" << *LLGraph_;
326
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
352
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";
369
370 if( verbose_ ) cout << "LP_SC::fwd() : Exporting to New RHS\n";
374
381
382 Epetra_Vector LLDiag( *LMap_ );
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_ );
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";
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.