Kokkos Node API and Local Linear Algebra Kernels Version of the Day
Kokkos_SeqSparseOps.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_MySparseOps_hpp
00043 #define __Kokkos_MySparseOps_hpp
00044 
00045 #include <Teuchos_DataAccess.hpp>
00046 #include <Teuchos_CompileTimeAssert.hpp>
00047 #include <iterator>
00048 #include <stdexcept>
00049 
00050 #include "Kokkos_ConfigDefs.hpp"
00051 #include "Kokkos_CrsMatrixBase.hpp"
00052 #include "Kokkos_CrsGraphBase.hpp"
00053 
00054 #include "Kokkos_MultiVector.hpp"
00055 #include "Kokkos_DefaultArithmetic.hpp"
00056 
00057 #include "Kokkos_Raw_SparseMatVec_decl.hpp"
00058 #include "Kokkos_Raw_SparseMatVec_def.hpp"
00059 #include "Kokkos_Raw_SparseTriangularSolve_decl.hpp"
00060 #include "Kokkos_Raw_SparseTriangularSolve_def.hpp"
00061 
00062 
00063 namespace Kokkos {
00064 
00068   template <class Ordinal,class Node>
00069   class SeqCrsGraph {
00070   public:
00071     typedef Ordinal ordinal_type;
00072     typedef Node node_type;
00073 
00074     SeqCrsGraph (size_t numRows,
00075                  const Teuchos::RCP<Node>& node,
00076                  const Teuchos::RCP<Teuchos::ParameterList>& params);
00077 
00078     Teuchos::RCP<Node> getNode() const {
00079       return node_;
00080     }
00081 
00082     size_t getNumRows () const {
00083       return numRows_;
00084     }
00085 
00086     Teuchos::ArrayRCP<const size_t> getPointers() const {
00087       return ptr_;
00088     }
00089 
00090     Teuchos::ArrayRCP<const Ordinal> getIndices() const {
00091       return ind_;
00092     }
00093 
00094     bool isInitialized() const {
00095       return isInitialized_;
00096     }
00097 
00102     bool isEmpty() const {
00103       return isEmpty_;
00104     }
00105 
00106     void setStructure (const Teuchos::ArrayRCP<const size_t>& ptr,
00107                        const Teuchos::ArrayRCP<const Ordinal>& ind);
00108     void setMatDesc (Teuchos::EUplo uplo, Teuchos::EDiag diag);
00109     void getMatDesc (Teuchos::EUplo &uplo, Teuchos::EDiag &diag) const;
00110 
00111   private:
00112     Teuchos::RCP<Node> node_;
00113     size_t numRows_;
00114     //Teuchos::RCP<ParameterList> params_;
00115     bool isInitialized_;
00116     bool isEmpty_;
00117     Teuchos::EUplo tri_uplo_;
00118     Teuchos::EDiag unit_diag_;
00119 
00120     Teuchos::ArrayRCP<const size_t> ptr_;
00121     Teuchos::ArrayRCP<const Ordinal> ind_;
00122   };
00123 
00124 
00131   template <class Scalar,
00132             class Ordinal,
00133             class Node>
00134   class SeqCrsMatrix {
00135   public:
00136     typedef Scalar scalar_type;
00137     typedef Ordinal ordinal_type;
00138     typedef Node node_type;
00139     typedef SeqCrsGraph<Ordinal,Node> graph_type;
00140 
00141     SeqCrsMatrix (const Teuchos::RCP<const SeqCrsGraph<Ordinal,Node> > &graph,
00142                   const Teuchos::RCP<Teuchos::ParameterList> &params);
00143 
00144     void setValues (const Teuchos::ArrayRCP<const Scalar>& val);
00145 
00146     Teuchos::ArrayRCP<const Scalar> getValues() const {
00147       return val_;
00148     }
00149 
00150     bool isInitialized() const {
00151       return isInitialized_;
00152     }
00153 
00154   private:
00155     Teuchos::RCP<const graph_type> graph_;
00156     Teuchos::ArrayRCP<const Scalar> val_;
00157     bool isInitialized_;
00158   };
00159 
00160   template <class Ordinal, class Node>
00161   SeqCrsGraph<Ordinal,Node>::
00162   SeqCrsGraph (size_t numRows,
00163                const Teuchos::RCP<Node> &node,
00164                const Teuchos::RCP<Teuchos::ParameterList> &params) :
00165     node_ (node),
00166     numRows_ (numRows),
00167     isInitialized_ (false),
00168     isEmpty_ (numRows == 0), // provisional; a matrix with numRows > 0
00169                              // may still have zero entries.
00170     tri_uplo_ (Teuchos::UNDEF_TRI),
00171     unit_diag_ (Teuchos::NON_UNIT_DIAG)
00172   {
00173     // Make sure that users only specialize for Kokkos Node types that
00174     // are host Nodes (vs. device Nodes, such as GPU Nodes).
00175     Teuchos::CompileTimeAssert<Node::isHostNode == false> cta; (void) cta;
00176 
00177     // We don't use params currently.
00178     (void) params;
00179   }
00180 
00181   template <class Ordinal, class Node>
00182   void
00183   SeqCrsGraph<Ordinal,Node>::
00184   setStructure (const Teuchos::ArrayRCP<const size_t> &ptr,
00185                 const Teuchos::ArrayRCP<const Ordinal> &ind)
00186   {
00187     std::string tfecfFuncName("setStructure(ptr,ind)");
00188     const size_t numRows = this->getNumRows();
00189 
00190     // mfh 19 June 2012: The tests expect std::runtime_error rather
00191     // than the arguably more appropriate std::invalid_argument, so
00192     // I'll throw std::runtime_error here.  Ditto for the other checks
00193     // below.
00194     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00195       ptr.is_null (),
00196       std::runtime_error,
00197       ": The input array 'ptr' must be nonnull, even for a matrix with zero "
00198       "rows.");
00199 
00200     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00201       (size_t) ptr.size() != numRows+1,
00202       std::runtime_error,
00203       ": Graph input data are not coherent:\n"
00204       "-- ptr.size() = " << ptr.size() << " != numRows+1 = "
00205       << (numRows+1) << ".");
00206 
00207     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00208       ptr[0] != 0,
00209       std::runtime_error,
00210       ": Graph input data are not coherent:\n"
00211       "-- ptr[0] = " << ptr[0] << " != 0.");
00212 
00213     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00214       (size_t) ind.size() != ptr[numRows],
00215       std::runtime_error,
00216       ": Graph input data are not coherent:\n"
00217       "-- ind.size() = " << ind.size() << " != ptr[numRows="
00218       << numRows << "] = " << ptr[numRows] << ".");
00219 
00220     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00221       isInitialized_,
00222       std::runtime_error,
00223       ": Graph has already been initialized."
00224     )
00225 
00226     const size_t numEntries = ptr[numRows];
00227     if (numRows == 0 || numEntries == 0) {
00228       isEmpty_ = true;
00229     }
00230     else {
00231       isEmpty_ = false;
00232     }
00233     ptr_ = ptr;
00234     ind_ = ind;
00235     isInitialized_ = true;
00236   }
00237 
00238   template <class Ordinal, class Node>
00239   void
00240   SeqCrsGraph<Ordinal,Node>::
00241   setMatDesc (Teuchos::EUplo uplo, Teuchos::EDiag diag)
00242   {
00243     tri_uplo_ = uplo;
00244     unit_diag_ = diag;
00245   }
00246 
00247   template <class Ordinal, class Node>
00248   void
00249   SeqCrsGraph<Ordinal,Node>::
00250   getMatDesc (Teuchos::EUplo &uplo, Teuchos::EDiag &diag) const
00251   {
00252     uplo = tri_uplo_;
00253     diag = unit_diag_;
00254   }
00255 
00256   template <class Scalar, class Ordinal, class Node>
00257   SeqCrsMatrix<Scalar,Ordinal,Node>::
00258   SeqCrsMatrix (const Teuchos::RCP<const SeqCrsGraph<Ordinal,Node> >& graph,
00259                 const Teuchos::RCP<Teuchos::ParameterList>& params) :
00260     graph_ (graph),
00261     isInitialized_ (false)
00262   {
00263     // Make sure that users only specialize for Kokkos Node types that
00264     // are host Nodes (vs. device Nodes, such as GPU Nodes).
00265     Teuchos::CompileTimeAssert<Node::isHostNode == false> cta; (void) cta;
00266 
00267     // We don't use params currently.
00268     (void) params;
00269   }
00270 
00271   template <class Scalar, class Ordinal, class Node>
00272   void
00273   SeqCrsMatrix<Scalar,Ordinal,Node>::
00274   setValues (const Teuchos::ArrayRCP<const Scalar> &val)
00275   {
00276     std::string tfecfFuncName("setValues(val)");
00277     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00278       isInitialized_, std::runtime_error, ": The matrix is already initialized."
00279     );
00280     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00281       graph_.is_null() && ! val.is_null(),
00282       std::runtime_error,
00283       ": The matrix has a null graph, but you're trying to give it a nonnull "
00284       "array of values."
00285     );
00286     val_ = val;
00287     if (val_.is_null ()) {
00288       isInitialized_ = false;
00289     }
00290     else {
00291       isInitialized_ = true;
00292     }
00293   }
00294 
00303   template <class Scalar, class Ordinal, class Node>
00304   class SeqSparseOps {
00305   public:
00307 
00308 
00310     typedef Scalar  scalar_type;
00312     typedef Ordinal ordinal_type;
00314     typedef Node    node_type;
00316     typedef SeqSparseOps<Scalar, Ordinal, Node> sparse_ops_type;
00317 
00319     template <class O, class N>
00320     struct graph {
00321       typedef SeqCrsGraph<O,N> graph_type;
00322     };
00323 
00325     template <class S, class O, class N>
00326     struct matrix {
00327       typedef SeqCrsMatrix<S,O,N> matrix_type;
00328     };
00329 
00343     template <class S2>
00344     struct bind_scalar {
00345       typedef SeqSparseOps<S2, Ordinal, Node> other_type;
00346     };
00347 
00349 
00350 
00351 
00355     SeqSparseOps (const Teuchos::RCP<Node>& node);
00356 
00358     ~SeqSparseOps();
00359 
00361 
00362 
00363 
00365     Teuchos::RCP<Node> getNode () const;
00366 
00368 
00369 
00370 
00372     static Teuchos::ArrayRCP<size_t>
00373     allocRowPtrs (const ArrayView<const size_t>& numEntriesPerRow);
00374 
00376     template <class T>
00377     static Teuchos::ArrayRCP<T>
00378     allocStorage (const Teuchos::ArrayView<const size_t>& rowPtrs);
00379 
00381     static void
00382     finalizeGraph (Teuchos::EUplo uplo,
00383                    Teuchos::EDiag diag,
00384                    SeqCrsGraph<Ordinal, Node>& graph,
00385                    const Teuchos::RCP<Teuchos::ParameterList> &params);
00386 
00388     static void
00389     finalizeMatrix (const SeqCrsGraph<Ordinal, Node>& graph,
00390                     SeqCrsMatrix<Scalar, Ordinal, Node>& matrix,
00391                     const Teuchos::RCP<Teuchos::ParameterList>& params);
00392 
00394     static void
00395     finalizeGraphAndMatrix (Teuchos::EUplo uplo,
00396                             Teuchos::EDiag diag,
00397                             SeqCrsGraph<Ordinal, Node>& graph,
00398                             SeqCrsMatrix<Scalar, Ordinal, Node>& matrix,
00399                             const Teuchos::RCP<Teuchos::ParameterList>& params);
00400 
00402     void
00403     setGraphAndMatrix (const Teuchos::RCP<const SeqCrsGraph<Ordinal,Node> > &graph,
00404                        const Teuchos::RCP<const SeqCrsMatrix<Scalar,Ordinal,Node> > &matrix);
00405 
00407 
00408 
00409 
00435     template <class DomainScalar, class RangeScalar>
00436     void
00437     multiply (Teuchos::ETransp trans,
00438               RangeScalar alpha,
00439               const MultiVector<DomainScalar,Node> &X,
00440               MultiVector<RangeScalar,Node> &Y) const;
00441 
00471     template <class DomainScalar, class RangeScalar>
00472     void
00473     multiply (Teuchos::ETransp trans,
00474               RangeScalar alpha,
00475               const MultiVector<DomainScalar,Node> &X,
00476               RangeScalar beta,
00477               MultiVector<RangeScalar,Node> &Y) const;
00478 
00499     template <class DomainScalar, class RangeScalar>
00500     void
00501     solve (Teuchos::ETransp trans,
00502            const MultiVector<DomainScalar,Node> &Y,
00503            MultiVector<RangeScalar,Node> &X) const;
00504 
00506 
00507   private:
00509     SeqSparseOps (const SeqSparseOps& source);
00510 
00512     Teuchos::RCP<Node> node_;
00513 
00514     // we do this one of two ways:
00515     // packed CRS: array of row pointers, array of indices, array of values.
00516 
00517     ArrayRCP<const Ordinal> ptr_;
00518     ArrayRCP<const Ordinal> ind_;
00519     ArrayRCP<const Scalar>  val_;
00520 
00521     Teuchos::EUplo  tri_uplo_;
00522     Teuchos::EDiag unit_diag_;
00523 
00524     size_t numRows_;
00525     bool isInitialized_;
00526     bool isEmpty_;
00527   };
00528 
00529   template <class Scalar, class Ordinal, class Node>
00530   void
00531   SeqSparseOps<Scalar,Ordinal,Node>::
00532   finalizeGraph (Teuchos::EUplo uplo,
00533                  Teuchos::EDiag diag,
00534                  SeqCrsGraph<Ordinal,Node>& graph,
00535                  const Teuchos::RCP<Teuchos::ParameterList>& params)
00536   {
00537     TEUCHOS_TEST_FOR_EXCEPTION(
00538       ! graph.isInitialized(), std::runtime_error,
00539       "Kokkos::SeqSparseOps::finalizeGraph: Graph has not yet been initialized."
00540     );
00541     graph.setMatDesc (uplo, diag);
00542   }
00543 
00544   template <class Scalar, class Ordinal, class Node>
00545   void
00546   SeqSparseOps<Scalar,Ordinal,Node>::
00547   finalizeMatrix (const SeqCrsGraph<Ordinal,Node> &graph,
00548                   SeqCrsMatrix<Scalar,Ordinal,Node> &matrix,
00549                   const Teuchos::RCP<Teuchos::ParameterList> &params)
00550   {
00551     TEUCHOS_TEST_FOR_EXCEPTION(
00552       ! matrix.isInitialized(), std::runtime_error,
00553       "Kokkos::SeqSparseOps::finalizeMatrix(graph,matrix,params): "
00554       "Matrix has not yet been initialized."
00555     );
00556   }
00557 
00558   template <class Scalar, class Ordinal, class Node>
00559   void
00560   SeqSparseOps<Scalar,Ordinal,Node>::
00561   finalizeGraphAndMatrix (Teuchos::EUplo uplo,
00562                           Teuchos::EDiag diag,
00563                           SeqCrsGraph<Ordinal,Node>& graph,
00564                           SeqCrsMatrix<Scalar,Ordinal,Node>& matrix,
00565                           const Teuchos::RCP<Teuchos::ParameterList>& params)
00566   {
00567     TEUCHOS_TEST_FOR_EXCEPTION(
00568       ! graph.isInitialized(), std::runtime_error,
00569       "Kokkos::SeqSparseOps::finalizeGraphAndMatrix(graph,matrix,params): "
00570       "Graph has not yet been initialized."
00571     );
00572     TEUCHOS_TEST_FOR_EXCEPTION(
00573       ! matrix.isInitialized(), std::runtime_error,
00574       "Kokkos::SeqSparseOps::finalizeGraphAndMatrix(graph,matrix,params): "
00575       "Matrix has not yet been initialized."
00576     );
00577     graph.setMatDesc (uplo, diag);
00578   }
00579 
00580 
00581   template<class Scalar, class Ordinal, class Node>
00582   SeqSparseOps<Scalar,Ordinal,Node>::
00583   SeqSparseOps (const Teuchos::RCP<Node>& node) :
00584     node_ (node),
00585     numRows_ (0),
00586     isInitialized_ (false),
00587     isEmpty_ (true) // Provisionally...
00588   {
00589     // Make sure that users only specialize SeqSparseOps for Kokkos
00590     // Node types that are host Nodes (vs. device Nodes, such as GPU
00591     // Nodes).
00592     Teuchos::CompileTimeAssert<Node::isHostNode == false> cta; (void)cta;
00593   }
00594 
00595   template<class Scalar, class Ordinal, class Node>
00596   SeqSparseOps<Scalar,Ordinal,Node>::~SeqSparseOps() {
00597   }
00598 
00599   template <class Scalar, class Ordinal, class Node>
00600   RCP<Node> SeqSparseOps<Scalar,Ordinal,Node>::getNode() const {
00601     return node_;
00602   }
00603 
00604   template <class Scalar, class Ordinal, class Node>
00605   void
00606   SeqSparseOps<Scalar,Ordinal,Node>::
00607   setGraphAndMatrix (const Teuchos::RCP<const SeqCrsGraph<Ordinal,Node> > &opgraph,
00608                      const Teuchos::RCP<const SeqCrsMatrix<Scalar,Ordinal,Node> > &opmatrix)
00609   {
00610     using Teuchos::ArrayRCP;
00611     using Teuchos::arcp;
00612     using Teuchos::arcp_const_cast;
00613     using Teuchos::null;
00614     // using std::cerr;
00615     // using std::endl;
00616 
00617     std::string tfecfFuncName("setGraphAndMatrix(uplo,diag,graph,matrix)");
00618     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00619       isInitialized_, std::runtime_error, " operators already initialized."
00620     );
00621 
00622     ArrayRCP<const size_t> ptr = opgraph->getPointers ();
00623     ArrayRCP<const Ordinal> ind = opgraph->getIndices ();
00624     ArrayRCP<const Scalar> val = opmatrix->getValues ();
00625     const size_t numRows = opgraph->getNumRows ();
00626 
00627     // cerr << "SeqSparseOps::setGraphAndMatrix: on entry to routine:" << endl
00628     //      << "ptr = ";
00629     // std::copy (ptr.begin(), ptr.end(), std::ostream_iterator<size_t> (cerr, " "));
00630     // cerr << endl << "ind = ";
00631     // std::copy (ind.begin(), ind.end(), std::ostream_iterator<Ordinal> (cerr, " "));
00632     // cerr << endl << "val = ";
00633     // std::copy (val.begin(), val.end(), std::ostream_iterator<Scalar> (cerr, " "));
00634     // cerr << endl << "numRows = " << numRows << endl;
00635 
00636     // Verify the input data before setting internal state.
00637     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00638       (size_t) ptr.size() != numRows + 1,
00639       std::invalid_argument,
00640       ": ptr.size() = " << ptr.size() << " != numRows+1 = " << (numRows + 1) << ".");
00641     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00642       ind.size() != val.size(),
00643       std::invalid_argument,
00644       ": ind.size() = " << ind.size() << " != val.size() = " << val.size()
00645       << ", for ptr = opgraph->getPointers() and ind = opgraph->getIndices().");
00646     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00647       ptr[numRows] != (size_t) ind.size(),
00648       std::invalid_argument,
00649       ": ptr[numRows = " << numRows << "] = " << ptr[numRows]
00650       << " != ind.size() = " << ind.size() << ", for ptr = "
00651       "opgraph->getPointers() and ind = opgraph->getIndices().");
00652 
00653     numRows_ = numRows;
00654     if (opgraph->isEmpty () || numRows_ == 0) {
00655       isEmpty_ = true;
00656       // We have to go through a little trouble because ptr_ is an
00657       // array of const Ordinal, but we need to set its first entry
00658       // here.
00659       ArrayRCP<Ordinal> myPtr = arcp<Ordinal> (1);
00660       myPtr[0] = 0;
00661       ptr_ = arcp_const_cast<const Ordinal> (myPtr);
00662       myPtr = null;
00663 
00664       // The matrix is empty, so ind_ and val_ have zero entries.
00665       ind_ = null;
00666       val_ = null;
00667     }
00668     else {
00669       isEmpty_ = false;
00670       // ptr is an array of size_t, ptr_ is an array of Ordinal.  We
00671       // compute with an array of Ordinal because Ordinal is generally
00672       // a 32-bit integer type, so it uses less bandwidth and storage.
00673       // We have to go through a little trouble because ptr_ is an
00674       // array of const Ordinal, but we need to set its entries here.
00675       //
00676       // FIXME (mfh 25 Jun 2012) In debug mode, we should check the
00677       // conversions from size_t to Ordinal for possible overflow.
00678       ArrayRCP<Ordinal> myPtr = arcp<Ordinal> (ptr.size ());
00679       std::copy (ptr.getRawPtr (),
00680                  ptr.getRawPtr () + ptr.size (),
00681                  myPtr.getRawPtr ());
00682       ptr_ = arcp_const_cast<const Ordinal> (myPtr);
00683       myPtr = null;
00684 
00685       // We can just use these arrays directly.
00686       ind_ = ind;
00687       val_ = val;
00688     }
00689     opgraph->getMatDesc (tri_uplo_, unit_diag_);
00690     isInitialized_ = true;
00691 
00692     // cerr << "SeqSparseOps::setGraphAndMatrix: on exit:" << endl
00693     //      << "ptr_ = ";
00694     // std::copy (ptr_.begin(), ptr_.end(), std::ostream_iterator<Ordinal> (cerr, " "));
00695     // cerr << endl << "ind_ = ";
00696     // std::copy (ind_.begin(), ind_.end(), std::ostream_iterator<Ordinal> (cerr, " "));
00697     // cerr << endl << "val_ = ";
00698     // std::copy (val_.begin(), val_.end(), std::ostream_iterator<Scalar> (cerr, " "));
00699 
00700     std::string triUplo;
00701     if (tri_uplo_ == Teuchos::UNDEF_TRI) {
00702       triUplo = "UNDEF_TRI";
00703     }
00704     else if (tri_uplo_ == Teuchos::LOWER_TRI) {
00705       triUplo = "LOWER_TRI";
00706     }
00707     else if (tri_uplo_ == Teuchos::UPPER_TRI) {
00708       triUplo = "UPPER_TRI";
00709     }
00710     std::string unitDiag;
00711     if (unit_diag_ == Teuchos::NON_UNIT_DIAG) {
00712       unitDiag = "NON_UNIT_DIAG";
00713     }
00714     else if (unit_diag_ == Teuchos::UNIT_DIAG) {
00715       unitDiag = "UNIT_DIAG";
00716     }
00717     // cerr << endl << "numRows_ = " << numRows_ << endl
00718     //      << "isEmpty_ = " << isEmpty_ << endl
00719     //      << "tri_uplo_ = " << triUplo << endl
00720     //      << "unit_diag_ = " << unitDiag << endl;
00721   }
00722 
00723   template <class Scalar, class Ordinal, class Node>
00724   template <class DomainScalar, class RangeScalar>
00725   void
00726   SeqSparseOps<Scalar,Ordinal,Node>::
00727   solve (Teuchos::ETransp trans,
00728          const MultiVector<DomainScalar, Node>& Y,
00729          MultiVector<RangeScalar, Node>& X) const
00730   {
00731     std::string tfecfFuncName("solve(trans,Y,X)");
00732     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00733       ! isInitialized_,
00734       std::runtime_error,
00735       ": The solve was not fully initialized."
00736     );
00737     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00738       X.getNumCols() != Y.getNumCols(),
00739       std::runtime_error,
00740       ": Input and output multivectors have different numbers of vectors (columns)."
00741     );
00742     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00743       X.getNumRows() < numRows_,
00744       std::runtime_error,
00745       ": Output multivector X does not have enough rows.  X.getNumRows() == "
00746       << X.getNumRows() << ", but the matrix has " << numRows_ << " rows."
00747     );
00748 
00749     if (numRows_ == 0) {
00750       return; // Nothing to do
00751     }
00752     else if (isEmpty_) {
00753       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00754         unit_diag_ != Teuchos::UNIT_DIAG,
00755         std::runtime_error,
00756         ": Solve with empty matrix is only valid for an implicit unit diagonal.");
00757       // solve I * X = Y for X = Y
00758       DefaultArithmetic<MultiVector<RangeScalar,Node> >::Assign (X, Y);
00759     }
00760     else {
00761       typedef ordinal_type OT;
00762       typedef scalar_type MST; // matrix scalar type
00763       typedef DomainScalar DST;
00764       typedef RangeScalar RST;
00765 
00766       RST* const X_raw = X.getValuesNonConst ().getRawPtr ();
00767       const Ordinal X_stride = (Ordinal) X.getStride ();
00768       const DST* const Y_raw = Y.getValues ().getRawPtr ();
00769       const Ordinal Y_stride = (Ordinal) Y.getStride ();
00770 
00771       const Ordinal* const ptr = ptr_.getRawPtr ();
00772       const Ordinal* const ind = ind_.getRawPtr ();
00773       const Scalar*  const val = val_.getRawPtr ();
00774       const Ordinal numRows = X.getNumRows ();
00775       const Ordinal numCols = Y.getNumRows ();
00776       const Ordinal numVecs = X.getNumCols ();
00777 
00778       if (trans == Teuchos::NO_TRANS) {
00779         if (tri_uplo_ == Teuchos::LOWER_TRI) {
00780           if (unit_diag_ == Teuchos::UNIT_DIAG) {
00781             using Kokkos::Raw::lowerTriSolveCsrColMajorUnitDiag;
00782             lowerTriSolveCsrColMajorUnitDiag<OT, MST, DST, RST> (numRows, numCols,
00783                                                                  numVecs,
00784                                                                  X_raw, X_stride,
00785                                                                  ptr, ind, val,
00786                                                                  Y_raw, Y_stride);
00787           }
00788           else { // non unit diagonal
00789             using Kokkos::Raw::lowerTriSolveCsrColMajor;
00790             lowerTriSolveCsrColMajor<OT, MST, DST, RST> (numRows, numCols,
00791                                                          numVecs,
00792                                                          X_raw, X_stride,
00793                                                          ptr, ind, val,
00794                                                          Y_raw, Y_stride);
00795           }
00796         }
00797         else { // upper triangular
00798           if (unit_diag_ == Teuchos::UNIT_DIAG) {
00799             using Kokkos::Raw::upperTriSolveCsrColMajorUnitDiag;
00800             upperTriSolveCsrColMajorUnitDiag<OT, MST, DST, RST> (numRows, numCols,
00801                                                                  numVecs,
00802                                                                  X_raw, X_stride,
00803                                                                  ptr, ind, val,
00804                                                                  Y_raw, Y_stride);
00805           }
00806           else { // non unit diagonal
00807             using Kokkos::Raw::upperTriSolveCsrColMajor;
00808             upperTriSolveCsrColMajor<OT, MST, DST, RST> (numRows, numCols,
00809                                                          numVecs,
00810                                                          X_raw, X_stride,
00811                                                          ptr, ind, val,
00812                                                          Y_raw, Y_stride);
00813           }
00814         }
00815       }
00816       else if (trans == Teuchos::TRANS) {
00817         if (tri_uplo_ == Teuchos::LOWER_TRI) {
00818           if (unit_diag_ == Teuchos::UNIT_DIAG) {
00819             using Kokkos::Raw::lowerTriSolveCscColMajorUnitDiag;
00820             // numRows resp. numCols come from the number of rows in Y
00821             // resp. X, so they still appear in the same order as
00822             // in the not transposed cases above.
00823             lowerTriSolveCscColMajorUnitDiag<OT, MST, DST, RST> (numRows, numCols,
00824                                                                  numVecs,
00825                                                                  X_raw, X_stride,
00826                                                                  ptr, ind, val,
00827                                                                  Y_raw, Y_stride);
00828           }
00829           else {
00830             using Kokkos::Raw::lowerTriSolveCscColMajor;
00831             lowerTriSolveCscColMajor<OT, MST, DST, RST> (numRows, numCols,
00832                                                          numVecs,
00833                                                          X_raw, X_stride,
00834                                                          ptr, ind, val,
00835                                                          Y_raw, Y_stride);
00836           }
00837         }
00838         else { // upper triangular
00839           if (unit_diag_ == Teuchos::UNIT_DIAG) {
00840             using Kokkos::Raw::upperTriSolveCscColMajorUnitDiag;
00841             upperTriSolveCscColMajorUnitDiag<OT, MST, DST, RST> (numRows, numCols,
00842                                                                  numVecs,
00843                                                                  X_raw, X_stride,
00844                                                                  ptr, ind, val,
00845                                                                  Y_raw, Y_stride);
00846           }
00847           else {
00848             using Kokkos::Raw::upperTriSolveCscColMajor;
00849             upperTriSolveCscColMajor<OT, MST, DST, RST> (numRows, numCols,
00850                                                          numVecs,
00851                                                          X_raw, X_stride,
00852                                                          ptr, ind, val,
00853                                                          Y_raw, Y_stride);
00854           }
00855         }
00856       }
00857       else if (trans == Teuchos::CONJ_TRANS) {
00858         if (tri_uplo_ == Teuchos::LOWER_TRI) {
00859           if (unit_diag_ == Teuchos::UNIT_DIAG) {
00860             using Kokkos::Raw::lowerTriSolveCscColMajorUnitDiagConj;
00861             lowerTriSolveCscColMajorUnitDiagConj<OT, MST, DST, RST> (numRows, numCols,
00862                                                                      numVecs,
00863                                                                      X_raw, X_stride,
00864                                                                      ptr, ind, val,
00865                                                                      Y_raw, Y_stride);
00866           }
00867           else {
00868             using Kokkos::Raw::lowerTriSolveCscColMajorConj;
00869             lowerTriSolveCscColMajorConj<OT, MST, DST, RST> (numRows, numCols,
00870                                                              numVecs,
00871                                                              X_raw, X_stride,
00872                                                              ptr, ind, val,
00873                                                              Y_raw, Y_stride);
00874           }
00875         }
00876         else { // upper triangular
00877           if (unit_diag_ == Teuchos::UNIT_DIAG) {
00878             using Kokkos::Raw::upperTriSolveCscColMajorUnitDiagConj;
00879             upperTriSolveCscColMajorUnitDiagConj<OT, MST, DST, RST> (numRows, numCols,
00880                                                                      numVecs,
00881                                                                      X_raw, X_stride,
00882                                                                      ptr, ind, val,
00883                                                                      Y_raw, Y_stride);
00884           }
00885           else {
00886             using Kokkos::Raw::upperTriSolveCscColMajorConj;
00887             upperTriSolveCscColMajorConj<OT, MST, DST, RST> (numRows, numCols,
00888                                                              numVecs,
00889                                                              X_raw, X_stride,
00890                                                              ptr, ind, val,
00891                                                              Y_raw, Y_stride);
00892           }
00893         }
00894       }
00895     }
00896   }
00897 
00898   template <class Scalar, class Ordinal, class Node>
00899   template <class DomainScalar, class RangeScalar>
00900   void
00901   SeqSparseOps<Scalar,Ordinal,Node>::
00902   multiply (Teuchos::ETransp trans,
00903             RangeScalar alpha,
00904             const MultiVector<DomainScalar, Node>& X,
00905             MultiVector<RangeScalar, Node>& Y) const
00906   {
00907     typedef DomainScalar DST;
00908     typedef RangeScalar RST;
00909     const RST beta = Teuchos::ScalarTraits<RST>::zero ();
00910     this->template multiply<DST, RST> (trans, alpha, X, beta, Y);
00911   }
00912 
00913   template <class Scalar, class Ordinal, class Node>
00914   template <class DomainScalar, class RangeScalar>
00915   void
00916   SeqSparseOps<Scalar,Ordinal,Node>::
00917   multiply (Teuchos::ETransp trans,
00918             RangeScalar alpha,
00919             const MultiVector<DomainScalar, Node> &X,
00920             RangeScalar beta,
00921             MultiVector<RangeScalar, Node> &Y) const
00922   {
00923     using std::cerr;
00924     using std::endl;
00925 
00926     std::string tfecfFuncName("multiply(trans,alpha,X,beta,Y)");
00927     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00928       ! isInitialized_,
00929       std::runtime_error,
00930       ": Sparse ops not initialized.");
00931     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00932       X.getNumCols() != Y.getNumCols(),
00933       std::runtime_error,
00934       ": X and Y do not have the same number of columns.");
00935 
00936     typedef ordinal_type OT;
00937     typedef scalar_type MST; // matrix scalar type
00938     typedef DomainScalar DST;
00939     typedef RangeScalar RST;
00940 
00941     // These dimensions come from the input and output multivectors,
00942     // so they apply for the transpose case as well.
00943     const Ordinal numRows = Y.getNumRows ();
00944     const Ordinal numCols = X.getNumRows ();
00945     const Ordinal numVecs = X.getNumCols ();
00946     RST* const Y_raw = Y.getValuesNonConst ().getRawPtr ();
00947     const Ordinal Y_stride = Y.getStride ();
00948     const DST* const X_raw = X.getValues ().getRawPtr ();
00949     const Ordinal X_stride = X.getStride ();
00950 
00951     const Ordinal* const ptr = ptr_.getRawPtr ();
00952     const Ordinal* const ind = ind_.getRawPtr ();
00953     const Scalar*  const val = val_.getRawPtr ();
00954 
00955     // The switch statement selects one of the hard-coded-numVecs
00956     // routines for certain values of numVecs.  Otherwise, it picks a
00957     // general routine.
00958     //
00959     // Note that we're taking numRows and numCols from Y resp. X.
00960     // Assuming that the dimensions of X and Y are correct, then
00961     // whether or not we're applying the transpose, the (transposed,
00962     // if applicable) matrix has dimensions numRows by numCols.
00963     // That's why we don't switch the order of numRows, numCols in the
00964     // invocations below.
00965     switch (numVecs) {
00966     case 1:
00967       if (trans == Teuchos::NO_TRANS) {
00968         using Kokkos::Raw::matVecCsrColMajor1Vec;
00969         matVecCsrColMajor1Vec<OT, MST, DST, RST> (numRows, numCols,
00970                                                   beta, Y_raw, Y_stride,
00971                                                   alpha, ptr, ind, val,
00972                                                   X_raw, X_stride);
00973       }
00974       else if (trans == Teuchos::TRANS) {
00975         using Kokkos::Raw::matVecCscColMajor1Vec;
00976         matVecCscColMajor1Vec<OT, MST, DST, RST> (numRows, numCols,
00977                                                   beta, Y_raw, Y_stride,
00978                                                   alpha, ptr, ind, val,
00979                                                   X_raw, X_stride);
00980       }
00981       else { // if (trans == Teuchos::CONJ_TRANS) {
00982         using Kokkos::Raw::matVecCscColMajorConj1Vec;
00983         matVecCscColMajorConj1Vec<OT, MST, DST, RST> (numRows, numCols,
00984                                                       beta, Y_raw, Y_stride,
00985                                                       alpha, ptr, ind, val,
00986                                                       X_raw, X_stride);
00987       }
00988       break;
00989     case 2:
00990       if (trans == Teuchos::NO_TRANS) {
00991         using Kokkos::Raw::matVecCsrColMajor2Vec;
00992         matVecCsrColMajor2Vec<OT, MST, DST, RST> (numRows, numCols,
00993                                                   beta, Y_raw, Y_stride,
00994                                                   alpha, ptr, ind, val,
00995                                                   X_raw, X_stride);
00996       }
00997       else if (trans == Teuchos::TRANS) {
00998         using Kokkos::Raw::matVecCscColMajor2Vec;
00999         matVecCscColMajor2Vec<OT, MST, DST, RST> (numRows, numCols,
01000                                                   beta, Y_raw, Y_stride,
01001                                                   alpha, ptr, ind, val,
01002                                                   X_raw, X_stride);
01003       }
01004       else { // if (trans == Teuchos::CONJ_TRANS) {
01005         using Kokkos::Raw::matVecCscColMajorConj2Vec;
01006         matVecCscColMajorConj2Vec<OT, MST, DST, RST> (numRows, numCols,
01007                                                       beta, Y_raw, Y_stride,
01008                                                       alpha, ptr, ind, val,
01009                                                       X_raw, X_stride);
01010       }
01011       break;
01012     case 3:
01013       if (trans == Teuchos::NO_TRANS) {
01014         using Kokkos::Raw::matVecCsrColMajor3Vec;
01015         matVecCsrColMajor3Vec<OT, MST, DST, RST> (numRows, numCols,
01016                                                   beta, Y_raw, Y_stride,
01017                                                   alpha, ptr, ind, val,
01018                                                   X_raw, X_stride);
01019       }
01020       else if (trans == Teuchos::TRANS) {
01021         using Kokkos::Raw::matVecCscColMajor3Vec;
01022         matVecCscColMajor3Vec<OT, MST, DST, RST> (numRows, numCols,
01023                                                   beta, Y_raw, Y_stride,
01024                                                   alpha, ptr, ind, val,
01025                                                   X_raw, X_stride);
01026       }
01027       else { // if (trans == Teuchos::CONJ_TRANS) {
01028         using Kokkos::Raw::matVecCscColMajorConj3Vec;
01029         matVecCscColMajorConj3Vec<OT, MST, DST, RST> (numRows, numCols,
01030                                                       beta, Y_raw, Y_stride,
01031                                                       alpha, ptr, ind, val,
01032                                                       X_raw, X_stride);
01033       }
01034       break;
01035     case 4:
01036       if (trans == Teuchos::NO_TRANS) {
01037         using Kokkos::Raw::matVecCsrColMajor4Vec;
01038         matVecCsrColMajor4Vec<OT, MST, DST, RST> (numRows, numCols,
01039                                                   beta, Y_raw, Y_stride,
01040                                                   alpha, ptr, ind, val,
01041                                                   X_raw, X_stride);
01042       }
01043       else if (trans == Teuchos::TRANS) {
01044         using Kokkos::Raw::matVecCscColMajor4Vec;
01045         matVecCscColMajor4Vec<OT, MST, DST, RST> (numRows, numCols,
01046                                                   beta, Y_raw, Y_stride,
01047                                                   alpha, ptr, ind, val,
01048                                                   X_raw, X_stride);
01049       }
01050       else { // if (trans == Teuchos::CONJ_TRANS) {
01051         using Kokkos::Raw::matVecCscColMajorConj4Vec;
01052         matVecCscColMajorConj4Vec<OT, MST, DST, RST> (numRows, numCols,
01053                                                       beta, Y_raw, Y_stride,
01054                                                       alpha, ptr, ind, val,
01055                                                       X_raw, X_stride);
01056       }
01057       break;
01058     default: // The "general case"
01059       if (trans == Teuchos::NO_TRANS) {
01060         using Kokkos::Raw::matVecCsrColMajor;
01061         matVecCsrColMajor<OT, MST, DST, RST> (numRows, numCols, numVecs,
01062                                               beta, Y_raw, Y_stride,
01063                                               alpha, ptr, ind, val,
01064                                               X_raw, X_stride);
01065       }
01066       else if (trans == Teuchos::TRANS) {
01067         using Kokkos::Raw::matVecCscColMajor;
01068         matVecCscColMajor<OT, MST, DST, RST> (numRows, numCols, numVecs,
01069                                               beta, Y_raw, Y_stride,
01070                                               alpha, ptr, ind, val,
01071                                               X_raw, X_stride);
01072       }
01073       else { // if (trans == Teuchos::CONJ_TRANS) {
01074         using Kokkos::Raw::matVecCscColMajorConj;
01075         matVecCscColMajorConj<OT, MST, DST, RST> (numRows, numCols, numVecs,
01076                                                   beta, Y_raw, Y_stride,
01077                                                   alpha, ptr, ind, val,
01078                                                   X_raw, X_stride);
01079       }
01080     }
01081   }
01082 
01083   template <class Scalar, class Ordinal, class Node>
01084   Teuchos::ArrayRCP<size_t>
01085   SeqSparseOps<Scalar, Ordinal, Node>::
01086   allocRowPtrs (const Teuchos::ArrayView<const size_t>& numEntriesPerRow)
01087   {
01088     Teuchos::ArrayRCP<size_t> ptr (numEntriesPerRow.size() + 1);
01089     ptr[0] = 0;
01090     std::partial_sum (numEntriesPerRow.getRawPtr (),
01091                       numEntriesPerRow.getRawPtr () + numEntriesPerRow.size (),
01092                       ptr.getRawPtr () + 1);
01093     return ptr;
01094   }
01095 
01096   template <class Scalar, class Ordinal, class Node>
01097   template <class T>
01098   Teuchos::ArrayRCP<T>
01099   SeqSparseOps<Scalar, Ordinal, Node>::
01100   allocStorage (const Teuchos::ArrayView<const size_t>& rowPtrs)
01101   {
01102     TEUCHOS_TEST_FOR_EXCEPTION(rowPtrs.size() == 0, std::invalid_argument,
01103       "SeqSparseOps::allocStorage: The input rowPtrs array must have length at "
01104       "least one, but rowPtrs.size() = " << rowPtrs.size() << ".");
01105 
01106     const size_t totalNumEntries = rowPtrs[rowPtrs.size() - 1];
01107     Teuchos::ArrayRCP<T> val (totalNumEntries);
01108     std::fill (val.getRawPtr (),
01109                val.getRawPtr () + totalNumEntries,
01110                Teuchos::ScalarTraits<T>::zero ());
01111     return val;
01112   }
01113 
01114 } // namespace Kokkos
01115 
01116 #endif // #ifndef __Kokkos_MySparseOps_hpp
01117 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends