WIP (Do Not Merge!!): Belos Small Dense Matrix Abstraction

packages/belos/epetra/src/Belos_Details_Epetra_registerSolverFactory.cpp

@@ -37,20 +37,21 @@ void registerSolverFactory() {
   typedef double ST;
   typedef Epetra_MultiVector MV;
   typedef Epetra_Operator OP;
+  typedef Teuchos::SerialDenseMatrix<int,ST> DM;
 
   Impl::registerSolverSubclassForTypes<BiCGStabSolMgr<ST,MV,OP>, ST, MV, OP> ("BICGSTAB");
-  Impl::registerSolverSubclassForTypes<BlockCGSolMgr<ST,MV,OP>, ST, MV, OP> ("BLOCK CG");
-  Impl::registerSolverSubclassForTypes<BlockGmresSolMgr<ST,MV,OP>, ST, MV, OP> ("BLOCK GMRES");
+  Impl::registerSolverSubclassForTypes<BlockCGSolMgr<ST,MV,OP,DM>, ST, MV, OP> ("BLOCK CG");
+  Impl::registerSolverSubclassForTypes<BlockGmresSolMgr<ST,MV,OP,DM>, ST, MV, OP> ("BLOCK GMRES");
   Impl::registerSolverSubclassForTypes<FixedPointSolMgr<ST,MV,OP>, ST, MV, OP> ("FIXED POINT");
-  Impl::registerSolverSubclassForTypes<GCRODRSolMgr<ST,MV,OP>, ST, MV, OP> ("GCRODR");
-  Impl::registerSolverSubclassForTypes<GmresPolySolMgr<ST,MV,OP>, ST, MV, OP> ("HYBRID BLOCK GMRES");
-  Impl::registerSolverSubclassForTypes<LSQRSolMgr<ST,MV,OP>, ST, MV, OP> ("LSQR");
-  Impl::registerSolverSubclassForTypes<MinresSolMgr<ST,MV,OP>, ST, MV, OP> ("MINRES");
-  Impl::registerSolverSubclassForTypes<PCPGSolMgr<ST,MV,OP>, ST, MV, OP> ("PCPG");
-  Impl::registerSolverSubclassForTypes<PseudoBlockCGSolMgr<ST,MV,OP>, ST, MV, OP> ("PSEUDOBLOCK CG");
-  Impl::registerSolverSubclassForTypes<PseudoBlockGmresSolMgr<ST,MV,OP>, ST, MV, OP> ("PSEUDOBLOCK GMRES");
+  Impl::registerSolverSubclassForTypes<GCRODRSolMgr<ST,MV,OP,DM>, ST, MV, OP> ("GCRODR");
+  Impl::registerSolverSubclassForTypes<GmresPolySolMgr<ST,MV,OP,DM>, ST, MV, OP> ("HYBRID BLOCK GMRES");
+  Impl::registerSolverSubclassForTypes<LSQRSolMgr<ST,MV,OP,DM>, ST, MV, OP> ("LSQR");
+  Impl::registerSolverSubclassForTypes<MinresSolMgr<ST,MV,OP,DM>, ST, MV, OP> ("MINRES");
+  Impl::registerSolverSubclassForTypes<PCPGSolMgr<ST,MV,OP,DM>, ST, MV, OP> ("PCPG");
+  Impl::registerSolverSubclassForTypes<PseudoBlockCGSolMgr<ST,MV,OP,DM>, ST, MV, OP> ("PSEUDOBLOCK CG");
+  Impl::registerSolverSubclassForTypes<PseudoBlockGmresSolMgr<ST,MV,OP,DM>, ST, MV, OP> ("PSEUDOBLOCK GMRES");
   Impl::registerSolverSubclassForTypes<PseudoBlockTFQMRSolMgr<ST,MV,OP>, ST, MV, OP> ("PSEUDOBLOCK TFQMR");
-  Impl::registerSolverSubclassForTypes<RCGSolMgr<ST,MV,OP>, ST, MV, OP> ("RCG");
+  Impl::registerSolverSubclassForTypes<RCGSolMgr<ST,MV,OP,DM>, ST, MV, OP> ("RCG");
   Impl::registerSolverSubclassForTypes<TFQMRSolMgr<ST,MV,OP>, ST, MV, OP> ("TFQMR");
 }
 

packages/belos/epetra/test/MVOPTester/cxx_main.cpp

