Actual source code: ex63.cxx
1: //
2: // ***********************************************************************
3: //
4: // Amesos2: Templated Direct Sparse Solver Package
5: // Copyright 2011 Sandia Corporation
6: //
7: // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8: // the U.S. Government retains certain rights in this software.
9: //
10: // Redistribution and use in source and binary forms, with or without
11: // modification, are permitted provided that the following conditions are
12: // met:
13: //
14: // 1. Redistributions of source code must retain the above copyright
15: // notice, this list of conditions and the following disclaimer.
16: //
17: // 2. Redistributions in binary form must reproduce the above copyright
18: // notice, this list of conditions and the following disclaimer in the
19: // documentation and/or other materials provided with the distribution.
20: //
21: // 3. Neither the name of the Corporation nor the names of the
22: // contributors may be used to endorse or promote products derived from
23: // this software without specific prior written permission.
24: //
25: // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26: // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27: // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28: // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29: // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30: // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31: // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32: // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33: // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34: // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35: // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36: //
37: // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38: //
39: // ***********************************************************************
40: //
42: /*
43: \file quick_solve.cpp
44: \author Eric Bavier <etbavie@sandia.gov>
45: \date Thu Jul 14 16:24:46 MDT 2011
47: \brief Intended to quickly check a solver interface
49: This example solves a simple sparse system of linear equations
50: using a given Amesos2 solver interface.
51: */
53: #include <Teuchos_ScalarTraits.hpp>
54: #include <Teuchos_RCP.hpp>
55: #include <Teuchos_GlobalMPISession.hpp>
56: #include <Teuchos_VerboseObject.hpp>
57: #include <Teuchos_CommandLineProcessor.hpp>
59: #include <Tpetra_DefaultPlatform.hpp>
60: #include <Tpetra_Map.hpp>
61: #include <Tpetra_MultiVector.hpp>
62: #include <Tpetra_CrsMatrix.hpp>
64: #include <MatrixMarket_Tpetra.hpp>
66: #include "Amesos2.hpp"
67: #include "Amesos2_Version.hpp"
69: #include <petsc.h>
71: int main(int argc, char *argv[])
72: {
73: Teuchos::GlobalMPISession mpiSession(&argc, &argv);
75: typedef double Scalar;
76: typedef Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
77: typedef int LO;
78: typedef int GO;
79: typedef Tpetra::DefaultPlatform::DefaultPlatformType Platform;
80: typedef Tpetra::DefaultPlatform::DefaultPlatformType::NodeType Node;
82: typedef Tpetra::CrsMatrix<Scalar, LO, GO> MAT;
83: typedef Tpetra::MultiVector<Scalar, LO, GO> MV;
85: using Teuchos::RCP;
86: using Teuchos::rcp;
87: using Teuchos::tuple;
88: using Tpetra::global_size_t;
90: Platform &platform = Tpetra::DefaultPlatform::getDefaultPlatform();
91: Teuchos::RCP<const Teuchos::Comm<int>> comm = platform.getComm();
92: Teuchos::RCP<Node> node = platform.getNode();
93: size_t myRank = comm->getRank();
95: RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
97: *fos << Amesos2::version() << std::endl << std::endl;
99: bool printMatrix = false;
100: bool printSolution = false;
101: bool printTiming = false;
102: bool printResidual = false;
103: bool printLUStats = false;
104: bool verbose = false;
105: std::string solver_name;
106: std::string filedir;
107: std::string filename;
108: Teuchos::CommandLineProcessor cmdp(false, false); // set last argument to false so Trilinos won't generate exceptionif it sees unrecognized option
109: // note it still prints a warning to stderr which cannot be stopped.
110: cmdp.setOption("verbose", "quiet", &verbose, "Print messages and results.");
111: cmdp.setOption("filedir", &filedir, "Directory where matrix-market files are located");
112: cmdp.setOption("filename", &filename, "Filename for Matrix-Market test matrix.");
113: cmdp.setOption("print-matrix", "no-print-matrix", &printMatrix, "Print the full matrix after reading it.");
114: cmdp.setOption("print-solution", "no-print-solution", &printSolution, "Print solution vector after solve.");
115: cmdp.setOption("print-timing", "no-print-timing", &printTiming, "Print solver timing statistics");
116: cmdp.setOption("print-residual", "no-print-residual", &printResidual, "Print solution residual");
117: cmdp.setOption("print-lu-stats", "no-print-lu-stats", &printLUStats, "Print nnz in L and U factors");
118: cmdp.setOption("solver", &solver_name, "Which TPL solver library to use.");
119: if (cmdp.parse(argc, argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) std::cerr << "Options unknown to Trilinos in command line" << std::endl;
121: // Before we do anything, check that the solver is enabled
122: if (!Amesos2::query(solver_name)) {
123: std::cerr << solver_name << " not enabled. Exiting..." << std::endl;
124: return EXIT_SUCCESS; // Otherwise CTest will pick it up as
125: // failure, which it isn't really
126: }
128: const size_t numVectors = 1;
130: // create a Map
131: global_size_t nrows = 6;
132: RCP<Tpetra::Map<LO, GO>> map = rcp(new Tpetra::Map<LO, GO>(nrows, 0, comm));
134: std::string mat_pathname = filedir + filename;
135: RCP<MAT> A = Tpetra::MatrixMarket::Reader<MAT>::readSparseFile(mat_pathname, comm, node);
137: if (printMatrix) {
138: A->describe(*fos, Teuchos::VERB_EXTREME);
139: } else if (verbose && myRank == 0) {
140: *fos << std::endl << A->description() << std::endl << std::endl;
141: }
143: // get the maps
144: RCP<const Tpetra::Map<LO, GO, Node>> dmnmap = A->getDomainMap();
145: RCP<const Tpetra::Map<LO, GO, Node>> rngmap = A->getRangeMap();
147: // Create random X
148: RCP<MV> Xhat = rcp(new MV(dmnmap, numVectors));
149: RCP<MV> X = rcp(new MV(dmnmap, numVectors));
150: X->setObjectLabel("X");
151: Xhat->setObjectLabel("Xhat");
152: X->randomize();
154: RCP<MV> B = rcp(new MV(rngmap, numVectors));
155: A->apply(*X, *B);
157: // Constructor from Factory
158: RCP<Amesos2::Solver<MAT, MV>> solver;
159: try {
160: solver = Amesos2::create<MAT, MV>(solver_name, A, Xhat, B);
161: } catch (std::invalid_argument e) {
162: *fos << e.what() << std::endl;
163: return 0;
164: }
166: solver->numericFactorization();
168: if (printLUStats && myRank == 0) {
169: Amesos2::Status solver_status = solver->getStatus();
170: *fos << "L+U nnz = " << solver_status.getNnzLU() << std::endl;
171: }
173: solver->solve();
175: if (printSolution) {
176: // Print the solution
177: Xhat->describe(*fos, Teuchos::VERB_EXTREME);
178: X->describe(*fos, Teuchos::VERB_EXTREME);
179: }
181: if (printTiming) {
182: // Print some timing statistics
183: solver->printTiming(*fos);
184: }
186: if (printResidual) {
187: Teuchos::Array<Magnitude> xhatnorms(numVectors);
188: Xhat->update(-1.0, *X, 1.0);
189: Xhat->norm2(xhatnorms());
190: if (myRank == 0) *fos << "Norm2 of Ax - b = " << xhatnorms << std::endl;
191: }
193: PetscFunctionBeginUser;
194: PetscCall(PetscInitialize(&argc, &argv, NULL, NULL));
195: PetscCall(PetscOptionsSetValue(NULL, "-options_left", "false"));
196: KSP ksp;
197: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
198: Mat Apetsc;
199: PetscCall(MatCreate(PETSC_COMM_WORLD, &Apetsc));
200: PetscViewer viewer;
201: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "${PETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int32-float64", FILE_MODE_READ, &viewer));
202: PetscCall(MatLoad(Apetsc, viewer));
203: Vec x, b;
204: PetscCall(MatCreateVecs(Apetsc, &x, &b));
205: PetscCall(PetscViewerDestroy(&viewer));
206: PetscCall(KSPSetOperators(ksp, Apetsc, Apetsc));
207: PetscCall(KSPSetFromOptions(ksp));
208: PetscCall(KSPSolve(ksp, x, b));
209: PetscCall(VecDestroy(&x));
210: PetscCall(VecDestroy(&b));
211: PetscCall(MatDestroy(&Apetsc));
212: PetscCall(KSPDestroy(&ksp));
213: PetscCall(PetscFinalize());
214: return 0;
215: }
217: /*TEST
219: build:
220: requires: trilinos
222: test:
223: requires: superlu
224: args: --filedir=${wPETSC_DIR}/share/petsc/datafiles/matrices/ --filename=amesos2_test_mat0.mtx --solver=SuperLU --print-residual=true -ksp_monitor -pc_type lu -pc_factor_mat_solver_type superlu -ksp_view -ksp_converged_reason
225: filter: grep -E -v "(Teuchos|Amesos2)"
227: test:
228: suffix: 2
229: requires: superlu_dist
230: args: --filedir=${wPETSC_DIR}/share/petsc/datafiles/matrices/ --filename=amesos2_test_mat0.mtx --solver=SuperLUDist --print-residual=true -ksp_monitor -pc_type lu -pc_factor_mat_solver_type superlu_dist -ksp_view -ksp_converged_reason
231: filter: grep -E -v "(Teuchos|Amesos2)"
233: TEST*/