Actual source code: ex63.cxx

petsc-3.6.4 2016-04-12
Report Typos and Errors
  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:   PetscErrorCode PetscOptionsSetValue("-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: }