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*/