Teko Version of the Day
Loading...
Searching...
No Matches
Teko_TpetraHelpers.cpp
1/*
2// @HEADER
3//
4// ***********************************************************************
5//
6// Teko: A package for block and physics based preconditioning
7// Copyright 2010 Sandia Corporation
8//
9// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10// the U.S. Government retains certain rights in this software.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
40//
41// ***********************************************************************
42//
43// @HEADER
44
45*/
46
47#include "Teko_TpetraHelpers.hpp"
48#include "Teko_ConfigDefs.hpp"
49
50#ifdef TEKO_HAVE_EPETRA
51#include "Thyra_EpetraLinearOp.hpp"
52#include "Thyra_EpetraThyraWrappers.hpp"
53
54// Epetra includes
55#include "Epetra_Vector.h"
56
57// EpetraExt includes
58#include "EpetraExt_ProductOperator.h"
59#include "EpetraExt_MatrixMatrix.h"
60
61#include "Teko_EpetraOperatorWrapper.hpp"
62#endif
63
64// Thyra Includes
65#include "Thyra_BlockedLinearOpBase.hpp"
66#include "Thyra_DefaultMultipliedLinearOp.hpp"
67#include "Thyra_DefaultDiagonalLinearOp.hpp"
68#include "Thyra_DefaultZeroLinearOp.hpp"
69#include "Thyra_DefaultBlockedLinearOp.hpp"
70
71#include "Thyra_SpmdVectorBase.hpp"
72#include "Thyra_SpmdVectorSpaceBase.hpp"
73#include "Thyra_ScalarProdVectorSpaceBase.hpp"
74
75// Teko includes
76#include "Teko_Utilities.hpp"
77
78// Tpetra
79#include "Thyra_TpetraLinearOp.hpp"
80#include "Thyra_TpetraMultiVector.hpp"
81#include "Tpetra_CrsMatrix.hpp"
82#include "Tpetra_Vector.hpp"
83#include "Thyra_TpetraThyraWrappers.hpp"
84#include "TpetraExt_MatrixMatrix.hpp"
85
86using Teuchos::RCP;
87using Teuchos::rcp;
88using Teuchos::rcpFromRef;
89using Teuchos::rcp_dynamic_cast;
90using Teuchos::null;
91
92namespace Teko {
93namespace TpetraHelpers {
94
105const Teuchos::RCP<const Thyra::LinearOpBase<ST> > thyraDiagOp(const RCP<const Tpetra::Vector<ST,LO,GO,NT> > & tv,const Tpetra::Map<LO,GO,NT> & map,
106 const std::string & lbl)
107{
108 const RCP<const Thyra::VectorBase<ST> > thyraVec // need a Thyra::VectorBase object
109 = Thyra::createConstVector<ST,LO,GO,NT>(tv,Thyra::createVectorSpace<ST,LO,GO,NT>(rcpFromRef(map)));
110 Teuchos::RCP<Thyra::LinearOpBase<ST> > op
111 = Teuchos::rcp(new Thyra::DefaultDiagonalLinearOp<ST>(thyraVec));
112 op->setObjectLabel(lbl);
113 return op;
114}
115
126const Teuchos::RCP<Thyra::LinearOpBase<ST> > thyraDiagOp(const RCP<Tpetra::Vector<ST,LO,GO,NT> > & tv,const Tpetra::Map<LO,GO,NT> & map,
127 const std::string & lbl)
128{
129 const RCP<Thyra::VectorBase<ST> > thyraVec // need a Thyra::VectorBase object
130 = Thyra::createVector<ST,LO,GO,NT>(tv,Thyra::createVectorSpace<ST,LO,GO,NT>(rcpFromRef(map)));
131 Teuchos::RCP<Thyra::LinearOpBase<ST> > op
132 = Teuchos::rcp(new Thyra::DefaultDiagonalLinearOp<ST>(thyraVec));
133 op->setObjectLabel(lbl);
134 return op;
135}
136
146void fillDefaultSpmdMultiVector(Teuchos::RCP<Thyra::TpetraMultiVector<ST,LO,GO,NT> > & spmdMV,
147 Teuchos::RCP<Tpetra::MultiVector<ST,LO,GO,NT> > & tpetraMV)
148{
149 // first get desired range and domain
150 //const RCP<const Thyra::SpmdVectorSpaceBase<ST> > range = spmdMV->spmdSpace();
151 const RCP<Thyra::TpetraVectorSpace<ST,LO,GO,NT> > range = Thyra::tpetraVectorSpace<ST,LO,GO,NT>(tpetraMV->getMap());
152 const RCP<const Thyra::ScalarProdVectorSpaceBase<ST> > domain
153 = rcp_dynamic_cast<const Thyra::ScalarProdVectorSpaceBase<ST> >(spmdMV->domain());
154
155 TEUCHOS_ASSERT((size_t) domain->dim()==tpetraMV->getNumVectors());
156
157 // New local view of raw data
158 if(!tpetraMV->isConstantStride())
159 TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement views of non-contiguous mult-vectors!
160
161 // Build the MultiVector
162 spmdMV->initialize(range, domain, tpetraMV);
163
164 // make sure the Epetra_MultiVector doesn't disappear prematurely
165 Teuchos::set_extra_data<RCP<Tpetra::MultiVector<ST,LO,GO,NT> > >(tpetraMV,"Tpetra::MultiVector",Teuchos::outArg(spmdMV));
166}
167
177void identityRowIndices(const Tpetra::Map<LO,GO,NT> & rowMap, const Tpetra::CrsMatrix<ST,LO,GO,NT> & mat,std::vector<GO> & outIndices)
178{
179 // loop over elements owned by this processor
180 for(size_t i=0;i<rowMap.getLocalNumElements();i++) {
181 bool rowIsIdentity = true;
182 GO rowGID = rowMap.getGlobalElement(i);
183
184 size_t numEntries = mat.getNumEntriesInGlobalRow (i);
185 auto indices = typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_global_inds_host_view_type(Kokkos::ViewAllocateWithoutInitializing("rowIndices"),numEntries);
186 auto values = typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_values_host_view_type(Kokkos::ViewAllocateWithoutInitializing("rowIndices"),numEntries);
187
188 mat.getGlobalRowCopy(rowGID,indices,values,numEntries);
189
190 // loop over the columns of this row
191 for(size_t j=0;j<numEntries;j++) {
192 GO colGID = indices(j);
193
194 // look at row entries
195 if(colGID==rowGID) rowIsIdentity &= values(j)==1.0;
196 else rowIsIdentity &= values(j)==0.0;
197
198 // not a dirchlet row...quit
199 if(not rowIsIdentity)
200 break;
201 }
202
203 // save a row that is dirchlet
204 if(rowIsIdentity)
205 outIndices.push_back(rowGID);
206 }
207}
208
219void zeroMultiVectorRowIndices(Tpetra::MultiVector<ST,LO,GO,NT> & mv,const std::vector<GO> & zeroIndices)
220{
221 LO colCnt = mv.getNumVectors();
222 std::vector<GO>::const_iterator itr;
223
224 // loop over the indices to zero
225 for(itr=zeroIndices.begin();itr!=zeroIndices.end();++itr) {
226
227 // loop over columns
228 for(int j=0;j<colCnt;j++)
229 mv.replaceGlobalValue(*itr,j,0.0);
230 }
231}
232
242ZeroedOperator::ZeroedOperator(const std::vector<GO> & zeroIndices,
243 const Teuchos::RCP<const Tpetra::Operator<ST,LO,GO,NT> > & op)
244 : zeroIndices_(zeroIndices), tpetraOp_(op)
245{ }
246
248void ZeroedOperator::apply(const Tpetra::MultiVector<ST,LO,GO,NT> & X, Tpetra::MultiVector<ST,LO,GO,NT> & Y, Teuchos::ETransp mode, ST alpha, ST beta) const
249{
250/*
251 Epetra_MultiVector temp(X);
252 zeroMultiVectorRowIndices(temp,zeroIndices_);
253 int result = epetraOp_->Apply(temp,Y);
254*/
255
256 tpetraOp_->apply(X,Y,mode,alpha,beta);
257
258 // zero a few of the rows
259 zeroMultiVectorRowIndices(Y,zeroIndices_);
260}
261
262bool isTpetraLinearOp(const LinearOp & op)
263{
264 // See if the operator is a TpetraLinearOp
265 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
266 if (!tOp.is_null())
267 return true;
268
269 // See if the operator is a wrapped TpetraLinearOp
270 ST scalar = 0.0;
271 Thyra::EOpTransp transp = Thyra::NOTRANS;
272 RCP<const Thyra::LinearOpBase<ST> > wrapped_op;
273 Thyra::unwrap(op, &scalar, &transp, &wrapped_op);
274 tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(wrapped_op);
275 if (!tOp.is_null())
276 return true;
277
278 return false;
279}
280
281RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > getTpetraCrsMatrix(const LinearOp & op, ST *scalar, bool *transp)
282{
283 // If the operator is a TpetraLinearOp
284 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
285 if(!tOp.is_null()){
286 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > matrix = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),true);
287 *scalar = 1.0;
288 *transp = false;
289 return matrix;
290 }
291
292 // If the operator is a wrapped TpetraLinearOp
293 RCP<const Thyra::LinearOpBase<ST> > wrapped_op;
294 Thyra::EOpTransp eTransp = Thyra::NOTRANS;
295 Thyra::unwrap(op, scalar, &eTransp, &wrapped_op);
296 tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(wrapped_op,true);
297 if(!tOp.is_null()){
298 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > matrix = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),true);
299 *transp = true;
300 if(eTransp == Thyra::NOTRANS)
301 *transp = false;
302 return matrix;
303 }
304
305 return Teuchos::null;
306}
307
308#ifdef TEKO_HAVE_EPETRA
309RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > epetraCrsMatrixToTpetra(const RCP<const Epetra_CrsMatrix> A_e, const RCP<const Teuchos::Comm<int> > comm)
310{
311 int* ptr;
312 int* ind;
313 double* val;
314
315 int info = A_e->ExtractCrsDataPointers (ptr, ind, val);
316 TEUCHOS_TEST_FOR_EXCEPTION(info!=0,std::logic_error, "Could not extract data from Epetra_CrsMatrix");
317 const LO numRows = A_e->Graph ().NumMyRows ();
318 const LO nnz = A_e->Graph ().NumMyEntries ();
319
320 Teuchos::ArrayRCP<size_t> ptr2 (numRows+1);
321 Teuchos::ArrayRCP<int> ind2 (nnz);
322 Teuchos::ArrayRCP<double> val2 (nnz);
323
324 std::copy (ptr, ptr + numRows + 1, ptr2.begin ());
325 std::copy (ind, ind + nnz, ind2.begin ());
326 std::copy (val, val + nnz, val2.begin ());
327
328 RCP<const Tpetra::Map<LO,GO,NT> > rowMap = epetraMapToTpetra(A_e->RowMap(),comm);
329 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > A_t = Tpetra::createCrsMatrix<ST,LO,GO,NT>(rowMap, A_e->GlobalMaxNumEntries());
330
331 RCP<const Tpetra::Map<LO,GO,NT> > domainMap = epetraMapToTpetra(A_e->OperatorDomainMap(),comm);
332 RCP<const Tpetra::Map<LO,GO,NT> > rangeMap = epetraMapToTpetra(A_e->OperatorRangeMap(),comm);
333 RCP<const Tpetra::Map<LO,GO,NT> > colMap = epetraMapToTpetra(A_e->ColMap(),comm);
334
335 A_t->replaceColMap(colMap);
336 A_t->setAllValues (ptr2, ind2, val2);
337 A_t->fillComplete(domainMap,rangeMap);
338 return A_t;
339
340}
341
342RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > nonConstEpetraCrsMatrixToTpetra(const RCP<Epetra_CrsMatrix> A_e, const RCP<const Teuchos::Comm<int> > comm)
343{
344 int* ptr;
345 int* ind;
346 double* val;
347
348 int info = A_e->ExtractCrsDataPointers (ptr, ind, val);
349 TEUCHOS_TEST_FOR_EXCEPTION(info!=0,std::logic_error, "Could not extract data from Epetra_CrsMatrix");
350 const LO numRows = A_e->Graph ().NumMyRows ();
351 const LO nnz = A_e->Graph ().NumMyEntries ();
352
353 Teuchos::ArrayRCP<size_t> ptr2 (numRows+1);
354 Teuchos::ArrayRCP<int> ind2 (nnz);
355 Teuchos::ArrayRCP<double> val2 (nnz);
356
357 std::copy (ptr, ptr + numRows + 1, ptr2.begin ());
358 std::copy (ind, ind + nnz, ind2.begin ());
359 std::copy (val, val + nnz, val2.begin ());
360
361 RCP<const Tpetra::Map<LO,GO,NT> > rowMap = epetraMapToTpetra(A_e->RowMap(),comm);
362 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > A_t = Tpetra::createCrsMatrix<ST,LO,GO,NT>(rowMap, A_e->GlobalMaxNumEntries());
363
364 RCP<const Tpetra::Map<LO,GO,NT> > domainMap = epetraMapToTpetra(A_e->OperatorDomainMap(),comm);
365 RCP<const Tpetra::Map<LO,GO,NT> > rangeMap = epetraMapToTpetra(A_e->OperatorRangeMap(),comm);
366 RCP<const Tpetra::Map<LO,GO,NT> > colMap = epetraMapToTpetra(A_e->ColMap(),comm);
367
368 A_t->replaceColMap(colMap);
369 A_t->setAllValues (ptr2, ind2, val2);
370 A_t->fillComplete(domainMap,rangeMap);
371 return A_t;
372
373}
374
375RCP<const Tpetra::Map<LO,GO,NT> > epetraMapToTpetra(const Epetra_Map eMap, const RCP<const Teuchos::Comm<int> > comm)
376{
377 std::vector<int> intGIDs(eMap.NumMyElements());
378 eMap.MyGlobalElements(&intGIDs[0]);
379
380 std::vector<GO> myGIDs(eMap.NumMyElements());
381 for(int k = 0; k < eMap.NumMyElements(); k++)
382 myGIDs[k] = (GO) intGIDs[k];
383
384 return rcp(new const Tpetra::Map<LO,GO,NT>(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(),Teuchos::ArrayView<GO>(myGIDs),0,comm));
385}
386#endif // TEKO_HAVE_EPETRA
387
388} // end namespace TpetraHelpers
389} // end namespace Teko
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.
ZeroedOperator(const std::vector< GO > &zeroIndices, const Teuchos::RCP< const Tpetra::Operator< ST, LO, GO, NT > > &op)
Constructor for a ZeroedOperator.