Thyra Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Thyra_EpetraExtAddTransformer.cpp
Go to the documentation of this file.
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
44#include "Thyra_AddedLinearOpBase.hpp"
45#include "Thyra_ScaledAdjointLinearOpBase.hpp"
49#include "Thyra_DiagonalLinearOpBase.hpp"
50#include "Thyra_DefaultDiagonalLinearOp.hpp"
51#include "Thyra_IdentityLinearOpBase.hpp"
52#include "Thyra_VectorStdOps.hpp"
53#include "Epetra_Map.h"
54#include "Epetra_LocalMap.h"
55#include "Epetra_SerialComm.h"
56#include "Epetra_Vector.h"
57#include "Epetra_CrsMatrix.h"
58#include "Teuchos_Assert.hpp"
59#include "EpetraExt_ConfigDefs.h"
60#include "EpetraExt_MatrixMatrix.h"
61#include "EpetraExt_MMHelpers.h"
62#include "EpetraExt_Transpose_RowMatrix.h"
63
64
65#include "EpetraExt_RowMatrixOut.h"
66
67
68namespace Thyra {
69
70
71// Overridden from LinearOpTransformerBase
72
73
80
81
84{
85 return nonconstEpetraLinearOp();
86}
87
88
91 const Ptr<LinearOpBase<double> > &op_inout) const
92{
93 using Thyra::unwrap;
94 using EpetraExt::MatrixMatrix;
95 using Teuchos::rcp;
96 using Teuchos::rcp_dynamic_cast;
98
99 //
100 // A) Get the component Thyra objects
101 //
102
105
106#ifdef TEUCHOS_DEBUG
107 TEUCHOS_ASSERT_EQUALITY( add_op.numOps(), 2 );
108#endif
109
110
111 // get properties of first operator: Transpose, scaler multiply...etc
112 const RCP<const LinearOpBase<double> > op_A = add_op.getOp(0);
113 double A_scalar = 0.0;
114 EOpTransp A_transp = NOTRANS;
116 unwrap( op_A, &A_scalar, &A_transp, &A );
117 TEUCHOS_ASSERT(A_transp==NOTRANS || A_transp==CONJTRANS); // sanity check
118
119 // get properties of third operator: Transpose, scaler multiply...etc
120 const RCP<const LinearOpBase<double> > op_B = add_op.getOp(1);
121 double B_scalar = 0.0;
122 EOpTransp B_transp = NOTRANS;
124 unwrap( op_B, &B_scalar, &B_transp, &B );
125 TEUCHOS_ASSERT(B_transp==NOTRANS || B_transp==CONJTRANS); // sanity check
126
127 //
128 // B) Extract out the Epetra_CrsMatrix objects and the vector
129 //
130
131 // first makre sure identity operators are represented as digaonal vectors
132 if(rcp_dynamic_cast<const Thyra::IdentityLinearOpBase<double> >(A)!=Teuchos::null) {
133 RCP<Thyra::VectorBase<double> > d = Thyra::createMember(A->domain(), "d");
134 Thyra::V_S( d.ptr(), 1.0 ); // ToDo: Set ton != 1.0 and generalize
135 A = Thyra::diagonal(d);
136 }
137 if(rcp_dynamic_cast<const Thyra::IdentityLinearOpBase<double> >(B)!=Teuchos::null) {
138 RCP<Thyra::VectorBase<double> > d = Thyra::createMember(B->domain(), "d");
139 Thyra::V_S( d.ptr(), 1.0 ); // ToDo: Set ton != 1.0 and generalize
140 B = Thyra::diagonal(d);
141 }
142
143
144 // see if exactly one operator is a diagonal linear op
149
150 // convert operators to Epetra_CrsMatrix
153 if(dA==Teuchos::null)
155 if(dB==Teuchos::null)
157
158 //
159 // C) Do the explicit addition
160 //
161
162 if(epetra_A!=Teuchos::null && epetra_B!=Teuchos::null) {
163
164 // allocate space for final addition: 3 steps
165 // 1. Get destination EpetraLinearOp
166 // 2. Extract RCP to destination Epetra_CrsMatrix
167 // 3. If neccessary, allocate new Epetra_CrsMatrix
169 RCP<Epetra_CrsMatrix> epetra_op =
171 Epetra_CrsMatrix * epetra_op_raw = epetra_op.get();
172
173 // perform addition
174 const int add_epetra_B_err
175 = EpetraExt::MatrixMatrix::Add(*epetra_A,A_transp==CONJTRANS,A_scalar,*epetra_B,B_transp==CONJTRANS,B_scalar,epetra_op_raw);
176 if(epetra_op==Teuchos::null)
177 epetra_op = Teuchos::rcp(epetra_op_raw);
178
180
181 epetra_op->FillComplete(epetra_A->DomainMap(),epetra_A->RangeMap());
182
183 // set output operator to use newly create epetra_op
184 thyra_epetra_op_inout.initialize(epetra_op);
185 }
186 else if((dA!=Teuchos::null && epetra_B!=Teuchos::null) ||
187 (dB!=Teuchos::null && epetra_A!=Teuchos::null)) {
188
189 // get unique addition values
191 double matScalar = (dA!=Teuchos::null) ? B_scalar : A_scalar;
192 RCP<const DiagonalLinearOpBase<double> > diag = (dA!=Teuchos::null) ? dA : dB;
193 double diagScalar = (dA!=Teuchos::null) ? A_scalar : B_scalar;
194
195 TEUCHOS_ASSERT(crsMat!=Teuchos::null);
196 TEUCHOS_ASSERT(diag!=Teuchos::null);
197
198 // get or allocate an object to use as the destination
200 RCP<Epetra_CrsMatrix> epetra_op =
202
203 if(epetra_op==Teuchos::null)
204 epetra_op = Teuchos::rcp(new Epetra_CrsMatrix(*crsMat));
205 else
206 *epetra_op = *crsMat;
207
208 // grab vector to add to diagonal
209 RCP<const Epetra_Vector> v = get_Epetra_Vector(epetra_op->OperatorDomainMap(),diag->getDiag());
210
211 if(matScalar!=1.0)
212 epetra_op->Scale(matScalar);
213
214 // grab digaonal from matrix, do summation, then replace the values
215 RCP<Epetra_Vector> diagonal = rcp(new Epetra_Vector(epetra_op->OperatorDomainMap()));
216 TEUCHOS_TEST_FOR_EXCEPTION(epetra_op->ExtractDiagonalCopy(*diagonal),std::runtime_error,
217 "Thyra::EpetraExtractAddTransformer::transform ExtractDiagonalCopy failed!");;
218 diagonal->Update(diagScalar,*v,1.0); // no need to scale matrix, already scaled
219 TEUCHOS_TEST_FOR_EXCEPTION(epetra_op->ReplaceDiagonalValues(*diagonal),std::runtime_error,
220 "Thyra::EpetraExtractAddTransformer::transform ReplaceDiagonalValues failed!");;
221
222 // set output operator to use newly create epetra_op
223 thyra_epetra_op_inout.initialize(epetra_op);
224 }
225 else {
226 TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,
227 "Your case of adding Epetra operators is not yet implemented! Contact the Thyra developers.");
228 }
229}
230
231
232} // namespace Thyra
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
Ptr< T > ptr() const
T * get() const
virtual void transform(const LinearOpBase< double > &op_in, const Ptr< LinearOpBase< double > > &op_inout) const
virtual RCP< LinearOpBase< double > > createOutputOp() const
virtual bool isCompatible(const LinearOpBase< double > &op_in) const
Concrete LinearOpBase adapter subclass for Epetra_Operator object.
RCP< Epetra_Vector > get_Epetra_Vector(const Epetra_Map &map, const RCP< VectorBase< double > > &v)
Get a non-const Epetra_Vector view from a non-const VectorBase object if possible.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
T_To & dyn_cast(T_From &from)
Teuchos::RCP< Epetra_Operator > get_Epetra_Operator(LinearOpBase< double > &op)
Full specialization for Scalar=double.