EpetraExt Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
58 namespace 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();
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 
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 
234  cout << "Permuted Graph:\n" << *NewGraph_;
235 
236  NewMatrix_ = new Epetra_CrsMatrix( Copy, *NewGraph_ );
237  NewMatrix_->Export( *OldMatrix_, *Exporter_, Insert );
239  cout << "Permuted Matrix:\n" << *NewMatrix_;
240 
241  UExporter_ = new Epetra_Export( *OldRowMap_, *UMap_ );
244 cout << "UExporter:\n" << *UExporter_;
245 cout << "RExporter:\n" << *RExporter_;
246 cout << "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 );
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 );
290  cout << "URMatrix:\n" << *URMatrix_;
291 
292  ULGraph_ = new Epetra_CrsGraph( Copy, *UMap_, *UMap_, 0 );
293  ULGraph_->Export( *OldGraph_, *UExporter_, Insert );
295  cout << "ULGraph:\n" << *ULGraph_;
296 
297  ULMatrix_ = new Epetra_CrsMatrix( Copy, *ULGraph_ );
298  ULMatrix_->Export( *OldMatrix_, *UExporter_, Insert );
300  cout << "ULMatrix:\n" << *ULMatrix_;
301 
302  RRGraph_ = new Epetra_CrsGraph( Copy, *RMap_, *RMap_, 0 );
303  RRGraph_->Export( *OldGraph_, *RExporter_, Insert );
305  cout << "RRGraph:\n" << *RRGraph_;
306 
307  RRMatrix_ = new Epetra_CrsMatrix( Copy, *RRGraph_ );
308  RRMatrix_->Export( *OldMatrix_, *RExporter_, Insert );
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 );
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 );
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 
354 
355  cout << "done with SC\n";
356 
357  return *NewProblem_;
358 }
359 
360 bool
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 
405 bool
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 
bool DistributedGlobal() const
int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
int NumMyRows() const
int GlobalMaxNumEntries() const
int FillComplete(bool OptimizeDataStorage=true)
int IndexBase() const
bool fwd()
Forward transfer of data from orig object input in the operator() method call to the new object creat...
const Epetra_Map & RowMap() const
int NumMyElements() const
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
NewTypeRef operator()(OriginalTypeRef orig)
Analysis of transform operation on original object and construction of new object.
int GID(int LID) const
int NumGlobalNonzeros() const
const Epetra_Comm & Comm() const
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
const Epetra_CrsGraph & Graph() const
bool rvs()
Reverse transfer of data from new object created in the operator() method call to the orig object inp...