Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
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
60namespace Tpetra {
61namespace 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.
67std::shared_ptr<Epetra_Comm>
68makeEpetraCommFromTeuchosComm (const Teuchos::Comm<int>& teuchosComm);
69
70} // namespace Details
71} // namespace Tpetra
72
73namespace { // (anonymous)
74
75template<class EpetraGlobalOrdinalType, class TpetraMapType>
76Epetra_Map
77tpetraToEpetraMapTmpl (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.getLocalNumElements ());
86 ArrayView<const TGO> global_index_list = tpetraMap.getLocalElementList ();
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
119namespace Tpetra {
120
122template<class TpetraMatrixType>
123class EpetraRowMatrix : public Epetra_BasicRowMatrix {
124public:
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
138private:
139 Teuchos::RCP<TpetraMatrixType> tpetra_matrix_;
140};//class EpetraRowMatrix
141
142template<class TpetraMatrixType>
143EpetraRowMatrix<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
181template<class TpetraMatrixType>
182int EpetraRowMatrix<TpetraMatrixType>::ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, double *Values, int * Indices) const
183{
184 using inds_view = typename TpetraMatrixType::nonconst_local_inds_host_view_type;
185 using vals_view = typename TpetraMatrixType::nonconst_values_host_view_type;
186 static_assert (std::is_same<typename TpetraMatrixType::scalar_type, double>::value,
187 "This code assumes that Tpetra::CrsMatrix's scalar_type is int.");
188 static_assert (std::is_same<typename TpetraMatrixType::local_ordinal_type, int>::value,
189 "This code assumes that Tpetra::CrsMatrix's local_ordinal_type is int.");
190 inds_view IndicesView(Indices, Length);
191 vals_view ValuesView(Values, Length);
192 size_t num_entries = NumEntries;
193 tpetra_matrix_->getLocalRowCopy(MyRow, IndicesView, ValuesView, num_entries);
194 NumEntries = num_entries;
195 return 0;
196}
197
198template<class TpetraMatrixType>
199int EpetraRowMatrix<TpetraMatrixType>::ExtractMyEntryView(int CurEntry, double * & Value, int & RowIndex, int & ColIndex)
200{
201 //not implemented
202 return -1;
203}
204
205template<class TpetraMatrixType>
206int EpetraRowMatrix<TpetraMatrixType>::ExtractMyEntryView(int CurEntry, double const * & Value, int & RowIndex, int & ColIndex) const
207{
208 //not implemented
209 return -1;
210}
211
212template<class TpetraMatrixType>
213int EpetraRowMatrix<TpetraMatrixType>::NumMyRowEntries(int MyRow, int & NumEntries) const
214{
215 NumEntries = tpetra_matrix_->getNumEntriesInLocalRow(MyRow);
216 return 0;
217}
218
219}//namespace Tpetra
220
221#endif // defined(HAVE_TPETRA_EPETRA)
222
223//here is the include-guard #endif:
224
225#endif
Implementation details of Tpetra.
Namespace Tpetra contains the class and methods constituting the Tpetra library.