Anasazi  Version of the Day
AnasaziTsqrOrthoManager.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright (2010) 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 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
31 
32 #ifndef __AnasaziTsqrOrthoManager_hpp
33 #define __AnasaziTsqrOrthoManager_hpp
34 
35 #include "AnasaziTsqrOrthoManagerImpl.hpp"
37 
38 
39 namespace Anasazi {
40 
62  template<class Scalar, class MV>
64  public:
65  typedef Scalar scalar_type;
66  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
69  typedef MV multivector_type;
70 
71  typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
72  typedef Teuchos::RCP<mat_type> mat_ptr;
73 
82  virtual int
83  normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const = 0;
84 
103  virtual int
105  MV& X_out,
106  Teuchos::Array<mat_ptr> C,
107  mat_ptr B,
108  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const = 0;
109 
112  };
113 
120  template<class Scalar, class MV>
122  public OrthoManager<Scalar, MV>,
123  public OutOfPlaceNormalizerMixin<Scalar, MV>,
124  public Teuchos::ParameterListAcceptor
125  {
126  public:
127  typedef Scalar scalar_type;
128  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
130  typedef MV multivector_type;
131 
132  typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
133  typedef Teuchos::RCP<mat_type> mat_ptr;
134 
135  void setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& params) {
136  impl_.setParameterList (params);
137  }
138 
139  Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList () {
140  return impl_.getNonconstParameterList ();
141  }
142 
143  Teuchos::RCP<Teuchos::ParameterList> unsetParameterList () {
144  return impl_.unsetParameterList ();
145  }
146 
154  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
155  return impl_.getValidParameters();
156  }
157 
167  Teuchos::RCP<const Teuchos::ParameterList> getFastParameters() const {
168  return impl_.getFastParameters();
169  }
170 
187  TsqrOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& params,
188  const std::string& label = "Anasazi") :
189  impl_ (params, label)
190  {}
191 
196  TsqrOrthoManager (const std::string& label) :
197  impl_ (label)
198  {}
199 
201  virtual ~TsqrOrthoManager() {}
202 
203  void innerProd (const MV &X, const MV& Y, mat_type& Z) const {
204  return impl_.innerProd (X, Y, Z);
205  }
206 
207  void norm (const MV& X, std::vector<magnitude_type>& normVec) const {
208  return impl_.norm (X, normVec);
209  }
210 
211  void
212  project (MV &X,
213  Teuchos::Array<Teuchos::RCP<const MV> > Q,
214  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> > > C
215  = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,Scalar> >(Teuchos::null))) const
216  {
217  return impl_.project (X, C, Q);
218  }
219 
220  int
221  normalize (MV &X, mat_ptr B = Teuchos::null) const
222  {
223  return impl_.normalize (X, B);
224  }
225 
226  int
228  Teuchos::Array<Teuchos::RCP<const MV> > Q,
229  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> > > C
230  = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,Scalar> >(Teuchos::null)),
231  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> > B = Teuchos::null) const
232  {
233  return impl_.projectAndNormalize (X, C, B, Q);
234  }
235 
252  int
253  normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const
254  {
255  return impl_.normalizeOutOfPlace (X, Q, B);
256  }
257 
278  int
280  MV& X_out,
281  Teuchos::Array<mat_ptr> C,
282  mat_ptr B,
283  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
284  {
285  return impl_.projectAndNormalizeOutOfPlace (X_in, X_out, C, B, Q);
286  }
287 
288  magnitude_type orthonormError (const MV &X) const {
289  return impl_.orthonormError (X);
290  }
291 
292  magnitude_type orthogError (const MV &X1, const MV &X2) const {
293  return impl_.orthogError (X1, X2);
294  }
295 
296  private:
302  mutable TsqrOrthoManagerImpl<Scalar, MV> impl_;
303  };
304 
305 
318  template<class Scalar, class MV, class OP>
320  public MatOrthoManager<Scalar, MV, OP>,
321  public OutOfPlaceNormalizerMixin<Scalar, MV>,
322  public Teuchos::ParameterListAcceptorDefaultBase
323  {
324  public:
325  typedef Scalar scalar_type;
326  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
328  typedef MV multivector_type;
330  typedef OP operator_type;
331 
332  typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
333  typedef Teuchos::RCP<mat_type> mat_ptr;
334 
335  private:
345 
349 
353 
357 
358  public:
379  TsqrMatOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& params,
380  const std::string& label = "Belos",
381  Teuchos::RCP<const OP> Op = Teuchos::null) :
382  MatOrthoManager<Scalar, MV, OP> (Op),
383  tsqr_ (params, label),
384  pSvqb_ (Teuchos::null) // Initialized lazily
385  {}
386 
395  TsqrMatOrthoManager (const std::string& label = "Belos",
396  Teuchos::RCP<const OP> Op = Teuchos::null) :
397  MatOrthoManager<Scalar, MV, OP> (Op),
398  tsqr_ (label),
399  pSvqb_ (Teuchos::null) // Initialized lazily
400  {}
401 
403  virtual ~TsqrMatOrthoManager() {}
404 
412  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
413  return tsqr_.getValidParameters ();
414  }
415 
425  Teuchos::RCP<const Teuchos::ParameterList> getFastParameters() {
426  return tsqr_.getFastParameters ();
427  }
428 
429  void setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& params) {
430  tsqr_.setParameterList (params);
431  }
432 
433  virtual void
434  setOp (Teuchos::RCP< const OP > Op)
435  {
436  // We override the base class' setOp() so that the
437  // SVQBOrthoManager instance gets the new op.
438  //
439  // The base_type notation helps C++ resolve the name for a
440  // member function of a templated base class.
441  base_type::setOp (Op); // base class gets a copy of the Op too
442 
443  if (! Op.is_null()) {
444  ensureSvqbInit (); // Make sure the SVQB object has been initialized
445  pSvqb_->setOp (Op);
446  }
447  }
448 
449  Teuchos::RCP<const OP> getOp () const {
450  // The base_type notation helps C++ resolve the name for a
451  // member function of a templated base class.
452  return base_type::getOp();
453  }
454 
455  void
456  projectMat (MV &X,
457  Teuchos::Array<Teuchos::RCP<const MV> > Q,
458  Teuchos::Array<Teuchos::RCP<mat_type> > C =
459  Teuchos::tuple (Teuchos::RCP<mat_type> (Teuchos::null)),
460  Teuchos::RCP<MV> MX = Teuchos::null,
461  Teuchos::Array<Teuchos::RCP<const MV> > MQ =
462  Teuchos::tuple (Teuchos::null)) const
463  {
464  if (getOp().is_null()) {
465  // FIXME (mfh 12 Jan 2011): Do we need to check if C is null
466  // and allocate space if so?
467  tsqr_.project (X, C, Q);
468  // FIXME (mfh 20 Jul 2010) What about MX and MQ?
469  //
470  // FIXME (mfh 12 Jan 2011) Need to port MVT::Assign from Belos to Anasazi
471 #if 0
472  if (! MX.is_null()) {
473  // MX gets a copy of X; M is the identity operator.
474  MVT::Assign (X, *MX);
475  }
476 #endif // 0
477  } else {
478  ensureSvqbInit ();
479  pSvqb_->projectMat (X, Q, C, MX, MQ);
480  }
481  }
482 
483  int
484  normalizeMat (MV &X,
485  mat_ptr B = Teuchos::null,
486  Teuchos::RCP<MV> MX = Teuchos::null) const
487  {
488  if (getOp().is_null()) {
489  // FIXME (mfh 12 Jan 2011): Do we need to check if B is null
490  // and allocate space if so?
491  const int rank = tsqr_.normalize (X, B);
492  // FIXME (mfh 20 Jul 2010) What about MX?
493  //
494  // FIXME (mfh 12 Jan 2011) Need to port MVT::Assign from Belos to Anasazi
495 #if 0
496  if (! MX.is_null()) {
497  // MX gets a copy of X; M is the identity operator.
498  MVT::Assign (X, *MX);
499  }
500 #endif // 0
501  return rank;
502  } else {
503  ensureSvqbInit ();
504  return pSvqb_->normalizeMat (X, B, MX);
505  }
506  }
507 
508  int
510  Teuchos::Array<Teuchos::RCP<const MV> > Q,
511  Teuchos::Array<Teuchos::RCP<mat_type> > C =
512  Teuchos::tuple (Teuchos::RCP<mat_type> (Teuchos::null)),
513  Teuchos::RCP<mat_type> B = Teuchos::null,
514  Teuchos::RCP<MV> MX = Teuchos::null,
515  Teuchos::Array<Teuchos::RCP<const MV> > MQ =
516  Teuchos::tuple (Teuchos::RCP<const MV> (Teuchos::null))) const
517  {
518  if (getOp().is_null()) {
519  // FIXME (mfh 12 Jan 2011): Do we need to check if C and B
520  // are null and allocate space if so?
521  const int rank = tsqr_.projectAndNormalize (X, C, B, Q);
522  // FIXME (mfh 20 Jul 2010) What about MX and MQ?
523  // mfh 12 Jan 2011: Ignore MQ (assume Q == MQ).
524  //
525  // FIXME (mfh 12 Jan 2011) Need to port MVT::Assign from Belos to Anasazi
526 #if 0
527  if (! MX.is_null()) {
528  // MX gets a copy of X; M is the identity operator.
529  MVT::Assign (X, *MX);
530  }
531 #endif // 0
532  return rank;
533  } else {
534  ensureSvqbInit ();
535  return pSvqb_->projectAndNormalizeMat (X, Q, C, B, MX, MQ);
536  }
537  }
538 
539  int
540  normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const
541  {
542  if (getOp().is_null()) {
543  return tsqr_.normalizeOutOfPlace (X, Q, B);
544  } else {
545  // SVQB normalizes in place, so we have to copy.
546  ensureSvqbInit ();
547  const int rank = pSvqb_->normalize (X, B);
548  MVT::Assign (X, Q);
549  return rank;
550  }
551  }
552 
553  int
555  MV& X_out,
556  Teuchos::Array<mat_ptr> C,
557  mat_ptr B,
558  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
559  {
560  using Teuchos::null;
561 
562  if (getOp().is_null()) {
563  return tsqr_.projectAndNormalizeOutOfPlace (X_in, X_out, C, B, Q);
564  } else {
565  ensureSvqbInit ();
566  // SVQB's projectAndNormalize wants an Array, not an ArrayView.
567  // Copy into a temporary Array and copy back afterwards.
568  Teuchos::Array<Teuchos::RCP<const MV> > Q_array (Q);
569  const int rank = pSvqb_->projectAndNormalize (X_in, Q_array, C, B);
570  Q.assign (Q_array);
571  // SVQB normalizes in place, so we have to copy X_in to X_out.
572  MVT::Assign (X_in, X_out);
573  return rank;
574  }
575  }
576 
577  magnitude_type
578  orthonormErrorMat (const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const
579  {
580  if (getOp().is_null()) {
581  return tsqr_.orthonormError (X);
582  // FIXME (mfh 20 Jul 2010) What about MX?
583  } else {
584  ensureSvqbInit ();
585  return pSvqb_->orthonormErrorMat (X, MX);
586  }
587  }
588 
589  magnitude_type
590  orthogErrorMat (const MV &X,
591  const MV &Y,
592  Teuchos::RCP<const MV> MX = Teuchos::null,
593  Teuchos::RCP<const MV> MY = Teuchos::null) const
594  {
595  if (getOp().is_null()) {
596  return tsqr_.orthogError (X, Y);
597  // FIXME (mfh 20 Jul 2010) What about MX and MY?
598  } else {
599  ensureSvqbInit ();
600  return pSvqb_->orthogErrorMat (X, Y, MX, MY);
601  }
602  }
603 
604  private:
606  void
607  ensureSvqbInit () const
608  {
609  // NOTE (mfh 12 Oct 2011) For now, we rely on the default values
610  // of SVQB parameters.
611  if (pSvqb_.is_null()) {
612  pSvqb_ = Teuchos::rcp (new svqb_type (getOp()));
613  }
614  }
615 
623  mutable tsqr_type tsqr_;
624 
629  mutable Teuchos::RCP<svqb_type> pSvqb_;
630  };
631 
632 } // namespace Anasazi
633 
634 #endif // __AnasaziTsqrOrthoManager_hpp
635 
magnitude_type orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector.
MatOrthoManager subclass using TSQR or SVQB.
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.
TSQR-based OrthoManager subclass implementation.
MV multivector_type
Multivector type with which this class was specialized.
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters()
Get "fast" parameters for TsqrMatOrthoManager.
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.
Mixin for out-of-place orthogonalization.
TsqrMatOrthoManager(const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null)
Constructor (that sets default parameters).
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.
MV multivector_type
Multivector type with which this class was specialized.
virtual ~TsqrMatOrthoManager()
Destructor (declared virtual for memory safety of derived classes).
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Default valid parameter list.
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).
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
virtual int normalizeOutOfPlace(MV &X, MV &Q, mat_ptr B) const =0
Normalize X into Q*B.
int normalizeOutOfPlace(MV &X, MV &Q, mat_ptr B) const
Normalize X into Q*B, overwriting X with invalid values.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
virtual ~OutOfPlaceNormalizerMixin()
Trivial virtual destructor, to silence compiler warnings.
int normalizeOutOfPlace(MV &X, MV &Q, mat_ptr B) const
Normalize X into Q*B.
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.
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
magnitude_type orthonormErrorMat(const MV &X, Teuchos::RCP< const MV > MX=Teuchos::null) const
This method computes the error in orthonormality of a multivector.
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.
Traits class which defines basic operations on multivectors.
Teuchos::RCP< const OP > getOp() const
Get operator used for inner product.
magnitude_type orthogError(const MV &X1, const MV &X2) const
This method computes the error in orthogonality of two multivectors.
void innerProd(const MV &X, const MV &Y, mat_type &Z) const
Provides the inner product defining the orthogonality concepts.
Orthogonalization manager based on the SVQB technique described in "A Block Orthogonalization Procedu...
TsqrOrthoManager(const Teuchos::RCP< Teuchos::ParameterList > &params, const std::string &label="Anasazi")
Constructor (that sets user-specified parameters).
TSQR-based OrthoManager subclass.
virtual ~TsqrOrthoManager()
Destructor, declared virtual for safe inheritance.
TsqrOrthoManager(const std::string &label)
Constructor (that sets default parameters).
OP operator_type
Operator type with which this class was specialized.
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get default parameters for TsqrMatOrthoManager.
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.
virtual void setOp(Teuchos::RCP< const OP > Op)
Set operator used for inner product.
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...
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.