diff --git a/packages/belos/src/BelosPCPGIter.hpp b/packages/belos/src/BelosPCPGIter.hpp
index 4ab3034ff1e5..50f085c89c61 100644
--- a/packages/belos/src/BelosPCPGIter.hpp
+++ b/packages/belos/src/BelosPCPGIter.hpp
@@ -55,9 +55,9 @@
 #include "BelosStatusTest.hpp"
 #include "BelosOperatorTraits.hpp"
 #include "BelosMultiVecTraits.hpp"
+#include "BelosDenseMatTraits.hpp"
 #include "BelosCGIteration.hpp"
 
-#include "Teuchos_SerialDenseMatrix.hpp"
 #include "Teuchos_ScalarTraits.hpp"
 #include "Teuchos_ParameterList.hpp"
 #include "Teuchos_TimeMonitor.hpp"
@@ -66,7 +66,7 @@
   \class Belos::PCPGIter
   
   \brief This class implements the PCPG iteration, where a
-  single-std::vector Krylov subspace is constructed.  The documentation
+  single-vector Krylov subspace is constructed.  The documentation
   refers to blocks, but note that at this point, all blocks have unit
   dimension.
  
@@ -82,7 +82,7 @@ namespace Belos {
    *
    * The structure is utilized by initialize() and getState().
    */
-  template <class ScalarType, class MV>
+  template <class ScalarType, class MV, class DM>
   struct PCPGIterState {
     /*! \brief The current dimension of the reduction.
      *
@@ -115,7 +115,7 @@ namespace Belos {
      *
      * The \c curDim by \c curDim D = diag(P'*AP) = U' * C
      */
-    Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > D;
+    Teuchos::RCP<const DM> D;
 
     PCPGIterState() : curDim(0), 
                       prevUdim(0), 
@@ -128,16 +128,17 @@ namespace Belos {
   
   //@}
   
-  template<class ScalarType, class MV, class OP>
-  class PCPGIter : virtual public Iteration<ScalarType,MV,OP> {
+  template<class ScalarType, class MV, class OP, class DM>
+  class PCPGIter : 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 DenseMatTraits<ScalarType,DM>    DMT;
     typedef Teuchos::ScalarTraits<ScalarType> SCT;
     typedef typename SCT::magnitudeType MagnitudeType;
     
@@ -153,8 +154,8 @@ namespace Belos {
      */
     PCPGIter( 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<MatOrthoManager<ScalarType,MV,OP> > &ortho,
+		const Teuchos::RCP<StatusTest<ScalarType,MV,OP,DM> > &tester,
+		const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP,DM> > &ortho,
 		Teuchos::ParameterList &params );
     
     //! Destructor.
@@ -203,14 +204,14 @@ namespace Belos {
      * \note For any pointer in \c newstate which directly points to the multivectors in 
      * the solver, the data is not (supposed to be) copied.
      */
-    void initialize(PCPGIterState<ScalarType,MV>& newstate);
+    void initialize(PCPGIterState<ScalarType,MV,DM>& newstate);
     
     /*! \brief Initialize the solver with the initial vectors from the linear problem.
      *  An exception is thrown if initialzed is called and newstate.R is null.
      */
     void initialize()
     {
-      PCPGIterState<ScalarType,MV> empty;
+      PCPGIterState<ScalarType,MV,DM> empty;
       initialize(empty);
     }
     
@@ -221,8 +222,8 @@ namespace Belos {
      * \returns A PCPGIterState object containing const pointers to the current
      * solver state.
      */
-    PCPGIterState<ScalarType,MV> getState() const {
-      PCPGIterState<ScalarType,MV> state;
+    PCPGIterState<ScalarType,MV,DM> getState() const {
+      PCPGIterState<ScalarType,MV,DM> state;
       state.Z = Z_;         // CG state
       state.P = P_;
       state.AP = AP_;
@@ -316,8 +317,8 @@ 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<OrthoManager<ScalarType,MV> >        ortho_;
+    const Teuchos::RCP<StatusTest<ScalarType,MV,OP,DM> >    stest_;
+    const Teuchos::RCP<OrthoManager<ScalarType,MV,DM> >     ortho_;
 
     //
     // Algorithmic parameters
@@ -375,16 +376,16 @@ namespace Belos {
     //
     // Projected matrices
     // D_ : Diagonal matrix of pivots D = P'AP 
-    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
+    Teuchos::RCP<DM> D_;
   };
   
   //////////////////////////////////////////////////////////////////////////////////////////////////
   // Constructor.
-  template<class ScalarType, class MV, class OP>
-  PCPGIter<ScalarType,MV,OP>::PCPGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
+  template<class ScalarType, class MV, class OP, class DM>
+  PCPGIter<ScalarType,MV,OP,DM>::PCPGIter(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<MatOrthoManager<ScalarType,MV,OP> > &ortho,
+					   const Teuchos::RCP<StatusTest<ScalarType,MV,OP,DM> > &tester,
+					   const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP,DM> > &ortho,
 					   Teuchos::ParameterList &params ):
     lp_(problem),
     om_(printer),
@@ -417,8 +418,8 @@ namespace Belos {
   
   //////////////////////////////////////////////////////////////////////////////////////////////////
   // Set the block size and adjust as necessary
-  template <class ScalarType, class MV, class OP>
-  void PCPGIter<ScalarType,MV,OP>::setSize( int savedBlocks )
+  template<class ScalarType, class MV, class OP, class DM>
+  void PCPGIter<ScalarType,MV,OP,DM>::setSize( int savedBlocks )
   {
     // allocate space only; perform no computation
     // Any change in size invalidates the state of the solver as implemented here.
@@ -437,8 +438,8 @@ namespace Belos {
 
   //////////////////////////////////////////////////////////////////////////////////////////////////
   // Enable the reuse of a single solver object for completely different linear systems
-  template <class ScalarType, class MV, class OP>
-  void PCPGIter<ScalarType,MV,OP>::resetState()
+  template<class ScalarType, class MV, class OP, class DM>
+  void PCPGIter<ScalarType,MV,OP,DM>::resetState()
   {
       stateStorageInitialized_ = false;
       initialized_ = false;
@@ -449,8 +450,8 @@ namespace Belos {
 
   //////////////////////////////////////////////////////////////////////////////////////////////////
   // Setup the state storage.  Called by either initialize or, if savedBlocks_ changes, setSize.
-  template <class ScalarType, class MV, class OP>
-  void PCPGIter<ScalarType,MV,OP>::setStateSize ()
+  template<class ScalarType, class MV, class OP, class DM> 
+  void PCPGIter<ScalarType,MV,OP,DM>::setStateSize ()
   {
     if (!stateStorageInitialized_) {
 
@@ -514,14 +515,14 @@ namespace Belos {
 	}
         if (keepDiagonal_) {
           if (D_ == Teuchos::null) {
-            D_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
+            D_ = DMT::Create();
           }
           if (initDiagonal_) {
-            D_->shape( newsd, newsd );
+            DMT::Reshape( *D_, newsd, newsd, true );
           }
           else {
             if (D_->numRows() < newsd || D_->numCols() < newsd) {
-              D_->shapeUninitialized( newsd, newsd );
+              DMT::Reshape( *D_, newsd, newsd, false );
             }
           }
         }
@@ -533,8 +534,8 @@ namespace Belos {
 
   //////////////////////////////////////////////////////////////////////////////////////////////////
   // Initialize the iteration object
-  template <class ScalarType, class MV, class OP>
-  void PCPGIter<ScalarType,MV,OP>::initialize(PCPGIterState<ScalarType,MV>& newstate)
+  template<class ScalarType, class MV, class OP, class DM>
+  void PCPGIter<ScalarType,MV,OP,DM>::initialize(PCPGIterState<ScalarType,MV,DM>& newstate)
   {
 
     TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
@@ -621,8 +622,8 @@ namespace Belos {
 
   //////////////////////////////////////////////////////////////////////////////////////////////////
   // Iterate until the status test informs us we should stop.
-  template <class ScalarType, class MV, class OP>
-  void PCPGIter<ScalarType,MV,OP>::iterate()
+  template<class ScalarType, class MV, class OP, class DM>
+  void PCPGIter<ScalarType,MV,OP,DM>::iterate()
   {
     //
     // Allocate/initialize data structures
@@ -633,10 +634,9 @@ namespace Belos {
     const bool debug = false;
 
     // Allocate memory for scalars.
-    Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
-    Teuchos::SerialDenseMatrix<int,ScalarType> pAp( 1, 1 );
-    Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
-    Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 );
+    ScalarType alpha, beta, rHz_old;
+    Teuchos::RCP<DM> pAp = DMT::Create(1,1);
+    Teuchos::RCP<DM> rHz = DMT::Create(1,1);
 
     if( iter_ != 0 )
       std::cout << " Iterate Warning: begin from nonzero iter_ ?" << std::endl;  //DMD
@@ -645,12 +645,12 @@ namespace Belos {
     std::vector<int> prevInd;
     Teuchos::RCP<const MV> Uprev;
     Teuchos::RCP<const MV> Cprev;
-    Teuchos::SerialDenseMatrix<int,ScalarType> CZ;
+    Teuchos::RCP<DM> CZ;
 
     if( prevUdim_ ){
       prevInd.resize( prevUdim_ );
       for( int i=0; i<prevUdim_ ; i++) prevInd[i] = i;
-      CZ.reshape( prevUdim_ , 1 );
+      CZ = DMT::Create( prevUdim_ , 1 );
       Uprev = MVT::CloneView(*U_, prevInd);
       Cprev = MVT::CloneView(*C_, prevInd);
     }
@@ -677,14 +677,9 @@ namespace Belos {
       Teuchos::RCP<MV> P; 
       curind[0] = curDim_ - 1;          // column = dimension - 1 
       P = MVT::CloneViewNonConst(*U_,curind); 
-      MVT::MvTransMv( one, *Cprev, *P, CZ );
-      MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P );       // P -= U*(C'Z)
+      MVT::MvTransMv( one, *Cprev, *P, *CZ );
+      MVT::MvTimesMatAddMv( -one, *Uprev, *CZ, one, *P );       // P -= U*(C'Z)
 
-      if( debug ){
-        MVT::MvTransMv( one, *Cprev, *P, CZ );
-        std::cout << " Input CZ post ortho " << std::endl;
-        CZ.print( std::cout );
-      }
       if( curDim_ == savedBlocks_ ){
         std::vector<int> zero_index(1);
         zero_index[0] = 0;
@@ -694,7 +689,8 @@ namespace Belos {
     }
 
     // Compute first <r,z> a.k.a. rHz
-    MVT::MvTransMv( one, *R_, *Z_, rHz );
+    MVT::MvTransMv( one, *R_, *Z_, *rHz );
+    DMT::SyncDeviceToHost( *rHz );
 
     ////////////////////////////////////////////////////////////////
     // iterate until the status test is satisfied
@@ -713,51 +709,52 @@ namespace Belos {
         P = MVT::CloneView(*U_,curind); 
         AP = MVT::CloneViewNonConst(*C_,curind); 
         lp_->applyOp( *P, *AP );
-        MVT::MvTransMv( one, *P, *AP, pAp );
+        MVT::MvTransMv( one, *P, *AP, *pAp );
       }else{
         if( prevUdim_ + iter_ == savedBlocks_ ){
           AP = MVT::CloneViewNonConst(*C_,curind); 
           lp_->applyOp( *P_, *AP );
-          MVT::MvTransMv( one, *P_, *AP, pAp );
+          MVT::MvTransMv( one, *P_, *AP, *pAp );
         }else{
           lp_->applyOp( *P_, *AP_ );
-          MVT::MvTransMv( one, *P_, *AP_, pAp );
+          MVT::MvTransMv( one, *P_, *AP_, *pAp );
         }
       }
+      DMT::SyncDeviceToHost( *pAp );
 
       if( keepDiagonal_  && prevUdim_ + iter_ <= savedBlocks_ )
-        (*D_)(iter_ -1 ,iter_ -1 ) = pAp(0,0);
+        DMT::Value( *D_, iter_-1 ,iter_-1 ) = DMT::ValueConst(*pAp,0,0);
 
       // positive pAp required 
-      TEUCHOS_TEST_FOR_EXCEPTION( pAp(0,0) <= zero, CGPositiveDefiniteFailure,
+      TEUCHOS_TEST_FOR_EXCEPTION( DMT::ValueConst(*pAp,0,0) <= zero, CGPositiveDefiniteFailure,
                           "Belos::PCPGIter::iterate(): non-positive value for p^H*A*p encountered!" );
 
       // alpha := <R_,Z_> / <P,AP>
-      alpha(0,0) = rHz(0,0) / pAp(0,0);
+      alpha = DMT::ValueConst(*rHz,0,0) / DMT::ValueConst(*pAp,0,0);
 
       // positive alpha required 
-      TEUCHOS_TEST_FOR_EXCEPTION( alpha(0,0) <= zero, CGPositiveDefiniteFailure,
+      TEUCHOS_TEST_FOR_EXCEPTION( alpha <= zero, CGPositiveDefiniteFailure,
                           "Belos::PCPGIter::iterate(): non-positive value for alpha encountered!" );
 
       // solution update  x += alpha * P
       if( curDim_ < savedBlocks_ ){
-         MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P, *cur_soln_vec );
+         MVT::MvAddMv( one, *cur_soln_vec, alpha, *P, *cur_soln_vec );
       }else{
-         MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
+         MVT::MvAddMv( one, *cur_soln_vec, alpha, *P_, *cur_soln_vec );
       }
       //lp_->updateSolution(); ... does nothing.
       //
       // The denominator of beta is saved before residual is updated [ old <R_, Z_> ].
       //
-      rHz_old(0,0) = rHz(0,0);
+      rHz_old = DMT::ValueConst(*rHz,0,0);
       //
       // residual update R_ := R_ - alpha * AP
       //
       if( prevUdim_ + iter_ <= savedBlocks_ ){
-         MVT::MvAddMv( one, *R_, -alpha(0,0), *AP, *R_ );
+         MVT::MvAddMv( one, *R_, -alpha, *AP, *R_ );
          AP = Teuchos::null;
       }else{
-         MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
+         MVT::MvAddMv( one, *R_, -alpha, *AP_, *R_ );
       }
       //
       // update beta := [ new <R_, Z_> ] / [ old <R_, Z_> ] and the search direction p.
@@ -768,23 +765,19 @@ namespace Belos {
         Z_ = R_;
       }
       //
-      MVT::MvTransMv( one, *R_, *Z_, rHz );
+      MVT::MvTransMv( one, *R_, *Z_, *rHz );
+      DMT::SyncDeviceToHost( *rHz );
       //
-      beta(0,0) = rHz(0,0) / rHz_old(0,0);
+      beta = DMT::ValueConst(*rHz,0,0) / rHz_old;
       //
       if( curDim_ < savedBlocks_ ){
          curDim_++;                                                         // update basis dim
          curind[0] = curDim_ - 1;
          Teuchos::RCP<MV> Pnext = MVT::CloneViewNonConst(*U_,curind);
-         MVT::MvAddMv( one, *Z_, beta(0,0), *P, *Pnext );
+         MVT::MvAddMv( one, *Z_, beta, *P, *Pnext );
          if( prevUdim_ ){ // Deflate seed space 
-             MVT::MvTransMv( one, *Cprev, *Z_, CZ );
-             MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *Pnext ); // Pnext -= U*(C'Z)
-             if( debug ){
-               std::cout << " Check CZ " << std::endl;
-               MVT::MvTransMv( one, *Cprev, *Pnext, CZ );
-               CZ.print( std::cout );
-             }
+             MVT::MvTransMv( one, *Cprev, *Z_, *CZ );
+             MVT::MvTimesMatAddMv( -one, *Uprev, *CZ, one, *Pnext ); // Pnext -= U*(C'Z)
          }
          P = Teuchos::null;
          if( curDim_ == savedBlocks_ ){
@@ -794,16 +787,10 @@ namespace Belos {
          }
          Pnext = Teuchos::null;
       }else{
-         MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
+         MVT::MvAddMv( one, *Z_, beta, *P_, *P_ );
          if( prevUdim_ ){ // Deflate seed space
-             MVT::MvTransMv( one, *Cprev, *Z_, CZ );
-             MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P_ );       // P_ -= U*(C'Z)
-
-             if( debug ){
-               std::cout << " Check CZ " << std::endl;
-               MVT::MvTransMv( one, *Cprev, *P_, CZ );
-               CZ.print( std::cout );
-             }
+             MVT::MvTransMv( one, *Cprev, *Z_, *CZ );
+             MVT::MvTimesMatAddMv( -one, *Uprev, *CZ, one, *P_ );       // P_ -= U*(C'Z)
          }
       }
       // CGB: 5/26/2010
diff --git a/packages/belos/src/BelosPCPGSolMgr.hpp b/packages/belos/src/BelosPCPGSolMgr.hpp
index 29e9d7112e9c..786518907b89 100644
--- a/packages/belos/src/BelosPCPGSolMgr.hpp
+++ b/packages/belos/src/BelosPCPGSolMgr.hpp
@@ -52,6 +52,7 @@
 #include "BelosSolverManager.hpp"
 
 #include "BelosPCPGIter.hpp"
+#include "BelosTeuchosDenseAdapter.hpp"
 
 #include "BelosOrthoManagerFactory.hpp"
 #include "BelosStatusTestMaxIters.hpp"
@@ -607,7 +608,7 @@ void PCPGSolMgr<ScalarType,MV,OP,DM,true>::setParameters( const Teuchos::RCP<Teu
     params_->set("Orthogonalization Constant",orthoKappa_);
     if (orthoType_=="DGKS") {
       if (orthoKappa_ > 0 && ortho_ != Teuchos::null && !changedOrthoType) {
-        Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
+        Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP,DM> >(ortho_)->setDepTol( orthoKappa_ );
       }
     }
   }
@@ -760,8 +761,8 @@ ReturnType PCPGSolMgr<ScalarType,MV,OP,DM,true>::solve() {
   //////////////////////////////////////////////////////////////////////////////////////
   // PCPG solver
 
-  Teuchos::RCP<PCPGIter<ScalarType,MV,OP> > pcpg_iter;
-  pcpg_iter = Teuchos::rcp( new PCPGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
+  Teuchos::RCP<PCPGIter<ScalarType,MV,OP,DM> > pcpg_iter;
+  pcpg_iter = Teuchos::rcp( new PCPGIter<ScalarType,MV,OP,DM>(problem_,printer_,outputTest_,ortho_,plist) );
   // Number of iterations required to generate initial recycle space (if needed)
 
   // Enter solve() iterations
@@ -835,7 +836,7 @@ ReturnType PCPGSolMgr<ScalarType,MV,OP,DM,true>::solve() {
 
 
       // Set the new state and initialize the solver.
-      PCPGIterState<ScalarType,MV> pcpgState; // fails if R == null.
+      PCPGIterState<ScalarType,MV,DM> pcpgState; // fails if R == null.
 
       pcpgState.R = R_;
       if( U_ != Teuchos::null ) pcpgState.U = U_;
@@ -908,7 +909,7 @@ ReturnType PCPGSolMgr<ScalarType,MV,OP,DM,true>::solve() {
       problem_->setCurrLS();
 
       // Get the state.   How did pcpgState die?
-      PCPGIterState<ScalarType,MV> oldState = pcpg_iter->getState();
+      PCPGIterState<ScalarType,MV,DM> oldState = pcpg_iter->getState();
 
       dimU_ = oldState.curDim;
       int q = oldState.prevUdim;
@@ -1059,7 +1060,7 @@ ReturnType PCPGSolMgr<ScalarType,MV,OP,DM,true>::solve() {
   // Save the convergence test value ("achieved tolerance") for this solve.
   {
     using Teuchos::rcp_dynamic_cast;
-    typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
+    typedef StatusTestGenResNorm<ScalarType,MV,OP,DM> conv_test_type;
     // testValues is nonnull and not persistent.
     const std::vector<MagnitudeType>* pTestValues =
       rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();