Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Tpetra_EpetraRowMatrix.hpp
1 /*
2 // @HEADER
3 // ***********************************************************************
4 //
5 // Tpetra: Templated Linear Algebra Services Package
6 // Copyright (2008) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ************************************************************************
41 // @HEADER
42 */
43 
44 #ifndef TPETRA_EPETRAROWMATRIX_HPP
45 #define TPETRA_EPETRAROWMATRIX_HPP
46 
47 #include "TpetraCore_config.h"
48 
49 #if defined(HAVE_TPETRA_EPETRA)
50 
51 #include <Epetra_Comm.h>
52 #include <Epetra_BasicRowMatrix.h>
53 #include <Tpetra_CrsMatrix.hpp>
54 #include <Teuchos_TestForException.hpp>
55 #include <memory> // std::shared_ptr
56 #include <stdexcept>
57 #include <type_traits>
58 #include <vector>
59 
60 namespace Tpetra {
61 namespace Details {
62 
63 // Epetra_MpiComm actually has reference-counting value semantics,
64 // just like std::shared_ptr. We only return
65 // std::shared_ptr<Epetra_Comm> because Epetra_Comm is an abstract
66 // base class, so we must return it by pointer.
67 std::shared_ptr<Epetra_Comm>
68 makeEpetraCommFromTeuchosComm (const Teuchos::Comm<int>& teuchosComm);
69 
70 } // namespace Details
71 } // namespace Tpetra
72 
73 namespace { // (anonymous)
74 
75 template<class EpetraGlobalOrdinalType, class TpetraMapType>
76 Epetra_Map
77 tpetraToEpetraMapTmpl (const TpetraMapType& tpetraMap)
78 {
79  using Teuchos::ArrayView;
80  typedef typename TpetraMapType::global_ordinal_type TGO;
81  typedef typename TpetraMapType::local_ordinal_type LO;
82  typedef EpetraGlobalOrdinalType EGO;
83 
84  const TGO gblNumInds = static_cast<TGO> (tpetraMap.getGlobalNumElements ());
85  const LO lclNumInds = static_cast<LO> (tpetraMap.getNodeNumElements ());
86  ArrayView<const TGO> global_index_list = tpetraMap.getNodeElementList ();
87 
88  std::vector<EGO> global_index_list_epetra;
89  const EGO* global_index_list_epetra_ptr = NULL;
90  if (std::is_same<TGO, EGO>::value) {
91  global_index_list_epetra_ptr =
92  reinterpret_cast<const EGO*> (global_index_list.getRawPtr ());
93  }
94  else {
95  global_index_list_epetra.resize (lclNumInds);
96  for (LO k = 0; k < lclNumInds; ++k) {
97  // TODO (mfh 11 Oct 2017) Detect overflow from TGO to EGO.
98  global_index_list_epetra[k] = static_cast<EGO> (global_index_list[k]);
99  }
100  global_index_list_epetra_ptr = global_index_list_epetra.data ();
101  }
102  const EGO indexBase = tpetraMap.getIndexBase ();
103  std::shared_ptr<Epetra_Comm> epetraComm =
104  Tpetra::Details::makeEpetraCommFromTeuchosComm (* (tpetraMap.getComm ()));
105  // It's OK for the shared_ptr to fall out of scope. Subclasses of
106  // Epetra_Comm have reference-counted value semantics, so passing
107  // the Epetra_Comm by const reference into Epetra_Map's constructor
108  // is safe.
109  //
110  // TODO (mfh 11 Oct 2017) Detect overflow from TGO to EGO, or from
111  // LO to int.
112  return Epetra_Map (static_cast<EGO> (gblNumInds),
113  static_cast<int> (lclNumInds),
114  global_index_list_epetra_ptr, indexBase, *epetraComm);
115 }
116 
117 } // namespace (anonymous)
118 
119 namespace Tpetra {
120 
122 template<class TpetraMatrixType>
123 class EpetraRowMatrix : public Epetra_BasicRowMatrix {
124 public:
125  EpetraRowMatrix(const Teuchos::RCP<TpetraMatrixType> &mat, const Epetra_Comm &comm);
126  virtual ~EpetraRowMatrix() {};
127 
128  int ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, double *Values, int * Indices) const;
129 
130  //not implemented
131  int ExtractMyEntryView(int CurEntry, double * & Value, int & RowIndex, int & ColIndex);
132 
133  //not implemented
134  int ExtractMyEntryView(int CurEntry, double const * & Value, int & RowIndex, int & ColIndex) const;
135 
136  int NumMyRowEntries(int MyRow, int & NumEntries) const;
137 
138 private:
139  Teuchos::RCP<TpetraMatrixType> tpetra_matrix_;
140 };//class EpetraRowMatrix
141 
142 template<class TpetraMatrixType>
143 EpetraRowMatrix<TpetraMatrixType>::EpetraRowMatrix(
144  const Teuchos::RCP<TpetraMatrixType> &mat, const Epetra_Comm &comm
145  )
146  : Epetra_BasicRowMatrix(comm),
147  tpetra_matrix_(mat)
148 {
149  using Teuchos::RCP;
150  typedef typename TpetraMatrixType::map_type tpetra_map_type;
151 #if ! defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
152  // Prefer using Epetra64 if it is enabled.
153  typedef long long EGO;
154 #elif ! defined(EPETRA_NO_32BIT_GLOBAL_INDICES)
155  // We don't have Epetra64, but we do have 32-bit indices.
156  typedef int EGO;
157 #else
158 # error "Epetra was not configured correctly. Neither 64-bit nor 32-bit indices are enabled."
159 #endif
160  const char tfecfFuncName[] = "EpetraRowMatrix: ";
161 
162  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
163  (mat.is_null (), std::invalid_argument,
164  "The input Tpetra matrix is null.");
165  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
166  (mat->getRowMap ().is_null (), std::invalid_argument,
167  "The input Tpetra matrix's row Map is null.");
168  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
169  (mat->getColMap ().is_null (), std::invalid_argument,
170  "The input Tpetra matrix's column Map is null.");
171 
172  RCP<const tpetra_map_type> tpetraRowMap = mat->getRowMap ();
173  Epetra_Map epetraRowMap =
174  tpetraToEpetraMapTmpl<EGO, tpetra_map_type> (*tpetraRowMap);
175  RCP<const tpetra_map_type> tpetraColMap = mat->getColMap ();
176  Epetra_Map epetraColMap =
177  tpetraToEpetraMapTmpl<EGO, tpetra_map_type> (*tpetraColMap);
178  this->SetMaps (epetraRowMap, epetraColMap);
179 }
180 
181 template<class TpetraMatrixType>
182 int EpetraRowMatrix<TpetraMatrixType>::ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, double *Values, int * Indices) const
183 {
184  static_assert (std::is_same<typename TpetraMatrixType::scalar_type, double>::value,
185  "This code assumes that Tpetra::CrsMatrix's scalar_type is int.");
186  static_assert (std::is_same<typename TpetraMatrixType::local_ordinal_type, int>::value,
187  "This code assumes that Tpetra::CrsMatrix's local_ordinal_type is int.");
188  Teuchos::ArrayView<int> inds(Indices, Length);
189  Teuchos::ArrayView<double> vals(Values, Length);
190  size_t num_entries = NumEntries;
191  tpetra_matrix_->getLocalRowCopy(MyRow, inds, vals, num_entries);
192  NumEntries = num_entries;
193  return 0;
194 }
195 
196 template<class TpetraMatrixType>
197 int EpetraRowMatrix<TpetraMatrixType>::ExtractMyEntryView(int CurEntry, double * & Value, int & RowIndex, int & ColIndex)
198 {
199  //not implemented
200  return -1;
201 }
202 
203 template<class TpetraMatrixType>
204 int EpetraRowMatrix<TpetraMatrixType>::ExtractMyEntryView(int CurEntry, double const * & Value, int & RowIndex, int & ColIndex) const
205 {
206  //not implemented
207  return -1;
208 }
209 
210 template<class TpetraMatrixType>
211 int EpetraRowMatrix<TpetraMatrixType>::NumMyRowEntries(int MyRow, int & NumEntries) const
212 {
213  NumEntries = tpetra_matrix_->getNumEntriesInLocalRow(MyRow);
214  return 0;
215 }
216 
217 }//namespace Tpetra
218 
219 #endif // defined(HAVE_TPETRA_EPETRA)
220 
221 //here is the include-guard #endif:
222 
223 #endif