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);
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;
302 cout <<
"NewGraph\n";
312 cout <<
"New CrsMatrix\n";
313 cout << *NewMatrix_ << endl;