Thyra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Thyra_EpetraExtAddTransformer.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 
43 #include "Thyra_EpetraExtAddTransformer.hpp"
44 #include "Thyra_AddedLinearOpBase.hpp"
45 #include "Thyra_ScaledAdjointLinearOpBase.hpp"
46 #include "Thyra_EpetraLinearOp.hpp"
47 #include "Thyra_get_Epetra_Operator.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 
68 namespace Thyra {
69 
70 
71 // Overridden from LinearOpTransformerBase
72 
73 
75  const LinearOpBase<double> &/* op_in */) const
76 {
79 }
80 
81 
84 {
85  return nonconstEpetraLinearOp();
86 }
87 
88 
90  const LinearOpBase<double> &op_in,
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;
97  using Teuchos::dyn_cast;
98 
99  //
100  // A) Get the component Thyra objects
101  //
102 
103  const AddedLinearOpBase<double> & add_op =
104  dyn_cast<const AddedLinearOpBase<double> >(op_in);
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
146  = rcp_dynamic_cast<const DiagonalLinearOpBase<double> >(A);
148  = rcp_dynamic_cast<const DiagonalLinearOpBase<double> >(B);
149 
150  // convert operators to Epetra_CrsMatrix
153  if(dA==Teuchos::null)
154  epetra_A = rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*A), true);
155  if(dB==Teuchos::null)
156  epetra_B = rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*B), true);
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
168  EpetraLinearOp &thyra_epetra_op_inout = dyn_cast<EpetraLinearOp>(*op_inout);
169  RCP<Epetra_CrsMatrix> epetra_op =
170  rcp_dynamic_cast<Epetra_CrsMatrix>(thyra_epetra_op_inout.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 
179  TEUCHOS_ASSERT_EQUALITY( add_epetra_B_err, 0 );
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
190  RCP<const Epetra_CrsMatrix> crsMat = (dA!=Teuchos::null) ? epetra_B : epetra_A;
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
199  EpetraLinearOp & thyra_epetra_op_inout = dyn_cast<EpetraLinearOp>(*op_inout);
200  RCP<Epetra_CrsMatrix> epetra_op =
201  rcp_dynamic_cast<Epetra_CrsMatrix>(thyra_epetra_op_inout.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
int ReplaceDiagonalValues(const Epetra_Vector &Diagonal)
Concrete LinearOpBase adapter subclass for Epetra_Operator object.
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
virtual void transform(const LinearOpBase< double > &op_in, const Ptr< LinearOpBase< double > > &op_inout) const
void unwrap(const LinearOpBase< Scalar > &Op, Scalar *scalar, EOpTransp *transp, const LinearOpBase< Scalar > **origOp)
Extract the overallScalar, overallTransp and const origOp from a const LinearOpBase object...
int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
Interface class for implicitly added linear operators.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Use the non-transposed operator.
T_To & dyn_cast(T_From &from)
Interface class for identity linear operators.
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types)...
T * get() const
const Epetra_Map & OperatorDomainMap() const
virtual int numOps() const =0
Returns the number of constituent operators.
int FillComplete(bool OptimizeDataStorage=true)
void initialize(const RCP< Epetra_Operator > &op, EOpTransp opTrans=NOTRANS, EApplyEpetraOpAs applyAs=EPETRA_OP_APPLY_APPLY, EAdjointEpetraOp adjointSupport=EPETRA_OP_ADJOINT_SUPPORTED, const RCP< const VectorSpaceBase< double > > &range=Teuchos::null, const RCP< const VectorSpaceBase< double > > &domain=Teuchos::null)
Fully initialize.
virtual bool isCompatible(const LinearOpBase< double > &op_in) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Ptr< T > ptr() const
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_UNREACHABLE_RETURN(dummyReturnVal)
Teuchos::RCP< Epetra_Operator > get_Epetra_Operator(LinearOpBase< Scalar > &)
Get smart pointer to non-const Epetra_Operator object from reference to a non-const EpetraLinearOp ac...
virtual Teuchos::RCP< const LinearOpBase< Scalar > > getOp(const int k) const =0
Return the kth constant constituent operator.
#define TEUCHOS_ASSERT(assertion_test)
int Scale(double ScalarConstant)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
virtual RCP< LinearOpBase< double > > createOutputOp() const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Interface class for for diagonal linear operators.
RCP< Epetra_Operator > epetra_op()