50 #include <Teuchos_ParameterList.hpp> 51 #include <Teuchos_RCP.hpp> 52 #include <Teuchos_CommandLineProcessor.hpp> 53 #include <Tpetra_DefaultPlatform.hpp> 54 #include <Tpetra_CrsMatrix.hpp> 55 #include <Tpetra_Vector.hpp> 56 #include <Galeri_XpetraMaps.hpp> 57 #include <Galeri_XpetraProblemFactory.hpp> 67 int main(
int narg,
char** arg)
70 Teuchos::GlobalMPISession mpiSession(&narg, &arg, NULL);
71 RCP<const Teuchos::Comm<int> > comm =
72 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
73 int me = comm->getRank();
77 typedef double scalar_t;
78 typedef Tpetra::Map<> Map_t;
79 typedef Map_t::local_ordinal_type localId_t;
80 typedef Map_t::global_ordinal_type globalId_t;
81 typedef Tpetra::CrsMatrix<scalar_t, localId_t, globalId_t> Matrix_t;
82 typedef Tpetra::MultiVector<scalar_t, localId_t, globalId_t> MultiVector_t;
83 typedef Tpetra::Vector<scalar_t, localId_t, globalId_t> Vector_t;
90 std::string method =
"scotch";
91 globalId_t nx = 50, ny = 40, nz = 30;
95 Teuchos::CommandLineProcessor cmdp (
false,
false);
96 cmdp.setOption(
"method", &method,
97 "Partitioning method to use: scotch or parmetis.");
98 cmdp.setOption(
"nx", &nx,
99 "number of gridpoints in X dimension for " 100 "mesh used to generate matrix; must be >= 1.");
101 cmdp.setOption(
"ny", &ny,
102 "number of gridpoints in Y dimension for " 103 "mesh used to generate matrix; must be >= 1.");
104 cmdp.setOption(
"nz", &nz,
105 "number of gridpoints in Z dimension for " 106 "mesh used to generate matrix; must be >= 1.");
107 cmdp.parse(narg, arg);
109 if ((nx < 1) || (ny < 1) || (nz < 1)) {
110 cout <<
"Input error: nx, ny and nz must be >= 1" << endl;
118 Teuchos::ParameterList galeriList;
119 galeriList.set(
"nx", nx);
120 galeriList.set(
"ny", ny);
121 galeriList.set(
"nz", nz);
124 RCP<Matrix_t> origMatrix;
127 RCP<const Map_t> map = rcp(
new Map_t(nGlobalElements, 0, comm));
129 typedef Galeri::Xpetra::Problem<Map_t,Matrix_t,MultiVector_t> Galeri_t;
130 RCP<Galeri_t> galeriProblem =
131 Galeri::Xpetra::BuildProblem<scalar_t, localId_t, globalId_t,
132 Map_t, Matrix_t, MultiVector_t>
133 (
"Laplace3D", map, galeriList);
134 origMatrix = galeriProblem->BuildMatrix();
136 catch (std::exception &e) {
137 cout <<
"Exception in Galeri matrix generation. " << e.what() << endl;
142 cout <<
"NumRows = " << origMatrix->getGlobalNumRows() << endl
143 <<
"NumNonzeros = " << origMatrix->getGlobalNumEntries() << endl
144 <<
"NumProcs = " << comm->getSize() << endl;
147 RCP<Vector_t> origVector, origProd;
148 origProd = Tpetra::createVector<scalar_t,localId_t,globalId_t>(
149 origMatrix->getRangeMap());
150 origVector = Tpetra::createVector<scalar_t,localId_t,globalId_t>(
151 origMatrix->getDomainMap());
152 origVector->randomize();
155 Teuchos::ParameterList param;
156 param.set(
"partitioning_approach",
"partition");
157 param.set(
"algorithm", method);
160 MatrixAdapter_t adapter(origMatrix);
168 catch (std::exception &e) {
169 cout <<
"Exception returned from solve(). " << e.what() << endl;
176 if (me == 0) cout <<
"Redistributing matrix..." << endl;
177 RCP<Matrix_t> redistribMatrix;
178 adapter.applyPartitioningSolution(*origMatrix, redistribMatrix,
181 if (me == 0) cout <<
"Redistributing vectors..." << endl;
182 RCP<Vector_t> redistribVector;
183 MultiVectorAdapter_t adapterVector(origVector);
184 adapterVector.applyPartitioningSolution(*origVector, redistribVector,
188 RCP<Vector_t> redistribProd;
189 redistribProd = Tpetra::createVector<scalar_t,localId_t,globalId_t>(
190 redistribMatrix->getRangeMap());
194 if (origMatrix->getGlobalNumRows() <= 50) {
195 cout << me <<
" ORIGINAL: ";
196 for (
size_t i = 0; i < origVector->getLocalLength(); i++)
197 cout << origVector->getMap()->getGlobalElement(i) <<
" ";
199 cout << me <<
" REDISTRIB: ";
200 for (
size_t i = 0; i < redistribVector->getLocalLength(); i++)
201 cout << redistribVector->getMap()->getGlobalElement(i) <<
" ";
209 if (me == 0) cout <<
"Matvec original..." << endl;
210 origMatrix->apply(*origVector, *origProd);
211 scalar_t origNorm = origProd->norm2();
213 cout <<
"Norm of Original matvec prod: " << origNorm << endl;
215 if (me == 0) cout <<
"Matvec redistributed..." << endl;
216 redistribMatrix->apply(*redistribVector, *redistribProd);
217 scalar_t redistribNorm = redistribProd->norm2();
219 cout <<
"Norm of Redistributed matvec prod: " << redistribNorm << endl;
222 const double epsilon = 0.00000001;
223 if (redistribNorm > origNorm+
epsilon || redistribNorm < origNorm-
epsilon)
224 cout <<
"Mat-Vec product changed; FAIL" << std::endl;
226 cout <<
"PASS" << endl;
Provides access for Zoltan2 to Xpetra::CrsMatrix data.
Defines the XpetraMultiVectorAdapter.
Defines the XpetraCrsMatrixAdapter class.
int main(int narg, char **arg)
An adapter for Xpetra::MultiVector.
const PartitioningSolution< Adapter > & getSolution()
Get the solution to the problem.
PartitioningProblem sets up partitioning problems for the user.
Defines the PartitioningProblem class.
void solve(bool updateInputData=true)
Direct the problem to create a solution.