EpetraExt Package Browser (Single Doxygen Collection)  Development
EpetraExt_BTF_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_CrsMatrix.h>
45 #include <Epetra_VbrMatrix.h>
46 #include <Epetra_CrsGraph.h>
47 #include <Epetra_Map.h>
48 #include <Epetra_BlockMap.h>
49 #include <Epetra_MultiVector.h>
50 #include <Epetra_LinearProblem.h>
51 
52 #include <set>
53 
54 using std::vector;
55 using std::map;
56 using std::set;
57 using std::cout;
58 using std::endl;
59 
60 #define MATTRANS_F77 F77_FUNC(mattrans,MATTRANS)
61 #define GENBTF_F77 F77_FUNC(genbtf,GENBTF)
62 
63 extern "C" {
64 extern void MATTRANS_F77( int*, int*, int*, int*, int*, int* );
65 extern void GENBTF_F77( int*, int*, int*, int*, int*, int*, int*, int*, int*,
66  int*, int*, int*, int*, int*, int*, int*, int*, int*,
67  int*, int*, int* );
68 }
69 
70 namespace EpetraExt {
71 
74 {
76 }
77 
78 void
81 {
82  if( NewProblem_ ) delete NewProblem_;
83 
84  if( NewMatrix_ ) delete NewMatrix_;
85 
86  if( NewLHS_ ) delete NewLHS_;
87  if( NewRHS_ ) delete NewRHS_;
88 
89  if( NewMap_ ) delete NewMap_;
90 
91  for( int i = 0; i < Blocks_.size(); ++i )
92  for( int j = 0; j < Blocks_[i].size(); ++j )
93  delete Blocks_[i][j];
94 }
95 
98 operator()( OriginalTypeRef orig )
99 {
100  changedLP_ = false;
101 
102  //check if there is a previous analysis and if it is valid
103  if( &orig == origObj_ && NewProblem_ )
104  {
105  int * indices;
106  double * values;
107  int currCnt;
108  int numRows = OrigMatrix_->NumMyRows();
109 
110  for( int i = 0; i < numRows; ++i )
111  if( ZeroElements_[i].size() )
112  {
113  int loc = 0;
114  OrigMatrix_->ExtractMyRowView( i, currCnt, values, indices );
115  for( int j = 0; j < currCnt; ++j )
116  if( ZeroElements_[i].count( indices[j] ) )
117  {
118  if( values[j] > threshold_ ) changedLP_ = true;
119  ++loc;
120  }
121  }
122 
123 // changedLP_ = true;
124 
125  //if changed, dump all the old stuff and start over
126  if( changedLP_ )
127  deleteNewObjs_();
128  else
129  return *newObj_;
130  }
131 
132  origObj_ = &orig;
133 
134  OrigMatrix_ = dynamic_cast<Epetra_CrsMatrix*>(orig.GetMatrix());
135  OrigLHS_ = orig.GetLHS();
136  OrigRHS_ = orig.GetRHS();
137 
139  { cout << "FAIL for Global!\n"; abort(); }
141  { cout << "FAIL for Global Indices!\n"; abort(); }
142 
143  int n = OrigMatrix_->NumMyRows();
144  int nnz = OrigMatrix_->NumMyNonzeros();
145 
146 // cout << "Orig Matrix:\n";
147 // cout << *OrigMatrix_ << endl;
148 
149  //create std CRS format
150  //also create graph without zero elements
151  vector<int> ia(n+1,0);
152  int maxEntries = OrigMatrix_->MaxNumEntries();
153  vector<int> ja_tmp(maxEntries);
154  vector<double> jva_tmp(maxEntries);
155  vector<int> ja(nnz);
156  int cnt;
157 
158  const Epetra_BlockMap & OldRowMap = OrigMatrix_->RowMap();
159  const Epetra_BlockMap & OldColMap = OrigMatrix_->ColMap();
160  Epetra_CrsGraph strippedGraph( Copy, OldRowMap, OldColMap, 0 );
161  ZeroElements_.resize(n);
162 
163  for( int i = 0; i < n; ++i )
164  {
165  ZeroElements_[i].clear();
166  OrigMatrix_->ExtractMyRowCopy( i, maxEntries, cnt, &jva_tmp[0], &ja_tmp[0] );
167  ia[i+1] = ia[i];
168  for( int j = 0; j < cnt; ++j )
169  {
170  if( fabs(jva_tmp[j]) > threshold_ )
171  ja[ ia[i+1]++ ] = ja_tmp[j];
172  else
173  ZeroElements_[i].insert( ja_tmp[j] );
174  }
175 
176  int new_cnt = ia[i+1] - ia[i];
177  strippedGraph.InsertMyIndices( i, new_cnt, &ja[ ia[i] ] );
178  }
179  nnz = ia[n];
180  strippedGraph.FillComplete();
181 
182  if( verbose_ > 2 )
183  {
184  cout << "Stripped Graph\n";
185  cout << strippedGraph;
186  }
187 
188  vector<int> iat(n+1,0);
189  vector<int> jat(nnz);
190  for( int i = 0; i < n; ++i )
191  for( int j = ia[i]; j < ia[i+1]; ++j )
192  ++iat[ ja[j]+1 ];
193  for( int i = 0; i < n; ++i )
194  iat[i+1] += iat[i];
195  for( int i = 0; i < n; ++i )
196  for( int j = ia[i]; j < ia[i+1]; ++j )
197  jat[ iat[ ja[j] ]++ ] = i;
198  for( int i = 0; i < n; ++i )
199  iat[n-i] = iat[n-i-1];
200  iat[0] = 0;
201 
202  //convert to Fortran indexing
203  for( int i = 0; i < n+1; ++i ) ++ia[i];
204  for( int i = 0; i < nnz; ++i ) ++ja[i];
205  for( int i = 0; i < n+1; ++i ) ++iat[i];
206  for( int i = 0; i < nnz; ++i ) ++jat[i];
207 
208  if( verbose_ > 2 )
209  {
210  cout << "-----------------------------------------\n";
211  cout << "CRS Format Graph\n";
212  cout << "-----------------------------------------\n";
213  for( int i = 0; i < n; ++i )
214  {
215  cout << i+1 << ": " << ia[i+1] << ": ";
216  for( int j = ia[i]-1; j<ia[i+1]-1; ++j )
217  cout << " " << ja[j];
218  cout << endl;
219  }
220  cout << "-----------------------------------------\n";
221  }
222 
223  if( verbose_ > 2 )
224  {
225  cout << "-----------------------------------------\n";
226  cout << "CCS Format Graph\n";
227  cout << "-----------------------------------------\n";
228  for( int i = 0; i < n; ++i )
229  {
230  cout << i+1 << ": " << iat[i+1] << ": ";
231  for( int j = iat[i]-1; j<iat[i+1]-1; ++j )
232  cout << " " << jat[j];
233  cout << endl;
234  }
235  cout << "-----------------------------------------\n";
236  }
237 
238  vector<int> w(10*n);
239 
240  vector<int> rowperm(n);
241  vector<int> colperm(n);
242 
243  //horizontal block
244  int nhrows, nhcols, hrzcmp;
245  //square block
246  int nsrows, sqcmpn;
247  //vertial block
248  int nvrows, nvcols, vrtcmp;
249 
250  vector<int> rcmstr(n+1);
251  vector<int> ccmstr(n+1);
252 
253  int msglvl = 0;
254  int output = 6;
255 
256  GENBTF_F77( &n, &n, &iat[0], &jat[0], &ia[0], &ja[0], &w[0],
257  &rowperm[0], &colperm[0], &nhrows, &nhcols,
258  &hrzcmp, &nsrows, &sqcmpn, &nvrows, &nvcols, &vrtcmp,
259  &rcmstr[0], &ccmstr[0], &msglvl, &output );
260 
261  //convert back to C indexing
262  for( int i = 0; i < n; ++i )
263  {
264  --rowperm[i];
265  --colperm[i];
266  }
267  for( int i = 0; (i<n+1) && (rcmstr[i]!=n+1); ++i )
268  {
269  --rcmstr[i];
270  --ccmstr[i];
271  }
272 
273  if( verbose_ > 0 )
274  {
275  cout << "-----------------------------------------\n";
276  cout << "BTF Output\n";
277  cout << "-----------------------------------------\n";
278 // cout << "RowPerm and ColPerm\n";
279 // for( int i = 0; i<n; ++i )
280 // cout << rowperm[i] << "\t" << colperm[i] << endl;
281  if( hrzcmp )
282  {
283  cout << "Num Horizontal: Rows, Cols, Comps\n";
284  cout << nhrows << "\t" << nhcols << "\t" << hrzcmp << endl;
285  }
286  cout << "Num Square: Rows, Comps\n";
287  cout << nsrows << "\t" << sqcmpn << endl;
288  if( vrtcmp )
289  {
290  cout << "Num Vertical: Rows, Cols, Comps\n";
291  cout << nvrows << "\t" << nvcols << "\t" << vrtcmp << endl;
292  }
293 // cout << "Row, Col of upper left pt in blocks\n";
294 // for( int i = 0; (i<n+1) && (rcmstr[i]!=n+1); ++i )
295 // cout << i << " " << rcmstr[i] << " " << ccmstr[i] << endl;
296 // cout << "-----------------------------------------\n";
297  }
298 
299  if( hrzcmp || vrtcmp )
300  {
301  cout << "FAILED! hrz cmp's:" << hrzcmp << " vrtcmp: " << vrtcmp << endl;
302  exit(0);
303  }
304 
305  rcmstr[sqcmpn] = n;
306 
307  //convert rowperm to OLD->NEW
308  //reverse ordering of permutation to get upper triangular
309  vector<int> rowperm_t( n );
310  vector<int> colperm_t( n );
311  for( int i = 0; i < n; ++i )
312  {
313  rowperm_t[i] = rowperm[i];
314  colperm_t[i] = colperm[i];
315  }
316 
317  //Generate New Domain and Range Maps
318  //for now, assume they start out as identical
319  OldGlobalElements_.resize(n);
320  OldRowMap.MyGlobalElements( &OldGlobalElements_[0] );
321 
322  vector<int> newDomainElements( n );
323  vector<int> newRangeElements( n );
324  for( int i = 0; i < n; ++i )
325  {
326  newDomainElements[ i ] = OldGlobalElements_[ rowperm_t[i] ];
327  newRangeElements[ i ] = OldGlobalElements_[ colperm_t[i] ];
328  }
329 
330  //Setup New Block Map Info
331  Blocks_.resize( sqcmpn );
332  BlockDim_.resize( sqcmpn );
333  for( int i = 0; i < sqcmpn; ++i )
334  {
335  BlockDim_[i] = rcmstr[i+1]-rcmstr[i];
336  for( int j = rcmstr[i]; j < rcmstr[i+1]; ++j )
337  {
338  BlockRowMap_[ newDomainElements[j] ] = i;
339  SubBlockRowMap_[ newDomainElements[j] ] = j-rcmstr[i];
340  BlockColMap_[ newRangeElements[j] ] = i;
341  SubBlockColMap_[ newRangeElements[j] ] = j-rcmstr[i];
342  }
343  }
344 
345  if( verbose_ > 2 )
346  {
347 /*
348  cout << "Block Mapping!\n";
349  cout << "--------------------------\n";
350  for( int i = 0; i < n; ++i )
351  {
352  cout << "Row: " << newDomainElements[i] << " " << BlockRowMap_[newDomainElements[i]] << " " <<
353  SubBlockRowMap_[newDomainElements[i]] << "\t" << "Col: " << newRangeElements[i] << " " <<
354  BlockColMap_[newRangeElements[i]] << " " << SubBlockColMap_[newRangeElements[i]] << endl;
355  }
356  for( int i = 0; i < sqcmpn; ++i )
357  cout << "BlockDim: " << i << " " << BlockDim_[i] << endl;
358  cout << "--------------------------\n";
359 */
360  int MinSize = 1000000000;
361  int MaxSize = 0;
362  for( int i = 0; i < sqcmpn; ++i )
363  {
364  if( MinSize > BlockDim_[i] ) MinSize = BlockDim_[i];
365  if( MaxSize < BlockDim_[i] ) MaxSize = BlockDim_[i];
366  }
367  cout << "Min Block Size: " << MinSize << " " << "Max Block Size: " << MaxSize << endl;
368  }
369 
370  vector<int> myBlockElements( sqcmpn );
371  for( int i = 0; i < sqcmpn; ++i ) myBlockElements[i] = i;
372  NewMap_ = new Epetra_BlockMap( sqcmpn, sqcmpn, &myBlockElements[0], &BlockDim_[0], 0, OldRowMap.Comm() );
373 
374  if( verbose_ > 2 )
375  {
376  cout << "New Block Map!\n";
377  cout << *NewMap_;
378  }
379 
380  //setup new graph
381  vector< set<int> > crsBlocks( sqcmpn );
382  BlockCnt_.resize( sqcmpn );
383  int maxLength = strippedGraph.MaxNumIndices();
384  vector<int> sIndices( maxLength );
385  int currLength;
386  for( int i = 0; i < n; ++i )
387  {
388  strippedGraph.ExtractGlobalRowCopy( OldGlobalElements_[i], maxLength, currLength, &sIndices[0] );
389  for( int j = 0; j < currLength; ++j )
390  crsBlocks[ BlockRowMap_[ OldGlobalElements_[i] ] ].insert( BlockColMap_[ sIndices[j] ] );
391  }
392 
393  for( int i = 0; i < sqcmpn; ++i )
394  {
395  BlockCnt_[i] = crsBlocks[i].size();
396  Blocks_[i].resize( BlockCnt_[i] );
397  }
398 
399  NewBlockRows_.resize( sqcmpn );
400  for( int i = 0; i < sqcmpn; ++i )
401  {
402  NewBlockRows_[i] = vector<int>( crsBlocks[i].begin(), crsBlocks[i].end() );
403  for( int j = 0; j < BlockCnt_[i]; ++j )
404  {
405  Blocks_[i][j] = new Epetra_SerialDenseMatrix();
406  Blocks_[i][j]->Shape( BlockDim_[i], BlockDim_[ NewBlockRows_[i][j] ] );
407  }
408  }
409 
410  //put data in Blocks_ and new LHS and RHS
411  NewLHS_ = new Epetra_MultiVector( *NewMap_, 1 );
412  NewRHS_ = new Epetra_MultiVector( *NewMap_, 1 );
413 
414  maxLength = OrigMatrix_->MaxNumEntries();
415  vector<int> indices( maxLength );
416  vector<double> values( maxLength );
417  for( int i = 0; i < n; ++i )
418  {
419  int BlockRow = BlockRowMap_[ OldGlobalElements_[i] ];
420  int SubBlockRow = SubBlockRowMap_[ OldGlobalElements_[i] ];
421  OrigMatrix_->ExtractGlobalRowCopy( OldGlobalElements_[i], maxLength, currLength, &values[0], &indices[0] );
422  for( int j = 0; j < currLength; ++j )
423  {
424  int BlockCol = BlockColMap_[ indices[j] ];
425  int SubBlockCol = SubBlockColMap_[ indices[j] ];
426  for( int k = 0; k < BlockCnt_[BlockRow]; ++k )
427  if( BlockCol == NewBlockRows_[BlockRow][k] )
428  {
429  if( values[j] > threshold_ )
430  {
431 // cout << BlockRow << " " << SubBlockRow << " " << BlockCol << " " << SubBlockCol << endl;
432 // cout << k << endl;
433  (*(Blocks_[BlockRow][k]))(SubBlockRow,SubBlockCol) = values[j];
434  break;
435  }
436  else
437  ZeroElements_[i].erase( OrigMatrix_->RowMap().LID( indices[j] ) );
438  }
439  }
440 
441 // NewLHS_->ReplaceGlobalValue( BlockCol, SubBlockCol, 0, (*OrigLHS_)[0][i] );
442  NewRHS_->ReplaceGlobalValue( BlockRow, SubBlockRow, 0, (*OrigRHS_)[0][i] );
443  }
444 
445  if( verbose_ > 2 )
446  {
447  cout << "Zero Elements: \n";
448  cout << "--------------\n";
449  int cnt = 0;
450  for( int i = 0; i < n; ++i )
451  {
452  set<int>::iterator iterSI = ZeroElements_[i].begin();
453  set<int>::iterator endSI = ZeroElements_[i].end();
454  for( ; iterSI != endSI; ++iterSI )
455  {
456  cout << " " << *iterSI;
457  ++cnt;
458  }
459  cout << endl;
460  }
461  cout << "ZE Cnt: " << cnt << endl;
462  cout << "--------------\n";
463  }
464 
465  //setup new matrix
467  for( int i = 0; i < sqcmpn; ++i )
468  {
470  for( int j = 0; j < BlockCnt_[i]; ++j )
471  NewMatrix_->SubmitBlockEntry( *(Blocks_[i][j]) );
473  }
475 
476  if( verbose_ > 2 )
477  {
478  cout << "New Block Matrix!\n";
479  cout << *NewMatrix_;
480  cout << "New Block LHS!\n";
481  cout << *NewLHS_;
482  cout << "New Block RHS!\n";
483  cout << *NewRHS_;
484  }
485 
486  //create new LP
489 
490  if( verbose_ ) cout << "-----------------------------------------\n";
491 
492  return *newObj_;
493 }
494 
495 bool
498 {
499  //zero out matrix
500  int NumBlockRows = BlockDim_.size();
501  for( int i = 0; i < NumBlockRows; ++i )
502  {
503  int NumBlocks = BlockCnt_[i];
504  for( int j = 0; j < NumBlocks; ++j )
505  {
506  int Size = BlockDim_[i] * BlockDim_[ NewBlockRows_[i][j] ];
507  double * p = Blocks_[i][j]->A();
508  for( int k = 0; k < Size; ++k ) *p++ = 0.0;
509  }
510  }
511 
512  int maxLength = OrigMatrix_->MaxNumEntries();
513  int n = OldGlobalElements_.size();
514  int currLength;
515  vector<int> indices( maxLength );
516  vector<double> values( maxLength );
517  for( int i = 0; i < n; ++i )
518  {
519  int BlockRow = BlockRowMap_[ OldGlobalElements_[i] ];
520  int SubBlockRow = SubBlockRowMap_[ OldGlobalElements_[i] ];
521  OrigMatrix_->ExtractGlobalRowCopy( OldGlobalElements_[i], maxLength, currLength, &values[0], &indices[0] );
522  for( int j = 0; j < currLength; ++j )
523  {
524  if( fabs(values[j]) > threshold_ )
525  {
526  int BlockCol = BlockColMap_[ indices[j] ];
527  int SubBlockCol = SubBlockColMap_[ indices[j] ];
528  for( int k = 0; k < BlockCnt_[BlockRow]; ++k )
529  if( BlockCol == NewBlockRows_[BlockRow][k] )
530  {
531  (*(Blocks_[BlockRow][k]))(SubBlockRow,SubBlockCol) = values[j];
532  break;
533  }
534  }
535  }
536 
537  NewRHS_->ReplaceGlobalValue( BlockRow, SubBlockRow, 0, (*OrigRHS_)[0][i] );
538  }
539 
540 /*
541  //fill matrix
542  int sqcmpn = BlockDim_.size();
543  for( int i = 0; i < sqcmpn; ++i )
544  {
545  NewMatrix_->BeginReplaceGlobalValues( i, NewBlockRows_[i].size(), &(NewBlockRows_[i])[0] );
546  for( int j = 0; j < NewBlockRows_[i].size(); ++j )
547  NewMatrix_->SubmitBlockEntry( Blocks_[i][j]->A(), Blocks_[i][j]->LDA(), Blocks_[i][j]->M(), Blocks_[i][j]->N() );
548  NewMatrix_->EndSubmitEntries();
549  }
550 */
551 
552  return true;
553 }
554 
555 bool
558 {
559  //copy data from NewLHS_ to OldLHS_;
560  int rowCnt = OrigLHS_->MyLength();
561  for( int i = 0; i < rowCnt; ++i )
562  {
563  int BlockCol = BlockColMap_[ OldGlobalElements_[i] ];
564  int SubBlockCol = SubBlockColMap_[ OldGlobalElements_[i] ];
565  (*OrigLHS_)[0][i] = (*NewLHS_)[0][ NewMap_->FirstPointInElement(BlockCol) + SubBlockCol ];
566  }
567 
568  return true;
569 }
570 
571 } //namespace EpetraExt
const Epetra_Map & RowMap() const
View
int MaxNumEntries() const
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
int MyGlobalElements(int *MyGlobalElementList) const
int MyLength() const
Copy
bool DistributedGlobal() const
int NumMyRows() const
int InsertMyIndices(int LocalRow, int NumIndices, int *Indices)
int ExtractGlobalRowCopy(int_type Row, int Length, int &NumEntries, double *values, int_type *Indices) const
std::vector< std::vector< int > > NewBlockRows_
int FirstPointInElement(int LID) const
bool fwd()
Forward transfer of data from orig object input in the operator() method call to the new object creat...
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
std::vector< std::set< int > > ZeroElements_
std::vector< std::vector< Epetra_SerialDenseMatrix * > > Blocks_
int ReplaceGlobalValue(int GlobalRow, int VectorIndex, double ScalarValue)
const Epetra_Comm & Comm() const
#define GENBTF_F77
NewTypeRef operator()(OriginalTypeRef orig)
const Epetra_Map & ColMap() const
int LID(int GID) const
int NumMyNonzeros() const
bool rvs()
Reverse transfer of data from new object created in the operator() method call to the orig object inp...
bool IndicesAreGlobal() const
#define MATTRANS_F77
int MaxNumIndices() const
int BeginInsertGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
int SubmitBlockEntry(double *Values, int LDA, int NumRows, int NumCols)
int ExtractGlobalRowCopy(int_type Row, int LenOfIndices, int &NumIndices, int_type *Indices) const