@@ -27,6 +27,7 @@
 #include "BelosMVOPTester.hpp"
 #include "BelosEpetraAdapter.hpp"
 #include "BelosOutputManager.hpp"
+#include "BelosTeuchosDenseAdapter.hpp"
 
 #include "Teuchos_StandardCatchMacros.hpp"
 
@@ -121,7 +122,7 @@ int main(int argc, char *argv[])
 
     // Issue several useful typedefs;
     typedef Belos::MultiVec<double> EMV;
-    //typedef Belos::Operator<double> EOP; // unused
+    typedef Teuchos::SerialDenseMatrix<int, double> TDM;
 
     // Create an Epetra_MultiVector for an initial std::vector to start the solver.
     // Note that this needs to have the same number of columns as the blocksize.
@@ -135,7 +136,7 @@ int main(int argc, char *argv[])
     }
 
     // test the Epetra adapter multivector
-    ierr = Belos::TestMultiVecTraits<double,EMV>(MyOM,ivec);
+    ierr = Belos::TestMultiVecTraits<double,EMV,TDM>(MyOM,ivec);
     gerr &= ierr;
     if (ierr) {
       MyOM->print(Belos::Warnings,"*** EpetraAdapter PASSED TestMultiVecTraits()\n");

packages/belos/epetra/test/OrthoManager/belos_orthomanager_epetra.cpp

@@ -315,7 +315,7 @@ main (int argc, char *argv[])
   // subclass, given a name for the subclass.  The name is not the
   // same as the class' syntactic name: e.g., "DKGS" is the name of
   // DkgsOrthoManager.
-  Belos::OrthoManagerFactory< scalar_type, MV, OP > factory;
+  Belos::OrthoManagerFactory< scalar_type, MV, OP, serial_matrix_type > factory;
   // The name of the (Mat)OrthoManager subclass to instantiate.
   std::string ortho (factory.defaultName());
 
@@ -528,7 +528,7 @@ main (int argc, char *argv[])
     // should return zero).
     int numFailed = 0;
     {
-      typedef Belos::Test::OrthoManagerTester<scalar_type, MV> tester_type;
+      typedef Belos::Test::OrthoManagerTester<scalar_type, MV, serial_matrix_type> tester_type;
       debugOut << "Running OrthoManager tests..." << endl;
       numFailed = tester_type::runTests (OM, isRankRevealing, S,
           sizeX1, sizeX2, outMan);

packages/belos/kokkos/example/GCRODRComplexKokkosExFile.cpp

@@ -17,6 +17,7 @@
 #include "BelosGCRODRSolMgr.hpp"
 #include "BelosOutputManager.hpp"
 
+#include "Teuchos_SerialDenseMatrix.hpp"
 #include "Teuchos_CommandLineProcessor.hpp"
 #include "Teuchos_ParameterList.hpp"
 #include "Teuchos_StandardCatchMacros.hpp"
@@ -41,11 +42,12 @@ bool success = true;
   typedef Kokkos::complex<double>           ST;
   typedef int                               OT;
   typedef Kokkos::DefaultExecutionSpace     EXSP;
-  typedef Teuchos::ScalarTraits<ST>        SCT;
+  typedef Teuchos::ScalarTraits<ST>         SCT;
+  typedef Teuchos::SerialDenseMatrix<int,ST>DM;
   typedef SCT::magnitudeType                MT;
-  typedef Belos::KokkosMultiVec<ST, EXSP>         MV;
-  typedef Belos::MultiVec<ST> KMV;
-  typedef Belos::Operator<ST> KOP; 
+  typedef Belos::KokkosMultiVec<ST, EXSP>   MV;
+  typedef Belos::MultiVec<ST>               KMV;
+  typedef Belos::Operator<ST>               KOP; 
   typedef Belos::MultiVecTraits<ST,KMV>     MVT;
   typedef Belos::OperatorTraits<ST,KMV,KOP>  OPT;
 
@@ -136,8 +138,8 @@ try {
   // *******************************************************************
   //
   // Create an iterative solver manager.
-  RCP< Belos::SolverManager<ST,KMV,KOP> > newSolver
-    = rcp( new Belos::GCRODRSolMgr<ST,KMV,KOP,true>(rcpFromRef(problem), rcpFromRef(belosList)) );
+  RCP< Belos::SolverManager<ST,KMV,KOP,DM> > newSolver
+    = rcp( new Belos::GCRODRSolMgr<ST,KMV,KOP,DM,true>(rcpFromRef(problem), rcpFromRef(belosList)) );
 
   //
   // **********Print out information about problem*******************

packages/belos/kokkos/example/GCRODRKokkosExFile.cpp

@@ -137,7 +137,7 @@ try {
   //
   // Create an iterative solver manager.
   RCP< Belos::SolverManager<ST,KMV,KOP> > newSolver
-    = rcp( new Belos::GCRODRSolMgr<ST,KMV,KOP,true>(rcpFromRef(problem), rcpFromRef(belosList)) );
+    = rcp( new Belos::GCRODRSolMgr<ST,KMV,KOP>(rcpFromRef(problem), rcpFromRef(belosList)) );
 
   //
   // **********Print out information about problem*******************

packages/belos/src/BelosBiCGStabIter.hpp

@@ -25,8 +25,6 @@
 #include "BelosOperatorTraits.hpp"
 #include "BelosMultiVecTraits.hpp"
 
-#include "Teuchos_SerialDenseMatrix.hpp"
-#include "Teuchos_SerialDenseVector.hpp"
 #include "Teuchos_ScalarTraits.hpp"
 #include "Teuchos_ParameterList.hpp"
 #include "Teuchos_TimeMonitor.hpp"
@@ -77,15 +75,15 @@ namespace Belos {
     }
   };
 
-  template<class ScalarType, class MV, class OP>
-  class BiCGStabIter : virtual public Iteration<ScalarType,MV,OP> {
+  template<class ScalarType, class MV, class OP, class DM>
+  class BiCGStabIter : virtual public Iteration<ScalarType,MV,OP,DM> {
 
   public:
 
     //
     // Convenience typedefs
     //
-    typedef MultiVecTraits<ScalarType,MV> MVT;
+    typedef MultiVecTraits<ScalarType,MV,DM> MVT;
     typedef OperatorTraits<ScalarType,MV,OP> OPT;
     typedef Teuchos::ScalarTraits<ScalarType> SCT;
     typedef typename SCT::magnitudeType MagnitudeType;
@@ -101,7 +99,7 @@ namespace Belos {
      */
     BiCGStabIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
                           const Teuchos::RCP<OutputManager<ScalarType> > &printer,
-                          const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
+                          const Teuchos::RCP<StatusTest<ScalarType,MV,OP,DM> > &tester,
                           Teuchos::ParameterList &params );
 
     //! Destructor.
@@ -235,7 +233,7 @@ namespace Belos {
     //
     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
-    const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
+    const Teuchos::RCP<StatusTest<ScalarType,MV,OP,DM> >       stest_;
 
     //
     // Algorithmic parameters
@@ -277,10 +275,10 @@ namespace Belos {
 
   //////////////////////////////////////////////////////////////////////////////////////////////////
   // Constructor.
-  template<class ScalarType, class MV, class OP>
-  BiCGStabIter<ScalarType,MV,OP>::BiCGStabIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
+  template<class ScalarType, class MV, class OP, class DM>
+  BiCGStabIter<ScalarType,MV,OP,DM>::BiCGStabIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
                                                                const Teuchos::RCP<OutputManager<ScalarType> > &printer,
-                                                               const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
+                                                               const Teuchos::RCP<StatusTest<ScalarType,MV,OP,DM> > &tester,
                                                                Teuchos::ParameterList &/* params */ ):
     lp_(problem),
     om_(printer),
@@ -295,8 +293,8 @@ namespace Belos {
 
   //////////////////////////////////////////////////////////////////////////////////////////////////
   // Initialize this iteration object
-  template <class ScalarType, class MV, class OP>
-  void BiCGStabIter<ScalarType,MV,OP>::initializeBiCGStab(BiCGStabIterationState<ScalarType,MV>& newstate)
+  template<class ScalarType, class MV, class OP, class DM>
+  void BiCGStabIter<ScalarType,MV,OP,DM>::initializeBiCGStab(BiCGStabIterationState<ScalarType,MV>& newstate)
   {
     // Check if there is any multivector to clone from.
     Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
@@ -425,8 +423,8 @@ namespace Belos {
 
   //////////////////////////////////////////////////////////////////////////////////////////////////
   // Iterate until the status test informs us we should stop.
-  template <class ScalarType, class MV, class OP>
-  void BiCGStabIter<ScalarType,MV,OP>::iterate()
+  template<class ScalarType, class MV, class OP, class DM>
+  void BiCGStabIter<ScalarType,MV,OP,DM>::iterate()
   {
     using Teuchos::RCP;
 
@@ -588,8 +586,8 @@ namespace Belos {
 
   //////////////////////////////////////////////////////////////////////////////////////////////////
   // Iterate until the status test informs us we should stop.
-  template <class ScalarType, class MV, class OP>
-  void BiCGStabIter<ScalarType,MV,OP>::axpy(const ScalarType alpha, const MV & A,
+  template<class ScalarType, class MV, class OP, class DM>
+  void BiCGStabIter<ScalarType,MV,OP,DM>::axpy(const ScalarType alpha, const MV & A,
                                             const std::vector<ScalarType> beta, const MV& B, MV& mv, bool minus)
   {
     Teuchos::RCP<const MV> A1, B1;

packages/belos/src/BelosBiCGStabSolMgr.hpp

@@ -26,6 +26,8 @@
 #include "BelosStatusTestCombo.hpp"
 #include "BelosStatusTestOutputFactory.hpp"
 #include "BelosOutputManager.hpp"
+
+#include "BelosTeuchosDenseAdapter.hpp"
 #ifdef BELOS_TEUCHOS_TIME_MONITOR
 #include "Teuchos_TimeMonitor.hpp"
 #endif
@@ -61,11 +63,11 @@ namespace Belos {
     BiCGStabSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
     {}};
 
-  template<class ScalarType, class MV, class OP>
-  class BiCGStabSolMgr : public SolverManager<ScalarType,MV,OP> {
+  template<class ScalarType, class MV, class OP, class DM = Teuchos::SerialDenseMatrix<int, ScalarType>>
+  class BiCGStabSolMgr : public SolverManager<ScalarType,MV,OP,DM> {
 
   private:
-    typedef MultiVecTraits<ScalarType,MV> MVT;
+    typedef MultiVecTraits<ScalarType,MV,DM> MVT;
     typedef OperatorTraits<ScalarType,MV,OP> OPT;
     typedef Teuchos::ScalarTraits<ScalarType> SCT;
     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
@@ -99,8 +101,8 @@ namespace Belos {
     virtual ~BiCGStabSolMgr() {};
 
     //! clone for Inverted Injection (DII)
-    Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
-      return Teuchos::rcp(new BiCGStabSolMgr<ScalarType,MV,OP>);
+    Teuchos::RCP<SolverManager<ScalarType, MV, OP, DM> > clone () const override {
+      return Teuchos::rcp(new BiCGStabSolMgr<ScalarType,MV,OP,DM>);
     }
     //@}
 
@@ -216,10 +218,10 @@ namespace Belos {
     Teuchos::RCP<std::ostream> outputStream_;
 
     // Status test.
-    Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
-    Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
-    Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
-    Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
+    Teuchos::RCP<StatusTest<ScalarType,MV,OP,DM> > sTest_;
+    Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP,DM> > maxIterTest_;
+    Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP,DM> > convTest_;
+    Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP,DM> > outputTest_;
 
     // Current parameter list.
     Teuchos::RCP<Teuchos::ParameterList> params_;
@@ -257,8 +259,8 @@ namespace Belos {
   };
 
 // Empty Constructor
-template<class ScalarType, class MV, class OP>
-BiCGStabSolMgr<ScalarType,MV,OP>::BiCGStabSolMgr() :
+template<class ScalarType, class MV, class OP, class DM>
+BiCGStabSolMgr<ScalarType,MV,OP,DM>::BiCGStabSolMgr() :
   outputStream_(Teuchos::rcpFromRef(std::cout)),
   convtol_(DefaultSolverParameters::convTol),
   maxIters_(maxIters_default_),
@@ -274,8 +276,8 @@ BiCGStabSolMgr<ScalarType,MV,OP>::BiCGStabSolMgr() :
 {}
 
 // Basic Constructor
-template<class ScalarType, class MV, class OP>
-BiCGStabSolMgr<ScalarType,MV,OP>::
+template<class ScalarType, class MV, class OP, class DM>
+BiCGStabSolMgr<ScalarType,MV,OP,DM>::
 BiCGStabSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
                      const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
   problem_(problem),
@@ -304,8 +306,8 @@ BiCGStabSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
   }
 }
 
-template<class ScalarType, class MV, class OP>
-void BiCGStabSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
+template<class ScalarType, class MV, class OP, class DM>
+void BiCGStabSolMgr<ScalarType,MV,OP,DM>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
 {
   using Teuchos::ParameterList;
   using Teuchos::parameterList;
@@ -400,8 +402,8 @@ void BiCGStabSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos
   }
 
   // Convergence
-  typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
-  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
+  typedef Belos::StatusTestCombo<ScalarType,MV,OP,DM>  StatusTestCombo_t;
+  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP,DM>  StatusTestResNorm_t;
 
   // Check for convergence tolerance
   if (params->isParameter("Convergence Tolerance")) {
@@ -482,7 +484,7 @@ void BiCGStabSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos
 
   // Basic test checks maximum iterations and native residual.
   if (maxIterTest_ == Teuchos::null)
-    maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
+    maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP,DM>( maxIters_ ) );
 
   // Implicit residual test, using the native residual to determine if convergence was achieved.
   if (convTest_ == Teuchos::null || newResTest) {
@@ -497,7 +499,7 @@ void BiCGStabSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos
 
     // Create the status test output class.
     // This class manages and formats the output from the status test.
-    StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
+    StatusTestOutputFactory<ScalarType,MV,OP,DM> stoFactory( outputStyle_ );
     outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
 
     // Set the solver string for the output test
@@ -519,9 +521,9 @@ void BiCGStabSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos
 }
 
 
-template<class ScalarType, class MV, class OP>
+template<class ScalarType, class MV, class OP, class DM>
 Teuchos::RCP<const Teuchos::ParameterList>
-BiCGStabSolMgr<ScalarType,MV,OP>::getValidParameters() const
+BiCGStabSolMgr<ScalarType,MV,OP,DM>::getValidParameters() const
 {
   using Teuchos::ParameterList;
   using Teuchos::parameterList;
@@ -575,8 +577,8 @@ BiCGStabSolMgr<ScalarType,MV,OP>::getValidParameters() const
 }
 
 
-template<class ScalarType, class MV, class OP>
-ReturnType BiCGStabSolMgr<ScalarType,MV,OP>::solve ()
+template<class ScalarType, class MV, class OP, class DM>
+ReturnType BiCGStabSolMgr<ScalarType,MV,OP,DM>::solve ()
 {
   // Set the current parameters if they were not set before.
   // NOTE:  This may occur if the user generated the solver manager with the default constructor and
@@ -624,8 +626,8 @@ ReturnType BiCGStabSolMgr<ScalarType,MV,OP>::solve ()
   //////////////////////////////////////////////////////////////////////////////////////
   // Pseudo-Block BiCGStab solver
 
-  Teuchos::RCP<BiCGStabIter<ScalarType,MV,OP> > bicgstab_iter
-    = Teuchos::rcp( new BiCGStabIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) );
+  Teuchos::RCP<BiCGStabIter<ScalarType,MV,OP,DM> > bicgstab_iter
+    = Teuchos::rcp( new BiCGStabIter<ScalarType,MV,OP,DM>(problem_,printer_,outputTest_,plist) );
 
   // Enter solve() iterations
   {
@@ -669,7 +671,7 @@ ReturnType BiCGStabSolMgr<ScalarType,MV,OP>::solve ()
           if ( convTest_->getStatus() == Passed ) {
 
             // Figure out which linear systems converged.
-            std::vector<int> convIdx = Teuchos::rcp_dynamic_cast<StatusTestGenResNorm<ScalarType,MV,OP> >(convTest_)->convIndices();
+            std::vector<int> convIdx = Teuchos::rcp_dynamic_cast<StatusTestGenResNorm<ScalarType,MV,OP,DM> >(convTest_)->convIndices();
 
             // If the number of converged linear systems is equal to the
             // number of current linear systems, then we are done with this block.
@@ -823,8 +825,8 @@ ReturnType BiCGStabSolMgr<ScalarType,MV,OP>::solve ()
 }
 
 //  This method requires the solver manager to return a std::string that describes itself.
-template<class ScalarType, class MV, class OP>
-std::string BiCGStabSolMgr<ScalarType,MV,OP>::description() const
+template<class ScalarType, class MV, class OP, class DM>
+std::string BiCGStabSolMgr<ScalarType,MV,OP,DM>::description() const
 {
   std::ostringstream oss;
   oss << "Belos::BiCGStabSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";