EpetraExt Package Browser (Single Doxygen Collection)  Development
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 
105 operator()( OriginalTypeRef orig )
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_ );
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 );
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 
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 
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...
bool rvs()
Reverse transfer of data from new object created in the operator() method call to the orig object inp...