Thyra Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Thyra_EpetraLinearOp.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
44#include "Thyra_SpmdMultiVectorBase.hpp"
45#include "Thyra_MultiVectorStdOps.hpp"
46#include "Thyra_AssertOp.hpp"
47#include "Teuchos_dyn_cast.hpp"
48#include "Teuchos_Assert.hpp"
49#include "Teuchos_getConst.hpp"
50#include "Teuchos_as.hpp"
52
53#include "Epetra_Map.h"
54#include "Epetra_Vector.h"
55#include "Epetra_Operator.h"
56#include "Epetra_CrsMatrix.h" // Printing and absolute row sums only!
57
58
59namespace Thyra {
60
61
62// Constructors / initializers / accessors
63
64
66 :isFullyInitialized_(false),
67 opTrans_(NOTRANS),
68 applyAs_(EPETRA_OP_APPLY_APPLY),
69 adjointSupport_(EPETRA_OP_ADJOINT_UNSUPPORTED)
70{}
71
72
74 const RCP<Epetra_Operator> &op,
75 EOpTransp opTrans,
76 EApplyEpetraOpAs applyAs,
77 EAdjointEpetraOp adjointSupport,
78 const RCP< const VectorSpaceBase<double> > &range_in,
79 const RCP< const VectorSpaceBase<double> > &domain_in
80 )
81{
82
83 using Teuchos::rcp_dynamic_cast;
84 typedef SpmdVectorSpaceBase<double> SPMDVSB;
85
86 // Validate input, allocate spaces, validate ...
87#ifdef TEUCHOS_DEBUG
88 TEUCHOS_TEST_FOR_EXCEPTION( is_null(op), std::invalid_argument,
89 "Thyra::EpetraLinearOp::initialize(...): Error!" );
90 // ToDo: Validate spmdRange, spmdDomain against op maps!
91#endif
92
93 RCP<const SPMDVSB> l_spmdRange;
94 if(!is_null(range_in))
95 l_spmdRange = rcp_dynamic_cast<const SPMDVSB>(range_in,true);
96 else
97 l_spmdRange = ( applyAs==EPETRA_OP_APPLY_APPLY
98 ? allocateRange(op,opTrans) : allocateDomain(op,opTrans) );
99
100 RCP<const SPMDVSB> l_spmdDomain;
101 if(!is_null(domain_in))
102 l_spmdDomain = rcp_dynamic_cast<const SPMDVSB>(domain_in,true);
103 else
104 l_spmdDomain = ( applyAs==EPETRA_OP_APPLY_APPLY
105 ? allocateDomain(op,opTrans) : allocateRange(op,opTrans) );
106
107 // Set data (no exceptions should be thrown now)
108 isFullyInitialized_ = true;
109 op_ = op;
110 rowMatrix_ = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(op_);
111 opTrans_ = opTrans;
112 applyAs_ = applyAs;
113 adjointSupport_ = adjointSupport;
114 range_ = l_spmdRange;
115 domain_ = l_spmdDomain;
116
117}
118
119
121 const RCP<const VectorSpaceBase<double> > &range_in,
122 const RCP<const VectorSpaceBase<double> > &domain_in,
123 const RCP<Epetra_Operator> &op,
124 EOpTransp opTrans,
125 EApplyEpetraOpAs applyAs,
126 EAdjointEpetraOp adjointSupport
127 )
128{
129
130 using Teuchos::rcp_dynamic_cast;
131 typedef SpmdVectorSpaceBase<double> SPMDVSB;
132
133 // Validate input, allocate spaces, validate ...
134#ifdef TEUCHOS_DEBUG
135 TEUCHOS_TEST_FOR_EXCEPTION( is_null(range_in), std::invalid_argument,
136 "Thyra::EpetraLinearOp::partiallyInitialize(...): Error!" );
137 TEUCHOS_TEST_FOR_EXCEPTION( is_null(domain_in), std::invalid_argument,
138 "Thyra::EpetraLinearOp::partiallyInitialize(...): Error!" );
139 TEUCHOS_TEST_FOR_EXCEPTION( is_null(op), std::invalid_argument,
140 "Thyra::EpetraLinearOp::partiallyInitialize(...): Error!" );
141#endif
142
144 l_spmdRange = rcp_dynamic_cast<const SPMDVSB>(range_in,true);
146 l_spmdDomain = rcp_dynamic_cast<const SPMDVSB>(domain_in,true);
147
148 // Set data (no exceptions should be thrown now)
149 isFullyInitialized_ = false;
150 op_ = op;
151 rowMatrix_ = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(op_);
152 opTrans_ = opTrans;
153 applyAs_ = applyAs;
154 adjointSupport_ = adjointSupport;
155 range_ = l_spmdRange;
156 domain_ = l_spmdDomain;
157
158}
159
160
161void EpetraLinearOp::setFullyInitialized(bool /* isFullyInitialized */)
162{
163 // ToDo: Validate that everything matches up!
164 isFullyInitialized_ = true;
165}
166
167
170 EOpTransp *opTrans,
171 EApplyEpetraOpAs *applyAs,
172 EAdjointEpetraOp *adjointSupport,
173 RCP<const VectorSpaceBase<double> > *range_out,
174 RCP<const VectorSpaceBase<double> > *domain_out
175 )
176{
177
178 if(op) *op = op_;
179 if(opTrans) *opTrans = opTrans_;
180 if(applyAs) *applyAs = applyAs_;
181 if(adjointSupport) *adjointSupport = adjointSupport_;
182 if(range_out) *range_out = range_;
183 if(domain_out) *domain_out = domain_;
184
185 isFullyInitialized_ = false;
186 op_ = Teuchos::null;
187 rowMatrix_ = Teuchos::null;
191 range_ = Teuchos::null;
192 domain_ = Teuchos::null;
193
194}
195
196
199{
200 return range_;
201}
202
203
206{
207 return domain_;
208}
209
210
213{
214 return op_;
215}
216
217
220{
221 return op_;
222}
223
224
225// Overridden from EpetraLinearOpBase
226
227
229 const Ptr<RCP<Epetra_Operator> > &epetraOp,
230 const Ptr<EOpTransp> &epetraOpTransp,
231 const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
232 const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
233 )
234{
235 *epetraOp = op_;
236 *epetraOpTransp = opTrans_;
237 *epetraOpApplyAs = applyAs_;
238 *epetraOpAdjointSupport = adjointSupport_;
239}
240
241
243 const Ptr<RCP<const Epetra_Operator> > &epetraOp,
244 const Ptr<EOpTransp> &epetraOpTransp,
245 const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
246 const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
247 ) const
248{
249 *epetraOp = op_;
250 *epetraOpTransp = opTrans_;
251 *epetraOpApplyAs = applyAs_;
252 *epetraOpAdjointSupport = adjointSupport_;
253}
254
255
256// Overridden from LinearOpBase
257
258
261{
262 return range_;
263}
264
265
268{
269 return domain_;
270}
271
272
275{
276 assert(0); // ToDo: Implement when needed
277 return Teuchos::null;
278}
279
280
281// Overridden from Teuchos::Describable
282
283
285{
286 std::ostringstream oss;
287 oss << Teuchos::Describable::description() << "{";
288 if(op_.get()) {
289 oss << "op=\'"<<typeName(*op_)<<"\'";
290 oss << ",rangeDim="<<this->range()->dim();
291 oss << ",domainDim="<<this->domain()->dim();
292 }
293 else {
294 oss << "op=NULL";
295 }
296 oss << "}";
297 return oss.str();
298}
299
300
302 FancyOStream &out,
303 const Teuchos::EVerbosityLevel verbLevel
304 ) const
305{
307 using Teuchos::as;
308 using Teuchos::rcp_dynamic_cast;
309 using Teuchos::OSTab;
310 using Teuchos::describe;
311 OSTab tab(out);
312 if ( as<int>(verbLevel) == as<int>(Teuchos::VERB_LOW) || is_null(op_)) {
313 out << this->description() << std::endl;
314 }
315 else if (includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM)) {
316 out
318 << "{"
319 << "rangeDim=" << this->range()->dim()
320 << ",domainDim=" << this->domain()->dim()
321 << "}\n";
322 OSTab tab2(out);
323 if (op_.get()) {
324 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
325 out << "opTrans="<<toString(opTrans_)<<"\n";
326 out << "applyAs="<<toString(applyAs_)<<"\n";
327 out << "adjointSupport="<<toString(adjointSupport_)<<"\n";
328 out << "op="<<typeName(*op_)<<"\n";
329 }
330 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
331 OSTab tab3(out);
333 csr_op = rcp_dynamic_cast<const Epetra_CrsMatrix>(op_);
334 if (!is_null(csr_op)) {
335 csr_op->Print(out);
336 }
337 }
338 }
339 else {
340 out << "op=NULL"<<"\n";
341 }
342 }
343}
344
345
346// protected
347
348
349// Protected member functions overridden from LinearOpBase
350
351
352bool EpetraLinearOp::opSupportedImpl(EOpTransp M_trans) const
353{
355 return false;
356 return ( M_trans == NOTRANS
358}
359
360
362 const EOpTransp M_trans,
363 const MultiVectorBase<double> &X_in,
364 const Ptr<MultiVectorBase<double> > &Y_inout,
365 const double alpha,
366 const double beta
367 ) const
368{
369
370 THYRA_FUNC_TIME_MONITOR("Thyra::EpetraLinearOp::euclideanApply");
371
372 const EOpTransp real_M_trans = real_trans(M_trans);
373
374#ifdef TEUCHOS_DEBUG
376 THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
377 "EpetraLinearOp::euclideanApply(...)", *this, M_trans, X_in, Y_inout
378 );
381 Exceptions::OpNotSupported,
382 "EpetraLinearOp::apply(...): *this was informed that adjoints "
383 "are not supported when initialized."
384 );
385#endif
386
387 const RCP<const VectorSpaceBase<double> > XY_domain = X_in.domain();
388 const int numCols = XY_domain->dim();
389
390 //
391 // Get Epetra_MultiVector objects for the arguments
392 //
393 // 2007/08/18: rabartl: Note: After profiling, I found that calling the more
394 // general functions get_Epetra_MultiVector(...) was too slow. These
395 // functions must ensure that memory is being remembered efficiently and the
396 // use of extra data with the RCP and other things is slow.
397 //
400 {
401 THYRA_FUNC_TIME_MONITOR_DIFF(
402 "Thyra::EpetraLinearOp::euclideanApply: Convert MultiVectors", MultiVectors);
403 // X
405 real_M_trans==NOTRANS ? getDomainMap() : getRangeMap(), X_in );
406 // Y
407 if( beta == 0 ) {
409 real_M_trans==NOTRANS ? getRangeMap() : getDomainMap(), *Y_inout );
410 }
411 }
412
413 //
414 // Set the operator mode
415 //
416
417 /* We need to save the transpose state here, and then reset it after
418 * application. The reason for this is that if we later apply the
419 * operator outside Thyra (in Aztec, for instance), it will remember
420 * the transpose flag set here. */
421 bool oldState = op_->UseTranspose();
422 op_->SetUseTranspose(
423 real_trans(trans_trans(opTrans_,M_trans)) == NOTRANS ? false : true );
424
425 //
426 // Perform the apply operation
427 //
428 {
429 THYRA_FUNC_TIME_MONITOR_DIFF(
430 "Thyra::EpetraLinearOp::euclideanApply: Apply", Apply);
431 if( beta == 0.0 ) {
432 // Y = M * X
434 THYRA_FUNC_TIME_MONITOR_DIFF(
435 "Thyra::EpetraLinearOp::euclideanApply: Apply(beta==0): Apply",
436 ApplyApply);
437 op_->Apply( *X, *Y );
438 }
440 THYRA_FUNC_TIME_MONITOR_DIFF(
441 "Thyra::EpetraLinearOp::euclideanApply: Apply(beta==0): ApplyInverse",
442 ApplyApplyInverse);
443 op_->ApplyInverse( *X, *Y );
444 }
445 else {
446#ifdef TEUCHOS_DEBUG
448#endif
449 }
450 // Y = alpha * Y
451 if( alpha != 1.0 ) {
452 THYRA_FUNC_TIME_MONITOR_DIFF(
453 "Thyra::EpetraLinearOp::euclideanApply: Apply(beta==0): Scale Y",
454 Scale);
455 Y->Scale(alpha);
456 }
457 }
458 else { // beta != 0.0
459 // Y_inout = beta * Y_inout
460 if(beta != 0.0) {
461 THYRA_FUNC_TIME_MONITOR_DIFF(
462 "Thyra::EpetraLinearOp::euclideanApply: Apply(beta!=0): Scale Y",
463 Scale);
464 scale( beta, Y_inout );
465 }
466 else {
467 THYRA_FUNC_TIME_MONITOR_DIFF(
468 "Thyra::EpetraLinearOp::euclideanApply: Apply(beta!=0): Y=0",
469 Apply2);
470 assign( Y_inout, 0.0 );
471 }
472 // T = M * X
473 Epetra_MultiVector T(op_->OperatorRangeMap(), numCols, false);
474 // NOTE: Above, op_->OperatorRange() will be right for either
475 // non-transpose or transpose because we have already set the
476 // UseTranspose flag correctly.
478 THYRA_FUNC_TIME_MONITOR_DIFF(
479 "Thyra::EpetraLinearOp::euclideanApply: Apply(beta!=0): Apply",
480 Apply2);
481 op_->Apply( *X, T );
482 }
484 THYRA_FUNC_TIME_MONITOR_DIFF(
485 "Thyra::EpetraLinearOp::euclideanApply: Apply(beta!=0): ApplyInverse",
486 ApplyInverse);
487 op_->ApplyInverse( *X, T );
488 }
489 else {
490#ifdef TEUCHOS_DEBUG
492#endif
493 }
494 // Y_inout += alpha * T
495 {
496 THYRA_FUNC_TIME_MONITOR_DIFF(
497 "Thyra::EpetraLinearOp::euclideanApply: Apply(beta!=0): Update Y",
498 Update);
499 update(
500 alpha,
503 Y_inout->range(),
504 XY_domain
505 ),
506 Y_inout
507 );
508 }
509 }
510 }
511
512 // Reset the transpose state
513 op_->SetUseTranspose(oldState);
514
515 // 2009/04/14: ToDo: This will not reset the transpose flag correctly if an
516 // exception is thrown!
517
518}
519
520
521// Protected member functions overridden from ScaledLinearOpBase
522
523
525{
526 return nonnull(rowMatrix_);
527}
528
529
531{
532 return nonnull(rowMatrix_);
533}
534
535
536void EpetraLinearOp::scaleLeftImpl(const VectorBase<double> &row_scaling_in)
537{
538 using Teuchos::rcpFromRef;
539 const RCP<const Epetra_Vector> row_scaling =
540 get_Epetra_Vector(getRangeMap(), rcpFromRef(row_scaling_in));
541 rowMatrix_->LeftScale(*row_scaling);
542}
543
544
545void EpetraLinearOp::scaleRightImpl(const VectorBase<double> &col_scaling_in)
546{
547 using Teuchos::rcpFromRef;
548 const RCP<const Epetra_Vector> col_scaling =
549 get_Epetra_Vector(getDomainMap(), rcpFromRef(col_scaling_in));
550 rowMatrix_->RightScale(*col_scaling);
551}
552
553
554// Protected member functions overridden from RowStatLinearOpBase
555
556
558 const RowStatLinearOpBaseUtils::ERowStat rowStat) const
559{
560 if (is_null(rowMatrix_)) {
561 return false;
562 }
563 switch (rowStat) {
564 case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
565 case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
566 return true;
567 default:
569 }
571}
572
573
575 const RowStatLinearOpBaseUtils::ERowStat rowStat,
576 const Ptr<VectorBase<double> > &rowStatVec_in
577 ) const
578{
579 using Teuchos::rcpFromPtr;
580 const RCP<Epetra_Vector> rowStatVec =
581 get_Epetra_Vector(getRangeMap(), rcpFromPtr(rowStatVec_in));
582 switch (rowStat) {
583 case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
584 rowMatrix_->InvRowSums(*rowStatVec);
585 break;
586 case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
587 // compute absolute row sum
588 computeAbsRowSum(*rowStatVec);
589 break;
590 default:
592 }
593}
594
595
596// Allocators for domain and range spaces
597
598
601 const RCP<Epetra_Operator> &op,
602 EOpTransp /* op_trans */
603 ) const
604{
605 return Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
606 create_VectorSpace(Teuchos::rcp(&op->OperatorDomainMap(),false))
607 );
608 // ToDo: What about the transpose argument???, test this!!!
609}
610
611
614 const RCP<Epetra_Operator> &op,
615 EOpTransp /* op_trans */
616 ) const
617{
618 return Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
619 create_VectorSpace(Teuchos::rcp(&op->OperatorRangeMap(),false))
620 );
621 // ToDo: What about the transpose argument???, test this!!!
622}
623
624// private
625
626
627const Epetra_Map& EpetraLinearOp::getRangeMap() const
628{
630 ? op_->OperatorRangeMap() : op_->OperatorDomainMap() );
631 // ToDo: What about the transpose argument???, test this!!!
632}
633
634
635const Epetra_Map& EpetraLinearOp::getDomainMap() const
636{
638 ? op_->OperatorDomainMap() : op_->OperatorRangeMap() );
639 // ToDo: What about the transpose argument???, test this!!!
640}
641
642void EpetraLinearOp::computeAbsRowSum(Epetra_Vector & x) const
643{
645
646 RCP<Epetra_CrsMatrix> crsMatrix
647 = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(rowMatrix_);
648
650 Exceptions::OpNotSupported,
651 "EpetraLinearOp::computeAbsRowSum(...): wrapped matrix must be of type "
652 "Epetra_CrsMatrix for this method. Other operator types are not supported."
653 );
654
655 //
656 // Put inverse of the sum of absolute values of the ith row of A in x[i].
657 // (this is a modified copy of Epetra_CrsMatrix::InvRowSums)
658 //
659
660 if (crsMatrix->Filled()) {
662 std::invalid_argument,
663 "EpetraLinearOp::computeAbsRowSum(...): Epetra_CrsMatrix must be filled"
664 );
665 }
666 int i, j;
667 x.PutScalar(0.0); // Make sure we sum into a vector of zeros.
668 double * xp = (double*)x.Values();
669 if (crsMatrix->Graph().RangeMap().SameAs(x.Map()) && crsMatrix->Exporter() != 0) {
670 Epetra_Vector x_tmp(crsMatrix->RowMap());
671 x_tmp.PutScalar(0.0);
672 double * x_tmp_p = (double*)x_tmp.Values();
673 for (i=0; i < crsMatrix->NumMyRows(); i++) {
674 int NumEntries = 0;
675 double * RowValues = 0;
676 crsMatrix->ExtractMyRowView(i,NumEntries,RowValues);
677 for (j=0; j < NumEntries; j++) x_tmp_p[i] += std::abs(RowValues[j]);
678 }
679 TEUCHOS_TEST_FOR_EXCEPT(0!=x.Export(x_tmp, *crsMatrix->Exporter(), Add)); //Export partial row sums to x.
680 }
681 else if (crsMatrix->Graph().RowMap().SameAs(x.Map())) {
682 for (i=0; i < crsMatrix->NumMyRows(); i++) {
683 int NumEntries = 0;
684 double * RowValues = 0;
685 crsMatrix->ExtractMyRowView(i,NumEntries,RowValues);
686 double scale = 0.0;
687 for (j=0; j < NumEntries; j++) scale += std::abs(RowValues[j]);
688 xp[i] = scale;
689 }
690 }
691 else { // x.Map different than both crsMatrix->Graph().RowMap() and crsMatrix->Graph().RangeMap()
692 TEUCHOS_TEST_FOR_EXCEPT(true); // The map of x must be the RowMap or RangeMap of A.
693 }
694}
695
696} // end namespace Thyra
697
698
699// Nonmembers
700
701
703Thyra::nonconstEpetraLinearOp()
704{
705 return Teuchos::rcp(new EpetraLinearOp());
706}
707
708
710Thyra::partialNonconstEpetraLinearOp(
711 const RCP<const VectorSpaceBase<double> > &range,
712 const RCP<const VectorSpaceBase<double> > &domain,
713 const RCP<Epetra_Operator> &op,
714 EOpTransp opTrans,
715 EApplyEpetraOpAs applyAs,
716 EAdjointEpetraOp adjointSupport
717 )
718{
719 RCP<EpetraLinearOp> thyraEpetraOp = Teuchos::rcp(new EpetraLinearOp());
720 thyraEpetraOp->partiallyInitialize(
721 range, domain,op,opTrans, applyAs, adjointSupport
722 );
723 return thyraEpetraOp;
724}
725
726
728Thyra::nonconstEpetraLinearOp(
729 const RCP<Epetra_Operator> &op,
730 EOpTransp opTrans,
731 EApplyEpetraOpAs applyAs,
732 EAdjointEpetraOp adjointSupport,
733 const RCP< const VectorSpaceBase<double> > &range,
734 const RCP< const VectorSpaceBase<double> > &domain
735 )
736{
737 RCP<EpetraLinearOp> thyraEpetraOp = Teuchos::rcp(new EpetraLinearOp());
738 thyraEpetraOp->initialize(
739 op,opTrans, applyAs, adjointSupport, range, domain
740 );
741 return thyraEpetraOp;
742}
743
744
746Thyra::epetraLinearOp(
747 const RCP<const Epetra_Operator> &op,
748 EOpTransp opTrans,
749 EApplyEpetraOpAs applyAs,
750 EAdjointEpetraOp adjointSupport,
751 const RCP<const VectorSpaceBase<double> > &range,
752 const RCP<const VectorSpaceBase<double> > &domain
753 )
754{
755 RCP<EpetraLinearOp> thyraEpetraOp = Teuchos::rcp(new EpetraLinearOp());
756 thyraEpetraOp->initialize(
757 Teuchos::rcp_const_cast<Epetra_Operator>(op), // Safe cast due to return type!
758 opTrans, applyAs, adjointSupport, range, domain
759 );
760 return thyraEpetraOp;
761}
762
763
765Thyra::nonconstEpetraLinearOp(
766 const RCP<Epetra_Operator> &op,
767 const std::string &label,
768 EOpTransp opTrans,
769 EApplyEpetraOpAs applyAs,
770 EAdjointEpetraOp adjointSupport,
771 const RCP<const VectorSpaceBase<double> > &range,
772 const RCP<const VectorSpaceBase<double> > &domain
773 )
774{
775 RCP<EpetraLinearOp> thyraEpetraOp = Teuchos::rcp(new EpetraLinearOp());
776 thyraEpetraOp->initialize(
777 op,opTrans, applyAs, adjointSupport, range, domain
778 );
779 thyraEpetraOp->setObjectLabel(label);
780 return thyraEpetraOp;
781}
782
783
785Thyra::epetraLinearOp(
786 const RCP<const Epetra_Operator> &op,
787 const std::string &label,
788 EOpTransp opTrans,
789 EApplyEpetraOpAs applyAs,
790 EAdjointEpetraOp adjointSupport,
791 const RCP< const SpmdVectorSpaceBase<double> > &range,
792 const RCP< const SpmdVectorSpaceBase<double> > &domain
793 )
794{
795 RCP<EpetraLinearOp> thyraEpetraOp = Teuchos::rcp(new EpetraLinearOp());
796 thyraEpetraOp->initialize(
797 Teuchos::rcp_const_cast<Epetra_Operator>(op), // Safe cast due to return type!
798 opTrans, applyAs, adjointSupport, range, domain
799 );
800 thyraEpetraOp->setObjectLabel(label);
801 return thyraEpetraOp;
802}
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
virtual std::string description() const
T * get() const
RCP< const SpmdVectorSpaceBase< double > > spmdDomain() const
Return a smart pointer to the SpmdVectorSpaceBase object for the domain.
std::string description() const
virtual void scaleLeftImpl(const VectorBase< double > &row_scaling)
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< double > &X, const Ptr< MultiVectorBase< double > > &Y, const double alpha, const double beta) const
virtual void scaleRightImpl(const VectorBase< double > &col_scaling)
RCP< Epetra_Operator > op_
EAdjointEpetraOp adjointSupport_
void getNonconstEpetraOpView(const Ptr< RCP< Epetra_Operator > > &epetraOp, const Ptr< EOpTransp > &epetraOpTransp, const Ptr< EApplyEpetraOpAs > &epetraOpApplyAs, const Ptr< EAdjointEpetraOp > &epetraOpAdjointSupport)
void getEpetraOpView(const Ptr< RCP< const Epetra_Operator > > &epetraOp, const Ptr< EOpTransp > &epetraOpTransp, const Ptr< EApplyEpetraOpAs > &epetraOpApplyAs, const Ptr< EAdjointEpetraOp > &epetraOpAdjointSupport) const
RCP< Epetra_RowMatrix > rowMatrix_
void uninitialize(RCP< Epetra_Operator > *op=NULL, EOpTransp *opTrans=NULL, EApplyEpetraOpAs *applyAs=NULL, EAdjointEpetraOp *adjointSupport=NULL, RCP< const VectorSpaceBase< double > > *range=NULL, RCP< const VectorSpaceBase< double > > *domain=NULL)
Set to uninitialized and optionally return the current state.
EpetraLinearOp()
Construct to uninitialized.
virtual bool supportsScaleRightImpl() const
void computeAbsRowSum(Epetra_Vector &rowStatVec_in) const
Compute the absolute row sum for this matrix.
RCP< const SpmdVectorSpaceBase< double > > domain_
void setFullyInitialized(bool isFullyInitialized=true)
Set to fully initialized.
const Epetra_Map & getDomainMap() const
virtual void getRowStatImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat, const Ptr< VectorBase< double > > &rowStatVec) const
void describe(FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual RCP< const SpmdVectorSpaceBase< double > > allocateRange(const RCP< Epetra_Operator > &op, EOpTransp op_trans) const
Allocate the range space of the operator.
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.
RCP< const LinearOpBase< double > > clone() const
void partiallyInitialize(const RCP< const VectorSpaceBase< double > > &range, const RCP< const VectorSpaceBase< double > > &domain, const RCP< Epetra_Operator > &op, EOpTransp opTrans=NOTRANS, EApplyEpetraOpAs applyAs=EPETRA_OP_APPLY_APPLY, EAdjointEpetraOp adjointSupport=EPETRA_OP_ADJOINT_SUPPORTED)
Partially initialize.
RCP< const SpmdVectorSpaceBase< double > > spmdRange() const
Return a smart pointer to the SpmdVectorSpaceBase object for the range.
virtual bool rowStatIsSupportedImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat) const
RCP< Epetra_Operator > epetra_op()
virtual RCP< const SpmdVectorSpaceBase< double > > allocateDomain(const RCP< Epetra_Operator > &op, EOpTransp op_trans) const
Allocate the domain space of the operator.
bool opSupportedImpl(EOpTransp M_trans) const
RCP< const VectorSpaceBase< double > > range() const
const Epetra_Map & getRangeMap() const
RCP< const SpmdVectorSpaceBase< double > > range_
virtual bool supportsScaleLeftImpl() const
RCP< const VectorSpaceBase< double > > domain() const
EApplyEpetraOpAs
Determine how the apply an Epetra_Operator as a linear operator.
RCP< MultiVectorBase< double > > create_MultiVector(const RCP< Epetra_MultiVector > &epetra_mv, const RCP< const VectorSpaceBase< double > > &range=Teuchos::null, const RCP< const VectorSpaceBase< double > > &domain=Teuchos::null)
Create a non-const MultiVectorBase object from a non-const Epetra_MultiVector object.
RCP< const VectorSpaceBase< double > > create_VectorSpace(const RCP< const Epetra_Map > &epetra_map)
Create an VectorSpaceBase object given an Epetra_Map object.
const std::string toString(const EAdjointEpetraOp adjointEpetraOp)
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.
RCP< Epetra_MultiVector > get_Epetra_MultiVector(const Epetra_Map &map, const RCP< MultiVectorBase< double > > &mv)
Get a non-const Epetra_MultiVector view from a non-const MultiVectorBase object if possible.
EAdjointEpetraOp
Determine if adjoints are supported on Epetra_Opeator or not.
@ EPETRA_OP_APPLY_APPLY
Apply using Epetra_Operator::Apply(...)
@ EPETRA_OP_APPLY_APPLY_INVERSE
Apply using Epetra_Operator::ApplyInverse(...)
@ EPETRA_OP_ADJOINT_UNSUPPORTED
Adjoint not supported.
@ EPETRA_OP_ADJOINT_SUPPORTED
Adjoint supported.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
bool is_null(const boost::shared_ptr< T > &p)
bool nonnull(const boost::shared_ptr< T > &p)
NOTRANS
TypeTo as(const TypeFrom &t)
std::string typeName(const T &t)
basic_FancyOStream< char > FancyOStream
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
basic_OSTab< char > OSTab
TRANS
const T & getConst(T &t)
TEUCHOSCORE_LIB_DLL_EXPORT bool includesVerbLevel(const EVerbosityLevel verbLevel, const EVerbosityLevel requestedVerbLevel, const bool isDefaultLevel=false)