Kokkos Node API and Local Linear Algebra Kernels Version of the Day
Kokkos_CuspSparseOps.hpp
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //          Kokkos: Node API and Parallel Node Kernels
00005 //              Copyright (2009) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ************************************************************************
00027 //@HEADER
00028 
00029 #ifndef KOKKOS_CUSPSPARSEOPS_HPP
00030 #define KOKKOS_CUSPSPARSEOPS_HPP
00031 
00032 #include <Teuchos_ArrayRCP.hpp>
00033 #include <Teuchos_TestForException.hpp>
00034 #include <Teuchos_TypeNameTraits.hpp>
00035 #include <Teuchos_CompileTimeAssert.hpp>
00036 #include <stdexcept>
00037 
00038 #include "Kokkos_ConfigDefs.hpp"
00039 #include "Kokkos_CrsMatrix.hpp" 
00040 #include "Kokkos_CrsGraph.hpp" 
00041 #include "Kokkos_MultiVector.hpp"
00042 #include "Kokkos_DefaultArithmetic.hpp"
00043 
00044 #include <cusp/hyb_matrix.h>
00045 #include <cusp/csr_matrix.h>
00046 #include <cusp/device/detail/spmv.h>
00047 
00048 namespace Kokkos {
00049 
00053   template <class Scalar, class Ordinal, class Node>
00054   class CUSPSparseOps {
00055   public:
00057 
00058 
00060     typedef Scalar  ScalarType;
00062     typedef Ordinal OrdinalType;
00064     typedef Node    NodeType;
00065 
00070     template <class S2>
00071     struct rebind {
00072       typedef CUSPSparseOps<S2,Ordinal,Node> other;
00073     };
00074 
00076 
00077 
00078 
00080     CUSPSparseOps(const RCP<Node> &node);
00081 
00083     ~CUSPSparseOps();
00084 
00086 
00087 
00088     
00090     RCP<Node> getNode() const;
00091 
00093 
00094 
00095 
00097     void initializeStructure(const CrsGraphHostCompute<Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > &graph);
00098 
00100     void initializeValues(const CrsMatrixHostCompute<Scalar,Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > &matrix);
00101 
00103     void clear();
00104 
00106 
00107 
00108 
00110     template <class DomainScalar, class RangeScalar>
00111     void multiply(Teuchos::ETransp trans, RangeScalar alpha, const MultiVector<DomainScalar,Node> &X, MultiVector<RangeScalar,Node> &Y) const;
00112 
00114     template <class DomainScalar, class RangeScalar>
00115     void multiply(Teuchos::ETransp trans, 
00116                   RangeScalar alpha, const MultiVector<DomainScalar,Node> &X, RangeScalar beta, MultiVector<RangeScalar,Node> &Y) const;
00117 
00119     template <class DomainScalar, class RangeScalar>
00120     void solve(Teuchos::ETransp trans, Teuchos::EUplo uplo, Teuchos::EDiag diag, 
00121                const MultiVector<DomainScalar,Node> &Y, MultiVector<RangeScalar,Node> &X) const;
00122 
00124 
00125   protected:
00127     CUSPSparseOps(const CUSPSparseOps& source);
00128 
00129     RCP<Node> node_;
00130     RCP<cusp::csr_matrix<Ordinal,Scalar,cusp::host_memory  > > hostCSR_;
00131     RCP<cusp::hyb_matrix<Ordinal,Scalar,cusp::device_memory> > devcHYB_;
00132 
00133     // we do this one of two ways: 
00134     // 1D/packed: arrays of offsets, array of ordinals, array of values.
00135     ArrayRCP<size_t>  begs1D_, ends1D_;
00136     ArrayRCP<Ordinal> inds1D_;
00137     // 2D: array of pointers
00138     ArrayRCP<ArrayRCP<Ordinal> > inds2D_;
00139 
00140     size_t numRows_;
00141     size_t numNZ_; 
00142     bool isEmpty_;
00143     bool indsInit_, valsInit_;
00144   };
00145 
00146   template<class Scalar, class Ordinal, class Node>
00147   CUSPSparseOps<Scalar,Ordinal,Node>::CUSPSparseOps(const RCP<Node> &node)
00148   : node_(node) 
00149   {
00150     clear();
00151   }
00152 
00153   template<class Scalar, class Ordinal, class Node>
00154   CUSPSparseOps<Scalar,Ordinal,Node>::~CUSPSparseOps() {
00155   }
00156 
00157   template <class Scalar, class Ordinal, class Node>
00158   RCP<Node> CUSPSparseOps<Scalar,Ordinal,Node>::getNode() const {
00159     return node_; 
00160   }
00161 
00162   template <class Scalar, class Ordinal, class Node>
00163   void CUSPSparseOps<Scalar,Ordinal,Node>::clear() {
00164     isEmpty_ = false;
00165     numRows_ = 0;
00166     numNZ_   = 0;
00167     indsInit_ = 0;
00168     valsInit_ = 0;
00169     //
00170     hostCSR_ = null;
00171     devcHYB_ = null;
00172     //
00173     begs1D_ = null;
00174     ends1D_ = null;
00175     inds1D_ = null;
00176     inds2D_ = null;
00177   }
00178 
00179   template <class Scalar, class Ordinal, class Node>
00180   void CUSPSparseOps<Scalar,Ordinal,Node>::initializeStructure(const CrsGraphHostCompute<Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > &graph) {
00181     TEST_FOR_EXCEPTION(indsInit_ == true || valsInit_ == true, std::runtime_error, 
00182         Teuchos::typeName(*this) << "::initializeStructure(): structure already initialized.");
00183     numRows_ = graph.getNumRows();
00184     if (graph.isEmpty() || numRows_ == 0) {
00185       isEmpty_ = true;
00186     }
00187     else if (graph.is1DStructure()) {
00188       const_cast<CrsGraphHostviceCompute<Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > &>(graph).getDeviceBuffers(inds, offs);
00189     }
00190     else {
00191       hostCSR_ = rcp(new cusp::csr_matrix<Ordinal,Scalar,cusp::host_memory>(numRows_,numCols_,numNZ_));
00192       if (graph.is2DStructure()) {
00193       }
00194       else {
00195         TEST_FOR_EXCEPT(true);
00196       }
00197     }
00198       ArrayRCP<Ordinal> inds;
00199       ArrayRCP<size_t > offs;
00200       const_cast<CrsGraphHostCompute<Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > &>(graph).getDeviceBuffers(inds, offs);
00201       pbuf_inds1D_    = inds;
00202       pbuf_offsets1D_ = offs;
00203     }
00204     indsInit_ = true;
00205   }
00206 
00207   template <class Scalar, class Ordinal, class Node>
00208   void CUSPSparseOps<Scalar,Ordinal,Node>::initializeValues(const CrsMatrixHostCompute<Scalar,Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > &matrix) {
00209     TEST_FOR_EXCEPTION(indsInit_ == false, std::runtime_error, 
00210         Teuchos::typeName(*this) << "::initializeValues(): must initialize values after graph.");
00211     TEST_FOR_EXCEPTION(numRows_ != matrix.getNumRows() || isEmpty_ != matrix.isEmpty(), std::runtime_error,
00212         Teuchos::typeName(*this) << "::initializeValues(): matrix not compatible with previously supplied graph.");
00213     if (!isEmpty_) {
00214       ArrayRCP<Scalar> vals;
00215       const_cast<CrsMatrixHostCompute<Scalar,Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > &>(matrix).getDeviceBuffer(vals);
00216       pbuf_vals1D_ = vals;
00217     }
00218     valsInit_ = true;
00219   }
00220 
00221 
00222   template <class Scalar, class Ordinal, class Node>
00223   template <class DomainScalar, class RangeScalar>
00224   void CUSPSparseOps<Scalar,Ordinal,Node>::solve(
00225                       Teuchos::ETransp trans, Teuchos::EUplo uplo, Teuchos::EDiag diag, 
00226                       const MultiVector<DomainScalar,Node> &Y,
00227                             MultiVector<RangeScalar,Node> &X) const {
00228     TEST_FOR_EXCEPTION(true, std::logic_error, 
00229       Teuchos::typeName(*this) << "::solve(): this class does not provide support for transposed multipication. Consider manually transposing the matrix.");
00230     return;
00231   }
00232 
00233 
00234   template <class Scalar, class Ordinal, class Node>
00235   template <class DomainScalar, class RangeScalar>
00236   void CUSPSparseOps<Scalar,Ordinal,Node>::multiply(
00237                                 Teuchos::ETransp trans, 
00238                                 RangeScalar alpha,
00239                                 const MultiVector<DomainScalar,Node> &X, 
00240                                 MultiVector<RangeScalar,Node> &Y) const {
00241     // beta is not provided and the output multivector enjoys overwrite semantics
00242     TEST_FOR_EXCEPTION(trans != Teuchos::NO_TRANS, std::logic_error, 
00243       Teuchos::typeName(*this) << "::multiply(): this class does not provide support for transposed multipication. Consider manually transposing the matrix.");
00244     TEST_FOR_EXCEPTION(indsInit_ == false || valsInit_ == false, std::runtime_error,
00245         Teuchos::typeName(*this) << "::multiply(): operation not fully initialized.");
00246     TEST_FOR_EXCEPT(X.getNumCols() != Y.getNumCols());
00247     ReadyBufferHelper<Node> rbh(node_);
00248     if (isEmpty_ == true) {
00249       // Y <= 0 * X
00250       //   <= 0
00251       DefaultArithmetic<MultiVector<RangeScalar,Node> >::Init(Y,Teuchos::ScalarTraits<RangeScalar>::zero());
00252     }
00253     else {
00254     }
00255     return;
00256   }
00257 
00258   template <class Scalar, class Ordinal, class Node>
00259   template <class DomainScalar, class RangeScalar>
00260   void CUSPSparseOps<Scalar,Ordinal,Node>::multiply(
00261                                 Teuchos::ETransp trans, 
00262                                 RangeScalar alpha, const MultiVector<DomainScalar,Node> &X, 
00263                                 RangeScalar beta, MultiVector<RangeScalar,Node> &Y) const { 
00264     TEST_FOR_EXCEPTION(true, std::logic_error,
00265       Teuchos::typeName(*this) << "::multiply(): CUSP does not support multiple scalar types for sparse matrix-vector multiplication.");
00266   }
00267 
00268   template <class Scalar, class Ordinal, class Node>
00269   template <>
00270   void CUSPSparseOps<Scalar,Ordinal,Node>::multiply(
00271                                 Teuchos::ETransp trans, 
00272                                 Scalar alpha, const MultiVector<Scalar,Node> &X, 
00273                                 Scalar beta, MultiVector<Scalar,Node> &Y) const {
00274     // beta is provided and the output multivector enjoys accumulate semantics
00275     TEST_FOR_EXCEPTION(indsInit_ == false || valsInit_ == false, std::runtime_error,
00276         Teuchos::typeName(*this) << "::multiply(): operation not fully initialized.");
00277     TEST_FOR_EXCEPT(X.getNumCols() != Y.getNumCols());
00278     TEST_FOR_EXCEPTION(trans != Teuchos::NO_TRANS, std::logic_error, 
00279       Teuchos::typeName(*this) << "::multiply(): this class does not provide support for transposed multipication. Consider manually transposing the matrix.");
00280     const size_t numRHS = X.getNumCols(),
00281                  Xstride = X.getStride(),
00282                  Ystride = Y.getStride();
00283     ReadyBufferHelper<Node> rbh(node_);
00284     rbh.begin();
00285     const Scalar * X = rbh.template addConstBuffer<Scalar>(X.getValues());
00286     Scalar       * Y = rbh.template addNonConstBuffer<Scalar>(Y.getValuesNonConst());
00287     rbh.end();
00288     for (int v=0; v != numRHS; ++v) {
00289       cusp::detail::device<int,Scalar>(*devcHyb_, X, Y);
00290       X += Xstride;  
00291       Y += Ystride;  
00292     }
00293     return;
00294   }
00295 
00308   template <class Ordinal, 
00309             class Node>
00310   class CrsGraph<Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > : public CrsGraphHostCompute<Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > {
00311   public:
00312     CrsGraph(size_t numRows, const Teuchos::RCP<Node> &node) : CrsGraphHostCompute<Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> >(numRows,node) {}
00313   private:
00314     CrsGraph(const CrsGraph<Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > &graph); // not implemented
00315   };
00316 
00325   template <class Scalar,
00326             class Ordinal, 
00327             class Node>
00328   class CrsMatrix<Scalar,Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > : public CrsMatrixHostCompute<Scalar,Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > {
00329   public:
00330     CrsMatrix() : CrsMatrixHostCompute<Scalar,Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> >() {}
00331     CrsMatrix(CrsGraph<Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > &graph) : CrsMatrixHostCompute<Scalar,Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> >(graph) {}
00332     CrsMatrix(const CrsGraph<Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > &graph) : CrsMatrixHostCompute<Scalar,Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> >(graph) {}
00333     void setStaticGraph(const CrsGraph<Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > &graph) {
00334       const CrsGraphHostCompute<Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > &hgraph = dynamic_cast<const CrsGraphHostCompute<Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > &>(graph);
00335       CrsMatrixHostCompute<Scalar,Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> >::setStaticGraph(hgraph);
00336     }
00337     void setOwnedGraph(CrsGraph<Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > &graph) {
00338       CrsGraphHostCompute<Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > &hgraph = dynamic_cast<CrsGraphHostCompute<Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > &>(graph);
00339       CrsMatrixHostCompute<Scalar,Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> >::setOwnedGraph(hgraph);
00340     }
00341   private:
00342     CrsMatrix(const CrsMatrix<Scalar,Ordinal,Node,CUSPSparseOps<void,Ordinal,Node> > &mat); // not implemented
00343   };
00344 
00345 
00346 } // namespace Kokkos
00347 
00348 #endif /* KOKKOS_DEFAULTSPARSEOPS_HPP */
00349 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends