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