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 //
41
42 /*
43 \file quick_solve.cpp
44 \author Eric Bavier <etbavie@sandia.gov>
45 \date Thu Jul 14 16:24:46 MDT 2011
46
47 \brief Intended to quickly check a solver interface
48
49 This example solves a simple sparse system of linear equations
50 using a given Amesos2 solver interface.
51 */
52
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>
58
59 #include <Tpetra_DefaultPlatform.hpp>
60 #include <Tpetra_Map.hpp>
61 #include <Tpetra_MultiVector.hpp>
62 #include <Tpetra_CrsMatrix.hpp>
63
64 #include <MatrixMarket_Tpetra.hpp>
65
66 #include "Amesos2.hpp"
67 #include "Amesos2_Version.hpp"
68
69 #include <petsc.h>
70
main(int argc,char * argv[])71 int main(int argc, char *argv[])
72 {
73 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
74
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;
81
82 typedef Tpetra::CrsMatrix<Scalar, LO, GO> MAT;
83 typedef Tpetra::MultiVector<Scalar, LO, GO> MV;
84
85 using Teuchos::RCP;
86 using Teuchos::rcp;
87 using Teuchos::tuple;
88 using Tpetra::global_size_t;
89
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();
94
95 RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
96
97 *fos << Amesos2::version() << std::endl << std::endl;
98
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;
120
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 }
127
128 const size_t numVectors = 1;
129
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));
133
134 std::string mat_pathname = filedir + filename;
135 RCP<MAT> A = Tpetra::MatrixMarket::Reader<MAT>::readSparseFile(mat_pathname, comm, node);
136
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 }
142
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();
146
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();
153
154 RCP<MV> B = rcp(new MV(rngmap, numVectors));
155 A->apply(*X, *B);
156
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 }
165
166 solver->numericFactorization();
167
168 if (printLUStats && myRank == 0) {
169 Amesos2::Status solver_status = solver->getStatus();
170 *fos << "L+U nnz = " << solver_status.getNnzLU() << std::endl;
171 }
172
173 solver->solve();
174
175 if (printSolution) {
176 // Print the solution
177 Xhat->describe(*fos, Teuchos::VERB_EXTREME);
178 X->describe(*fos, Teuchos::VERB_EXTREME);
179 }
180
181 if (printTiming) {
182 // Print some timing statistics
183 solver->printTiming(*fos);
184 }
185
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 }
192
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 }
216
217 /*TEST
218
219 build:
220 requires: trilinos
221
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)"
226
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)"
232
233 TEST*/
234