EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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>
51
52#include <set>
53
54using std::vector;
55using std::map;
56using std::set;
57using std::cout;
58using std::endl;
59
60#define MATTRANS_F77 F77_FUNC(mattrans,MATTRANS)
61#define GENBTF_F77 F77_FUNC(genbtf,GENBTF)
62
63extern "C" {
64extern void MATTRANS_F77( int*, int*, int*, int*, int*, int* );
65extern 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
70namespace EpetraExt {
71
77
78void
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_ )
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(); }
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
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 )
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
495bool
497fwd()
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
555bool
557rvs()
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
#define GENBTF_F77
#define MATTRANS_F77
View
Copy
std::vector< std::vector< Epetra_SerialDenseMatrix * > > Blocks_
NewTypeRef operator()(OriginalTypeRef orig)
Analysis of transform operation on original object and construction of new object.
std::vector< std::vector< int > > NewBlockRows_
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...
std::vector< std::set< int > > ZeroElements_
int FirstPointInElement(int LID) const
const Epetra_Map & RowMap() const
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
int NumMyNonzeros() const
int ExtractGlobalRowCopy(int_type Row, int Length, int &NumEntries, double *values, int_type *Indices) const
const Epetra_Map & ColMap() const
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
int MaxNumEntries() const
int NumMyRows() const
bool IndicesAreGlobal() const
int MyLength() const
int ReplaceGlobalValue(int GlobalRow, int VectorIndex, double ScalarValue)
int SubmitBlockEntry(double *Values, int LDA, int NumRows, int NumCols)
int BeginInsertGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.