Thyra  Version of the Day
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 
52 namespace Thyra {
53 
54 
55 // public
56 
57 
58 // Constructors/Initializers/Accessors
59 
60 
61 template<class Scalar>
63  : blockFillIsActive_(false), numDiagBlocks_(0)
64 {}
65 
66 
67 template<class Scalar>
69  const RCP<PhysicallyBlockedLinearOpBase<Scalar> > &blocks
70  )
71 {
72  assertAndSetBlockStructure(*blocks);
73  blocks_.initialize(blocks);
74 }
75 
76 
77 template<class Scalar>
79  const RCP<const PhysicallyBlockedLinearOpBase<Scalar> > &blocks
80  )
81 {
82  assertAndSetBlockStructure(*blocks);
83  blocks_.initialize(blocks);
84 }
85 
86 
87 template<class Scalar>
88 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
90 {
91  return blocks_.getNonconstObj();
92 }
93 
94 
95 template<class Scalar>
96 RCP<const PhysicallyBlockedLinearOpBase<Scalar> >
98 {
99  return blocks_.getConstObj();
100 }
101 
102 
103 // Overridden from PhysicallyBlockedLinearOpWithSolveBase
104 
105 
106 template<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 
116 template<class Scalar>
118  const int i, const int j,
119  const Teuchos::RCP<LinearOpWithSolveBase<Scalar> > &block
120  )
121 {
122  setLOWSBlockImpl(i,j,block);
123 }
124 
125 
126 template<class Scalar>
128  const int i, const int j,
129  const Teuchos::RCP<const LinearOpWithSolveBase<Scalar> > &block
130  )
131 {
132  setLOWSBlockImpl(i,j,block);
133 }
134 
135 
136 // Overridden from PhysicallyBlockedLinearOpBase
137 
138 
139 template<class Scalar>
141 {
142  assertBlockFillIsActive(false);
143  TEUCHOS_TEST_FOR_EXCEPT("ToDo: Have not implemented flexible block fill yet!");
144 }
145 
146 
147 template<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 #endif
156  assertBlockFillIsActive(false);
157  numDiagBlocks_ = numRowBlocks;
158  diagonalBlocks_.resize(numDiagBlocks_);
159  productRange_ = null;
160  productDomain_ = null;
161  blockFillIsActive_ = true;
162 }
163 
164 
165 template<class Scalar>
167  const Teuchos::RCP<const ProductVectorSpaceBase<Scalar> > &productRange_in,
168  const Teuchos::RCP<const ProductVectorSpaceBase<Scalar> > &productDomain_in
169  )
170 {
171 #ifdef THYRA_DEBUG
172  TEUCHOS_TEST_FOR_EXCEPT( is_null(productRange_in) );
173  TEUCHOS_TEST_FOR_EXCEPT( is_null(productDomain_in) );
174  TEUCHOS_TEST_FOR_EXCEPT( productRange_in->numBlocks() != productDomain_in->numBlocks() );
175 #endif
176  assertBlockFillIsActive(false);
177  productRange_ = productRange_in;
178  productDomain_ = productDomain_in;
179  numDiagBlocks_ = productRange_in->numBlocks();
180  diagonalBlocks_.resize(numDiagBlocks_);
181  blockFillIsActive_ = true;
182 }
183 
184 
185 template<class Scalar>
187 {
188  return blockFillIsActive_;
189 }
190 
191 
192 template<class Scalar>
194  const int i, const int j
195  ) const
196 {
197  assertBlockFillIsActive(true);
198  assertBlockRowCol(i,j);
199  return false; // ToDo: Change this once we accept off-diagonal blocks
200 }
201 
202 
203 template<class Scalar>
205  const int i, const int j,
206  const Teuchos::RCP<LinearOpBase<Scalar> > &block
207  )
208 {
209  assertBlockFillIsActive(true);
210  TEUCHOS_TEST_FOR_EXCEPT("Error, we don't support off-diagonal LOB objects yet!");
211 }
212 
213 
214 template<class Scalar>
216  const int i, const int j,
217  const Teuchos::RCP<const LinearOpBase<Scalar> > &block
218  )
219 {
220  assertBlockFillIsActive(true);
221  TEUCHOS_TEST_FOR_EXCEPT("Error, we don't support off-diagonal LOB objects yet!");
222 }
223 
224 
225 template<class Scalar>
227 {
228  assertBlockFillIsActive(true);
229  Array<RCP<const VectorSpaceBase<Scalar> > > rangeSpaces;
230  Array<RCP<const VectorSpaceBase<Scalar> > > domainSpaces;
231  for ( int k = 0; k < numDiagBlocks_; ++k ) {
232  const RCP<const LinearOpWithSolveBase<Scalar> > lows_k =
233  diagonalBlocks_[k].getConstObj();
234  TEUCHOS_TEST_FOR_EXCEPTION(is_null(lows_k), std::logic_error,
235  "Error, the block diagonal k="<<k<<" can not be null when ending block fill!"
236  );
237  if (is_null(productRange_)) {
238  rangeSpaces.push_back(lows_k->range());
239  domainSpaces.push_back(lows_k->domain());
240  }
241  }
242  if (is_null(productRange_)) {
243  productRange_ = productVectorSpace<Scalar>(rangeSpaces());
244  productDomain_ = productVectorSpace<Scalar>(domainSpaces());
245  }
246  blockFillIsActive_ = false;
247 }
248 
249 
250 template<class Scalar>
252 {
253  assertBlockFillIsActive(false);
254  productRange_ = Teuchos::null;
255  productDomain_ = Teuchos::null;
256  numDiagBlocks_ = 0;
257  diagonalBlocks_.resize(0);
258 }
259 
260 
261 // Overridden from BlockedLinearOpWithSolveBase
262 
263 
264 template<class Scalar>
265 Teuchos::RCP<LinearOpWithSolveBase<Scalar> >
267  const int i, const int j
268  )
269 {
270  assertBlockFillIsActive(false);
271  assertBlockRowCol(i,j);
272  if (i!=j)
273  return Teuchos::null;
274  return diagonalBlocks_[i].getNonconstObj();
275 }
276 
277 
278 template<class Scalar>
279 Teuchos::RCP<const LinearOpWithSolveBase<Scalar> >
281  const int i, const int j
282  ) const
283 {
284  assertBlockFillIsActive(false);
285  assertBlockRowCol(i, j);
286  if (i != j)
287  return Teuchos::null;
288  return diagonalBlocks_[i].getConstObj();
289 }
290 
291 
292 // Overridden from BlockedLinearOpBase
293 
294 
295 template<class Scalar>
296 Teuchos::RCP<const ProductVectorSpaceBase<Scalar> >
298 {
299  return productRange_;
300 }
301 
302 
303 template<class Scalar>
304 Teuchos::RCP<const ProductVectorSpaceBase<Scalar> >
306 {
307  return productDomain_;
308 }
309 
310 
311 template<class Scalar>
313  const int i, const int j
314  ) const
315 {
316  assertBlockFillIsActive(false);
317  assertBlockRowCol(i,j);
318  if (i!=j)
319  return false; // ToDo: Update this when off-diagonals are supported!
320  return !is_null(diagonalBlocks_[i].getConstObj());
321 }
322 
323 
324 template<class Scalar>
326  const int i, const int j
327  ) const
328 {
329  assertBlockFillIsActive(true);
330  assertBlockRowCol(i,j);
331  return diagonalBlocks_[i].isConst();
332 }
333 
334 
335 template<class Scalar>
336 Teuchos::RCP<LinearOpBase<Scalar> >
338  const int i, const int j
339  )
340 {
341  assertBlockFillIsActive(true);
342  assertBlockRowCol(i,j);
343  if (i!=j)
344  return Teuchos::null; // ToDo: Update when off-diagonals are supported!
345  return this->getNonconstLOWSBlock(i,j);
346 }
347 
348 
349 template<class Scalar>
350 Teuchos::RCP<const LinearOpBase<Scalar> >
352  const int i, const int j
353  ) const
354 {
355  assertBlockFillIsActive(true);
356  assertBlockRowCol(i,j);
357  if (i!=j)
358  return Teuchos::null; // ToDo: Update when off-diagonals are supported!
359  return this->getLOWSBlock(i,j);
360 }
361 
362 
363 // Overridden from LinearOpBase
364 
365 
366 template<class Scalar>
367 Teuchos::RCP< const VectorSpaceBase<Scalar> >
369 {
370  return this->productRange();
371 }
372 
373 
374 template<class Scalar>
375 Teuchos::RCP< const VectorSpaceBase<Scalar> >
377 {
378  return this->productDomain();
379 }
380 
381 
382 template<class Scalar>
383 Teuchos::RCP<const LinearOpBase<Scalar> >
385 {
386  return Teuchos::null; // ToDo: Implement clone when needed!
387 }
388 
389 
390 // Overridden from Teuchos::Describable
391 
392 
393 template<class Scalar>
394 std::string
396 {
397  assertBlockFillIsActive(false);
398  std::ostringstream oss;
399  oss
400  << Teuchos::Describable::description() << "{"
401  << "numDiagBlocks="<<numDiagBlocks_
402  << "}";
403  return oss.str();
404 }
405 
406 
407 template<class Scalar>
409  Teuchos::FancyOStream &out,
410  const Teuchos::EVerbosityLevel verbLevel
411  ) const
412 {
413  assertBlockFillIsActive(false);
414  Teuchos::Describable::describe(out, verbLevel);
415  // ToDo: Fill in a better version of this!
416 }
417 
418 
419 // protected
420 
421 
422 // Overridden from LinearOpBase
423 
424 
425 template<class Scalar>
427  EOpTransp M_trans
428  ) const
429 {
430  using Thyra::opSupported;
431  assertBlockFillIsActive(false);
432  for ( int k = 0; k < numDiagBlocks_; ++k ) {
433  if ( !opSupported(*diagonalBlocks_[k].getConstObj(),M_trans) )
434  return false;
435  }
436  return true;
437  // ToDo: To be safe we really should do a collective reduction of all
438  // clusters of processes. However, for the typical use case, every block
439  // will return the same info and we should be safe!
440 }
441 
442 
443 template<class Scalar>
445  const EOpTransp M_trans,
446  const MultiVectorBase<Scalar> &X_in,
447  const Ptr<MultiVectorBase<Scalar> > &Y_inout,
448  const Scalar alpha,
449  const Scalar beta
450  ) const
451 {
452 
453  using Teuchos::RCP;
454  using Teuchos::dyn_cast;
455  using Thyra::apply;
456 
457 #ifdef THYRA_DEBUG
459  "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::apply(...)",
460  *this, M_trans, X_in, &*Y_inout
461  );
462 #endif // THYRA_DEBUG
463 
464  //
465  // Y = alpha * op(M) * X + beta*Y
466  //
467  // =>
468  //
469  // Y[i] = beta+Y[i] + alpha*op(Op)[i]*X[i], for i=0...numDiagBlocks-1
470  //
471  // ToDo: Update to handle off diagonal blocks when needed!
472  //
473 
475  &X = dyn_cast<const ProductMultiVectorBase<Scalar> >(X_in);
477  &Y = dyn_cast<ProductMultiVectorBase<Scalar> >(*Y_inout);
478 
479  for ( int i = 0; i < numDiagBlocks_; ++ i ) {
480  Thyra::apply( *diagonalBlocks_[i].getConstObj(), M_trans,
481  *X.getMultiVectorBlock(i), Y.getNonconstMultiVectorBlock(i).ptr(),
482  alpha, beta
483  );
484  }
485 
486 }
487 
488 
489 // Overridden from LinearOpWithSolveBase
490 
491 
492 template<class Scalar>
493 bool
495  EOpTransp M_trans
496  ) const
497 {
498  assertBlockFillIsActive(false);
499  for ( int k = 0; k < numDiagBlocks_; ++k ) {
500  if ( !Thyra::solveSupports( *diagonalBlocks_[k].getConstObj(), M_trans ) )
501  return false;
502  }
503  return true;
504  // ToDo: To be safe we really should do a collective reduction of all
505  // clusters of processes. However, for the typical use case, every block
506  // will return the same info and we should be safe!
507 }
508 
509 
510 template<class Scalar>
511 bool
513  EOpTransp M_trans, const SolveMeasureType& solveMeasureType
514  ) const
515 {
516  using Thyra::solveSupportsSolveMeasureType;
517  assertBlockFillIsActive(false);
518  for ( int k = 0; k < numDiagBlocks_; ++k ) {
519  if (
520  !solveSupportsSolveMeasureType(
521  *diagonalBlocks_[k].getConstObj(),
522  M_trans, solveMeasureType
523  )
524  )
525  {
526  return false;
527  }
528  }
529  return true;
530 }
531 
532 
533 template<class Scalar>
536  const EOpTransp M_trans,
537  const MultiVectorBase<Scalar> &B_in,
538  const Ptr<MultiVectorBase<Scalar> > &X_inout,
539  const Ptr<const SolveCriteria<Scalar> > solveCriteria
540  ) const
541 {
542 
543  using Teuchos::RCP;
544  using Teuchos::dyn_cast;
545  using Thyra::solve;
546 
547 #ifdef THYRA_DEBUG
549  "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::apply(...)",
550  *this, M_trans, *X_inout, &B_in
551  );
552  TEUCHOS_TEST_FOR_EXCEPT(!this->solveSupportsImpl(M_trans));
553  TEUCHOS_TEST_FOR_EXCEPTION(
554  nonnull(solveCriteria) && !solveCriteria->solveMeasureType.useDefault(),
555  std::logic_error,
556  "Error, we can't handle any non-default solve criteria yet!"
557  );
558  // ToDo: If solve criteria is to be handled, then we will have to be very
559  // carefull how it it interpreted in terms of the individual period solves!
560 #endif // THYRA_DEBUG
561 
562  //
563  // Y = alpha * inv(op(M)) * X + beta*Y
564  //
565  // =>
566  //
567  // X[i] = inv(op(Op[i]))*B[i], for i=0...numDiagBlocks-1
568  //
569  // ToDo: Update to handle off diagonal blocks when needed!
570  //
571 
573  &B = dyn_cast<const ProductMultiVectorBase<Scalar> >(B_in);
575  &X = dyn_cast<ProductMultiVectorBase<Scalar> >(*X_inout);
576 
577  for ( int i = 0; i < numDiagBlocks_; ++ i ) {
578  const RCP<const LinearOpWithSolveBase<Scalar> >
579  Op_k = diagonalBlocks_[i].getConstObj();
580  Op_k->setOStream(this->getOStream());
581  Op_k->setVerbLevel(this->getVerbLevel());
582  Thyra::solve( *Op_k, M_trans, *B.getMultiVectorBlock(i),
583  X.getNonconstMultiVectorBlock(i).ptr() );
584  // ToDo: Pass in solve criteria when needed!
585  }
586 
587  return SolveStatus<Scalar>();
588 
589 }
590 
591 
592 
593 // private
594 
595 
596 template<class Scalar>
598  bool blockFillIsActive_in
599  ) const
600 {
601 #ifdef THYRA_DEBUG
602  TEUCHOS_TEST_FOR_EXCEPT(!(blockFillIsActive_==blockFillIsActive_in));
603 #endif
604 }
605 
606 
607 template<class Scalar>
608 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertBlockRowCol(
609  const int i, const int j
610  ) const
611 {
612 #ifdef THYRA_DEBUG
613  TEUCHOS_TEST_FOR_EXCEPTION(
614  !( 0 <= i && i < numDiagBlocks_ ), std::logic_error,
615  "Error, i="<<i<<" does not fall in the range [0,"<<numDiagBlocks_-1<<"]!"
616  );
617  TEUCHOS_TEST_FOR_EXCEPTION(
618  !( 0 <= j && j < numDiagBlocks_ ), std::logic_error,
619  "Error, j="<<j<<" does not fall in the range [0,"<<numDiagBlocks_-1<<"]!"
620  );
621 #endif
622 }
623 
624 
625 template<class Scalar>
626 template<class LinearOpWithSolveType>
627 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::setLOWSBlockImpl(
628  const int i, const int j,
629  const Teuchos::RCP<LinearOpWithSolveType> &block
630  )
631 {
632  assertBlockFillIsActive(true);
633 #ifdef THYRA_DEBUG
634  TEUCHOS_ASSERT_INEQUALITY( i, >=, 0 );
635  TEUCHOS_ASSERT_INEQUALITY( j, >=, 0 );
636  TEUCHOS_ASSERT_INEQUALITY( i, <, numDiagBlocks_ );
637  TEUCHOS_ASSERT_INEQUALITY( j, <, numDiagBlocks_ );
638  TEUCHOS_TEST_FOR_EXCEPTION(
639  !this->acceptsLOWSBlock(i,j), std::logic_error,
640  "Error, this DefaultBlockedTriangularLinearOpWithSolve does not accept\n"
641  "LOWSB objects for the block i="<<i<<", j="<<j<<"!"
642  );
643 #endif
644  diagonalBlocks_[i] = block;
645 }
646 
647 
648 template<class Scalar>
649 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertAndSetBlockStructure(
650  const PhysicallyBlockedLinearOpBase<Scalar>& blocks
651  )
652 {
653 #ifdef THYRA_DEBUG
655  "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertAndSetBlockStructure(blocks)",
656  *blocks.range(), *this->range()
657  );
659  "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertAndSetBlockStructure(blocks)",
660  *blocks.domain(), *this->domain()
661  );
662  // ToDo: Make sure that all of the blocks are above or below the diagonal
663  // but not both!
664 #endif
665  // ToDo: Set if this is an upper or lower triangular block operator.
666 }
667 
668 
669 } // namespace Thyra
670 
671 
672 #endif // THYRA_DEFAULT_BLOCKED_TRIANGULAR_LINEAR_OP_WITH_SOLVE_DEF_HPP
Base interface for product multi-vectors.
RCP< LinearOpWithSolveBase< Scalar > > getNonconstLOWSBlock(const int i, const int j)
#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...
void setBlocks(const RCP< const PhysicallyBlockedLinearOpBase< Scalar > > &blocks)
Base class for all linear operators that can support a high-level solve operation.
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
void setNonconstBlock(const int i, const int j, const RCP< LinearOpBase< Scalar > > &block)
Concrete composite LinearOpWithSolveBase subclass that creates single upper or lower block triangular...
#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...
bool solveSupportsSolveMeasureTypeImpl(EOpTransp M_trans, const SolveMeasureType &solveMeasureType) const
void setBlock(const int i, const int j, const RCP< const LinearOpBase< Scalar > > &block)
SolveStatus< Scalar > solveImpl(const EOpTransp transp, const MultiVectorBase< Scalar > &B, const Ptr< MultiVectorBase< Scalar > > &X, const Ptr< const SolveCriteria< Scalar > > solveCriteria) const
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
RCP< const LinearOpWithSolveBase< Scalar > > getLOWSBlock(const int i, const int j) const
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.
RCP< const LinearOpBase< Scalar > > getBlock(const int i, const int j) const
RCP< LinearOpBase< Scalar > > getNonconstBlock(const int i, const int j)
Interface for a collection of column vectors called a multi-vector.
void setLOWSBlock(const int i, const int j, const RCP< const LinearOpWithSolveBase< Scalar > > &block)
Base interface for physically blocked linear operators.
Simple struct for the return status from a solve.
Base class for all linear operators.
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.
Simple struct that defines the requested solution criteria for a solve.
void setNonconstBlocks(const RCP< PhysicallyBlockedLinearOpBase< Scalar > > &blocks)