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: Teuchos::GlobalMPISession mpiSession(&argc,&argv);
74: typedef double Scalar;
75: typedef Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
76: typedef int LO;
77: typedef int GO;
78: typedef Tpetra::DefaultPlatform::DefaultPlatformType Platform;
79: typedef Tpetra::DefaultPlatform::DefaultPlatformType::NodeType Node;
81: typedef Tpetra::CrsMatrix<Scalar,LO,GO> MAT;
82: typedef Tpetra::MultiVector<Scalar,LO,GO> MV;
84: using Tpetra::global_size_t;
85: using Teuchos::tuple;
86: using Teuchos::RCP;
87: using Teuchos::rcp;
89: Platform &platform = Tpetra::DefaultPlatform::getDefaultPlatform();
90: Teuchos::RCP<const Teuchos::Comm<int> > comm = platform.getComm();
91: Teuchos::RCP<Node> node = platform.getNode();
92: size_t myRank = comm->getRank();
94: RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
96: *fos << Amesos2::version() << std::endl << std::endl;
98: bool printMatrix = false;
99: bool printSolution = false;
100: bool printTiming = false;
101: bool printResidual = false;
102: bool printLUStats = false;
103: bool verbose = false;
104: std::string solver_name;
105: std::string filedir;
106: std::string filename;
107: Teuchos::CommandLineProcessor cmdp(false,false); // set last argument to false so Trilinos won't generate exceptionif it sees unrecognized option
108: // note it still prints a warning to stderr which cannot be stopped.
109: cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
110: cmdp.setOption("filedir",&filedir,"Directory where matrix-market files are located");
111: cmdp.setOption("filename",&filename,"Filename for Matrix-Market test matrix.");
112: cmdp.setOption("print-matrix","no-print-matrix",&printMatrix,"Print the full matrix after reading it.");
113: cmdp.setOption("print-solution","no-print-solution",&printSolution,"Print solution vector after solve.");
114: cmdp.setOption("print-timing","no-print-timing",&printTiming,"Print solver timing statistics");
115: cmdp.setOption("print-residual","no-print-residual",&printResidual,"Print solution residual");
116: cmdp.setOption("print-lu-stats","no-print-lu-stats",&printLUStats,"Print nnz in L and U factors");
117: cmdp.setOption("solver", &solver_name, "Which TPL solver library to use.");
118: if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
119: std::cerr << "Options unknown to Trilinos in command line"<< std::endl;
120: }
122: // Before we do anything, check that the solver is enabled
123: if (!Amesos2::query(solver_name)){
124: std::cerr << solver_name << " not enabled. Exiting..." << std::endl;
125: return EXIT_SUCCESS; // Otherwise CTest will pick it up as
126: // failure, which it isn't really
127: }
129: const size_t numVectors = 1;
131: // create a Map
132: global_size_t nrows = 6;
133: RCP<Tpetra::Map<LO,GO> > map
134: = rcp( new Tpetra::Map<LO,GO>(nrows,0,comm));
136: std::string mat_pathname = filedir + filename;
137: RCP<MAT> A = Tpetra::MatrixMarket::Reader<MAT>::readSparseFile(mat_pathname,comm,node);
139: if (printMatrix){
140: A->describe(*fos, Teuchos::VERB_EXTREME);
141: }
142: else if (verbose && myRank==0){
143: *fos << std::endl << A->description() << std::endl << std::endl;
144: }
146: // get the maps
147: RCP<const Tpetra::Map<LO,GO,Node> > dmnmap = A->getDomainMap();
148: RCP<const Tpetra::Map<LO,GO,Node> > rngmap = A->getRangeMap();
150: // Create random X
151: RCP<MV> Xhat = rcp( new MV(dmnmap,numVectors));
152: RCP<MV> X = rcp( new MV(dmnmap,numVectors));
153: X->setObjectLabel("X");
154: Xhat->setObjectLabel("Xhat");
155: X->randomize();
157: RCP<MV> B = rcp(new MV(rngmap,numVectors));
158: A->apply(*X, *B);
160: // Constructor from Factory
161: RCP<Amesos2::Solver<MAT,MV> > solver;
162: try{
163: solver = Amesos2::create<MAT,MV>(solver_name, A, Xhat, B);
164: } catch (std::invalid_argument e){
165: *fos << e.what() << std::endl;
166: return 0;
167: }
169: solver->numericFactorization();
171: if (printLUStats && myRank == 0){
172: Amesos2::Status solver_status = solver->getStatus();
173: *fos << "L+U nnz = " << solver_status.getNnzLU() << std::endl;
174: }
176: solver->solve();
178: if (printSolution){
179: // Print the solution
180: Xhat->describe(*fos,Teuchos::VERB_EXTREME);
181: X->describe(*fos,Teuchos::VERB_EXTREME);
182: }
184: if (printTiming){
185: // Print some timing statistics
186: solver->printTiming(*fos);
187: }
189: if (printResidual){
190: Teuchos::Array<Magnitude> xhatnorms(numVectors);
191: Xhat->update(-1.0, *X, 1.0);
192: Xhat->norm2(xhatnorms());
193: if (myRank == 0){
194: *fos << "Norm2 of Ax - b = " << xhatnorms << std::endl;
195: }
196: }
198: PetscErrorCodePetscInitialize(&argc,&argv,NULL,NULL);if (ierr) return ierr;
199: PetscOptionsSetValue(NULL,"-options_left","false");
200: KSP ksp;
201: KSPCreate(PETSC_COMM_WORLD,&ksp);
202: Mat Apetsc;
203: MatCreate(PETSC_COMM_WORLD,&Apetsc);
204: PetscViewer viewer;
205: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"${PETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int32-float64",FILE_MODE_READ,&viewer);
206: MatLoad(Apetsc,viewer);
207: Vec x,b;
208: MatCreateVecs(Apetsc,&x,&b);
209: PetscViewerDestroy(&viewer);
210: KSPSetOperators(ksp,Apetsc,Apetsc);
211: KSPSetFromOptions(ksp);
212: KSPSolve(ksp,x,b);
213: VecDestroy(&x);
214: VecDestroy(&b);
215: MatDestroy(&Apetsc);
216: KSPDestroy(&ksp);
217: PetscFinalize();
218: return ierr;
219: }
222: /*TEST
224: build:
225: requires: trilinos
227: test:
228: requires: superlu
229: 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
230: filter: egrep -v "(Teuchos|Amesos2)"
232: test:
233: suffix: 2
234: requires: superlu_dist
235: 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
236: filter: egrep -v "(Teuchos|Amesos2)"
238: TEST*/