Kokkos Node API and Local Linear Algebra Kernels Version of the Day
Kokkos_DefaultRelaxation.hpp
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //          Kokkos: Node API and Parallel Node Kernels
00005 //              Copyright (2008) Sandia Corporation
00006 // 
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 // 
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00038 // 
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 #ifndef KOKKOS_DEFAULTRELAXATION_HPP
00043 #define KOKKOS_DEFAULTRELAXATION_HPP
00044 
00045 #include <stdio.h>
00046 #include <stdexcept>
00047 
00048 #include <Teuchos_ArrayRCP.hpp>
00049 #include <Teuchos_DataAccess.hpp>
00050 #include <Teuchos_Assert.hpp>
00051 #include <Teuchos_TypeNameTraits.hpp>
00052 #include <Teuchos_ScalarTraits.hpp>
00053 
00054 #include "Kokkos_ConfigDefs.hpp"
00055 #include "Kokkos_MultiVector.hpp"
00056 #include "Kokkos_NodeHelpers.hpp"
00057 #include "Kokkos_DefaultArithmetic.hpp"
00058 #include "Kokkos_DefaultRelaxationKernelOps.hpp"
00059 #include "Kokkos_DefaultSparseOps.hpp"
00060 
00061 namespace Kokkos {
00062 
00070   template <class Scalar, class Ordinal, class Node>
00071   class DefaultRelaxation {
00072   public:
00073     typedef Scalar  ScalarType;
00074     typedef Ordinal OrdinalType;
00075     typedef Node    NodeType;
00076 
00078 
00080 
00082     DefaultRelaxation(const RCP<Node> &node = DefaultNode::getDefaultNode());
00083 
00085     ~DefaultRelaxation();
00086 
00088 
00090 
00091     
00093     RCP<Node> getNode() const;
00094 
00096 
00098 
00100 
00105     template <class UNSUPPORTED_GRAPH, class UNSUPPORTED_MATRIX>
00106     void initializeData(const RCP<const UNSUPPORTED_GRAPH > &graph, 
00107                         const RCP<const UNSUPPORTED_MATRIX> &matrix);
00108 
00110     void initializeData(const RCP<const DefaultCrsGraph<Ordinal,Node>         > &graph, 
00111                         const RCP<const DefaultCrsMatrix<Scalar,Ordinal,Node> > &matrix);
00112 
00117     void setDiagonal(MultiVector<Scalar,Node> & diag);    
00118 
00120     void clear();
00121 
00123 
00125 
00127 
00129 
00131     void sweep_jacobi(Scalar dampingFactor_, MultiVector<Scalar,Node> &X, const MultiVector<Scalar,Node> &B) const;
00132 
00133 #ifdef ENABLE_ALL_OTHER_RELAXATION
00134 
00135     void sweep_fine_hybrid(Scalar dampingFactor_, MultiVector<Scalar,Node> &X, const MultiVector<Scalar,Node> &B) const;
00136 
00138     void sweep_coarse_hybrid(Scalar dampingFactor_, size_t num_chunks, MultiVector<Scalar,Node> &X, const MultiVector<Scalar,Node> &B) const;
00139 
00140 #endif //ifdef ENABLE_ALL_OTHER_RELAXATION
00141 
00142     void setup_chebyshev(const Scalar lambda_max, const Scalar lambda_min);
00143 
00145     void sweep_chebyshev(MultiVector<Scalar,Node> &X, const MultiVector<Scalar,Node> &B) const;
00146 
00148 
00149   protected:
00151     DefaultRelaxation(const DefaultRelaxation& source);
00152 
00154     // NTS: This should eventually disappear into some other Kokkos class...
00155     void ExtractDiagonal();    
00156 
00158     bool UpdateJacobiTemp(size_t num_vectors, size_t vec_leng) const;
00159 
00161     bool UpdateChebyTemp(size_t num_vectors, size_t vec_leng) const;
00162 
00164     RCP<Node> node_;
00165 
00166     // we do this one of two ways: 
00167     // packed: array of offsets, ordinals, values. obviously the smallest footprint.
00168     ArrayRCP<const Ordinal> inds_;
00169     ArrayRCP<const size_t>  ptrs_;
00170     ArrayRCP<const Scalar>  vals_;
00171 
00173     ArrayRCP<Scalar> diagonal_;
00174 
00176     mutable ArrayRCP<Scalar> tempJacobiVector_;
00177     mutable size_t lastNumJacobiVectors_;
00178 
00179     // Arrays containing temp vectors for Chebyshev
00180     mutable ArrayRCP<Scalar> tempChebyVectorX_;
00181     mutable ArrayRCP<Scalar> tempChebyVectorW_;
00182     mutable size_t lastNumChebyVectors_;
00183     mutable bool cheby_setup_done_;
00184     mutable bool first_cheby_iteration_;
00185 
00186     // Constants for Chebyshev
00187     mutable Scalar lmin_,lmax_,delta_,s1_,oneOverTheta_,rho_,rho_new_,dtemp1_,dtemp2_;
00188 
00189     size_t numRows_;
00190     bool isInit_, isEmpty_;
00191   }; //class DefaultRelaxation declaration
00192 
00193   /**********************************************************************/
00194   template<class Scalar, class Ordinal, class Node>
00195   DefaultRelaxation<Scalar,Ordinal,Node>::DefaultRelaxation(const RCP<Node> &node)
00196   : node_(node)
00197   , lastNumJacobiVectors_(0)
00198   , lastNumChebyVectors_(0)
00199   , cheby_setup_done_(false)
00200   , first_cheby_iteration_(false)  
00201   , isInit_(false)
00202   , isEmpty_(false)
00203   {
00204     lmin_=lmax_=delta_=s1_=oneOverTheta_=rho_=rho_new_=dtemp1_=dtemp2_=Teuchos::ScalarTraits<Scalar>::zero();
00205   }
00206 
00207   /**********************************************************************/
00208   template<class Scalar, class Ordinal, class Node>
00209   DefaultRelaxation<Scalar,Ordinal,Node>::~DefaultRelaxation() {
00210     clear();
00211   }
00212 
00213   /**********************************************************************/
00214   template <class Scalar, class Ordinal, class Node>
00215   template <class UNSUPPORTED_GRAPH, class UNSUPPORTED_MATRIX>
00216   void DefaultRelaxation<Scalar,Ordinal,Node>::initializeData(
00217                         const RCP<const UNSUPPORTED_GRAPH > &graph,
00218                         const RCP<const UNSUPPORTED_MATRIX> &matrix) 
00219   {
00220     // not implemented for general graphs and matrices
00221     TEUCHOS_TEST_FOR_EXCEPT(true);
00222   }
00223 
00224   /**********************************************************************/
00225   template <class Scalar, class Ordinal, class Node>
00226   void DefaultRelaxation<Scalar,Ordinal,Node>::initializeData(
00227                 const RCP<const DefaultCrsGraph<Ordinal,Node>         > &graph, 
00228                 const RCP<const DefaultCrsMatrix<Scalar,Ordinal,Node> > &matrix)
00229   {
00230     std::string tfecfFuncName("initializeData()");
00231     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00232         isInit_ == true,
00233         std::runtime_error, " structure already initialized."
00234     )
00235     numRows_ = graph->getNumRows();
00236     if (graph->isEmpty() || numRows_ == 0) {
00237       isEmpty_ = true;
00238     }
00239     else {
00240       isEmpty_ = false;
00241       ptrs_ = graph->getPointers();
00242       inds_ = graph->getIndices();
00243       vals_ = matrix->getValues();
00244       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00245           inds_.size() != vals_.size(),
00246           std::runtime_error, " graph and matrix are not congruent."
00247       )
00248     }
00249     isInit_ = true;
00250     ExtractDiagonal();
00251   } //initializeData()
00252 
00253   /**********************************************************************/
00254   template <class Scalar, class Ordinal, class Node>
00255   RCP<Node> DefaultRelaxation<Scalar,Ordinal,Node>::getNode() const { 
00256     return node_; 
00257   }
00258 
00259   /**********************************************************************/
00260   template <class Scalar, class Ordinal, class Node>
00261   void DefaultRelaxation<Scalar,Ordinal,Node>::clear() {
00262     ptrs_ = null;
00263     inds_ = null;
00264     vals_ = null;
00265     diagonal_         = null;
00266     tempJacobiVector_ = null;
00267     tempChebyVectorW_ = null;
00268     tempChebyVectorX_ = null;
00269     isInit_   = false;
00270     isEmpty_  = false;
00271     lastNumJacobiVectors_=0;
00272     lastNumChebyVectors_=0;
00273   }
00274 
00275   /**********************************************************************/
00276   // Sets the diagonal inverted for relaxation using a MultiVector
00277   template <class Scalar, class Ordinal, class Node>
00278   void DefaultRelaxation<Scalar,Ordinal,Node>::setDiagonal(MultiVector<Scalar,Node> & diag){
00279     // Allocate space for diagonal
00280     diagonal_ = arcp<Scalar>(numRows_);    
00281 
00282     // NTS: Copy diag over
00283 
00284     // Make it fail for now...
00285     TEUCHOS_TEST_FOR_EXCEPT(true);
00286   }
00287     
00288   /**********************************************************************/
00289   template <class Scalar, class Ordinal, class Node>
00290   void DefaultRelaxation<Scalar,Ordinal,Node>::ExtractDiagonal(){    
00291     TEUCHOS_TEST_FOR_EXCEPTION(isInit_ == false, std::runtime_error, 
00292         Teuchos::typeName(*this) << "::ExtractDiagonal(): initializeValues() hasn't been called.");
00293 
00294     // Allocate space for diagonal
00295     diagonal_ = arcp<Scalar>(numRows_);    
00296       
00297     // Extract the diagonal 
00298     typedef ExtractDiagonalOp<Scalar,Ordinal>  Op;
00299     ReadyBufferHelper<Node> rbh(node_);
00300     Op wdp;
00301     rbh.begin();
00302     wdp.numRows = numRows_;
00303     wdp.ptrs    = rbh.template addConstBuffer<size_t>(ptrs_);
00304     wdp.inds    = rbh.template addConstBuffer<Ordinal>(inds_);
00305     wdp.vals    = rbh.template addConstBuffer<Scalar>(vals_);
00306     wdp.diag    = rbh.template addNonConstBuffer<Scalar>(diagonal_);
00307     rbh.end();
00308     node_->template parallel_for<Op>(0,numRows_,wdp);
00309   } //ExtractDiagonal()
00310 
00311   /**********************************************************************/
00312   template <class Scalar, class Ordinal, class Node>
00313   bool DefaultRelaxation<Scalar,Ordinal,Node>::UpdateJacobiTemp(size_t num_vectors, size_t vec_leng) const{
00314     // Re-allocate memory if needed
00315     if(num_vectors > lastNumJacobiVectors_){
00316       tempJacobiVector_= null;
00317       lastNumJacobiVectors_=num_vectors;
00318       tempJacobiVector_ = arcp<Scalar>(vec_leng*lastNumJacobiVectors_); 
00319       return true;
00320     }
00321     return false;
00322   } //UpdateJacobiTemp()
00323 
00324   /**********************************************************************/
00325   template <class Scalar, class Ordinal, class Node>
00326   bool DefaultRelaxation<Scalar,Ordinal,Node>::UpdateChebyTemp(size_t num_vectors, size_t vec_leng) const{
00327     // Re-allocate memory if needed
00328     if(num_vectors > lastNumChebyVectors_){
00329       tempChebyVectorW_= null;
00330       tempChebyVectorX_= null;
00331       lastNumChebyVectors_ = num_vectors;
00332       tempChebyVectorW_ = arcp<Scalar>(numRows_*lastNumChebyVectors_);
00333       tempChebyVectorX_ = arcp<Scalar>(vec_leng*lastNumChebyVectors_); 
00334       return true;
00335     }
00336     return false;
00337   } //UpdateChebyTemp()
00338 
00339 
00340 #ifdef ENABLE_ALL_OTHER_RELAXATION
00341   /**********************************************************************/
00342   template <class Scalar, class Ordinal, class Node>
00343   void DefaultRelaxation<Scalar,Ordinal,Node>::sweep_fine_hybrid(Scalar dampingFactor_,
00344                          MultiVector<Scalar,Node> &X, const MultiVector<Scalar,Node> &B) const{
00345     typedef DefaultFineGrainHybridGaussSeidelOp<Scalar,Ordinal>  Op;
00346 
00347     TEUCHOS_TEST_FOR_EXCEPTION(isInit_ == false, std::runtime_error,
00348         Teuchos::typeName(*this) << "::sweep_fine_hybrid(): operation not fully initialized.");
00349     TEUCHOS_TEST_FOR_EXCEPT(X.getNumCols() != B.getNumCols());
00350     ReadyBufferHelper<Node> rbh(node_);
00351     if (isEmpty_ == true) {
00352       // This makes no sense to try to call ...
00353       TEUCHOS_TEST_FOR_EXCEPT(true);
00354     }
00355     else {
00356       Op wdp;
00357       rbh.begin();
00358       wdp.numRows = numRows_;
00359       wdp.ptrs    = rbh.template addConstBuffer<size_t>(ptrs_);
00360       wdp.inds    = rbh.template addConstBuffer<Ordinal>(inds_);
00361       wdp.vals    = rbh.template addConstBuffer<Scalar>(vals_);
00362       wdp.x       = rbh.template addNonConstBuffer<Scalar>(X.getValuesNonConst());
00363       wdp.b       = rbh.template addConstBuffer<Scalar>(B.getValues());
00364       wdp.diag    = rbh.template addConstBuffer<Scalar>(diagonal_);
00365       wdp.damping_factor = dampingFactor_;
00366       wdp.xstride = X.getStride();
00367       wdp.bstride = B.getStride();
00368       rbh.end();
00369       const size_t numRHS = X.getNumCols();
00370       node_->template parallel_for<Op>(0,numRows_*numRHS,wdp);
00371     }
00372     return;
00373   } //sweep_fine_hybrid()
00374 
00375   /**********************************************************************/
00376   template <class Scalar, class Ordinal, class Node>
00377   void DefaultRelaxation<Scalar,Ordinal,Node>::sweep_coarse_hybrid(Scalar dampingFactor_,size_t num_chunks,
00378                          MultiVector<Scalar,Node> &X, const MultiVector<Scalar,Node> &B) const{
00379     typedef DefaultCoarseGrainHybridGaussSeidelOp<Scalar,Ordinal>  Op;
00380 
00381     TEUCHOS_TEST_FOR_EXCEPTION(isInit_ == false, std::runtime_error,
00382         Teuchos::typeName(*this) << "::sweep_coarse_hybrid(): operation not fully initialized.");
00383     TEUCHOS_TEST_FOR_EXCEPT(X.getNumCols() != B.getNumCols());
00384     ReadyBufferHelper<Node> rbh(node_);   
00385 
00386     if (isEmpty_ == true) {
00387       // This makes no sense to try to call ...
00388       TEUCHOS_TEST_FOR_EXCEPT(true);
00389     }
00390     else {
00391       Op wdp;
00392       rbh.begin();
00393       wdp.numRows = numRows_;
00394       wdp.numChunks = num_chunks;
00395       wdp.inds    = rbh.template addConstBuffer<Ordinal>(inds_);
00396       wdp.vals    = rbh.template addConstBuffer<Scalar>(vals_);
00397       wdp.x       = rbh.template addNonConstBuffer<Scalar>(X.getValuesNonConst());
00398       wdp.b       = rbh.template addConstBuffer<Scalar>(B.getValues());
00399       wdp.diag    = rbh.template addConstBuffer<Scalar>(diagonal_);
00400       wdp.damping_factor = dampingFactor_;
00401       wdp.xstride = X.getStride();
00402       wdp.bstride = B.getStride();
00403       rbh.end();
00404       const size_t numRHS = X.getNumCols();
00405       node_->template parallel_for<Op>(0,num_chunks*numRHS,wdp);
00406     }
00407     return;
00408   } //sweep_coarse_hybrid()
00409 #endif //ifdef ENABLE_ALL_OTHER_RELAXATION
00410 
00411   
00412   /********************************************************************/
00413   template <class Scalar, class Ordinal, class Node>
00414   void DefaultRelaxation<Scalar,Ordinal,Node>::sweep_jacobi(Scalar dampingFactor_,
00415                          MultiVector<Scalar,Node> &X, const MultiVector<Scalar,Node> &B) const{
00416     typedef DefaultJacobiOp<Scalar,Ordinal>  Op;
00417 
00418     TEUCHOS_TEST_FOR_EXCEPTION(isInit_ == false, std::runtime_error,
00419         Teuchos::typeName(*this) << "::sweep_jacobi(): operation not fully initialized.");
00420     TEUCHOS_TEST_FOR_EXCEPT(X.getNumCols() != B.getNumCols());
00421     ReadyBufferHelper<Node> rbh(node_);
00422 
00423     // Copy x over to the temp vector
00424     // NTS: The MultiVector copy constructor is a View. We need to do this the hard way.
00425     MultiVector<Scalar,Node> X0(X.getNode());
00426     UpdateJacobiTemp(X.getNumCols(),X.getNumRows());
00427     X0.initializeValues(X.getNumRows(),B.getNumCols(),tempJacobiVector_,B.getStride());
00428     DefaultArithmetic<MultiVector<Scalar,Node> >::Assign(X0,X);
00429 
00430     if (isEmpty_ == true) {
00431       // This makes no sense to try to call ...
00432       TEUCHOS_TEST_FOR_EXCEPT(true);
00433     }
00434     else {
00435       Op wdp;
00436       rbh.begin();
00437       wdp.numRows = numRows_;
00438       wdp.ptrs    = rbh.template addConstBuffer<size_t>(ptrs_);
00439       wdp.inds    = rbh.template addConstBuffer<Ordinal>(inds_);
00440       wdp.vals    = rbh.template addConstBuffer<Scalar>(vals_);
00441       wdp.diag    = rbh.template addConstBuffer<Scalar>(diagonal_);
00442       wdp.x       = rbh.template addNonConstBuffer<Scalar>(X.getValuesNonConst());
00443       wdp.x0      = rbh.template addConstBuffer<Scalar>(X0.getValues());
00444       wdp.b       = rbh.template addConstBuffer<Scalar>(B.getValues());
00445       wdp.damping_factor = dampingFactor_;
00446       wdp.xstride = X.getStride();
00447       wdp.bstride = B.getStride();
00448       rbh.end();
00449       const size_t numRHS = X.getNumCols();
00450       node_->template parallel_for<Op>(0,numRows_*numRHS,wdp);
00451     }
00452     return;
00453   } //sweep_jacobi()
00454 
00455   /********************************************************************/
00456   template <class Scalar, class Ordinal, class Node>
00457   void DefaultRelaxation<Scalar,Ordinal,Node>::setup_chebyshev(const Scalar lambda_max, const Scalar lambda_min){
00458     lmax_=lambda_max;
00459     lmin_=lambda_min;
00460 
00461     Scalar alpha = lmin_;
00462     Scalar beta  = 1.1*lmax_;
00463     delta_ = 2.0 / (beta-alpha);
00464     Scalar theta_ = 0.5 * (beta+alpha);
00465     s1_    = theta_ * delta_;
00466     oneOverTheta_=1.0 / theta_;
00467     rho_=1.0/s1_;
00468 
00469     first_cheby_iteration_=true;
00470     cheby_setup_done_=true;
00471   } //setup_chebyshev()
00472 
00473 
00474   /********************************************************************/
00475   template <class Scalar, class Ordinal, class Node>
00476   void DefaultRelaxation<Scalar,Ordinal,Node>::sweep_chebyshev(MultiVector<Scalar,Node> &X, const MultiVector<Scalar,Node> &B) const{
00477     typedef DefaultChebyshevOp<Scalar,Ordinal>  Op;
00478 
00479     TEUCHOS_TEST_FOR_EXCEPTION(isInit_ == false, std::runtime_error,
00480                Teuchos::typeName(*this) << "::sweep_chebyshev(): operation not fully initialized.");
00481     TEUCHOS_TEST_FOR_EXCEPT(X.getNumCols() != B.getNumCols());
00482     //size_t xstride = X.getStride();
00483     //size_t bstride = B.getStride();
00484     //TEUCHOS_TEST_FOR_EXCEPT(xstride != bstride); //TODO JJH don't think we want this b/c X is imported, B is local
00485     TEUCHOS_TEST_FOR_EXCEPT(cheby_setup_done_ == false);
00486     
00487     ReadyBufferHelper<Node> rbh(node_);
00488     
00489     // If the number of multivectors has changed, we need to run Cheby from scratch
00490     first_cheby_iteration_= UpdateChebyTemp(X.getNumCols(),X.getNumRows()) || first_cheby_iteration_;
00491     
00492     // Copy X to X0 so we get the matvec correct
00493     MultiVector<Scalar,Node> X0(X.getNode());
00494     X0.initializeValues(X.getNumRows(),X.getNumCols(),tempChebyVectorX_,X.getStride());
00495     DefaultArithmetic<MultiVector<Scalar,Node> >::Assign(X0,X);
00496 
00497     // Update Scalars 
00498     if(!first_cheby_iteration_){
00499       rho_new_ = 1.0 / (2.0*s1_ - rho_);
00500       dtemp1_=rho_*rho_new_;
00501       dtemp2_=2*rho_new_*delta_;
00502       rho_=rho_new_;
00503     }    
00504     
00505     if (isEmpty_ == true) {
00506       // This makes no sense to try to call ...
00507       TEUCHOS_TEST_FOR_EXCEPT(true);
00508     }
00509     else {
00510       Op wdp;
00511       rbh.begin();
00512       wdp.numRows = numRows_;
00513       wdp.ptrs    = rbh.template addConstBuffer<size_t>(ptrs_);
00514       wdp.inds    = rbh.template addConstBuffer<Ordinal>(inds_);
00515       wdp.vals    = rbh.template addConstBuffer<Scalar>(vals_);
00516       wdp.x       = rbh.template addNonConstBuffer<Scalar>(X.getValuesNonConst());
00517       wdp.x0      = rbh.template addConstBuffer<Scalar>(tempChebyVectorX_);
00518       wdp.w       = rbh.template addNonConstBuffer<Scalar>(tempChebyVectorW_);
00519       wdp.b       = rbh.template addConstBuffer<Scalar>(B.getValues());
00520       wdp.diag    = rbh.template addConstBuffer<Scalar>(diagonal_);
00521       wdp.stride = X.getStride();
00522       wdp.first_step=first_cheby_iteration_;
00523       wdp.zero_initial_guess=false;//HAQ Fix me later
00524       wdp.oneOverTheta = oneOverTheta_;
00525       wdp.dtemp1 = dtemp1_;
00526       wdp.dtemp2 = dtemp2_;
00527       
00528       rbh.end();    
00529       const size_t numRHS = X.getNumCols();
00530       node_->template parallel_for<Op>(0,numRows_*numRHS,wdp);
00531     }
00532     
00533     first_cheby_iteration_=false;
00534 
00535   } //sweep_chebyshev()
00536 
00537 } // namespace Kokkos
00538 
00539 #endif // KOKKOS_DEFAULTRELAXATION_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends