Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_SolverCore_def.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 //
00003 // ***********************************************************************
00004 //
00005 //           Amesos2: Templated Direct Sparse Solver Package
00006 //                  Copyright 2011 Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ***********************************************************************
00041 //
00042 // @HEADER
00043 
00053 #ifndef AMESOS2_SOLVERCORE_DEF_HPP
00054 #define AMESOS2_SOLVERCORE_DEF_HPP
00055 
00056 #include "Amesos2_MatrixAdapter_def.hpp"
00057 #include "Amesos2_MultiVecAdapter_def.hpp"
00058 
00059 #include "Amesos2_Util.hpp"
00060 
00061 
00062 namespace Amesos2 {
00063 
00064 
00065 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00066 SolverCore<ConcreteSolver,Matrix,Vector>::SolverCore(
00067   Teuchos::RCP<const Matrix> A,
00068   Teuchos::RCP<Vector>       X,
00069   Teuchos::RCP<const Vector> B )
00070   : matrixA_(createConstMatrixAdapter<Matrix>(A))
00071   , multiVecX_(X) // may be null
00072   , multiVecB_(B) // may be null
00073   , globalNumRows_(matrixA_->getGlobalNumRows())
00074   , globalNumCols_(matrixA_->getGlobalNumCols())
00075   , globalNumNonZeros_(matrixA_->getGlobalNNZ())
00076   , rowIndexBase_(matrixA_->getRowIndexBase())
00077   , columnIndexBase_(matrixA_->getColumnIndexBase())
00078   , rank_(Teuchos::rank(*this->getComm()))
00079   , root_(rank_ == 0)
00080   , nprocs_(Teuchos::size(*this->getComm()))
00081 {
00082     TEUCHOS_TEST_FOR_EXCEPTION(
00083     !matrixShapeOK(),
00084     std::invalid_argument,
00085     "Matrix shape inappropriate for this solver");
00086 }
00087 
00088 
00090 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00091 SolverCore<ConcreteSolver,Matrix,Vector>::~SolverCore( )
00092 {
00093   // Nothing
00094 }
00095 
00096 
00097 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00098 Solver<Matrix,Vector>&
00099 SolverCore<ConcreteSolver,Matrix,Vector>::preOrdering()
00100 {
00101 #ifdef HAVE_AMESOS2_TIMERS
00102   Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
00103 #endif
00104 
00105   loadA(PREORDERING);
00106 
00107   static_cast<solver_type*>(this)->preOrdering_impl();
00108   ++status_.numPreOrder_;
00109   status_.last_phase_ = PREORDERING;
00110 
00111   return *this;
00112 }
00113 
00114 
00115 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00116 Solver<Matrix,Vector>&
00117 SolverCore<ConcreteSolver,Matrix,Vector>::symbolicFactorization()
00118 {
00119 #ifdef HAVE_AMESOS2_TIMERS
00120   Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
00121 #endif
00122 
00123   if( !status_.preOrderingDone() ){
00124     preOrdering();
00125     if( !matrix_loaded_ ) loadA(SYMBFACT);
00126   } else {
00127     loadA(SYMBFACT);
00128   }
00129 
00130   static_cast<solver_type*>(this)->symbolicFactorization_impl();
00131   ++status_.numSymbolicFact_;
00132   status_.last_phase_ = SYMBFACT;
00133 
00134   return *this;
00135 }
00136 
00137 
00138 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00139 Solver<Matrix,Vector>&
00140 SolverCore<ConcreteSolver,Matrix,Vector>::numericFactorization()
00141 {
00142 #ifdef HAVE_AMESOS2_TIMERS
00143   Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
00144 #endif
00145 
00146   if( !status_.symbolicFactorizationDone() ){
00147     symbolicFactorization();
00148     if( !matrix_loaded_ ) loadA(NUMFACT);
00149   } else {
00150     loadA(NUMFACT);
00151   }
00152 
00153   static_cast<solver_type*>(this)->numericFactorization_impl();
00154   ++status_.numNumericFact_;
00155   status_.last_phase_ = NUMFACT;
00156 
00157   return *this;
00158 }
00159 
00160 
00161 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00162 void
00163 SolverCore<ConcreteSolver,Matrix,Vector>::solve()
00164 {
00165   solve(multiVecX_.ptr(), multiVecB_.ptr());
00166 }
00167 
00168 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00169 void
00170 SolverCore<ConcreteSolver,Matrix,Vector>::solve(const Teuchos::Ptr<Vector> X,
00171                                                 const Teuchos::Ptr<const Vector> B) const
00172 {
00173 #ifdef HAVE_AMESOS2_TIMERS
00174   Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
00175 #endif
00176 
00177   X.assert_not_null();
00178   B.assert_not_null();
00179 
00180   const Teuchos::RCP<MultiVecAdapter<Vector> > x =
00181     createMultiVecAdapter<Vector>(Teuchos::rcpFromPtr(X));
00182   const Teuchos::RCP<const MultiVecAdapter<Vector> > b =
00183     createConstMultiVecAdapter<Vector>(Teuchos::rcpFromPtr(B));
00184 
00185 #ifdef HAVE_AMESOS2_DEBUG
00186   // Check some required properties of X and B
00187   TEUCHOS_TEST_FOR_EXCEPTION(x->getGlobalLength() != matrixA_->getGlobalNumCols(),
00188                      std::invalid_argument,
00189                      "MultiVector X must have length equal to the number of "
00190                      "global columns in A");
00191 
00192   TEUCHOS_TEST_FOR_EXCEPTION(b->getGlobalLength() != matrixA_->getGlobalNumRows(),
00193                      std::invalid_argument,
00194                      "MultiVector B must have length equal to the number of "
00195                      "global rows in A");
00196 
00197   TEUCHOS_TEST_FOR_EXCEPTION(x->getGlobalNumVectors() != b->getGlobalNumVectors(),
00198                      std::invalid_argument,
00199                      "X and B MultiVectors must have the same number of vectors");
00200 #endif  // HAVE_AMESOS2_DEBUG
00201 
00202   if( !status_.numericFactorizationDone() ){
00203     // This casting-away of constness is probably OK because this
00204     // function is meant to be "logically const"
00205     const_cast<type*>(this)->numericFactorization();
00206   }
00207 
00208   static_cast<const solver_type*>(this)->solve_impl(Teuchos::outArg(*x), Teuchos::ptrInArg(*b));
00209   ++status_.numSolve_;
00210   status_.last_phase_ = SOLVE;
00211 }
00212 
00213 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00214 void
00215 SolverCore<ConcreteSolver,Matrix,Vector>::solve(Vector* X, const Vector* B) const
00216 {
00217   solve(Teuchos::ptr(X), Teuchos::ptr(B));
00218 }
00219 
00220 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00221 bool
00222 SolverCore<ConcreteSolver,Matrix,Vector>::matrixShapeOK()
00223 {
00224 #ifdef HAVE_AMESOS2_TIMERS
00225   Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
00226 #endif
00227 
00228   return( static_cast<solver_type*>(this)->matrixShapeOK_impl() );
00229 }
00230 
00231   // The RCP should probably be to a const Matrix, but that throws a
00232   // wrench in some of the traits mechanisms that aren't expecting
00233   // const types
00234 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00235 void
00236 SolverCore<ConcreteSolver,Matrix,Vector>::setA( const Teuchos::RCP<const Matrix> a,
00237                                                 EPhase keep_phase )
00238 {
00239   matrixA_ = createConstMatrixAdapter(a);
00240 
00241 #ifdef HAVE_AMESOS2_DEBUG
00242   TEUCHOS_TEST_FOR_EXCEPTION( (keep_phase != CLEAN) &&
00243                       (globalNumRows_ != matrixA_->getGlobalNumRows() ||
00244                        globalNumCols_ != matrixA_->getGlobalNumCols()),
00245                       std::invalid_argument,
00246                       "Dimensions of new matrix be the same as the old matrix if "
00247                       "keeping any solver phase" );
00248 #endif
00249 
00250   status_.last_phase_ = keep_phase;
00251 
00252   // Reset phase counters
00253   switch( status_.last_phase_ ){
00254   case CLEAN:
00255     status_.numPreOrder_ = 0;
00256   case PREORDERING:
00257     status_.numSymbolicFact_ = 0;
00258   case SYMBFACT:
00259     status_.numNumericFact_ = 0;
00260   case NUMFACT:                 // probably won't ever happen by itself
00261     status_.numSolve_ = 0;
00262   case SOLVE:                   // probably won't ever happen
00263     break;
00264   }
00265 
00266   // Re-get the matrix dimensions in case they have changed
00267   globalNumNonZeros_ = matrixA_->getGlobalNNZ();
00268   globalNumCols_     = matrixA_->getGlobalNumCols();
00269   globalNumRows_     = matrixA_->getGlobalNumRows();
00270 }
00271 
00272 
00273 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00274 Solver<Matrix,Vector>&
00275 SolverCore<ConcreteSolver,Matrix,Vector>::setParameters(
00276   const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
00277 {
00278 #ifdef HAVE_AMESOS2_TIMERS
00279   Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
00280 #endif
00281 
00282   if( parameterList->name() == "Amesos2" ){
00283     // First, validate the parameterList
00284     Teuchos::RCP<const Teuchos::ParameterList> valid_params = getValidParameters();
00285     parameterList->validateParameters(*valid_params);
00286 
00287     // Do everything here that is for generic status and control parameters
00288     control_.setControlParameters(parameterList);
00289 
00290     // Finally, hook to the implementation's parameter list parser
00291     // First check if there is a dedicated sublist for this solver and use that if there is
00292     if( parameterList->isSublist(name()) ){
00293       // Have control look through the solver's parameter list to see if
00294       // there is anything it recognizes (mostly the "Transpose" parameter)
00295       control_.setControlParameters(Teuchos::sublist(parameterList, name()));
00296 
00297       static_cast<solver_type*>(this)->setParameters_impl(Teuchos::sublist(parameterList, name()));
00298     }
00299   }
00300 
00301   return *this;
00302 }
00303 
00304 
00305 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00306 Teuchos::RCP<const Teuchos::ParameterList>
00307 SolverCore<ConcreteSolver,Matrix,Vector>::getValidParameters() const
00308 {
00309 #ifdef HAVE_AMESOS2_TIMERS
00310   Teuchos::TimeMonitor LocalTimer1( timers_.totalTime_ );
00311 #endif
00312 
00313   using Teuchos::ParameterList;
00314   using Teuchos::RCP;
00315   using Teuchos::rcp;
00316 
00317   RCP<ParameterList> control_params = rcp(new ParameterList("Amesos2 Control"));
00318   control_params->set("Transpose", false, "Whether to solve with the matrix transpose");
00319   //  control_params->set("AddToDiag", "");
00320   //  control_params->set("AddZeroToDiag", false);
00321   //  control_params->set("MatrixProperty", "general");
00322   //  control_params->set("Reindex", false);
00323 
00324   RCP<const ParameterList>
00325     solver_params = static_cast<const solver_type*>(this)->getValidParameters_impl();
00326   // inject the "Transpose" parameter into the solver's valid parameters
00327   Teuchos::rcp_const_cast<ParameterList>(solver_params)->set("Transpose", false,
00328                    "Whether to solve with the "
00329                    "matrix transpose");
00330 
00331   RCP<ParameterList> amesos2_params = rcp(new ParameterList("Amesos2"));
00332   amesos2_params->setParameters(*control_params);
00333   amesos2_params->set(name(), *solver_params);
00334 
00335   return amesos2_params;
00336 }
00337 
00338 
00339 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00340 std::string
00341 SolverCore<ConcreteSolver,Matrix,Vector>::description() const
00342 {
00343   std::ostringstream oss;
00344   oss << name() << " solver interface";
00345   return oss.str();
00346 }
00347 
00348 
00349 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00350 void
00351 SolverCore<ConcreteSolver,Matrix,Vector>::describe(
00352   Teuchos::FancyOStream &out,
00353   const Teuchos::EVerbosityLevel verbLevel) const
00354 {
00355   if( matrixA_.is_null() || (rank_ != 0) ){ return; }
00356   using std::endl;
00357   using std::setw;
00358   using Teuchos::VERB_DEFAULT;
00359   using Teuchos::VERB_NONE;
00360   using Teuchos::VERB_LOW;
00361   using Teuchos::VERB_MEDIUM;
00362   using Teuchos::VERB_HIGH;
00363   using Teuchos::VERB_EXTREME;
00364   Teuchos::EVerbosityLevel vl = verbLevel;
00365   if (vl == VERB_DEFAULT) vl = VERB_LOW;
00366   Teuchos::RCP<const Teuchos::Comm<int> > comm = this->getComm();
00367   size_t width = 1;
00368   for( size_t dec = 10; dec < globalNumRows_; dec *= 10 ) {
00369     ++width;
00370   }
00371   width = std::max<size_t>(width,size_t(11)) + 2;
00372   Teuchos::OSTab tab(out);
00373   //    none: print nothing
00374   //     low: print O(1) info from node 0
00375   //  medium: print O(P) info, num entries per node
00376   //    high: print O(N) info, num entries per row
00377   // extreme: print O(NNZ) info: print indices and values
00378   //
00379   // for medium and higher, print constituent objects at specified verbLevel
00380   if( vl != VERB_NONE ) {
00381     std::string p = name();
00382     Util::printLine(out);
00383     out << this->description() << std::endl << std::endl;
00384 
00385     out << p << "Matrix has " << globalNumRows_ << "rows"
00386         << " and " << globalNumNonZeros_ << "nonzeros"
00387         << std::endl;
00388     if( vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME ){
00389       out << p << "Nonzero elements per row = "
00390           << globalNumNonZeros_ / globalNumRows_
00391           << std::endl;
00392       out << p << "Percentage of nonzero elements = "
00393           << 100.0 * globalNumNonZeros_ / (globalNumRows_ * globalNumCols_)
00394           << std::endl;
00395     }
00396     if( vl == VERB_HIGH || vl == VERB_EXTREME ){
00397       out << p << "Use transpose = " << control_.useTranspose_
00398           << std::endl;
00399     }
00400     if ( vl == VERB_EXTREME ){
00401       printTiming(out,vl);
00402     }
00403     Util::printLine(out);
00404   }
00405 }
00406 
00407 
00408 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00409 void
00410 SolverCore<ConcreteSolver,Matrix,Vector>::printTiming(
00411   Teuchos::FancyOStream &out,
00412   const Teuchos::EVerbosityLevel verbLevel) const
00413 {
00414   if( matrixA_.is_null() || (rank_ != 0) ){ return; }
00415 
00416   double preTime  = timers_.preOrderTime_.totalElapsedTime();
00417   double symTime  = timers_.symFactTime_.totalElapsedTime();
00418   double numTime  = timers_.numFactTime_.totalElapsedTime();
00419   double solTime  = timers_.solveTime_.totalElapsedTime();
00420   double totTime  = timers_.totalTime_.totalElapsedTime();
00421   double overhead = totTime - (preTime + symTime + numTime + solTime);
00422 
00423   std::string p = name() + " : ";
00424   Util::printLine(out);
00425 
00426   out << p << "Time to convert matrix to implementation format = "
00427       << timers_.mtxConvTime_.totalElapsedTime() << " (s)"
00428       << std::endl;
00429   out << p << "Time to redistribute matrix = "
00430       << timers_.mtxRedistTime_.totalElapsedTime() << " (s)"
00431       << std::endl;
00432 
00433   out << p << "Time to convert vectors to implementation format = "
00434       << timers_.vecConvTime_.totalElapsedTime() << " (s)"
00435       << std::endl;
00436   out << p << "Time to redistribute vectors = "
00437       << timers_.vecRedistTime_.totalElapsedTime() << " (s)"
00438       << std::endl;
00439 
00440   out << p << "Number of pre-orderings = "
00441       << status_.getNumPreOrder()
00442       << std::endl;
00443   out << p << "Time for pre-ordering = "
00444       << preTime << " (s), avg = "
00445       << preTime / status_.getNumPreOrder() << " (s)"
00446       << std::endl;
00447 
00448   out << p << "Number of symbolic factorizations = "
00449       << status_.getNumSymbolicFact()
00450       << std::endl;
00451   out << p << "Time for sym fact = "
00452       << symTime << " (s), avg = "
00453       << symTime / status_.getNumSymbolicFact() << " (s)"
00454       << std::endl;
00455 
00456   out << p << "Number of numeric factorizations = "
00457       << status_.getNumNumericFact()
00458       << std::endl;
00459   out << p << "Time for num fact = "
00460       << numTime << " (s), avg = "
00461       << numTime / status_.getNumNumericFact() << " (s)"
00462       << std::endl;
00463 
00464   out << p << "Number of solve phases = "
00465       << status_.getNumSolve()
00466       << std::endl;
00467   out << p << "Time for solve = "
00468       << solTime << " (s), avg = "
00469       << solTime / status_.getNumSolve() << " (s)"
00470       << std::endl;
00471 
00472   out << p << "Total time spent in Amesos2 = "
00473       << totTime << " (s)"
00474       << std::endl;
00475   out << p << "Total time spent in the Amesos2 interface = "
00476       << overhead << " (s)"
00477       << std::endl;
00478   out << p << "  (the above time does not include solver time)"
00479       << std::endl;
00480   out << p << "Amesos2 interface time / total time = "
00481       << overhead / totTime
00482       << std::endl;
00483   Util::printLine(out);
00484 }
00485 
00486 
00487 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00488 void
00489 SolverCore<ConcreteSolver,Matrix,Vector>::getTiming(
00490   Teuchos::ParameterList& timingParameterList) const
00491 {}
00492 
00493 
00494 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00495 std::string
00496 SolverCore<ConcreteSolver,Matrix,Vector>::name() const
00497 {
00498   std::string solverName = solver_type::name;
00499   return solverName;
00500 }
00501 
00502 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00503 void
00504 SolverCore<ConcreteSolver,Matrix,Vector>::loadA(EPhase current_phase)
00505 {
00506   matrix_loaded_ = static_cast<solver_type*>(this)->loadA_impl(current_phase);
00507 }
00508 
00509 
00510 } // end namespace Amesos2
00511 
00512 #endif  // AMESOS2_SOLVERCORE_DEF_HPP