Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DefaultBlockedTriangularLinearOpWithSolve_def.hpp
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#ifndef THYRA_DEFAULT_BLOCKED_TRIANGULAR_LINEAR_OP_WITH_SOLVE_DEF_HPP
43#define THYRA_DEFAULT_BLOCKED_TRIANGULAR_LINEAR_OP_WITH_SOLVE_DEF_HPP
44
45
46#include "Thyra_DefaultBlockedTriangularLinearOpWithSolve_decl.hpp"
47#include "Thyra_ProductMultiVectorBase.hpp"
48#include "Thyra_DefaultProductVectorSpace.hpp"
49#include "Thyra_AssertOp.hpp"
50
51
52namespace Thyra {
53
54
55// public
56
57
58// Constructors/Initializers/Accessors
59
60
61template<class Scalar>
65
66
67template<class Scalar>
70 )
71{
72 assertAndSetBlockStructure(*blocks);
73 blocks_.initialize(blocks);
74}
75
76
77template<class Scalar>
79 const RCP<const PhysicallyBlockedLinearOpBase<Scalar> > &blocks
80 )
81{
82 assertAndSetBlockStructure(*blocks);
83 blocks_.initialize(blocks);
84}
85
86
87template<class Scalar>
93
94
95template<class Scalar>
98{
99 return blocks_.getConstObj();
100}
101
102
103// Overridden from PhysicallyBlockedLinearOpWithSolveBase
104
105
106template<class Scalar>
108 const int i, const int j
109 ) const
110{
111 assertBlockFillIsActive(true);
112 assertBlockRowCol(i,j);
113 return i==j; // Only accept LOWS blocks along the diagonal!
114}
115
116template<class Scalar>
118 const int i, const int j,
120 )
121{
122 setLOWSBlockImpl(i,j,block);
123}
124
125
126template<class Scalar>
128 const int i, const int j,
130 )
131{
132 setLOWSBlockImpl(i,j,block);
133}
134
135
136// Overridden from PhysicallyBlockedLinearOpBase
137
138
139template<class Scalar>
141{
142 assertBlockFillIsActive(false);
143 TEUCHOS_TEST_FOR_EXCEPT("ToDo: Have not implemented flexible block fill yet!");
144}
145
146
147template<class Scalar>
149 const int numRowBlocks, const int numColBlocks
150 )
151{
152 using Teuchos::null;
153#ifdef THYRA_DEBUG
154 TEUCHOS_ASSERT_EQUALITY(numRowBlocks, numColBlocks);
155#else
156 (void)numColBlocks;
157#endif
158 assertBlockFillIsActive(false);
159 numDiagBlocks_ = numRowBlocks;
160 diagonalBlocks_.resize(numDiagBlocks_);
161 productRange_ = null;
162 productDomain_ = null;
163 blockFillIsActive_ = true;
164}
165
166
167template<class Scalar>
169 const Teuchos::RCP<const ProductVectorSpaceBase<Scalar> > &productRange_in,
170 const Teuchos::RCP<const ProductVectorSpaceBase<Scalar> > &productDomain_in
171 )
172{
173#ifdef THYRA_DEBUG
174 TEUCHOS_TEST_FOR_EXCEPT( is_null(productRange_in) );
175 TEUCHOS_TEST_FOR_EXCEPT( is_null(productDomain_in) );
176 TEUCHOS_TEST_FOR_EXCEPT( productRange_in->numBlocks() != productDomain_in->numBlocks() );
177#endif
178 assertBlockFillIsActive(false);
179 productRange_ = productRange_in;
180 productDomain_ = productDomain_in;
181 numDiagBlocks_ = productRange_in->numBlocks();
182 diagonalBlocks_.resize(numDiagBlocks_);
183 blockFillIsActive_ = true;
184}
185
186
187template<class Scalar>
189{
190 return blockFillIsActive_;
191}
192
193
194template<class Scalar>
196 const int i, const int j
197 ) const
198{
199 assertBlockFillIsActive(true);
200 assertBlockRowCol(i,j);
201 return false; // ToDo: Change this once we accept off-diagonal blocks
202}
203
204
205template<class Scalar>
207 const int /* i */, const int /* j */,
208 const Teuchos::RCP<LinearOpBase<Scalar> > &/* block */
209 )
210{
211 assertBlockFillIsActive(true);
212 TEUCHOS_TEST_FOR_EXCEPT("Error, we don't support off-diagonal LOB objects yet!");
213}
214
215
216template<class Scalar>
218 const int /* i */, const int /* j */,
219 const Teuchos::RCP<const LinearOpBase<Scalar> > &/* block */
220 )
221{
222 assertBlockFillIsActive(true);
223 TEUCHOS_TEST_FOR_EXCEPT("Error, we don't support off-diagonal LOB objects yet!");
224}
225
226
227template<class Scalar>
229{
230 assertBlockFillIsActive(true);
233 for ( int k = 0; k < numDiagBlocks_; ++k ) {
235 diagonalBlocks_[k].getConstObj();
236 TEUCHOS_TEST_FOR_EXCEPTION(is_null(lows_k), std::logic_error,
237 "Error, the block diagonal k="<<k<<" can not be null when ending block fill!"
238 );
239 if (is_null(productRange_)) {
240 rangeSpaces.push_back(lows_k->range());
241 domainSpaces.push_back(lows_k->domain());
242 }
243 }
244 if (is_null(productRange_)) {
245 productRange_ = productVectorSpace<Scalar>(rangeSpaces());
246 productDomain_ = productVectorSpace<Scalar>(domainSpaces());
247 }
248 blockFillIsActive_ = false;
249}
250
251
252template<class Scalar>
254{
255 assertBlockFillIsActive(false);
256 productRange_ = Teuchos::null;
257 productDomain_ = Teuchos::null;
258 numDiagBlocks_ = 0;
259 diagonalBlocks_.resize(0);
260}
261
262
263// Overridden from BlockedLinearOpWithSolveBase
264
265
266template<class Scalar>
269 const int i, const int j
270 )
271{
272 assertBlockFillIsActive(false);
273 assertBlockRowCol(i,j);
274 if (i!=j)
275 return Teuchos::null;
276 return diagonalBlocks_[i].getNonconstObj();
277}
278
279
280template<class Scalar>
283 const int i, const int j
284 ) const
285{
286 assertBlockFillIsActive(false);
287 assertBlockRowCol(i, j);
288 if (i != j)
289 return Teuchos::null;
290 return diagonalBlocks_[i].getConstObj();
291}
292
293
294// Overridden from BlockedLinearOpBase
295
296
297template<class Scalar>
300{
301 return productRange_;
302}
303
304
305template<class Scalar>
308{
309 return productDomain_;
310}
311
312
313template<class Scalar>
315 const int i, const int j
316 ) const
317{
318 assertBlockFillIsActive(false);
319 assertBlockRowCol(i,j);
320 if (i!=j)
321 return false; // ToDo: Update this when off-diagonals are supported!
322 return !is_null(diagonalBlocks_[i].getConstObj());
323}
324
325
326template<class Scalar>
328 const int i, const int j
329 ) const
330{
331 assertBlockFillIsActive(true);
332 assertBlockRowCol(i,j);
333 return diagonalBlocks_[i].isConst();
334}
335
336
337template<class Scalar>
340 const int i, const int j
341 )
342{
343 assertBlockFillIsActive(true);
344 assertBlockRowCol(i,j);
345 if (i!=j)
346 return Teuchos::null; // ToDo: Update when off-diagonals are supported!
347 return this->getNonconstLOWSBlock(i,j);
348}
349
350
351template<class Scalar>
354 const int i, const int j
355 ) const
356{
357 assertBlockFillIsActive(true);
358 assertBlockRowCol(i,j);
359 if (i!=j)
360 return Teuchos::null; // ToDo: Update when off-diagonals are supported!
361 return this->getLOWSBlock(i,j);
362}
363
364
365// Overridden from LinearOpBase
366
367
368template<class Scalar>
371{
372 return this->productRange();
373}
374
375
376template<class Scalar>
379{
380 return this->productDomain();
381}
382
383
384template<class Scalar>
387{
388 return Teuchos::null; // ToDo: Implement clone when needed!
389}
390
391
392// Overridden from Teuchos::Describable
393
394
395template<class Scalar>
396std::string
398{
399 assertBlockFillIsActive(false);
400 std::ostringstream oss;
401 oss
403 << "numDiagBlocks="<<numDiagBlocks_
404 << "}";
405 return oss.str();
406}
407
408
409template<class Scalar>
412 const Teuchos::EVerbosityLevel verbLevel
413 ) const
414{
415 assertBlockFillIsActive(false);
416 Teuchos::Describable::describe(out, verbLevel);
417 // ToDo: Fill in a better version of this!
418}
419
420
421// protected
422
423
424// Overridden from LinearOpBase
425
426
427template<class Scalar>
429 EOpTransp M_trans
430 ) const
431{
432 using Thyra::opSupported;
433 assertBlockFillIsActive(false);
434 for ( int k = 0; k < numDiagBlocks_; ++k ) {
435 if ( !opSupported(*diagonalBlocks_[k].getConstObj(),M_trans) )
436 return false;
437 }
438 return true;
439 // ToDo: To be safe we really should do a collective reduction of all
440 // clusters of processes. However, for the typical use case, every block
441 // will return the same info and we should be safe!
442}
443
444
445template<class Scalar>
447 const EOpTransp M_trans,
448 const MultiVectorBase<Scalar> &X_in,
449 const Ptr<MultiVectorBase<Scalar> > &Y_inout,
450 const Scalar alpha,
451 const Scalar beta
452 ) const
453{
454
455 using Teuchos::RCP;
456 using Teuchos::dyn_cast;
457 using Thyra::apply;
458
459#ifdef THYRA_DEBUG
461 "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::apply(...)",
462 *this, M_trans, X_in, &*Y_inout
463 );
464#endif // THYRA_DEBUG
465
466 //
467 // Y = alpha * op(M) * X + beta*Y
468 //
469 // =>
470 //
471 // Y[i] = beta+Y[i] + alpha*op(Op)[i]*X[i], for i=0...numDiagBlocks-1
472 //
473 // ToDo: Update to handle off diagonal blocks when needed!
474 //
475
477 &X = dyn_cast<const ProductMultiVectorBase<Scalar> >(X_in);
479 &Y = dyn_cast<ProductMultiVectorBase<Scalar> >(*Y_inout);
480
481 for ( int i = 0; i < numDiagBlocks_; ++ i ) {
482 Thyra::apply( *diagonalBlocks_[i].getConstObj(), M_trans,
483 *X.getMultiVectorBlock(i), Y.getNonconstMultiVectorBlock(i).ptr(),
484 alpha, beta
485 );
486 }
487
488}
489
490
491// Overridden from LinearOpWithSolveBase
492
493
494template<class Scalar>
495bool
497 EOpTransp M_trans
498 ) const
499{
500 assertBlockFillIsActive(false);
501 for ( int k = 0; k < numDiagBlocks_; ++k ) {
502 if ( !Thyra::solveSupports( *diagonalBlocks_[k].getConstObj(), M_trans ) )
503 return false;
504 }
505 return true;
506 // ToDo: To be safe we really should do a collective reduction of all
507 // clusters of processes. However, for the typical use case, every block
508 // will return the same info and we should be safe!
509}
510
511
512template<class Scalar>
513bool
515 EOpTransp M_trans, const SolveMeasureType& solveMeasureType
516 ) const
517{
518 using Thyra::solveSupportsSolveMeasureType;
519 assertBlockFillIsActive(false);
520 for ( int k = 0; k < numDiagBlocks_; ++k ) {
521 if (
522 !solveSupportsSolveMeasureType(
523 *diagonalBlocks_[k].getConstObj(),
524 M_trans, solveMeasureType
525 )
526 )
527 {
528 return false;
529 }
530 }
531 return true;
532}
533
534
535template<class Scalar>
538 const EOpTransp M_trans,
539 const MultiVectorBase<Scalar> &B_in,
540 const Ptr<MultiVectorBase<Scalar> > &X_inout,
541 const Ptr<const SolveCriteria<Scalar> > solveCriteria
542 ) const
543{
544
545 using Teuchos::RCP;
546 using Teuchos::dyn_cast;
547 using Thyra::solve;
548
549#ifdef THYRA_DEBUG
551 "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::apply(...)",
552 *this, M_trans, *X_inout, &B_in
553 );
554 TEUCHOS_TEST_FOR_EXCEPT(!this->solveSupportsImpl(M_trans));
555 // TEUCHOS_TEST_FOR_EXCEPTION(
556 // nonnull(solveCriteria) && !solveCriteria->solveMeasureType.useDefault(),
557 // std::logic_error,
558 // "Error, we can't handle any non-default solve criteria yet!"
559 // );
560 // ToDo: If solve criteria is to be handled, then we will have to be very
561 // carefull how it it interpreted in terms of the individual period solves!
562#endif // THYRA_DEBUG
563
564 //
565 // Y = alpha * inv(op(M)) * X + beta*Y
566 //
567 // =>
568 //
569 // X[i] = inv(op(Op[i]))*B[i], for i=0...numDiagBlocks-1
570 //
571 // ToDo: Update to handle off diagonal blocks when needed!
572 //
573
575 &B = dyn_cast<const ProductMultiVectorBase<Scalar> >(B_in);
577 &X = dyn_cast<ProductMultiVectorBase<Scalar> >(*X_inout);
578
579 bool converged = true;
580 for ( int i = 0; i < numDiagBlocks_; ++ i ) {
582 Op_k = diagonalBlocks_[i].getConstObj();
583 Op_k->setOStream(this->getOStream());
584 Op_k->setVerbLevel(this->getVerbLevel());
585 SolveStatus<Scalar> status =
586 Thyra::solve( *Op_k, M_trans, *B.getMultiVectorBlock(i),
587 X.getNonconstMultiVectorBlock(i).ptr(), solveCriteria );
589 converged = false;
590 // ToDo: Pass in solve criteria when needed!
591 }
592
593 SolveStatus<Scalar> solveStatus;
594 solveStatus.solveStatus =
596
597 return solveStatus;
598
599}
600
601
602
603// private
604
605
606template<class Scalar>
608 bool blockFillIsActive_in
609 ) const
610{
611#ifdef THYRA_DEBUG
612 TEUCHOS_TEST_FOR_EXCEPT(!(blockFillIsActive_==blockFillIsActive_in));
613#else
614 (void)blockFillIsActive_in;
615#endif
616}
617
618
619template<class Scalar>
620void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertBlockRowCol(
621 const int i, const int j
622 ) const
623{
624#ifdef THYRA_DEBUG
626 !( 0 <= i && i < numDiagBlocks_ ), std::logic_error,
627 "Error, i="<<i<<" does not fall in the range [0,"<<numDiagBlocks_-1<<"]!"
628 );
630 !( 0 <= j && j < numDiagBlocks_ ), std::logic_error,
631 "Error, j="<<j<<" does not fall in the range [0,"<<numDiagBlocks_-1<<"]!"
632 );
633#else
634 (void)i;
635 (void)j;
636#endif
637}
638
639
640template<class Scalar>
641template<class LinearOpWithSolveType>
642void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::setLOWSBlockImpl(
643 const int i, const int j,
645 )
646{
647 assertBlockFillIsActive(true);
648#ifdef THYRA_DEBUG
649 TEUCHOS_ASSERT_INEQUALITY( i, >=, 0 );
650 TEUCHOS_ASSERT_INEQUALITY( j, >=, 0 );
651 TEUCHOS_ASSERT_INEQUALITY( i, <, numDiagBlocks_ );
652 TEUCHOS_ASSERT_INEQUALITY( j, <, numDiagBlocks_ );
654 !this->acceptsLOWSBlock(i,j), std::logic_error,
655 "Error, this DefaultBlockedTriangularLinearOpWithSolve does not accept\n"
656 "LOWSB objects for the block i="<<i<<", j="<<j<<"!"
657 );
658#else
659 (void)j;
660#endif
661 diagonalBlocks_[i] = block;
662}
663
664
665template<class Scalar>
666void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertAndSetBlockStructure(
667 const PhysicallyBlockedLinearOpBase<Scalar>& blocks
668 )
669{
670#ifdef THYRA_DEBUG
672 "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertAndSetBlockStructure(blocks)",
673 *blocks.range(), *this->range()
674 );
676 "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertAndSetBlockStructure(blocks)",
677 *blocks.domain(), *this->domain()
678 );
679 // ToDo: Make sure that all of the blocks are above or below the diagonal
680 // but not both!
681#else
682 (void)blocks;
683#endif
684 // ToDo: Set if this is an upper or lower triangular block operator.
685}
686
687
688} // namespace Thyra
689
690
691#endif // THYRA_DEFAULT_BLOCKED_TRIANGULAR_LINEAR_OP_WITH_SOLVE_DEF_HPP
void push_back(const value_type &x)
virtual void describe(FancyOStream &out, const EVerbosityLevel verbLevel=verbLevel_default) const
virtual std::string description() const
Concrete composite LinearOpWithSolveBase subclass that creates single upper or lower block triangular...
void setNonconstBlocks(const RCP< PhysicallyBlockedLinearOpBase< Scalar > > &blocks)
SolveStatus< Scalar > solveImpl(const EOpTransp transp, const MultiVectorBase< Scalar > &B, const Ptr< MultiVectorBase< Scalar > > &X, const Ptr< const SolveCriteria< Scalar > > solveCriteria) const
void setLOWSBlock(const int i, const int j, const RCP< const LinearOpWithSolveBase< Scalar > > &block)
RCP< const LinearOpWithSolveBase< Scalar > > getLOWSBlock(const int i, const int j) const
RCP< const LinearOpBase< Scalar > > getBlock(const int i, const int j) const
void setNonconstBlock(const int i, const int j, const RCP< LinearOpBase< Scalar > > &block)
void setBlocks(const RCP< const PhysicallyBlockedLinearOpBase< Scalar > > &blocks)
bool solveSupportsSolveMeasureTypeImpl(EOpTransp M_trans, const SolveMeasureType &solveMeasureType) const
std::string description() const
Prints just the name DefaultBlockedTriangularLinearOpWithSolve along with the overall dimensions and ...
void setNonconstLOWSBlock(const int i, const int j, const RCP< LinearOpWithSolveBase< Scalar > > &block)
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Prints the details about the constituent linear operators.
void setBlock(const int i, const int j, const RCP< const LinearOpBase< Scalar > > &block)
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
RCP< LinearOpWithSolveBase< Scalar > > getNonconstLOWSBlock(const int i, const int j)
Base class for all linear operators.
Base class for all linear operators that can support a high-level solve operation.
Interface for a collection of column vectors called a multi-vector.
Base interface for physically blocked linear operators.
Base interface for product multi-vectors.
virtual Teuchos::RCP< const MultiVectorBase< Scalar > > getMultiVectorBlock(const int k) const =0
Returns a non-persisting const view of the (zero-based) kth block multi-vector.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEUCHOS_ASSERT_INEQUALITY(val1, comp, val2)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
@ SOLVE_STATUS_UNCONVERGED
The requested solution criteria has likely not been achieved.
@ SOLVE_STATUS_CONVERGED
The requested solution criteria has likely been achieved.
#define THYRA_ASSERT_VEC_SPACES(FUNC_NAME, VS1, VS2)
This is a very useful macro that should be used to validate that two vector spaces are compatible.
#define THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(FUNC_NAME, M, M_T, X, Y)
This is a very useful macro that should be used to validate that the spaces for the multi-vector vers...
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
T_To & dyn_cast(T_From &from)
Simple struct that defines the requested solution criteria for a solve.
Simple struct for the return status from a solve.
ESolveStatus solveStatus
The return status of the solve.