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   , rank_(Teuchos::rank(*this->getComm()))
00077   , root_(rank_ == 0)
00078   , nprocs_(Teuchos::size(*this->getComm()))
00079 {
00080     TEUCHOS_TEST_FOR_EXCEPTION(
00081     !matrixShapeOK(),
00082     std::invalid_argument,
00083     "Matrix shape inappropriate for this solver");
00084 }
00085 
00086 
00088 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00089 SolverCore<ConcreteSolver,Matrix,Vector>::~SolverCore( )
00090 {
00091   // Nothing
00092 }
00093 
00094 
00095 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00096 Solver<Matrix,Vector>&
00097 SolverCore<ConcreteSolver,Matrix,Vector>::preOrdering()
00098 {
00099 #ifdef HAVE_AMESOS2_TIMERS
00100   Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
00101 #endif
00102 
00103   loadA(PREORDERING);
00104 
00105   static_cast<solver_type*>(this)->preOrdering_impl();
00106   ++status_.numPreOrder_;
00107   status_.last_phase_ = PREORDERING;
00108 
00109   return *this;
00110 }
00111 
00112 
00113 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00114 Solver<Matrix,Vector>&
00115 SolverCore<ConcreteSolver,Matrix,Vector>::symbolicFactorization()
00116 {
00117 #ifdef HAVE_AMESOS2_TIMERS
00118   Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
00119 #endif
00120 
00121   if( !status_.preOrderingDone() ){
00122     preOrdering();
00123     if( !matrix_loaded_ ) loadA(SYMBFACT);
00124   } else {
00125     loadA(SYMBFACT);
00126   }
00127 
00128   static_cast<solver_type*>(this)->symbolicFactorization_impl();
00129   ++status_.numSymbolicFact_;
00130   status_.last_phase_ = SYMBFACT;
00131 
00132   return *this;
00133 }
00134 
00135 
00136 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00137 Solver<Matrix,Vector>&
00138 SolverCore<ConcreteSolver,Matrix,Vector>::numericFactorization()
00139 {
00140 #ifdef HAVE_AMESOS2_TIMERS
00141   Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
00142 #endif
00143 
00144   if( !status_.symbolicFactorizationDone() ){
00145     symbolicFactorization();
00146     if( !matrix_loaded_ ) loadA(NUMFACT);
00147   } else {
00148     loadA(NUMFACT);
00149   }
00150 
00151   static_cast<solver_type*>(this)->numericFactorization_impl();
00152   ++status_.numNumericFact_;
00153   status_.last_phase_ = NUMFACT;
00154 
00155   return *this;
00156 }
00157 
00158 
00159 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00160 void
00161 SolverCore<ConcreteSolver,Matrix,Vector>::solve()
00162 {
00163   solve(multiVecX_.ptr(), multiVecB_.ptr());
00164 }
00165 
00166 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00167 void
00168 SolverCore<ConcreteSolver,Matrix,Vector>::solve(const Teuchos::Ptr<Vector> X,
00169                                                 const Teuchos::Ptr<const Vector> B) const
00170 {
00171 #ifdef HAVE_AMESOS2_TIMERS
00172   Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
00173 #endif
00174 
00175   X.assert_not_null();
00176   B.assert_not_null();
00177 
00178   const Teuchos::RCP<MultiVecAdapter<Vector> > x =
00179     createMultiVecAdapter<Vector>(Teuchos::rcpFromPtr(X));
00180   const Teuchos::RCP<const MultiVecAdapter<Vector> > b =
00181     createConstMultiVecAdapter<Vector>(Teuchos::rcpFromPtr(B));
00182 
00183 #ifdef HAVE_AMESOS2_DEBUG
00184   // Check some required properties of X and B
00185   TEUCHOS_TEST_FOR_EXCEPTION(x->getGlobalLength() != matrixA_->getGlobalNumCols(),
00186                      std::invalid_argument,
00187                      "MultiVector X must have length equal to the number of "
00188                      "global columns in A");
00189 
00190   TEUCHOS_TEST_FOR_EXCEPTION(b->getGlobalLength() != matrixA_->getGlobalNumRows(),
00191                      std::invalid_argument,
00192                      "MultiVector B must have length equal to the number of "
00193                      "global rows in A");
00194 
00195   TEUCHOS_TEST_FOR_EXCEPTION(x->getGlobalNumVectors() != b->getGlobalNumVectors(),
00196                      std::invalid_argument,
00197                      "X and B MultiVectors must have the same number of vectors");
00198 #endif  // HAVE_AMESOS2_DEBUG
00199 
00200   if( !status_.numericFactorizationDone() ){
00201     // This casting-away of constness is probably OK because this
00202     // function is meant to be "logically const"
00203     const_cast<type*>(this)->numericFactorization();
00204   }
00205 
00206   static_cast<const solver_type*>(this)->solve_impl(Teuchos::outArg(*x), Teuchos::ptrInArg(*b));
00207   ++status_.numSolve_;
00208   status_.last_phase_ = SOLVE;
00209 }
00210 
00211 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00212 void
00213 SolverCore<ConcreteSolver,Matrix,Vector>::solve(Vector* X, const Vector* B) const
00214 {
00215   solve(Teuchos::ptr(X), Teuchos::ptr(B));
00216 }
00217 
00218 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00219 bool
00220 SolverCore<ConcreteSolver,Matrix,Vector>::matrixShapeOK()
00221 {
00222 #ifdef HAVE_AMESOS2_TIMERS
00223   Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
00224 #endif
00225 
00226   return( static_cast<solver_type*>(this)->matrixShapeOK_impl() );
00227 }
00228 
00229   // The RCP should probably be to a const Matrix, but that throws a
00230   // wrench in some of the traits mechanisms that aren't expecting
00231   // const types
00232 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00233 void
00234 SolverCore<ConcreteSolver,Matrix,Vector>::setA( const Teuchos::RCP<const Matrix> a,
00235                                                 EPhase keep_phase )
00236 {
00237   matrixA_ = createConstMatrixAdapter(a);
00238 
00239 #ifdef HAVE_AMESOS2_DEBUG
00240   TEUCHOS_TEST_FOR_EXCEPTION( (keep_phase != CLEAN) &&
00241                       (globalNumRows_ != matrixA_->getGlobalNumRows() ||
00242                        globalNumCols_ != matrixA_->getGlobalNumCols()),
00243                       std::invalid_argument,
00244                       "Dimensions of new matrix be the same as the old matrix if "
00245                       "keeping any solver phase" );
00246 #endif
00247 
00248   status_.last_phase_ = keep_phase;
00249 
00250   // Reset phase counters
00251   switch( status_.last_phase_ ){
00252   case CLEAN:
00253     status_.numPreOrder_ = 0;
00254   case PREORDERING:
00255     status_.numSymbolicFact_ = 0;
00256   case SYMBFACT:
00257     status_.numNumericFact_ = 0;
00258   case NUMFACT:                 // probably won't ever happen by itself
00259     status_.numSolve_ = 0;
00260   case SOLVE:                   // probably won't ever happen
00261     break;
00262   }
00263 
00264   // Re-get the matrix dimensions in case they have changed
00265   globalNumNonZeros_ = matrixA_->getGlobalNNZ();
00266   globalNumCols_     = matrixA_->getGlobalNumCols();
00267   globalNumRows_     = matrixA_->getGlobalNumRows();
00268 }
00269 
00270 
00271 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00272 Solver<Matrix,Vector>&
00273 SolverCore<ConcreteSolver,Matrix,Vector>::setParameters(
00274   const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
00275 {
00276 #ifdef HAVE_AMESOS2_TIMERS
00277   Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
00278 #endif
00279 
00280   if( parameterList->name() == "Amesos2" ){
00281     // First, validate the parameterList
00282     Teuchos::RCP<const Teuchos::ParameterList> valid_params = getValidParameters();
00283     parameterList->validateParameters(*valid_params);
00284 
00285     // Do everything here that is for generic status and control parameters
00286     control_.setControlParameters(parameterList);
00287 
00288     // Finally, hook to the implementation's parameter list parser
00289     // First check if there is a dedicated sublist for this solver and use that if there is
00290     if( parameterList->isSublist(name()) ){
00291       // Have control look through the solver's parameter list to see if
00292       // there is anything it recognizes (mostly the "Transpose" parameter)
00293       control_.setControlParameters(Teuchos::sublist(parameterList, name()));
00294 
00295       static_cast<solver_type*>(this)->setParameters_impl(Teuchos::sublist(parameterList, name()));
00296     }
00297   }
00298 
00299   return *this;
00300 }
00301 
00302 
00303 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00304 Teuchos::RCP<const Teuchos::ParameterList>
00305 SolverCore<ConcreteSolver,Matrix,Vector>::getValidParameters() const
00306 {
00307 #ifdef HAVE_AMESOS2_TIMERS
00308   Teuchos::TimeMonitor LocalTimer1( timers_.totalTime_ );
00309 #endif
00310 
00311   using Teuchos::ParameterList;
00312   using Teuchos::RCP;
00313   using Teuchos::rcp;
00314 
00315   RCP<ParameterList> control_params = rcp(new ParameterList("Amesos2 Control"));
00316   control_params->set("Transpose", false, "Whether to solve with the matrix transpose");
00317   //  control_params->set("AddToDiag", "");
00318   //  control_params->set("AddZeroToDiag", false);
00319   //  control_params->set("MatrixProperty", "general");
00320   //  control_params->set("Reindex", false);
00321 
00322   RCP<const ParameterList>
00323     solver_params = static_cast<const solver_type*>(this)->getValidParameters_impl();
00324   // inject the "Transpose" parameter into the solver's valid parameters
00325   Teuchos::rcp_const_cast<ParameterList>(solver_params)->set("Transpose", false,
00326                    "Whether to solve with the "
00327                    "matrix transpose");
00328 
00329   RCP<ParameterList> amesos2_params = rcp(new ParameterList("Amesos2"));
00330   amesos2_params->setParameters(*control_params);
00331   amesos2_params->set(name(), *solver_params);
00332 
00333   return amesos2_params;
00334 }
00335 
00336 
00337 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00338 std::string
00339 SolverCore<ConcreteSolver,Matrix,Vector>::description() const
00340 {
00341   std::ostringstream oss;
00342   oss << name() << " solver interface";
00343   return oss.str();
00344 }
00345 
00346 
00347 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00348 void
00349 SolverCore<ConcreteSolver,Matrix,Vector>::describe(
00350   Teuchos::FancyOStream &out,
00351   const Teuchos::EVerbosityLevel verbLevel) const
00352 {
00353   if( matrixA_.is_null() || (rank_ != 0) ){ return; }
00354   using std::endl;
00355   using std::setw;
00356   using Teuchos::VERB_DEFAULT;
00357   using Teuchos::VERB_NONE;
00358   using Teuchos::VERB_LOW;
00359   using Teuchos::VERB_MEDIUM;
00360   using Teuchos::VERB_HIGH;
00361   using Teuchos::VERB_EXTREME;
00362   Teuchos::EVerbosityLevel vl = verbLevel;
00363   if (vl == VERB_DEFAULT) vl = VERB_LOW;
00364   Teuchos::RCP<const Teuchos::Comm<int> > comm = this->getComm();
00365   size_t width = 1;
00366   for( size_t dec = 10; dec < globalNumRows_; dec *= 10 ) {
00367     ++width;
00368   }
00369   width = std::max<size_t>(width,11) + 2;
00370   Teuchos::OSTab tab(out);
00371   //    none: print nothing
00372   //     low: print O(1) info from node 0
00373   //  medium: print O(P) info, num entries per node
00374   //    high: print O(N) info, num entries per row
00375   // extreme: print O(NNZ) info: print indices and values
00376   //
00377   // for medium and higher, print constituent objects at specified verbLevel
00378   if( vl != VERB_NONE ) {
00379     std::string p = name();
00380     Util::printLine(out);
00381     out << this->description() << std::endl << std::endl;
00382 
00383     out << p << "Matrix has " << globalNumRows_ << "rows"
00384         << " and " << globalNumNonZeros_ << "nonzeros"
00385         << std::endl;
00386     if( vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME ){
00387       out << p << "Nonzero elements per row = "
00388           << globalNumNonZeros_ / globalNumRows_
00389           << std::endl;
00390       out << p << "Percentage of nonzero elements = "
00391           << 100.0 * globalNumNonZeros_ / (globalNumRows_ * globalNumCols_)
00392           << std::endl;
00393     }
00394     if( vl == VERB_HIGH || vl == VERB_EXTREME ){
00395       out << p << "Use transpose = " << control_.useTranspose_
00396           << std::endl;
00397     }
00398     if ( vl == VERB_EXTREME ){
00399       printTiming(out,vl);
00400     }
00401     Util::printLine(out);
00402   }
00403 }
00404 
00405 
00406 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00407 void
00408 SolverCore<ConcreteSolver,Matrix,Vector>::printTiming(
00409   Teuchos::FancyOStream &out,
00410   const Teuchos::EVerbosityLevel verbLevel) const
00411 {
00412   if( matrixA_.is_null() || (rank_ != 0) ){ return; }
00413 
00414   double preTime  = timers_.preOrderTime_.totalElapsedTime();
00415   double symTime  = timers_.symFactTime_.totalElapsedTime();
00416   double numTime  = timers_.numFactTime_.totalElapsedTime();
00417   double solTime  = timers_.solveTime_.totalElapsedTime();
00418   double totTime  = timers_.totalTime_.totalElapsedTime();
00419   double overhead = totTime - (preTime + symTime + numTime + solTime);
00420 
00421   std::string p = name() + " : ";
00422   Util::printLine(out);
00423 
00424   out << p << "Time to convert matrix to implementation format = "
00425       << timers_.mtxConvTime_.totalElapsedTime() << " (s)"
00426       << std::endl;
00427   out << p << "Time to redistribute matrix = "
00428       << timers_.mtxRedistTime_.totalElapsedTime() << " (s)"
00429       << std::endl;
00430 
00431   out << p << "Time to convert vectors to implementation format = "
00432       << timers_.vecConvTime_.totalElapsedTime() << " (s)"
00433       << std::endl;
00434   out << p << "Time to redistribute vectors = "
00435       << timers_.vecRedistTime_.totalElapsedTime() << " (s)"
00436       << std::endl;
00437 
00438   out << p << "Number of pre-orderings = "
00439       << status_.getNumPreOrder()
00440       << std::endl;
00441   out << p << "Time for pre-ordering = "
00442       << preTime << " (s), avg = "
00443       << preTime / status_.getNumPreOrder() << " (s)"
00444       << std::endl;
00445 
00446   out << p << "Number of symbolic factorizations = "
00447       << status_.getNumSymbolicFact()
00448       << std::endl;
00449   out << p << "Time for sym fact = "
00450       << symTime << " (s), avg = "
00451       << symTime / status_.getNumSymbolicFact() << " (s)"
00452       << std::endl;
00453 
00454   out << p << "Number of numeric factorizations = "
00455       << status_.getNumNumericFact()
00456       << std::endl;
00457   out << p << "Time for num fact = "
00458       << numTime << " (s), avg = "
00459       << numTime / status_.getNumNumericFact() << " (s)"
00460       << std::endl;
00461 
00462   out << p << "Number of solve phases = "
00463       << status_.getNumSolve()
00464       << std::endl;
00465   out << p << "Time for solve = "
00466       << solTime << " (s), avg = "
00467       << solTime / status_.getNumSolve() << " (s)"
00468       << std::endl;
00469 
00470   out << p << "Total time spent in Amesos2 = "
00471       << totTime << " (s)"
00472       << std::endl;
00473   out << p << "Total time spent in the Amesos2 interface = "
00474       << overhead << " (s)"
00475       << std::endl;
00476   out << p << "  (the above time does not include solver time)"
00477       << std::endl;
00478   out << p << "Amesos2 interface time / total time = "
00479       << overhead / totTime
00480       << std::endl;
00481   Util::printLine(out);
00482 }
00483 
00484 
00485 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00486 void
00487 SolverCore<ConcreteSolver,Matrix,Vector>::getTiming(
00488   Teuchos::ParameterList& timingParameterList) const
00489 {}
00490 
00491 
00492 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00493 std::string
00494 SolverCore<ConcreteSolver,Matrix,Vector>::name() const
00495 {
00496   std::string solverName = solver_type::name;
00497   return solverName;
00498 }
00499 
00500 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
00501 void
00502 SolverCore<ConcreteSolver,Matrix,Vector>::loadA(EPhase current_phase)
00503 {
00504   matrix_loaded_ = static_cast<solver_type*>(this)->loadA_impl(current_phase);
00505 }
00506 
00507 
00508 } // end namespace Amesos2
00509 
00510 #endif  // AMESOS2_SOLVERCORE_DEF_HPP