45 #include <Epetra_Import.h> 46 #include <Epetra_CrsMatrix.h> 47 #include <Epetra_CrsGraph.h> 48 #include <Epetra_Map.h> 56 #define MATTRANS_F77 F77_FUNC(mattrans,MATTRANS) 57 #define GENBTF_F77 F77_FUNC(genbtf,GENBTF) 60 extern void MATTRANS_F77(
int*,
int*,
int*,
int*,
int*,
int* );
61 extern void GENBTF_F77(
int*,
int*,
int*,
int*,
int*,
int*,
int*,
int*,
int*,
62 int*,
int*,
int*,
int*,
int*,
int*,
int*,
int*,
int*,
71 if( NewRowMap_ )
delete NewRowMap_;
72 if( NewColMap_ )
delete NewColMap_;
74 if( Importer_ )
delete Importer_;
76 if( NewMatrix_ )
delete NewMatrix_;
77 if( NewGraph_ )
delete NewGraph_;
86 if( orig.RowMap().DistributedGlobal() )
87 { cout <<
"FAIL for Global!\n"; abort(); }
88 if( orig.IndicesAreGlobal() )
89 { cout <<
"FAIL for Global Indices!\n"; abort(); }
91 int n = orig.NumMyRows();
92 int nnz = orig.NumMyNonzeros();
96 cout <<
"Orig Matrix:\n";
102 vector<int> ia(n+1,0);
103 int maxEntries = orig.MaxNumEntries();
104 vector<int> ja_tmp(maxEntries);
105 vector<double> jva_tmp(maxEntries);
109 const Epetra_BlockMap & OldRowMap = orig.RowMap();
110 const Epetra_BlockMap & OldColMap = orig.ColMap();
111 Epetra_CrsGraph strippedGraph( Copy, OldRowMap, OldColMap, 0 );
113 for(
int i = 0; i < n; ++i )
115 orig.ExtractMyRowCopy( i, maxEntries, cnt, &jva_tmp[0], &ja_tmp[0] );
117 for(
int j = 0; j < cnt; ++j )
118 if( fabs(jva_tmp[j]) > threshold_ )
119 ja[ ia[i+1]++ ] = ja_tmp[j];
121 int new_cnt = ia[i+1] - ia[i];
122 strippedGraph.InsertMyIndices( i, new_cnt, &ja[ ia[i] ] );
125 strippedGraph.FillComplete();
129 cout <<
"Stripped Graph\n";
130 cout << strippedGraph;
133 vector<int> iat(n+1,0);
134 vector<int> jat(nnz);
135 for(
int i = 0; i < n; ++i )
136 for(
int j = ia[i]; j < ia[i+1]; ++j )
138 for(
int i = 0; i < n; ++i )
140 for(
int i = 0; i < n; ++i )
141 for(
int j = ia[i]; j < ia[i+1]; ++j )
142 jat[ iat[ ja[j] ]++ ] = i;
143 for(
int i = 0; i < n; ++i )
144 iat[n-i] = iat[n-i-1];
148 for(
int i = 0; i < n+1; ++i ) ++ia[i];
149 for(
int i = 0; i < nnz; ++i ) ++ja[i];
150 for(
int i = 0; i < n+1; ++i ) ++iat[i];
151 for(
int i = 0; i < nnz; ++i ) ++jat[i];
155 cout <<
"-----------------------------------------\n";
156 cout <<
"CRS Format Graph\n";
157 cout <<
"-----------------------------------------\n";
158 for(
int i = 0; i < n; ++i )
160 cout << i+1 <<
": " << ia[i+1] <<
": ";
161 for(
int j = ia[i]-1; j<ia[i+1]-1; ++j )
162 cout <<
" " << ja[j];
165 cout <<
"-----------------------------------------\n";
180 cout <<
"-----------------------------------------\n";
181 cout <<
"CCS Format Graph\n";
182 cout <<
"-----------------------------------------\n";
183 for(
int i = 0; i < n; ++i )
185 cout << i+1 <<
": " << iat[i+1] <<
": ";
186 for(
int j = iat[i]-1; j<iat[i+1]-1; ++j )
187 cout <<
" " << jat[j];
190 cout <<
"-----------------------------------------\n";
195 vector<int> rowperm(n);
196 vector<int> colperm(n);
199 int nhrows, nhcols, hrzcmp;
203 int nvrows, nvcols, vrtcmp;
205 vector<int> rcmstr(n+1);
206 vector<int> ccmstr(n+1);
211 GENBTF_F77( &n, &n, &iat[0], &jat[0], &ia[0], &ja[0], &w[0],
212 &rowperm[0], &colperm[0], &nhrows, &nhcols,
213 &hrzcmp, &nsrows, &sqcmpn, &nvrows, &nvcols, &vrtcmp,
214 &rcmstr[0], &ccmstr[0], &msglvl, &output );
217 for(
int i = 0; i < n; ++i )
222 for(
int i = 0; (i<n+1) && (rcmstr[i]!=n+1); ++i )
230 cout <<
"-----------------------------------------\n";
231 cout <<
"BTF Output\n";
232 cout <<
"-----------------------------------------\n";
233 cout <<
"RowPerm and ColPerm\n";
234 for(
int i = 0; i<n; ++i )
235 cout << rowperm[i] <<
"\t" << colperm[i] << endl;
238 cout <<
"Num Horizontal: Rows, Cols, Comps\n";
239 cout << nhrows <<
"\t" << nhcols <<
"\t" << hrzcmp << endl;
241 cout <<
"Num Square: Rows, Comps\n";
242 cout << nsrows <<
"\t" << sqcmpn << endl;
245 cout <<
"Num Vertical: Rows, Cols, Comps\n";
246 cout << nvrows <<
"\t" << nvcols <<
"\t" << vrtcmp << endl;
248 cout <<
"Row, Col of upper left pt in blocks\n";
249 for(
int i = 0; (i<n+1) && (rcmstr[i]!=n+1); ++i )
250 cout << i <<
" " << rcmstr[i] <<
" " << ccmstr[i] << endl;
251 cout <<
"-----------------------------------------\n";
254 if( hrzcmp || vrtcmp )
256 cout <<
"FAILED! hrz cmp's:" << hrzcmp <<
" vrtcmp: " << vrtcmp << endl;
262 vector<int> rowperm_t( n );
263 vector<int> colperm_t( n );
264 for(
int i = 0; i < n; ++i )
267 rowperm_t[i] = rowperm[i];
268 colperm_t[i] = colperm[i];
273 vector<int> myElements( n );
274 OldRowMap.MyGlobalElements( &myElements[0] );
276 vector<int> newDomainElements( n );
277 vector<int> newRangeElements( n );
278 for(
int i = 0; i < n; ++i )
280 newDomainElements[ i ] = myElements[ rowperm_t[i] ];
281 newRangeElements[ i ] = myElements[ colperm_t[i] ];
284 NewRowMap_ =
new Epetra_Map( n, n, &newDomainElements[0], OldRowMap.IndexBase(), OldRowMap.Comm() );
285 NewColMap_ =
new Epetra_Map( n, n, &newRangeElements[0], OldColMap.IndexBase(), OldColMap.Comm() );
289 cout <<
"New Row Map\n";
290 cout << *NewRowMap_ << endl;
291 cout <<
"New ColMap\n";
292 cout << *NewColMap_ << endl;
296 NewGraph_ =
new Epetra_CrsGraph( Copy, *NewRowMap_, *NewColMap_, 0 );
297 Importer_ =
new Epetra_Import( *NewRowMap_, OldRowMap );
298 NewGraph_->Import( strippedGraph, *Importer_, Insert );
299 NewGraph_->FillComplete();
302 cout <<
"NewGraph\n";
306 NewMatrix_ =
new Epetra_CrsMatrix( Copy, *NewGraph_ );
307 NewMatrix_->Import( orig, *Importer_, Insert );
308 NewMatrix_->FillComplete();
312 cout <<
"New CrsMatrix\n";
313 cout << *NewMatrix_ << endl;
325 int ret = NewMatrix_->Import( *
origObj_, *Importer_, Insert );
326 if (ret<0)
return false;
334 int ret =
origObj_->Export( *NewMatrix_, *Importer_, Insert );
335 if (ret<0)
return false;
NewTypeRef operator()(OriginalTypeRef orig)
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...