56 #ifdef HAVE_SHYLUCORE_MPI 57 #include "Epetra_MpiComm.h" 59 #include "Epetra_SerialComm.h" 61 #include "Epetra_CrsMatrix.h" 64 #include "EpetraExt_RowMatrixOut.h" 70 #include "Isorropia_Epetra.hpp" 71 #include "Isorropia_EpetraRedistributor.hpp" 72 #include "Isorropia_EpetraPartitioner.hpp" 80 Epetra_CrsMatrix *balanceAndRedistribute(Epetra_CrsMatrix *A,
86 Epetra_Map ARowMap = A->RowMap();
93 Isorropia::Epetra::Partitioner *partitioner =
new 94 Isorropia::Epetra::Partitioner(A, isoList,
false);
95 partitioner->partition();
97 Isorropia::Epetra::Redistributor rd(partitioner);
98 Epetra_CrsMatrix *newA;
99 rd.redistribute(*A, newA);
101 EpetraExt::RowMatrixToMatlabFile(
"A.mat", *newA);
108 void checkMaps(Epetra_CrsMatrix *A)
111 Epetra_Map AColMap = A->ColMap();
116 Epetra_Map ADomainMap = A->DomainMap();
118 int nelems = ADomainMap.NumMyElements();
122 int *dom_cols = ADomainMap.MyGlobalElements();
126 Epetra_Map ARangeMap = A->RangeMap();
131 int *ran_cols = ARangeMap.MyGlobalElements();
135 Epetra_Map ARowMap = A->RowMap();
140 int *rows = ARowMap.MyGlobalElements();
146 for (
int i = 0; i < nelems ; i++)
149 assert(dom_cols[i] == ran_cols[i]);
150 assert(rows[i] == ran_cols[i]);
155 void findLocalColumns(Epetra_CrsMatrix *A,
int *gvals,
int &SNumGlobalCols)
158 int n = A->NumGlobalRows();
160 Epetra_Map AColMap = A->ColMap();
161 int ncols = AColMap.NumMyElements();
162 int *cols = AColMap.MyGlobalElements();
166 int *vals =
new int[n];
167 for (
int i = 0; i < n ; i++)
174 for (
int i = 0; i < ncols ; i++)
180 A->Comm().SumAll(vals, gvals, n);
183 for (
int i = 0; i < n ; i++)
202 void findNarrowSeparator(Epetra_CrsMatrix *A,
int *gvals)
207 int n = A->NumGlobalRows();
209 int myPID = A->Comm().MyPID();
212 Epetra_Map rMap = A->RowMap();
213 Epetra_Map cMap = A->ColMap();
214 int *rows = rMap.MyGlobalElements();
215 int relems = rMap.NumMyElements();
217 int *vals =
new int[n];
218 int *allGIDs =
new int[n];
219 for (
int i = 0; i < n ; i++)
225 for (
int i = 0; i < relems ; i++)
227 vals[rows[i]] = myPID;
234 A->Comm().SumAll(vals, allGIDs, n);
237 for (
int i = 0; i < n ; i++)
241 for (
int i = 0; i < relems; i++)
248 bool movetoBlockDiagonal =
false;
253 (void) A->ExtractMyRowView(i, nentries, values, indices);
256 assert(nentries != 0);
257 for (
int j = 0; j < nentries; j++)
259 cgid = cMap.GID(indices[j]);
261 if (gvals[cgid] == 1 || allGIDs[cgid] == myPID)
266 if (allGIDs[cgid] < myPID)
269 movetoBlockDiagonal =
true;
279 movetoBlockDiagonal =
false;
309 if (movetoBlockDiagonal)
326 for (
int i = 0; i < n ; i++)
329 A->Comm().SumAll(vals, allGIDs, n);
330 for (
int i = 0; i < n ; i++)
345 void findBlockElems(Epetra_CrsMatrix *A,
int nrows,
int *rows,
int *gvals,
346 int Lnr,
int *LeftElems,
347 int Rnr,
int *RightElems,
string s1,
string s2,
bool cols)
351 int rcnt = 0;
int lcnt = 0;
353 ostringstream ssmsg1;
354 ostringstream ssmsg2;
357 ostringstream fnamestr;
358 fnamestr << s1 <<
".mat";
359 string fname = fnamestr.str();
360 ofstream outD(fname.c_str());
362 ostringstream fnamestrR;
363 fnamestrR << s2 <<
".mat";
364 string fnameR = fnamestrR.str();
365 ofstream outR(fnameR.c_str());
370 for (
int i = 0; i < nrows; i++)
373 assert (gvals[gid] >= 1);
380 if (cols && A->LRID(gid) == -1)
continue;
382 LeftElems[lcnt++] = gid;
383 ssmsg1 << gid <<
" ";
391 RightElems[rcnt++] = gid;
392 ssmsg2 << gid <<
" ";
405 cout << ssmsg1.str() << endl;
406 cout << ssmsg2.str() << endl;
408 ssmsg1.clear(); ssmsg1.str(
"");
409 ssmsg2.clear(); ssmsg2.str(
"");