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> 60 #define MATTRANS_F77 F77_FUNC(mattrans,MATTRANS) 61 #define GENBTF_F77 F77_FUNC(genbtf,GENBTF) 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*,
82 if( NewProblem_ )
delete NewProblem_;
84 if( NewMatrix_ )
delete NewMatrix_;
86 if( NewLHS_ )
delete NewLHS_;
87 if( NewRHS_ )
delete NewRHS_;
89 if( NewMap_ )
delete NewMap_;
91 for(
int i = 0; i < Blocks_.size(); ++i )
92 for(
int j = 0; j < Blocks_[i].size(); ++j )
103 if( &orig ==
origObj_ && NewProblem_ )
108 int numRows = OrigMatrix_->NumMyRows();
110 for(
int i = 0; i < numRows; ++i )
111 if( ZeroElements_[i].size() )
114 OrigMatrix_->ExtractMyRowView( i, currCnt, values, indices );
115 for(
int j = 0; j < currCnt; ++j )
116 if( ZeroElements_[i].count( indices[j] ) )
118 if( values[j] > threshold_ ) changedLP_ =
true;
134 OrigMatrix_ =
dynamic_cast<Epetra_CrsMatrix*
>(orig.GetMatrix());
135 OrigLHS_ = orig.GetLHS();
136 OrigRHS_ = orig.GetRHS();
138 if( OrigMatrix_->RowMap().DistributedGlobal() )
139 { cout <<
"FAIL for Global!\n"; abort(); }
140 if( OrigMatrix_->IndicesAreGlobal() )
141 { cout <<
"FAIL for Global Indices!\n"; abort(); }
143 int n = OrigMatrix_->NumMyRows();
144 int nnz = OrigMatrix_->NumMyNonzeros();
151 vector<int> ia(n+1,0);
152 int maxEntries = OrigMatrix_->MaxNumEntries();
153 vector<int> ja_tmp(maxEntries);
154 vector<double> jva_tmp(maxEntries);
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);
163 for(
int i = 0; i < n; ++i )
165 ZeroElements_[i].clear();
166 OrigMatrix_->ExtractMyRowCopy( i, maxEntries, cnt, &jva_tmp[0], &ja_tmp[0] );
168 for(
int j = 0; j < cnt; ++j )
170 if( fabs(jva_tmp[j]) > threshold_ )
171 ja[ ia[i+1]++ ] = ja_tmp[j];
173 ZeroElements_[i].insert( ja_tmp[j] );
176 int new_cnt = ia[i+1] - ia[i];
177 strippedGraph.InsertMyIndices( i, new_cnt, &ja[ ia[i] ] );
180 strippedGraph.FillComplete();
184 cout <<
"Stripped Graph\n";
185 cout << strippedGraph;
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 )
193 for(
int i = 0; i < n; ++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];
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];
210 cout <<
"-----------------------------------------\n";
211 cout <<
"CRS Format Graph\n";
212 cout <<
"-----------------------------------------\n";
213 for(
int i = 0; i < n; ++i )
215 cout << i+1 <<
": " << ia[i+1] <<
": ";
216 for(
int j = ia[i]-1; j<ia[i+1]-1; ++j )
217 cout <<
" " << ja[j];
220 cout <<
"-----------------------------------------\n";
225 cout <<
"-----------------------------------------\n";
226 cout <<
"CCS Format Graph\n";
227 cout <<
"-----------------------------------------\n";
228 for(
int i = 0; i < n; ++i )
230 cout << i+1 <<
": " << iat[i+1] <<
": ";
231 for(
int j = iat[i]-1; j<iat[i+1]-1; ++j )
232 cout <<
" " << jat[j];
235 cout <<
"-----------------------------------------\n";
240 vector<int> rowperm(n);
241 vector<int> colperm(n);
244 int nhrows, nhcols, hrzcmp;
248 int nvrows, nvcols, vrtcmp;
250 vector<int> rcmstr(n+1);
251 vector<int> ccmstr(n+1);
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 );
262 for(
int i = 0; i < n; ++i )
267 for(
int i = 0; (i<n+1) && (rcmstr[i]!=n+1); ++i )
275 cout <<
"-----------------------------------------\n";
276 cout <<
"BTF Output\n";
277 cout <<
"-----------------------------------------\n";
283 cout <<
"Num Horizontal: Rows, Cols, Comps\n";
284 cout << nhrows <<
"\t" << nhcols <<
"\t" << hrzcmp << endl;
286 cout <<
"Num Square: Rows, Comps\n";
287 cout << nsrows <<
"\t" << sqcmpn << endl;
290 cout <<
"Num Vertical: Rows, Cols, Comps\n";
291 cout << nvrows <<
"\t" << nvcols <<
"\t" << vrtcmp << endl;
299 if( hrzcmp || vrtcmp )
301 cout <<
"FAILED! hrz cmp's:" << hrzcmp <<
" vrtcmp: " << vrtcmp << endl;
309 vector<int> rowperm_t( n );
310 vector<int> colperm_t( n );
311 for(
int i = 0; i < n; ++i )
313 rowperm_t[i] = rowperm[i];
314 colperm_t[i] = colperm[i];
319 OldGlobalElements_.resize(n);
320 OldRowMap.MyGlobalElements( &OldGlobalElements_[0] );
322 vector<int> newDomainElements( n );
323 vector<int> newRangeElements( n );
324 for(
int i = 0; i < n; ++i )
326 newDomainElements[ i ] = OldGlobalElements_[ rowperm_t[i] ];
327 newRangeElements[ i ] = OldGlobalElements_[ colperm_t[i] ];
331 Blocks_.resize( sqcmpn );
332 BlockDim_.resize( sqcmpn );
333 for(
int i = 0; i < sqcmpn; ++i )
335 BlockDim_[i] = rcmstr[i+1]-rcmstr[i];
336 for(
int j = rcmstr[i]; j < rcmstr[i+1]; ++j )
338 BlockRowMap_[ newDomainElements[j] ] = i;
339 SubBlockRowMap_[ newDomainElements[j] ] = j-rcmstr[i];
340 BlockColMap_[ newRangeElements[j] ] = i;
341 SubBlockColMap_[ newRangeElements[j] ] = j-rcmstr[i];
360 int MinSize = 1000000000;
362 for(
int i = 0; i < sqcmpn; ++i )
364 if( MinSize > BlockDim_[i] ) MinSize = BlockDim_[i];
365 if( MaxSize < BlockDim_[i] ) MaxSize = BlockDim_[i];
367 cout <<
"Min Block Size: " << MinSize <<
" " <<
"Max Block Size: " << MaxSize << endl;
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() );
376 cout <<
"New Block Map!\n";
381 vector< set<int> > crsBlocks( sqcmpn );
382 BlockCnt_.resize( sqcmpn );
383 int maxLength = strippedGraph.MaxNumIndices();
384 vector<int> sIndices( maxLength );
386 for(
int i = 0; i < n; ++i )
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] ] );
393 for(
int i = 0; i < sqcmpn; ++i )
395 BlockCnt_[i] = crsBlocks[i].size();
396 Blocks_[i].resize( BlockCnt_[i] );
399 NewBlockRows_.resize( sqcmpn );
400 for(
int i = 0; i < sqcmpn; ++i )
402 NewBlockRows_[i] = vector<int>( crsBlocks[i].begin(), crsBlocks[i].end() );
403 for(
int j = 0; j < BlockCnt_[i]; ++j )
405 Blocks_[i][j] =
new Epetra_SerialDenseMatrix();
406 Blocks_[i][j]->Shape( BlockDim_[i], BlockDim_[ NewBlockRows_[i][j] ] );
411 NewLHS_ =
new Epetra_MultiVector( *NewMap_, 1 );
412 NewRHS_ =
new Epetra_MultiVector( *NewMap_, 1 );
414 maxLength = OrigMatrix_->MaxNumEntries();
415 vector<int> indices( maxLength );
416 vector<double> values( maxLength );
417 for(
int i = 0; i < n; ++i )
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 )
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] )
429 if( values[j] > threshold_ )
433 (*(Blocks_[BlockRow][k]))(SubBlockRow,SubBlockCol) = values[j];
437 ZeroElements_[i].erase( OrigMatrix_->RowMap().LID( indices[j] ) );
442 NewRHS_->ReplaceGlobalValue( BlockRow, SubBlockRow, 0, (*OrigRHS_)[0][i] );
447 cout <<
"Zero Elements: \n";
448 cout <<
"--------------\n";
450 for(
int i = 0; i < n; ++i )
452 set<int>::iterator iterSI = ZeroElements_[i].begin();
453 set<int>::iterator endSI = ZeroElements_[i].end();
454 for( ; iterSI != endSI; ++iterSI )
456 cout <<
" " << *iterSI;
461 cout <<
"ZE Cnt: " << cnt << endl;
462 cout <<
"--------------\n";
466 NewMatrix_ =
new Epetra_VbrMatrix( View, *NewMap_, &BlockCnt_[0] );
467 for(
int i = 0; i < sqcmpn; ++i )
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();
474 NewMatrix_->FillComplete();
478 cout <<
"New Block Matrix!\n";
480 cout <<
"New Block LHS!\n";
482 cout <<
"New Block RHS!\n";
487 NewProblem_ =
new Epetra_LinearProblem( NewMatrix_, NewLHS_, NewRHS_ );
490 if( verbose_ ) cout <<
"-----------------------------------------\n";
500 int NumBlockRows = BlockDim_.size();
501 for(
int i = 0; i < NumBlockRows; ++i )
503 int NumBlocks = BlockCnt_[i];
504 for(
int j = 0; j < NumBlocks; ++j )
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;
512 int maxLength = OrigMatrix_->MaxNumEntries();
513 int n = OldGlobalElements_.size();
515 vector<int> indices( maxLength );
516 vector<double> values( maxLength );
517 for(
int i = 0; i < n; ++i )
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 )
524 if( fabs(values[j]) > threshold_ )
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] )
531 (*(Blocks_[BlockRow][k]))(SubBlockRow,SubBlockCol) = values[j];
537 NewRHS_->ReplaceGlobalValue( BlockRow, SubBlockRow, 0, (*OrigRHS_)[0][i] );
560 int rowCnt = OrigLHS_->MyLength();
561 for(
int i = 0; i < rowCnt; ++i )
563 int BlockCol = BlockColMap_[ OldGlobalElements_[i] ];
564 int SubBlockCol = SubBlockColMap_[ OldGlobalElements_[i] ];
565 (*OrigLHS_)[0][i] = (*NewLHS_)[0][ NewMap_->FirstPointInElement(BlockCol) + SubBlockCol ];
bool fwd()
Forward transfer of data from orig object input in the operator() method call to the new object creat...
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
NewTypeRef operator()(OriginalTypeRef orig)
bool rvs()
Reverse transfer of data from new object created in the operator() method call to the orig object inp...