Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
BelosTsqrOrthoManager.hpp
Go to the documentation of this file.
1//@HEADER
2// ************************************************************************
3//
4// Belos: Block Linear Solvers Package
5// Copyright 2004 Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
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 Michael A. Heroux (maherou@sandia.gov)
38//
39// ************************************************************************
40//@HEADER
41
44
45#ifndef __BelosTsqrOrthoManager_hpp
46#define __BelosTsqrOrthoManager_hpp
47
49// Belos doesn't have an SVQB orthomanager implemented yet, so we just
50// use a default orthomanager for the case where the inner product
51// matrix is nontrivial.
53
54
55namespace Belos {
56
78template<class Scalar, class MV>
80public:
81 typedef Scalar scalar_type;
82 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
85 typedef MV multivector_type;
86
87 typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
88 typedef Teuchos::RCP<mat_type> mat_ptr;
89
98 virtual int
99 normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const = 0;
100
119 virtual int
121 MV& X_out,
122 Teuchos::Array<mat_ptr> C,
123 mat_ptr B,
124 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const = 0;
125
128};
129
136template<class Scalar, class MV>
138 public OrthoManager<Scalar, MV>,
139 public OutOfPlaceNormalizerMixin<Scalar, MV>
140{
141public:
142 typedef Scalar scalar_type;
143 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
146
147 typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
148 typedef Teuchos::RCP<mat_type> mat_ptr;
149
150 void setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& params) {
151 impl_.setParameterList (params);
152 }
153
154 Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList () {
155 return impl_.getNonconstParameterList ();
156 }
157
158 Teuchos::RCP<Teuchos::ParameterList> unsetParameterList () {
159 return impl_.unsetParameterList ();
160 }
161
169 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
170 return impl_.getValidParameters();
171 }
172
182 Teuchos::RCP<const Teuchos::ParameterList> getFastParameters() const {
183 return impl_.getFastParameters();
184 }
185
201 TsqrOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& params,
202 const std::string& label = "Belos") :
203 impl_ (params, label)
204 {}
205
210 TsqrOrthoManager (const std::string& label) :
211 impl_ (label)
212 {}
213
215 virtual ~TsqrOrthoManager() {}
216
238 void
240 {
241 impl_.setReorthogonalizationCallback (callback);
242 }
243
244 void innerProd (const MV &X, const MV &Y, mat_type& Z) const {
245 return impl_.innerProd (X, Y, Z);
246 }
247
248 void norm (const MV& X, std::vector<magnitude_type>& normVec) const {
249 return impl_.norm (X, normVec);
250 }
251
252 void
253 project (MV &X,
254 Teuchos::Array<mat_ptr> C,
255 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
256 {
257 return impl_.project (X, C, Q);
258 }
259
260 int
261 normalize (MV &X, mat_ptr B) const
262 {
263 return impl_.normalize (X, B);
264 }
265
266protected:
267 virtual int
269 Teuchos::Array<mat_ptr> C,
270 mat_ptr B,
271 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
272 {
273 return impl_.projectAndNormalize (X, C, B, Q);
274 }
275
276public:
293 int
294 normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const
295 {
296 return impl_.normalizeOutOfPlace (X, Q, B);
297 }
298
319 int
321 MV& X_out,
322 Teuchos::Array<mat_ptr> C,
323 mat_ptr B,
324 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
325 {
326 return impl_.projectAndNormalizeOutOfPlace (X_in, X_out, C, B, Q);
327 }
328
329 magnitude_type orthonormError (const MV &X) const {
330 return impl_.orthonormError (X);
331 }
332
333 magnitude_type orthogError (const MV &X1, const MV &X2) const {
334 return impl_.orthogError (X1, X2);
335 }
336
344 void setLabel (const std::string& label) {
345 impl_.setLabel (label);
346 }
347
348 const std::string& getLabel() const { return impl_.getLabel(); }
349
350private:
357
359 std::string label_;
360};
361
362
376template<class Scalar, class MV, class OP>
378 public MatOrthoManager<Scalar, MV, OP>,
379 public OutOfPlaceNormalizerMixin<Scalar, MV>
380{
381public:
382 typedef Scalar scalar_type;
383 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
387 typedef OP operator_type;
388
389 typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
390 typedef Teuchos::RCP<mat_type> mat_ptr;
391
392private:
402
406
410
414
415public:
436 TsqrMatOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& params,
437 const std::string& label = "Belos",
438 Teuchos::RCP<const OP> Op = Teuchos::null) :
439 MatOrthoManager<Scalar, MV, OP> (Op),
440 tsqr_ (params, label),
441 pDgks_ (Teuchos::null) // Initialized lazily
442 {}
443
452 TsqrMatOrthoManager (const std::string& label = "Belos",
453 Teuchos::RCP<const OP> Op = Teuchos::null) :
454 MatOrthoManager<Scalar, MV, OP> (Op),
455 tsqr_ (label),
456 pDgks_ (Teuchos::null) // Initialized lazily
457 {}
458
461
469 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
470 return tsqr_.getValidParameters ();
471 }
472
482 Teuchos::RCP<const Teuchos::ParameterList> getFastParameters() {
483 return tsqr_.getFastParameters ();
484 }
485
486 void setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& params) {
487 tsqr_.setParameterList (params);
488 }
489
490 const std::string& getLabel() const { return tsqr_.getLabel (); }
491
492 void
493 setOp (Teuchos::RCP<const OP> Op)
494 {
495 // We override the base class' setOp() so that the
496 // DGKSOrthoManager instance gets the new op.
497 //
498 // The base_type notation helps C++ resolve the name for a
499 // member function of a templated base class.
500 base_type::setOp (Op); // base class gets a copy of the Op too
501
502 if (! Op.is_null()) {
503 ensureDgksInit (); // Make sure the DGKS object has been initialized
504 pDgks_->setOp (Op);
505 }
506 }
507
508 Teuchos::RCP<const OP> getOp () const {
509 // The base_type notation helps C++ resolve the name for a
510 // member function of a templated base class.
511 return base_type::getOp();
512 }
513
514 void
515 project (MV &X,
516 Teuchos::RCP<MV> MX,
517 Teuchos::Array<mat_ptr> C,
518 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
519 {
520 if (getOp().is_null()) {
521 tsqr_.project (X, C, Q);
522 if (! MX.is_null()) {
523 // MX gets a copy of X; M is the identity operator.
524 MVT::Assign (X, *MX);
525 }
526 } else {
528 pDgks_->project (X, MX, C, Q);
529 }
530 }
531
532 void
533 project (MV &X,
534 Teuchos::Array<mat_ptr> C,
535 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
536 {
537 project (X, Teuchos::null, C, Q);
538 }
539
540 int
541 normalize (MV& X, Teuchos::RCP<MV> MX, mat_ptr B) const
542 {
543 if (getOp().is_null()) {
544 const int rank = tsqr_.normalize (X, B);
545 if (! MX.is_null()) {
546 // MX gets a copy of X; M is the identity operator.
547 MVT::Assign (X, *MX);
548 }
549 return rank;
550 } else {
552 return pDgks_->normalize (X, MX, B);
553 }
554 }
555
556 int normalize (MV& X, mat_ptr B) const {
557 return normalize (X, Teuchos::null, B);
558 }
559
560 // Attempted fix for a warning about hiding
561 // OrthoManager::projectAndNormalize(). Unfortunately, this fix turns out
562 // to produce a compilation error with cray++, see bug #6129,
563 // <https://software.sandia.gov/bugzilla/show_bug.cgi?id=6129>.
564 //using Belos::OrthoManager<Scalar, MV>::projectAndNormalize;
565
566protected:
567 virtual int
569 Teuchos::RCP<MV> MX,
570 Teuchos::Array<mat_ptr> C,
571 mat_ptr B,
572 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
573 {
574 if (getOp().is_null()) {
575 const int rank = tsqr_.projectAndNormalize (X, C, B, Q);
576 if (! MX.is_null()) {
577 // MX gets a copy of X; M is the identity operator.
578 MVT::Assign (X, *MX);
579 }
580 return rank;
581 } else {
583 return pDgks_->projectAndNormalize (X, MX, C, B, Q);
584 }
585 }
586
587public:
588 int
589 normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const
590 {
591 if (getOp().is_null()) {
592 return tsqr_.normalizeOutOfPlace (X, Q, B);
593 } else {
594 // DGKS normalizes in place, so we have to copy.
596 const int rank = pDgks_->normalize (X, B);
597 MVT::Assign (X, Q);
598 return rank;
599 }
600 }
601
602 int
604 MV& X_out,
605 Teuchos::Array<mat_ptr> C,
606 mat_ptr B,
607 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
608 {
609 using Teuchos::null;
610
611 if (getOp().is_null()) {
612 return tsqr_.projectAndNormalizeOutOfPlace (X_in, X_out, C, B, Q);
613 } else {
614 // DGKS normalizes in place, so we have to copy.
616 const int rank = pDgks_->projectAndNormalize (X_in, null, C, B, Q);
617 MVT::Assign (X_in, X_out);
618 return rank;
619 }
620 }
621
623 orthonormError (const MV &X, Teuchos::RCP<const MV> MX) const
624 {
625 if (getOp().is_null()) {
626 return tsqr_.orthonormError (X); // Ignore MX
627 } else {
629 return pDgks_->orthonormError (X, MX);
630 }
631 }
632
633 magnitude_type orthonormError (const MV &X) const {
634 return orthonormError (X, Teuchos::null);
635 }
636
637 magnitude_type orthogError (const MV &X1, const MV &X2) const {
638 return orthogError (X1, Teuchos::null, X2);
639 }
640
642 orthogError (const MV &X1,
643 Teuchos::RCP<const MV> MX1,
644 const MV &X2) const
645 {
646 if (getOp ().is_null ()) {
647 // Ignore MX1, since we don't need to write to it.
648 return tsqr_.orthogError (X1, X2);
649 } else {
651 return pDgks_->orthogError (X1, MX1, X2);
652 }
653 }
654
655 void
656 setLabel (const std::string& label)
657 {
658 tsqr_.setLabel (label);
659
660 // Make sure DGKS gets the new label, if it's initialized.
661 // Otherwise, it will get the new label on initialization.
662 if (! pDgks_.is_null ()) {
663 pDgks_->setLabel (label);
664 }
665 }
666
667private:
669 void
671 {
672 // NOTE (mfh 11 Jan 2011) DGKS has a parameter that needs to be
673 // set. For now, we just use the default value of the parameter.
674 if (pDgks_.is_null ()) {
675 pDgks_ = Teuchos::rcp (new dgks_type (getLabel (), getOp ()));
676 }
677 }
678
687
693 mutable Teuchos::RCP<dgks_type> pDgks_;
694};
695
696} // namespace Belos
697
698#endif // __BelosTsqrOrthoManager_hpp
Classical Gram-Schmidt (with DGKS correction) implementation of the Belos::OrthoManager class.
Orthogonalization manager back end based on Tall Skinny QR (TSQR)
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
void setOp(Teuchos::RCP< const OP > Op)
Set operator.
Teuchos::RCP< const OP > getOp() const
Get operator.
Traits class which defines basic operations on multivectors.
static void Assign(const MV &A, MV &mv)
mv := A
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Mixin for out-of-place orthogonalization.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
MV multivector_type
Multivector type with which this class was specialized.
virtual int projectAndNormalizeOutOfPlace(MV &X_in, MV &X_out, Teuchos::Array< mat_ptr > C, mat_ptr B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const =0
Project and normalize X_in into X_out.
virtual int normalizeOutOfPlace(MV &X, MV &Q, mat_ptr B) const =0
Normalize X into Q*B.
Teuchos::SerialDenseMatrix< int, Scalar > mat_type
virtual ~OutOfPlaceNormalizerMixin()
Trivial virtual destructor, to silence compiler warnings.
Interface of callback invoked by TsqrOrthoManager on reorthogonalization.
MatOrthoManager subclass using TSQR or DGKS.
void ensureDgksInit() const
Ensure that the DGKSOrthoManager instance is initialized.
int normalize(MV &X, Teuchos::RCP< MV > MX, mat_ptr B) const
virtual ~TsqrMatOrthoManager()
Destructor (declared virtual for memory safety of derived classes).
MultiVecTraits< Scalar, MV > MVT
Traits class for the multivector type.
TsqrOrthoManagerImpl< Scalar, MV > tsqr_type
Implementation of TSQR-based orthogonalization.
magnitude_type orthogError(const MV &X1, Teuchos::RCP< const MV > MX1, const MV &X2) const
This method computes the error in orthogonality of two multivectors. The method has the option of exp...
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters()
Get "fast" parameters for TsqrMatOrthoManager.
int normalize(MV &X, mat_ptr B) const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &params)
void project(MV &X, Teuchos::Array< mat_ptr > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Teuchos::SerialDenseMatrix< int, Scalar > mat_type
TsqrMatOrthoManager(const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null)
Constructor (that sets default parameters).
magnitude_type orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector.
int normalizeOutOfPlace(MV &X, MV &Q, mat_ptr B) const
Normalize X into Q*B.
tsqr_type tsqr_
TSQR + BGS orthogonalization manager implementation.
void setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
Teuchos::RCP< const OP > getOp() const
MV multivector_type
Multivector type with which this class was specialized.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get default parameters for TsqrMatOrthoManager.
MatOrthoManager< Scalar, MV, OP > base_type
This will be used to help C++ resolve getOp().
magnitude_type orthogError(const MV &X1, const MV &X2) const
This method computes the error in orthogonality of two multivectors. This method.
DGKSOrthoManager< Scalar, MV, OP > dgks_type
Implementation of DGKS-based orthogonalization.
void setOp(Teuchos::RCP< const OP > Op)
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
magnitude_type orthonormError(const MV &X, Teuchos::RCP< const MV > MX) const
This method computes the error in orthonormality of a multivector. The method has the option of explo...
Teuchos::RCP< dgks_type > pDgks_
DGKS orthogonalization manager.
TsqrMatOrthoManager(const Teuchos::RCP< Teuchos::ParameterList > &params, const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null)
Constructor (that sets user-specified parameters).
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
void project(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< mat_ptr > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
int projectAndNormalizeOutOfPlace(MV &X_in, MV &X_out, Teuchos::Array< mat_ptr > C, mat_ptr B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Project and normalize X_in into X_out.
Teuchos::RCP< mat_type > mat_ptr
OP operator_type
Operator type with which this class was specialized.
virtual int projectAndNormalizeWithMxImpl(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< mat_ptr > C, mat_ptr B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
TSQR-based OrthoManager subclass implementation.
int projectAndNormalizeOutOfPlace(MV &X_in, MV &X_out, Teuchos::Array< mat_ptr > C, mat_ptr B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q)
Project and normalize X_in into X_out; overwrite X_in.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &params)
Set parameters from the given parameter list.
int normalizeOutOfPlace(MV &X, MV &Q, mat_ptr B)
Normalize X into Q*B, overwriting X.
void setLabel(const std::string &label)
Set the label for timers.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Default valid parameter list.
int projectAndNormalize(MV &X, Teuchos::Array< mat_ptr > C, mat_ptr B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q)
Project X against Q and normalize X.
magnitude_type orthogError(const MV &X1, const MV &X2) const
Return the Frobenius norm of the inner product of X1 with itself.
void project(MV &X, Teuchos::Array< mat_ptr > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q)
Compute and .
int normalize(MV &X, mat_ptr B)
Orthogonalize the columns of X in place.
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters()
Get "fast" parameters for TsqrOrthoManagerImpl.
const std::string & getLabel() const
Get the label for timers (if timers are enabled).
magnitude_type orthonormError(const MV &X) const
Return .
TSQR-based OrthoManager subclass.
TsqrOrthoManager(const std::string &label)
Constructor (that sets default parameters).
magnitude_type orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector.
void project(MV &X, Teuchos::Array< mat_ptr > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
void norm(const MV &X, std::vector< magnitude_type > &normVec) const
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters() const
Get "fast" parameters for TsqrOrthoManager.
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
void innerProd(const MV &X, const MV &Y, mat_type &Z) const
Provides the inner product defining the orthogonality concepts.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Default valid parameter list.
virtual int projectAndNormalizeImpl(MV &X, Teuchos::Array< mat_ptr > C, mat_ptr B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
magnitude_type orthogError(const MV &X1, const MV &X2) const
This method computes the error in orthogonality of two multivectors.
int normalize(MV &X, mat_ptr B) const
int normalizeOutOfPlace(MV &X, MV &Q, mat_ptr B) const
Normalize X into Q*B, overwriting X with invalid values.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &params)
TsqrOrthoManager(const Teuchos::RCP< Teuchos::ParameterList > &params, const std::string &label="Belos")
Constructor (that sets user-specified parameters).
TsqrOrthoManagerImpl< Scalar, MV > impl_
The implementation of TSQR.
std::string label_
Label for timers (if timers are enabled at compile time).
Teuchos::SerialDenseMatrix< int, Scalar > mat_type
virtual ~TsqrOrthoManager()
Destructor, declared virtual for safe inheritance.
int projectAndNormalizeOutOfPlace(MV &X_in, MV &X_out, Teuchos::Array< mat_ptr > C, mat_ptr B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Project and normalize X_in into X_out; overwrite X_in.
Teuchos::RCP< mat_type > mat_ptr
void setReorthogonalizationCallback(const Teuchos::RCP< ReorthogonalizationCallback< Scalar > > &callback)
Set callback to be invoked on reorthogonalization.
void setLabel(const std::string &label)
Set the label for (the timers for) this orthogonalization manager, and create new timers if the label...