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)
61extern 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);
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];
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 );
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] ];
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;
326 if (ret<0)
return false;
335 if (ret<0)
return false;
NewTypeRef operator()(OriginalTypeRef orig)
Analysis of transform operation on original object and construction of new object.
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...
int MyGlobalElements(int *MyGlobalElementList) const
const Epetra_Comm & Comm() const
int InsertMyIndices(int LocalRow, int NumIndices, int *Indices)
int FillComplete(bool OptimizeDataStorage=true)
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.