53 #include <Teuchos_ParameterList.hpp>
54 #include <Teuchos_RCP.hpp>
55 #include <Teuchos_FancyOStream.hpp>
56 #include <Teuchos_CommandLineProcessor.hpp>
57 #include <Tpetra_CrsMatrix.hpp>
58 #include <Tpetra_Vector.hpp>
59 #include <MatrixMarket_Tpetra.hpp>
83 typedef Tpetra::CrsMatrix<z2TestScalar, z2TestLO, z2TestGO>
SparseMatrix;
85 typedef Tpetra::Vector<z2TestScalar, z2TestLO, z2TestGO>
VectorT;
86 typedef VectorT::node_type
Node;
94 typedef Tpetra::Vector<int, z2TestLO, z2TestGO>
IntVector;
97 #define epsilon 0.00000001
103 std::string inputFile =
"";
104 std::string outputFile =
"";
106 bool verbose =
false;
107 bool distributeInput =
true;
108 bool haveFailure =
false;
113 bool isNormalized =
false;
114 bool isGeneralized =
false;
115 std::string precType =
"jacobi";
116 std::string initialGuess =
"random";
117 bool useFullOrtho =
true;
120 Tpetra::ScopeGuard tscope(&narg, &arg);
121 RCP<const Teuchos::Comm<int> > comm = Tpetra::getDefaultComm();
122 int me = comm->getRank();
123 int commsize = comm->getSize();
126 Teuchos::CommandLineProcessor cmdp (
false,
false);
127 cmdp.setOption(
"inputPath", &inputPath,
128 "Path to the MatrixMarket or Zoltan file to be read; "
129 "if not specified, a default path will be used.");
130 cmdp.setOption(
"inputFile", &inputFile,
131 "Name of the Matrix Market or Zoltan file to read; "
132 "if not specified, a matrix will be generated by MueLu.");
133 cmdp.setOption(
"outputFile", &outputFile,
134 "Name of the Matrix Market sparse matrix file to write, "
135 "echoing the input/generated matrix.");
136 cmdp.setOption(
"vertexWeights", &nVwgts,
137 "Number of weights to generate for each vertex");
138 cmdp.setOption(
"verbose",
"quiet", &verbose,
139 "Print messages and results.");
140 cmdp.setOption(
"distribute",
"no-distribute", &distributeInput,
141 "indicate whether or not to distribute "
142 "input across the communicator");
144 cmdp.setOption(
"normalized",
"combinatorial", &isNormalized,
145 "indicate whether or not to use a normalized Laplacian.");
146 cmdp.setOption(
"generalized",
"non-generalized", &isGeneralized,
147 "indicate whether or not to use a generalized Laplacian.");
148 cmdp.setOption(
"precond", &precType,
149 "indicate which preconditioner to use [muelu|jacobi|polynomial].");
150 cmdp.setOption(
"initialGuess", &initialGuess,
151 "initial guess for LOBPCG");
152 cmdp.setOption(
"useFullOrtho",
"partialOrtho", &useFullOrtho,
153 "use full orthogonalization.");
164 std::string matrixType(
"Laplace3D");
166 cmdp.setOption(
"x", &xdim,
167 "number of gridpoints in X dimension for "
168 "mesh used to generate matrix.");
169 cmdp.setOption(
"y", &ydim,
170 "number of gridpoints in Y dimension for "
171 "mesh used to generate matrix.");
172 cmdp.setOption(
"z", &zdim,
173 "number of gridpoints in Z dimension for "
174 "mesh used to generate matrix.");
175 cmdp.setOption(
"matrix", &matrixType,
176 "Matrix type: Laplace1D, Laplace2D, or Laplace3D");
179 cmdp.parse(narg, arg);
181 RCP<UserInputForTests> uinput;
185 true, distributeInput));
189 true, distributeInput));
191 RCP<SparseMatrix> origMatrix = uinput->getUITpetraCrsMatrix();
193 if (origMatrix->getGlobalNumRows() < 40) {
194 Teuchos::FancyOStream out(Teuchos::rcp(&std::cout,
false));
195 origMatrix->describe(out, Teuchos::VERB_EXTREME);
199 if (outputFile !=
"") {
201 Tpetra::MatrixMarket::Writer<SparseMatrix>::writeSparseFile(outputFile,
202 origMatrix, verbose);
206 std::cout <<
"NumRows = " << origMatrix->getGlobalNumRows() << std::endl
207 <<
"NumNonzeros = " << origMatrix->getGlobalNumEntries() << std::endl
208 <<
"NumProcs = " << comm->getSize() << std::endl
209 <<
"NumLocalRows (rank 0) = " << origMatrix->getLocalNumRows() << std::endl;
212 RCP<VectorT> origVector, origProd;
213 origProd = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
214 origMatrix->getRangeMap());
215 origVector = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
216 origMatrix->getDomainMap());
217 origVector->randomize();
220 Teuchos::RCP<Teuchos::ParameterList> params(
new Teuchos::ParameterList);
221 params->set(
"num_global_parts", commsize);
222 Teuchos::RCP<Teuchos::ParameterList> sphynxParams(
new Teuchos::ParameterList);
223 sphynxParams->set(
"sphynx_skip_preprocessing",
true);
224 sphynxParams->set(
"sphynx_preconditioner_type", precType);
225 sphynxParams->set(
"sphynx_verbosity", verbose ? 1 : 0);
226 sphynxParams->set(
"sphynx_initial_guess", initialGuess);
227 sphynxParams->set(
"sphynx_use_full_ortho", useFullOrtho);
228 std::string problemType =
"combinatorial";
230 problemType =
"normalized";
231 else if(isGeneralized)
232 problemType =
"generalized";
233 sphynxParams->set(
"sphynx_problem_type", problemType);
236 Teuchos::RCP<SparseGraphAdapter> adapter = Teuchos::rcp(
new SparseGraphAdapter(origMatrix->getCrsGraph(), nVwgts));
245 size_t nrows = origMatrix->getLocalNumRows();
248 for (
size_t i = 0; i < nrows; i++) {
249 size_t idx = i * nVwgts;
250 vwgts[idx] =
zscalar_t(origMatrix->getRowMap()->getGlobalElement(i));
251 for (
int j = 1; j < nVwgts; j++) vwgts[idx+j] = 1.;
253 for (
int j = 0; j < nVwgts; j++) {
254 if (j !=
NNZ_IDX) adapter->setVertexWeights(&vwgts[j], nVwgts, j);
255 else adapter->setVertexWeightIsDegree(
NNZ_IDX);
264 if (me == 0) std::cout <<
"Calling solve() " << std::endl;
268 if (me == 0) std::cout <<
"Done solve() " << std::endl;
270 catch (std::runtime_error &e) {
272 std::cout <<
"Runtime exception returned from solve(): " << e.what();
273 if (!strncmp(e.what(),
"BUILD ERROR", 11)) {
275 std::cout <<
" PASS" << std::endl;
280 std::cout <<
" FAIL" << std::endl;
284 catch (std::logic_error &e) {
286 std::cout <<
"Logic exception returned from solve(): " << e.what()
287 <<
" FAIL" << std::endl;
290 catch (std::bad_alloc &e) {
292 std::cout <<
"Bad_alloc exception returned from solve(): " << e.what()
293 <<
" FAIL" << std::endl;
296 catch (std::exception &e) {
298 std::cout <<
"Unknown exception returned from solve(). " << e.what()
299 <<
" FAIL" << std::endl;
305 size_t checkNparts = comm->getSize();
307 size_t checkLength = origMatrix->getLocalNumRows();
311 size_t *countPerPart =
new size_t[checkNparts];
312 size_t *globalCountPerPart =
new size_t[checkNparts];
315 for (
size_t i = 0; i < checkNparts; i++) countPerPart[i] = 0;
316 for (
size_t i = 0; i < checkNparts * nVwgts; i++) wtPerPart[i] = 0.;
318 for (
size_t i = 0; i < checkLength; i++) {
319 if (
size_t(checkParts[i]) >= checkNparts)
320 std::cout <<
"Invalid Part " << checkParts[i] <<
": FAIL" << std::endl;
321 countPerPart[checkParts[i]]++;
322 for (
int j = 0; j < nVwgts; j++) {
324 wtPerPart[checkParts[i]*nVwgts+j] += vwgts[i*nVwgts+j];
326 wtPerPart[checkParts[i]*nVwgts+j] += origMatrix->getNumEntriesInLocalRow(i);
329 Teuchos::reduceAll<int, size_t>(*comm, Teuchos::REDUCE_SUM, checkNparts,
330 countPerPart, globalCountPerPart);
331 Teuchos::reduceAll<int, zscalar_t>(*comm, Teuchos::REDUCE_SUM,
333 wtPerPart, globalWtPerPart);
335 size_t min = std::numeric_limits<std::size_t>::max();
338 size_t minrank = 0, maxrank = 0;
339 for (
size_t i = 0; i < checkNparts; i++) {
340 if (globalCountPerPart[i] < min) {min = globalCountPerPart[i]; minrank = i;}
341 if (globalCountPerPart[i] > max) {max = globalCountPerPart[i]; maxrank = i;}
342 sum += globalCountPerPart[i];
346 float avg = (float) sum / (
float) checkNparts;
347 std::cout <<
"Minimum count: " << min <<
" on rank " << minrank << std::endl;
348 std::cout <<
"Maximum count: " << max <<
" on rank " << maxrank << std::endl;
349 std::cout <<
"Average count: " << avg << std::endl;
350 std::cout <<
"Total count: " << sum
351 << (sum != origMatrix->getGlobalNumRows()
352 ?
"Work was lost; FAIL"
355 std::cout <<
"Imbalance: " << max / avg << std::endl;
357 std::vector<zscalar_t> minwt(nVwgts, std::numeric_limits<zscalar_t>::max());
358 std::vector<zscalar_t> maxwt(nVwgts, 0.);
359 std::vector<zscalar_t> sumwt(nVwgts, 0.);
360 for (
size_t i = 0; i < checkNparts; i++) {
361 for (
int j = 0; j < nVwgts; j++) {
362 size_t idx = i*nVwgts+j;
363 if (globalWtPerPart[idx] < minwt[j]) minwt[j] = globalWtPerPart[idx];
364 if (globalWtPerPart[idx] > maxwt[j]) maxwt[j] = globalWtPerPart[idx];
365 sumwt[j] += globalWtPerPart[idx];
368 for (
int j = 0; j < nVwgts; j++) {
369 float avgwt = (float) sumwt[j] / (
float) checkNparts;
370 std::cout << std::endl;
371 std::cout <<
"Minimum weight[" << j <<
"]: " << minwt[j] << std::endl;
372 std::cout <<
"Maximum weight[" << j <<
"]: " << maxwt[j] << std::endl;
373 std::cout <<
"Average weight[" << j <<
"]: " << avgwt << std::endl;
374 std::cout <<
"Imbalance: " << maxwt[j] / avgwt << std::endl;
379 delete [] countPerPart;
381 delete [] globalCountPerPart;
382 delete [] globalWtPerPart;
386 if (me == 0) std::cout <<
"Redistributing matrix..." << std::endl;
390 problem.getSolution());
391 if (redistribMatrix->getGlobalNumRows() < 40) {
392 Teuchos::FancyOStream out(Teuchos::rcp(&std::cout,
false));
393 redistribMatrix->describe(out, Teuchos::VERB_EXTREME);
396 if (me == 0) std::cout <<
"Redistributing vectors..." << std::endl;
400 problem.getSolution());
402 RCP<VectorT> redistribProd;
403 redistribProd = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
404 redistribMatrix->getRangeMap());
408 RCP<IntVector> origIntVec;
410 origIntVec = Tpetra::createVector<int,z2TestLO,z2TestGO>(
411 origMatrix->getRangeMap());
412 for (
size_t i = 0; i < origIntVec->getLocalLength(); i++)
413 origIntVec->replaceLocalValue(i, me);
417 problem.getSolution());
418 int origIntNorm = origIntVec->norm1();
419 int redistIntNorm = redistIntVec->norm1();
420 if (me == 0) std::cout <<
"IntegerVectorTest: " << origIntNorm <<
" == "
421 << redistIntNorm <<
" ?";
422 if (origIntNorm != redistIntNorm) {
423 if (me == 0) std::cout <<
" FAIL" << std::endl;
426 else if (me == 0) std::cout <<
" OK" << std::endl;
432 if (me == 0) std::cout <<
"Matvec original..." << std::endl;
433 origMatrix->apply(*origVector, *origProd);
436 std::cout <<
"Norm of Original matvec prod: " << origNorm << std::endl;
438 if (me == 0) std::cout <<
"Matvec redistributed..." << std::endl;
439 redistribMatrix->apply(*redistribVector, *redistribProd);
442 std::cout <<
"Norm of Redistributed matvec prod: " << redistribNorm << std::endl;
444 if (redistribNorm > origNorm+
epsilon || redistribNorm < origNorm-
epsilon) {
449 delete redistribVector;
450 delete redistribMatrix;
454 std::cout <<
"Mat-Vec product changed; FAIL" << std::endl;
458 std::cout <<
"PASS" << std::endl;
Provides access for Zoltan2 to Xpetra::CrsMatrix data.
Tpetra::CrsMatrix< z2TestScalar, z2TestLO, z2TestGO > SparseMatrix
Provides access for Zoltan2 to Xpetra::CrsGraph data.
Zoltan2::XpetraCrsGraphAdapter< SparseGraph > SparseGraphAdapter
int main(int narg, char **arg)
common code used by tests
typename InputTraits< User >::part_t part_t
Defines the XpetraMultiVectorAdapter.
Defines XpetraCrsGraphAdapter class.
Defines the XpetraCrsMatrixAdapter class.
void applyPartitioningSolution(const User &in, User *&out, const PartitioningSolution< Adapter > &solution) const
An adapter for Xpetra::MultiVector.
Tpetra::Map::local_ordinal_type zlno_t
Zoltan2::XpetraMultiVectorAdapter< IntVector > IntVectorAdapter
Tpetra::CrsGraph< z2TestLO, z2TestGO > SparseGraph
Zoltan2::XpetraCrsMatrixAdapter< SparseMatrix > SparseMatrixAdapter
Tpetra::Vector< int, z2TestLO, z2TestGO > IntVector
void applyPartitioningSolution(const User &in, User *&out, const PartitioningSolution< Adapter > &solution) const
Zoltan2::XpetraMultiVectorAdapter< Vector > MultiVectorAdapter
Tpetra::Map::global_ordinal_type zgno_t
Tpetra::Vector< z2TestScalar, z2TestLO, z2TestGO > VectorT
void solve(bool updateInputData=true)
Direct the problem to create a solution.
std::string testDataFilePath(".")