Teko Version of the Day
Loading...
Searching...
No Matches
Teko_TpetraHelpers.cpp
1// @HEADER
2// *****************************************************************************
3// Teko: A package for block and physics based preconditioning
4//
5// Copyright 2010 NTESS and the Teko contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#include "Teko_TpetraHelpers.hpp"
11#include "Teko_ConfigDefs.hpp"
12
13#ifdef TEKO_HAVE_EPETRA
14#include "Thyra_EpetraLinearOp.hpp"
15#include "Thyra_EpetraThyraWrappers.hpp"
16
17// Epetra includes
18#include "Epetra_Vector.h"
19
20// EpetraExt includes
21#include "EpetraExt_ProductOperator.h"
22#include "EpetraExt_MatrixMatrix.h"
23
24#include "Teko_EpetraOperatorWrapper.hpp"
25#endif
26
27// Thyra Includes
28#include "Thyra_BlockedLinearOpBase.hpp"
29#include "Thyra_DefaultMultipliedLinearOp.hpp"
30#include "Thyra_DefaultDiagonalLinearOp.hpp"
31#include "Thyra_DefaultZeroLinearOp.hpp"
32#include "Thyra_DefaultBlockedLinearOp.hpp"
33
34#include "Thyra_SpmdVectorBase.hpp"
35#include "Thyra_SpmdVectorSpaceBase.hpp"
36#include "Thyra_ScalarProdVectorSpaceBase.hpp"
37
38// Teko includes
39#include "Teko_Utilities.hpp"
40
41// Tpetra
42#include "Thyra_TpetraLinearOp.hpp"
43#include "Thyra_TpetraMultiVector.hpp"
44#include "Tpetra_CrsMatrix.hpp"
45#include "Tpetra_Vector.hpp"
46#include "Thyra_TpetraThyraWrappers.hpp"
47#include "TpetraExt_MatrixMatrix.hpp"
48
49using Teuchos::null;
50using Teuchos::RCP;
51using Teuchos::rcp;
52using Teuchos::rcp_dynamic_cast;
53using Teuchos::rcpFromRef;
54
55namespace Teko {
56namespace TpetraHelpers {
57
68const Teuchos::RCP<const Thyra::LinearOpBase<ST> > thyraDiagOp(
69 const RCP<const Tpetra::Vector<ST, LO, GO, NT> >& tv, const Tpetra::Map<LO, GO, NT>& map,
70 const std::string& lbl) {
71 const RCP<const Thyra::VectorBase<ST> > thyraVec // need a Thyra::VectorBase object
72 = Thyra::createConstVector<ST, LO, GO, NT>(
73 tv, Thyra::createVectorSpace<ST, LO, GO, NT>(rcpFromRef(map)));
74 Teuchos::RCP<Thyra::LinearOpBase<ST> > op =
75 Teuchos::rcp(new Thyra::DefaultDiagonalLinearOp<ST>(thyraVec));
76 op->setObjectLabel(lbl);
77 return op;
78}
79
90const Teuchos::RCP<Thyra::LinearOpBase<ST> > thyraDiagOp(
91 const RCP<Tpetra::Vector<ST, LO, GO, NT> >& tv, const Tpetra::Map<LO, GO, NT>& map,
92 const std::string& lbl) {
93 const RCP<Thyra::VectorBase<ST> > thyraVec // need a Thyra::VectorBase object
94 = Thyra::createVector<ST, LO, GO, NT>(
95 tv, Thyra::createVectorSpace<ST, LO, GO, NT>(rcpFromRef(map)));
96 Teuchos::RCP<Thyra::LinearOpBase<ST> > op =
97 Teuchos::rcp(new Thyra::DefaultDiagonalLinearOp<ST>(thyraVec));
98 op->setObjectLabel(lbl);
99 return op;
100}
101
111void fillDefaultSpmdMultiVector(Teuchos::RCP<Thyra::TpetraMultiVector<ST, LO, GO, NT> >& spmdMV,
112 Teuchos::RCP<Tpetra::MultiVector<ST, LO, GO, NT> >& tpetraMV) {
113 // first get desired range and domain
114 // const RCP<const Thyra::SpmdVectorSpaceBase<ST> > range = spmdMV->spmdSpace();
115 const RCP<Thyra::TpetraVectorSpace<ST, LO, GO, NT> > range =
116 Thyra::tpetraVectorSpace<ST, LO, GO, NT>(tpetraMV->getMap());
117 const RCP<const Thyra::ScalarProdVectorSpaceBase<ST> > domain =
118 rcp_dynamic_cast<const Thyra::ScalarProdVectorSpaceBase<ST> >(spmdMV->domain());
119
120 TEUCHOS_ASSERT((size_t)domain->dim() == tpetraMV->getNumVectors());
121
122 // New local view of raw data
123 if (!tpetraMV->isConstantStride())
124 TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement views of non-contiguous mult-vectors!
125
126 // Build the MultiVector
127 spmdMV->initialize(range, domain, tpetraMV);
128
129 // make sure the Epetra_MultiVector doesn't disappear prematurely
130 Teuchos::set_extra_data<RCP<Tpetra::MultiVector<ST, LO, GO, NT> > >(
131 tpetraMV, "Tpetra::MultiVector", Teuchos::outArg(spmdMV));
132}
133
143void identityRowIndices(const Tpetra::Map<LO, GO, NT>& rowMap,
144 const Tpetra::CrsMatrix<ST, LO, GO, NT>& mat, std::vector<GO>& outIndices) {
145 // loop over elements owned by this processor
146 for (size_t i = 0; i < rowMap.getLocalNumElements(); i++) {
147 bool rowIsIdentity = true;
148 GO rowGID = rowMap.getGlobalElement(i);
149
150 size_t numEntries = mat.getNumEntriesInGlobalRow(i);
151 auto indices = typename Tpetra::CrsMatrix<ST, LO, GO, NT>::nonconst_global_inds_host_view_type(
152 Kokkos::ViewAllocateWithoutInitializing("rowIndices"), numEntries);
153 auto values = typename Tpetra::CrsMatrix<ST, LO, GO, NT>::nonconst_values_host_view_type(
154 Kokkos::ViewAllocateWithoutInitializing("rowIndices"), numEntries);
155
156 mat.getGlobalRowCopy(rowGID, indices, values, numEntries);
157
158 // loop over the columns of this row
159 for (size_t j = 0; j < numEntries; j++) {
160 GO colGID = indices(j);
161
162 // look at row entries
163 if (colGID == rowGID)
164 rowIsIdentity &= values(j) == 1.0;
165 else
166 rowIsIdentity &= values(j) == 0.0;
167
168 // not a dirchlet row...quit
169 if (not rowIsIdentity) break;
170 }
171
172 // save a row that is dirchlet
173 if (rowIsIdentity) outIndices.push_back(rowGID);
174 }
175}
176
187void zeroMultiVectorRowIndices(Tpetra::MultiVector<ST, LO, GO, NT>& mv,
188 const std::vector<GO>& zeroIndices) {
189 LO colCnt = mv.getNumVectors();
190 std::vector<GO>::const_iterator itr;
191
192 // loop over the indices to zero
193 for (itr = zeroIndices.begin(); itr != zeroIndices.end(); ++itr) {
194 // loop over columns
195 for (int j = 0; j < colCnt; j++) mv.replaceGlobalValue(*itr, j, 0.0);
196 }
197}
198
208ZeroedOperator::ZeroedOperator(const std::vector<GO>& zeroIndices,
209 const Teuchos::RCP<const Tpetra::Operator<ST, LO, GO, NT> >& op)
210 : zeroIndices_(zeroIndices), tpetraOp_(op) {}
211
213void ZeroedOperator::apply(const Tpetra::MultiVector<ST, LO, GO, NT>& X,
214 Tpetra::MultiVector<ST, LO, GO, NT>& Y, Teuchos::ETransp mode, ST alpha,
215 ST beta) const {
216 /*
217 Epetra_MultiVector temp(X);
218 zeroMultiVectorRowIndices(temp,zeroIndices_);
219 int result = epetraOp_->Apply(temp,Y);
220 */
221
222 tpetraOp_->apply(X, Y, mode, alpha, beta);
223
224 // zero a few of the rows
225 zeroMultiVectorRowIndices(Y, zeroIndices_);
226}
227
228bool isTpetraLinearOp(const LinearOp& op) {
229 // See if the operator is a TpetraLinearOp
230 RCP<const Thyra::TpetraLinearOp<ST, LO, GO, NT> > tOp =
231 rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(op);
232 if (!tOp.is_null()) return true;
233
234 // See if the operator is a wrapped TpetraLinearOp
235 ST scalar = 0.0;
236 Thyra::EOpTransp transp = Thyra::NOTRANS;
237 RCP<const Thyra::LinearOpBase<ST> > wrapped_op;
238 Thyra::unwrap(op, &scalar, &transp, &wrapped_op);
239 tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(wrapped_op);
240 if (!tOp.is_null()) return true;
241
242 return false;
243}
244
245RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> > getTpetraCrsMatrix(const LinearOp& op, ST* scalar,
246 bool* transp) {
247 // If the operator is a TpetraLinearOp
248 RCP<const Thyra::TpetraLinearOp<ST, LO, GO, NT> > tOp =
249 rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(op);
250 if (!tOp.is_null()) {
251 RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> > matrix =
252 rcp_dynamic_cast<const Tpetra::CrsMatrix<ST, LO, GO, NT> >(tOp->getConstTpetraOperator(),
253 true);
254 *scalar = 1.0;
255 *transp = false;
256 return matrix;
257 }
258
259 // If the operator is a wrapped TpetraLinearOp
260 RCP<const Thyra::LinearOpBase<ST> > wrapped_op;
261 Thyra::EOpTransp eTransp = Thyra::NOTRANS;
262 Thyra::unwrap(op, scalar, &eTransp, &wrapped_op);
263 tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(wrapped_op, true);
264 if (!tOp.is_null()) {
265 RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> > matrix =
266 rcp_dynamic_cast<const Tpetra::CrsMatrix<ST, LO, GO, NT> >(tOp->getConstTpetraOperator(),
267 true);
268 *transp = true;
269 if (eTransp == Thyra::NOTRANS) *transp = false;
270 return matrix;
271 }
272
273 return Teuchos::null;
274}
275
276#ifdef TEKO_HAVE_EPETRA
277RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> > epetraCrsMatrixToTpetra(
278 const RCP<const Epetra_CrsMatrix> A_e, const RCP<const Teuchos::Comm<int> > comm) {
279 int* ptr;
280 int* ind;
281 double* val;
282
283 int info = A_e->ExtractCrsDataPointers(ptr, ind, val);
284 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
285 "Could not extract data from Epetra_CrsMatrix");
286 const LO numRows = A_e->Graph().NumMyRows();
287 const LO nnz = A_e->Graph().NumMyEntries();
288
289 Teuchos::ArrayRCP<size_t> ptr2(numRows + 1);
290 Teuchos::ArrayRCP<int> ind2(nnz);
291 Teuchos::ArrayRCP<double> val2(nnz);
292
293 std::copy(ptr, ptr + numRows + 1, ptr2.begin());
294 std::copy(ind, ind + nnz, ind2.begin());
295 std::copy(val, val + nnz, val2.begin());
296
297 RCP<const Tpetra::Map<LO, GO, NT> > rowMap = epetraMapToTpetra(A_e->RowMap(), comm);
298 RCP<Tpetra::CrsMatrix<ST, LO, GO, NT> > A_t =
299 Tpetra::createCrsMatrix<ST, LO, GO, NT>(rowMap, A_e->GlobalMaxNumEntries());
300
301 RCP<const Tpetra::Map<LO, GO, NT> > domainMap = epetraMapToTpetra(A_e->OperatorDomainMap(), comm);
302 RCP<const Tpetra::Map<LO, GO, NT> > rangeMap = epetraMapToTpetra(A_e->OperatorRangeMap(), comm);
303 RCP<const Tpetra::Map<LO, GO, NT> > colMap = epetraMapToTpetra(A_e->ColMap(), comm);
304
305 A_t->replaceColMap(colMap);
306 A_t->setAllValues(ptr2, ind2, val2);
307 A_t->fillComplete(domainMap, rangeMap);
308 return A_t;
309}
310
311RCP<Tpetra::CrsMatrix<ST, LO, GO, NT> > nonConstEpetraCrsMatrixToTpetra(
312 const RCP<Epetra_CrsMatrix> A_e, const RCP<const Teuchos::Comm<int> > comm) {
313 int* ptr;
314 int* ind;
315 double* val;
316
317 int info = A_e->ExtractCrsDataPointers(ptr, ind, val);
318 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
319 "Could not extract data from Epetra_CrsMatrix");
320 const LO numRows = A_e->Graph().NumMyRows();
321 const LO nnz = A_e->Graph().NumMyEntries();
322
323 Teuchos::ArrayRCP<size_t> ptr2(numRows + 1);
324 Teuchos::ArrayRCP<int> ind2(nnz);
325 Teuchos::ArrayRCP<double> val2(nnz);
326
327 std::copy(ptr, ptr + numRows + 1, ptr2.begin());
328 std::copy(ind, ind + nnz, ind2.begin());
329 std::copy(val, val + nnz, val2.begin());
330
331 RCP<const Tpetra::Map<LO, GO, NT> > rowMap = epetraMapToTpetra(A_e->RowMap(), comm);
332 RCP<Tpetra::CrsMatrix<ST, LO, GO, NT> > A_t =
333 Tpetra::createCrsMatrix<ST, LO, GO, NT>(rowMap, A_e->GlobalMaxNumEntries());
334
335 RCP<const Tpetra::Map<LO, GO, NT> > domainMap = epetraMapToTpetra(A_e->OperatorDomainMap(), comm);
336 RCP<const Tpetra::Map<LO, GO, NT> > rangeMap = epetraMapToTpetra(A_e->OperatorRangeMap(), comm);
337 RCP<const Tpetra::Map<LO, GO, NT> > colMap = epetraMapToTpetra(A_e->ColMap(), comm);
338
339 A_t->replaceColMap(colMap);
340 A_t->setAllValues(ptr2, ind2, val2);
341 A_t->fillComplete(domainMap, rangeMap);
342 return A_t;
343}
344
345RCP<const Tpetra::Map<LO, GO, NT> > epetraMapToTpetra(const Epetra_Map eMap,
346 const RCP<const Teuchos::Comm<int> > comm) {
347 std::vector<int> intGIDs(eMap.NumMyElements());
348 eMap.MyGlobalElements(&intGIDs[0]);
349
350 std::vector<GO> myGIDs(eMap.NumMyElements());
351 for (int k = 0; k < eMap.NumMyElements(); k++) myGIDs[k] = (GO)intGIDs[k];
352
353 return rcp(
354 new const Tpetra::Map<LO, GO, NT>(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(),
355 Teuchos::ArrayView<GO>(myGIDs), 0, comm));
356}
357#endif // TEKO_HAVE_EPETRA
358
359} // end namespace TpetraHelpers
360} // end namespace Teko
Teko_Utilities.hpp
Teko::TpetraHelpers::ZeroedOperator::apply
void apply(const Tpetra::MultiVector< ST, LO, GO, NT > &X, Tpetra::MultiVector< ST, LO, GO, NT > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, ST alpha=Teuchos::ScalarTraits< ST >::one(), ST beta=Teuchos::ScalarTraits< ST >::zero()) const
Perform a matrix-vector product with certain rows zeroed out.
Definition Teko_TpetraHelpers.cpp:213
Teko::TpetraHelpers::ZeroedOperator::ZeroedOperator
ZeroedOperator(const std::vector< GO > &zeroIndices, const Teuchos::RCP< const Tpetra::Operator< ST, LO, GO, NT > > &op)
Constructor for a ZeroedOperator.
Definition Teko_TpetraHelpers.cpp:208
Generated on for Teko by doxygen 1.15.0