Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziTsqrOrthoManager.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Anasazi: Block Eigensolvers Package
5// Copyright 2004 Sandia Corporation
6//
7// Under 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 __AnasaziTsqrOrthoManager_hpp
46#define __AnasaziTsqrOrthoManager_hpp
47
48#include "AnasaziTsqrOrthoManagerImpl.hpp"
50
51
52namespace Anasazi {
53
75 template<class Scalar, class MV>
77 public:
78 typedef Scalar scalar_type;
79 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
82 typedef MV multivector_type;
83
84 typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
85 typedef Teuchos::RCP<mat_type> mat_ptr;
86
95 virtual int
96 normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const = 0;
97
116 virtual int
118 MV& X_out,
119 Teuchos::Array<mat_ptr> C,
120 mat_ptr B,
121 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const = 0;
122
125 };
126
133 template<class Scalar, class MV>
135 public OrthoManager<Scalar, MV>,
136 public OutOfPlaceNormalizerMixin<Scalar, MV>,
137 public Teuchos::ParameterListAcceptor
138 {
139 public:
140 typedef Scalar scalar_type;
141 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
143 typedef MV multivector_type;
144
145 typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
146 typedef Teuchos::RCP<mat_type> mat_ptr;
147
148 void setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& params) {
149 impl_.setParameterList (params);
150 }
151
152 Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList () {
153 return impl_.getNonconstParameterList ();
154 }
155
156 Teuchos::RCP<Teuchos::ParameterList> unsetParameterList () {
157 return impl_.unsetParameterList ();
158 }
159
167 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
168 return impl_.getValidParameters();
169 }
170
180 Teuchos::RCP<const Teuchos::ParameterList> getFastParameters() const {
181 return impl_.getFastParameters();
182 }
183
200 TsqrOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& params,
201 const std::string& label = "Anasazi") :
202 impl_ (params, label)
203 {}
204
209 TsqrOrthoManager (const std::string& label) :
210 impl_ (label)
211 {}
212
214 virtual ~TsqrOrthoManager() {}
215
216 void innerProd (const MV &X, const MV& Y, mat_type& Z) const {
217 return impl_.innerProd (X, Y, Z);
218 }
219
220 void norm (const MV& X, std::vector<magnitude_type>& normVec) const {
221 return impl_.norm (X, normVec);
222 }
223
224 void
225 project (MV &X,
226 Teuchos::Array<Teuchos::RCP<const MV> > Q,
227 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> > > C
228 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,Scalar> >(Teuchos::null))) const
229 {
230 return impl_.project (X, C, Q);
231 }
232
233 int
234 normalize (MV &X, mat_ptr B = Teuchos::null) const
235 {
236 return impl_.normalize (X, B);
237 }
238
239 int
241 Teuchos::Array<Teuchos::RCP<const MV> > Q,
242 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> > > C
243 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,Scalar> >(Teuchos::null)),
244 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> > B = Teuchos::null) const
245 {
246 return impl_.projectAndNormalize (X, C, B, Q);
247 }
248
265 int
266 normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const
267 {
268 return impl_.normalizeOutOfPlace (X, Q, B);
269 }
270
291 int
293 MV& X_out,
294 Teuchos::Array<mat_ptr> C,
295 mat_ptr B,
296 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
297 {
298 return impl_.projectAndNormalizeOutOfPlace (X_in, X_out, C, B, Q);
299 }
300
301 magnitude_type orthonormError (const MV &X) const {
302 return impl_.orthonormError (X);
303 }
304
305 magnitude_type orthogError (const MV &X1, const MV &X2) const {
306 return impl_.orthogError (X1, X2);
307 }
308
309 private:
315 mutable TsqrOrthoManagerImpl<Scalar, MV> impl_;
316 };
317
318
331 template<class Scalar, class MV, class OP>
333 public MatOrthoManager<Scalar, MV, OP>,
334 public OutOfPlaceNormalizerMixin<Scalar, MV>,
335 public Teuchos::ParameterListAcceptorDefaultBase
336 {
337 public:
338 typedef Scalar scalar_type;
339 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
343 typedef OP operator_type;
344
345 typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
346 typedef Teuchos::RCP<mat_type> mat_ptr;
347
348 private:
357 typedef MatOrthoManager<Scalar, MV, OP> base_type;
358
361 typedef TsqrOrthoManagerImpl<Scalar, MV> tsqr_type;
362
365 typedef SVQBOrthoManager<Scalar, MV, OP> svqb_type;
366
369 typedef MultiVecTraits<Scalar, MV> MVT;
370
371 public:
392 TsqrMatOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& params,
393 const std::string& label = "Belos",
394 Teuchos::RCP<const OP> Op = Teuchos::null) :
395 MatOrthoManager<Scalar, MV, OP> (Op),
396 tsqr_ (params, label),
397 pSvqb_ (Teuchos::null) // Initialized lazily
398 {}
399
408 TsqrMatOrthoManager (const std::string& label = "Belos",
409 Teuchos::RCP<const OP> Op = Teuchos::null) :
410 MatOrthoManager<Scalar, MV, OP> (Op),
411 tsqr_ (label),
412 pSvqb_ (Teuchos::null) // Initialized lazily
413 {}
414
417
425 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
426 return tsqr_.getValidParameters ();
427 }
428
438 Teuchos::RCP<const Teuchos::ParameterList> getFastParameters() {
439 return tsqr_.getFastParameters ();
440 }
441
442 void setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& params) {
443 tsqr_.setParameterList (params);
444 }
445
446 virtual void
447 setOp (Teuchos::RCP< const OP > Op)
448 {
449 // We override the base class' setOp() so that the
450 // SVQBOrthoManager instance gets the new op.
451 //
452 // The base_type notation helps C++ resolve the name for a
453 // member function of a templated base class.
454 base_type::setOp (Op); // base class gets a copy of the Op too
455
456 if (! Op.is_null()) {
457 ensureSvqbInit (); // Make sure the SVQB object has been initialized
458 pSvqb_->setOp (Op);
459 }
460 }
461
462 Teuchos::RCP<const OP> getOp () const {
463 // The base_type notation helps C++ resolve the name for a
464 // member function of a templated base class.
465 return base_type::getOp();
466 }
467
468 void
469 projectMat (MV &X,
470 Teuchos::Array<Teuchos::RCP<const MV> > Q,
471 Teuchos::Array<Teuchos::RCP<mat_type> > C =
472 Teuchos::tuple (Teuchos::RCP<mat_type> (Teuchos::null)),
473 Teuchos::RCP<MV> MX = Teuchos::null,
474 Teuchos::Array<Teuchos::RCP<const MV> > MQ =
475 Teuchos::tuple (Teuchos::null)) const
476 {
477 if (getOp().is_null()) {
478 // FIXME (mfh 12 Jan 2011): Do we need to check if C is null
479 // and allocate space if so?
480 tsqr_.project (X, C, Q);
481 // FIXME (mfh 20 Jul 2010) What about MX and MQ?
482 //
483 // FIXME (mfh 12 Jan 2011) Need to port MVT::Assign from Belos to Anasazi
484#if 0
485 if (! MX.is_null()) {
486 // MX gets a copy of X; M is the identity operator.
487 MVT::Assign (X, *MX);
488 }
489#endif // 0
490 } else {
491 ensureSvqbInit ();
492 pSvqb_->projectMat (X, Q, C, MX, MQ);
493 }
494 }
495
496 int
497 normalizeMat (MV &X,
498 mat_ptr B = Teuchos::null,
499 Teuchos::RCP<MV> MX = Teuchos::null) const
500 {
501 if (getOp().is_null()) {
502 // FIXME (mfh 12 Jan 2011): Do we need to check if B is null
503 // and allocate space if so?
504 const int rank = tsqr_.normalize (X, B);
505 // FIXME (mfh 20 Jul 2010) What about MX?
506 //
507 // FIXME (mfh 12 Jan 2011) Need to port MVT::Assign from Belos to Anasazi
508#if 0
509 if (! MX.is_null()) {
510 // MX gets a copy of X; M is the identity operator.
511 MVT::Assign (X, *MX);
512 }
513#endif // 0
514 return rank;
515 } else {
516 ensureSvqbInit ();
517 return pSvqb_->normalizeMat (X, B, MX);
518 }
519 }
520
521 int
523 Teuchos::Array<Teuchos::RCP<const MV> > Q,
524 Teuchos::Array<Teuchos::RCP<mat_type> > C =
525 Teuchos::tuple (Teuchos::RCP<mat_type> (Teuchos::null)),
526 Teuchos::RCP<mat_type> B = Teuchos::null,
527 Teuchos::RCP<MV> MX = Teuchos::null,
528 Teuchos::Array<Teuchos::RCP<const MV> > MQ =
529 Teuchos::tuple (Teuchos::RCP<const MV> (Teuchos::null))) const
530 {
531 if (getOp().is_null()) {
532 // FIXME (mfh 12 Jan 2011): Do we need to check if C and B
533 // are null and allocate space if so?
534 const int rank = tsqr_.projectAndNormalize (X, C, B, Q);
535 // FIXME (mfh 20 Jul 2010) What about MX and MQ?
536 // mfh 12 Jan 2011: Ignore MQ (assume Q == MQ).
537 //
538 // FIXME (mfh 12 Jan 2011) Need to port MVT::Assign from Belos to Anasazi
539#if 0
540 if (! MX.is_null()) {
541 // MX gets a copy of X; M is the identity operator.
542 MVT::Assign (X, *MX);
543 }
544#endif // 0
545 return rank;
546 } else {
547 ensureSvqbInit ();
548 return pSvqb_->projectAndNormalizeMat (X, Q, C, B, MX, MQ);
549 }
550 }
551
552 int
553 normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const
554 {
555 if (getOp().is_null()) {
556 return tsqr_.normalizeOutOfPlace (X, Q, B);
557 } else {
558 // SVQB normalizes in place, so we have to copy.
559 ensureSvqbInit ();
560 const int rank = pSvqb_->normalize (X, B);
561 MVT::Assign (X, Q);
562 return rank;
563 }
564 }
565
566 int
568 MV& X_out,
569 Teuchos::Array<mat_ptr> C,
570 mat_ptr B,
571 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
572 {
573 using Teuchos::null;
574
575 if (getOp().is_null()) {
576 return tsqr_.projectAndNormalizeOutOfPlace (X_in, X_out, C, B, Q);
577 } else {
578 ensureSvqbInit ();
579 // SVQB's projectAndNormalize wants an Array, not an ArrayView.
580 // Copy into a temporary Array and copy back afterwards.
581 Teuchos::Array<Teuchos::RCP<const MV> > Q_array (Q);
582 const int rank = pSvqb_->projectAndNormalize (X_in, Q_array, C, B);
583 Q.assign (Q_array);
584 // SVQB normalizes in place, so we have to copy X_in to X_out.
585 MVT::Assign (X_in, X_out);
586 return rank;
587 }
588 }
589
590 magnitude_type
591 orthonormErrorMat (const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const
592 {
593 if (getOp().is_null()) {
594 return tsqr_.orthonormError (X);
595 // FIXME (mfh 20 Jul 2010) What about MX?
596 } else {
597 ensureSvqbInit ();
598 return pSvqb_->orthonormErrorMat (X, MX);
599 }
600 }
601
602 magnitude_type
603 orthogErrorMat (const MV &X,
604 const MV &Y,
605 Teuchos::RCP<const MV> MX = Teuchos::null,
606 Teuchos::RCP<const MV> MY = Teuchos::null) const
607 {
608 if (getOp().is_null()) {
609 return tsqr_.orthogError (X, Y);
610 // FIXME (mfh 20 Jul 2010) What about MX and MY?
611 } else {
612 ensureSvqbInit ();
613 return pSvqb_->orthogErrorMat (X, Y, MX, MY);
614 }
615 }
616
617 private:
619 void
620 ensureSvqbInit () const
621 {
622 // NOTE (mfh 12 Oct 2011) For now, we rely on the default values
623 // of SVQB parameters.
624 if (pSvqb_.is_null()) {
625 pSvqb_ = Teuchos::rcp (new svqb_type (getOp()));
626 }
627 }
628
636 mutable tsqr_type tsqr_;
637
642 mutable Teuchos::RCP<svqb_type> pSvqb_;
643 };
644
645} // namespace Anasazi
646
647#endif // __AnasaziTsqrOrthoManager_hpp
648
Orthogonalization manager based on the SVQB technique described in "A Block Orthogonalization Procedu...
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
virtual Teuchos::RCP< const OP > getOp() const
Get operator used for inner product.
virtual void setOp(Teuchos::RCP< const OP > Op)
Set operator used for inner product.
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
Mixin for out-of-place orthogonalization.
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 ~OutOfPlaceNormalizerMixin()
Trivial virtual destructor, to silence compiler warnings.
virtual int normalizeOutOfPlace(MV &X, MV &Q, mat_ptr B) const =0
Normalize X into Q*B.
MV multivector_type
Multivector type with which this class was specialized.
MatOrthoManager subclass using TSQR or SVQB.
int projectAndNormalizeMat(MV &X, Teuchos::Array< Teuchos::RCP< const MV > > Q, Teuchos::Array< Teuchos::RCP< mat_type > > C=Teuchos::tuple(Teuchos::RCP< mat_type >(Teuchos::null)), Teuchos::RCP< mat_type > B=Teuchos::null, Teuchos::RCP< MV > MX=Teuchos::null, Teuchos::Array< Teuchos::RCP< const MV > > MQ=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null))) const
Provides matrix-based projection/orthonormalization method.
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters()
Get "fast" parameters for TsqrMatOrthoManager.
OP operator_type
Operator type with which this class was specialized.
virtual void setOp(Teuchos::RCP< const OP > Op)
Set operator used for inner product.
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.
int normalizeOutOfPlace(MV &X, MV &Q, mat_ptr B) const
Normalize X into Q*B.
TsqrMatOrthoManager(const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null)
Constructor (that sets default parameters).
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).
magnitude_type orthonormErrorMat(const MV &X, Teuchos::RCP< const MV > MX=Teuchos::null) const
This method computes the error in orthonormality of a multivector.
magnitude_type orthogErrorMat(const MV &X, const MV &Y, Teuchos::RCP< const MV > MX=Teuchos::null, Teuchos::RCP< const MV > MY=Teuchos::null) const
This method computes the error in orthogonality of two multivectors.
Teuchos::RCP< const OP > getOp() const
Get operator used for inner product.
MV multivector_type
Multivector type with which this class was specialized.
void projectMat(MV &X, Teuchos::Array< Teuchos::RCP< const MV > > Q, Teuchos::Array< Teuchos::RCP< mat_type > > C=Teuchos::tuple(Teuchos::RCP< mat_type >(Teuchos::null)), Teuchos::RCP< MV > MX=Teuchos::null, Teuchos::Array< Teuchos::RCP< const MV > > MQ=Teuchos::tuple(Teuchos::null)) const
Provides matrix-based projection method.
virtual ~TsqrMatOrthoManager()
Destructor (declared virtual for memory safety of derived classes).
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get default parameters for TsqrMatOrthoManager.
TSQR-based OrthoManager subclass.
magnitude_type orthogError(const MV &X1, const MV &X2) const
This method computes the error in orthogonality of two multivectors.
virtual ~TsqrOrthoManager()
Destructor, declared virtual for safe inheritance.
int projectAndNormalize(MV &X, Teuchos::Array< Teuchos::RCP< const MV > > Q, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > > C=Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > >(Teuchos::null)), Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > B=Teuchos::null) const
Given a set of bases Q[i] and a multivector X, this method computes an orthonormal basis for .
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters() const
Get "fast" parameters for TsqrOrthoManager.
magnitude_type orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector.
TsqrOrthoManager(const Teuchos::RCP< Teuchos::ParameterList > &params, const std::string &label="Anasazi")
Constructor (that sets user-specified parameters).
TsqrOrthoManager(const std::string &label)
Constructor (that sets default parameters).
int normalizeOutOfPlace(MV &X, MV &Q, mat_ptr B) const
Normalize X into Q*B, overwriting X with invalid values.
void innerProd(const MV &X, const MV &Y, mat_type &Z) const
Provides the inner product defining the orthogonality concepts.
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< const Teuchos::ParameterList > getValidParameters() const
Default valid parameter list.
void project(MV &X, Teuchos::Array< Teuchos::RCP< const MV > > Q, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > > C=Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > >(Teuchos::null))) const
Given a list of mutually orthogonal and internally orthonormal bases Q, this method projects a multiv...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.