EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 {
75  deleteNewObjs_();
76 }
77 
78 void
79 LinearProblem_BTF::
80 deleteNewObjs_()
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 
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 
138  if( OrigMatrix_->RowMap().DistributedGlobal() )
139  { cout << "FAIL for Global!\n"; abort(); }
140  if( OrigMatrix_->IndicesAreGlobal() )
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
466  NewMatrix_ = new Epetra_VbrMatrix( View, *NewMap_, &BlockCnt_[0] );
467  for( int i = 0; i < sqcmpn; ++i )
468  {
469  NewMatrix_->BeginInsertGlobalValues( i, BlockCnt_[i], &(NewBlockRows_[i])[0] );
470  for( int j = 0; j < BlockCnt_[i]; ++j )
471  NewMatrix_->SubmitBlockEntry( *(Blocks_[i][j]) );
472  NewMatrix_->EndSubmitEntries();
473  }
474  NewMatrix_->FillComplete();
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
487  NewProblem_ = new Epetra_LinearProblem( NewMatrix_, NewLHS_, NewRHS_ );
488  newObj_ = NewProblem_;
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
bool DistributedGlobal() const
int MyGlobalElements(int *MyGlobalElementList) const
int InsertMyIndices(int LocalRow, int NumIndices, int *Indices)
const Epetra_Map & ColMap() 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
const Epetra_Map & RowMap() const
#define GENBTF_F77
int FirstPointInElement(int LID) const
int MaxNumEntries() const
int ExtractGlobalRowCopy(int GlobalRow, int Length, int &NumEntries, double *Values, int *Indices) const
int NumMyRows() const
int LID(int GID) const
int ExtractGlobalRowCopy(int GlobalRow, int LenOfIndices, int &NumIndices, int *Indices) const
NewTypeRef operator()(OriginalTypeRef orig)
Analysis of transform operation on original object and construction of new object.
const Epetra_Comm & Comm() const
bool IndicesAreGlobal() const
#define MATTRANS_F77
bool rvs()
Reverse transfer of data from new object created in the operator() method call to the orig object inp...
int MaxNumIndices() const
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
int NumMyNonzeros() const