Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_EpetraOperatorWrapper.cpp
1// @HEADER
2// ***********************************************************************
3//
4// Thyra: Interfaces and Support for Abstract Numerical Algorithms
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#include "Thyra_EpetraOperatorWrapper.hpp"
44#include "Thyra_DetachedSpmdVectorView.hpp"
45#include "Thyra_DefaultProductVector.hpp"
46#include "Thyra_ProductVectorSpaceBase.hpp"
47#include "Thyra_SpmdVectorBase.hpp"
48#include "Thyra_EpetraLinearOp.hpp"
49
50#ifdef HAVE_MPI
51# include "Epetra_MpiComm.h"
52#endif
53#include "Epetra_SerialComm.h"
54#include "Epetra_Vector.h"
55
56#ifdef HAVE_MPI
58#endif
59#include "Teuchos_DefaultSerialComm.hpp"
60
61
62namespace Thyra {
63
64
65// Constructor, utilties
66
67
69 const RCP<const LinearOpBase<double> > &thyraOp
70 )
71 : useTranspose_(false),
72 thyraOp_(thyraOp),
73 range_(thyraOp->range()),
74 domain_(thyraOp->domain()),
75 comm_(getEpetraComm(*thyraOp)),
76 rangeMap_(get_Epetra_Map(*range_, comm_)),
77 domainMap_(get_Epetra_Map(*domain_, comm_)),
78 label_(thyraOp->description())
79{;}
80
81
83 const Ptr<VectorBase<double> > &thyraVec) const
84{
85
86 using Teuchos::rcpFromPtr;
87 using Teuchos::rcp_dynamic_cast;
88
89 const int numVecs = x.NumVectors();
90
91 TEUCHOS_TEST_FOR_EXCEPTION(numVecs != 1, std::runtime_error,
92 "epetraToThyra does not work with MV dimension != 1");
93
94 const RCP<ProductVectorBase<double> > prodThyraVec =
95 castOrCreateNonconstProductVectorBase(rcpFromPtr(thyraVec));
96
97 const ArrayView<const double> epetraData(x[0], x.Map().NumMyElements());
98 // NOTE: I tried using Epetra_MultiVector::operator()(int) to return an
99 // Epetra_Vector object but it has a defect when Reset(...) is called which
100 // results in a memory access error (see bug 4700).
101
102 int offset = 0;
103 const int numBlocks = prodThyraVec->productSpace()->numBlocks();
104 for (int b = 0; b < numBlocks; ++b) {
105 const RCP<VectorBase<double> > vec_b = prodThyraVec->getNonconstVectorBlock(b);
106 const RCP<const SpmdVectorSpaceBase<double> > spmd_vs_b =
107 rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(vec_b->space(), true);
109 const ArrayRCP<double> thyraData = view.sv().values();
110 const int localNumElems = spmd_vs_b->localSubDim();
111 for (int i=0; i < localNumElems; ++i) {
112 thyraData[i] = epetraData[i+offset];
113 }
114 offset += localNumElems;
115 }
116
117}
118
119
121 const VectorBase<double>& thyraVec, Epetra_MultiVector& x) const
122{
123
124 using Teuchos::rcpFromRef;
125 using Teuchos::rcp_dynamic_cast;
126
127 const int numVecs = x.NumVectors();
128
129 TEUCHOS_TEST_FOR_EXCEPTION(numVecs != 1, std::runtime_error,
130 "epetraToThyra does not work with MV dimension != 1");
131
132 const RCP<const ProductVectorBase<double> > prodThyraVec =
133 castOrCreateProductVectorBase(rcpFromRef(thyraVec));
134
135 const ArrayView<double> epetraData(x[0], x.Map().NumMyElements());
136 // NOTE: See above!
137
138 int offset = 0;
139 const int numBlocks = prodThyraVec->productSpace()->numBlocks();
140 for (int b = 0; b < numBlocks; ++b) {
141 const RCP<const VectorBase<double> > vec_b = prodThyraVec->getVectorBlock(b);
142 const RCP<const SpmdVectorSpaceBase<double> > spmd_vs_b =
143 rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(vec_b->space(), true);
145 const ArrayRCP<const double> thyraData = view.sv().values();
146 const int localNumElems = spmd_vs_b->localSubDim();
147 for (int i=0; i < localNumElems; ++i) {
148 epetraData[i+offset] = thyraData[i];
149 }
150 offset += localNumElems;
151 }
152
153}
154
155
156// Overridden from Epetra_Operator
157
158
160 Epetra_MultiVector& Y) const
161{
162
164 opRange = ( !useTranspose_ ? range_ : domain_ ),
165 opDomain = ( !useTranspose_ ? domain_ : range_ );
166
167 const RCP<VectorBase<double> > tx = createMember(opDomain);
168 copyEpetraIntoThyra(X, tx.ptr());
169
170 const RCP<VectorBase<double> > ty = createMember(opRange);
171
172 Thyra::apply<double>( *thyraOp_, !useTranspose_ ? NOTRANS : CONJTRANS, *tx, ty.ptr());
173
174 copyThyraIntoEpetra(*ty, Y);
175
176 return 0;
177
178}
179
180
182 Epetra_MultiVector& /* Y */) const
183{
184 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
185 "EpetraOperatorWrapper::ApplyInverse not implemented");
187}
188
189
191{
192 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
193 "EpetraOperatorWrapper::NormInf not implemated");
195}
196
197
198// private
199
200
202EpetraOperatorWrapper::getEpetraComm(const LinearOpBase<double>& thyraOp)
203{
204
205 using Teuchos::rcp;
206 using Teuchos::rcp_dynamic_cast;
208#ifdef HAVE_MPI
209 using Teuchos::MpiComm;
210#endif
211
213
215 rcp_dynamic_cast<const ProductVectorSpaceBase<double> >(vs);
216
217 if (nonnull(prod_vs))
218 vs = prod_vs->getBlock(0);
219
221 rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(vs, true);
222
223 return get_Epetra_Comm(*spmd_vs->getComm());
224
225}
226
227
228} // namespace Thyra
229
230
232Thyra::makeEpetraWrapper(const RCP<const LinearOpBase<double> > &thyraOp)
233{
234 return epetraLinearOp(
235 Teuchos::rcp(new EpetraOperatorWrapper(thyraOp))
236 );
237}
int NumMyElements() const
const Epetra_BlockMap & Map() const
int NumVectors() const
Ptr< T > ptr() const
Create an explicit detached non-mutable (const) view of all of the local elements on this process of ...
Create an explicit detached mutable (non-const) view of all of the local elements on this process of ...
void copyEpetraIntoThyra(const Epetra_MultiVector &x, const Ptr< VectorBase< double > > &thyraVec) const
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
EpetraOperatorWrapper(const RCP< const LinearOpBase< double > > &thyraOp)
void copyThyraIntoEpetra(const VectorBase< double > &thyraVec, Epetra_MultiVector &x) const
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Base class for all linear operators.
virtual RCP< const VectorSpaceBase< Scalar > > range() const =0
Return a smart pointer for the range space for this operator.
Abstract interface for finite-dimensional dense vectors.
RCP< const Epetra_Comm > get_Epetra_Comm(const Teuchos::Comm< Ordinal > &comm)
Get (or create) and Epetra_Comm given a Teuchos::Comm object.
RCP< const Epetra_Map > get_Epetra_Map(const VectorSpaceBase< double > &vs, const RCP< const Epetra_Comm > &comm)
Get (or create) an Epetra_Map object given an VectorSpaceBase object an optionally an extra Epetra_Co...
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
@ NOTRANS
Use the non-transposed operator.
@ CONJTRANS
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types).
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
T_To & dyn_cast(T_From &from)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)