Belos  Version of the Day
BelosIMGSOrthoManager.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 
42 
47 #ifndef BELOS_IMGS_ORTHOMANAGER_HPP
48 #define BELOS_IMGS_ORTHOMANAGER_HPP
49 
57 // #define ORTHO_DEBUG
58 
59 #include "BelosConfigDefs.hpp"
60 #include "BelosMultiVecTraits.hpp"
61 #include "BelosOperatorTraits.hpp"
62 #include "BelosMatOrthoManager.hpp"
63 #include "Teuchos_SerialDenseMatrix.hpp"
64 #include "Teuchos_SerialDenseVector.hpp"
65 
66 #include "Teuchos_as.hpp"
67 #ifdef BELOS_TEUCHOS_TIME_MONITOR
68 #include "Teuchos_TimeMonitor.hpp"
69 #endif // BELOS_TEUCHOS_TIME_MONITOR
70 
71 namespace Belos {
72 
74  template<class ScalarType, class MV, class OP>
75  Teuchos::RCP<Teuchos::ParameterList> getIMGSDefaultParameters ();
76 
78  template<class ScalarType, class MV, class OP>
79  Teuchos::RCP<Teuchos::ParameterList> getIMGSFastParameters();
80 
81  template<class ScalarType, class MV, class OP>
82  class IMGSOrthoManager :
83  public MatOrthoManager<ScalarType,MV,OP>
84  {
85  private:
86  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
87  typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
88  typedef Teuchos::ScalarTraits<ScalarType> SCT;
89  typedef MultiVecTraits<ScalarType,MV> MVT;
90  typedef OperatorTraits<ScalarType,MV,OP> OPT;
91 
92  public:
94 
95 
97  IMGSOrthoManager( const std::string& label = "Belos",
98  Teuchos::RCP<const OP> Op = Teuchos::null,
99  const int max_ortho_steps = max_ortho_steps_default_,
100  const MagnitudeType blk_tol = blk_tol_default_,
101  const MagnitudeType sing_tol = sing_tol_default_ )
102  : MatOrthoManager<ScalarType,MV,OP>(Op),
103  max_ortho_steps_( max_ortho_steps ),
104  blk_tol_( blk_tol ),
105  sing_tol_( sing_tol ),
106  label_( label )
107  {
108 #ifdef BELOS_TEUCHOS_TIME_MONITOR
109  std::stringstream ss;
110  ss << label_ + ": IMGS[" << max_ortho_steps_ << "]";
111 
112  std::string orthoLabel = ss.str() + ": Orthogonalization";
113  timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
114 
115  std::string updateLabel = ss.str() + ": Ortho (Update)";
116  timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
117 
118  std::string normLabel = ss.str() + ": Ortho (Norm)";
119  timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
120 
121  std::string ipLabel = ss.str() + ": Ortho (Inner Product)";
122  timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
123 #endif
124  }
125 
127  IMGSOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& plist,
128  const std::string& label = "Belos",
129  Teuchos::RCP<const OP> Op = Teuchos::null) :
130  MatOrthoManager<ScalarType,MV,OP>(Op),
131  max_ortho_steps_ (max_ortho_steps_default_),
132  blk_tol_ (blk_tol_default_),
133  sing_tol_ (sing_tol_default_),
134  label_ (label)
135  {
136  setParameterList (plist);
137 
138 #ifdef BELOS_TEUCHOS_TIME_MONITOR
139  std::stringstream ss;
140  ss << label_ + ": IMGS[" << max_ortho_steps_ << "]";
141 
142  std::string orthoLabel = ss.str() + ": Orthogonalization";
143  timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
144 
145  std::string updateLabel = ss.str() + ": Ortho (Update)";
146  timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
147 
148  std::string normLabel = ss.str() + ": Ortho (Norm)";
149  timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
150 
151  std::string ipLabel = ss.str() + ": Ortho (Inner Product)";
152  timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
153 #endif
154  }
155 
157  ~IMGSOrthoManager() {}
159 
161 
162  void
163  setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
164  {
165  using Teuchos::Exceptions::InvalidParameterName;
166  using Teuchos::ParameterList;
167  using Teuchos::parameterList;
168  using Teuchos::RCP;
169 
170  RCP<const ParameterList> defaultParams = getValidParameters();
171  RCP<ParameterList> params;
172  if (plist.is_null()) {
173  params = parameterList (*defaultParams);
174  } else {
175  params = plist;
176  // Some users might want to specify "blkTol" as "depTol". Due
177  // to this case, we don't invoke
178  // validateParametersAndSetDefaults on params. Instead, we go
179  // through the parameter list one parameter at a time and look
180  // for alternatives.
181  }
182 
183  // Using temporary variables and fetching all values before
184  // setting the output arguments ensures the strong exception
185  // guarantee for this function: if an exception is thrown, no
186  // externally visible side effects (in this case, setting the
187  // output arguments) have taken place.
188  int maxNumOrthogPasses;
189  MagnitudeType blkTol;
190  MagnitudeType singTol;
191 
192  try {
193  maxNumOrthogPasses = params->get<int> ("maxNumOrthogPasses");
194  } catch (InvalidParameterName&) {
195  maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
196  params->set ("maxNumOrthogPasses", maxNumOrthogPasses);
197  }
198 
199  // Handling of the "blkTol" parameter is a special case. This
200  // is because some users may prefer to call this parameter
201  // "depTol" for consistency with DGKS. However, our default
202  // parameter list calls this "blkTol", and we don't want the
203  // default list's value to override the user's value. Thus, we
204  // first check the user's parameter list for both names, and
205  // only then access the default parameter list.
206  try {
207  blkTol = params->get<MagnitudeType> ("blkTol");
208  } catch (InvalidParameterName&) {
209  try {
210  blkTol = params->get<MagnitudeType> ("depTol");
211  // "depTol" is the wrong name, so remove it and replace with
212  // "blkTol". We'll set "blkTol" below.
213  params->remove ("depTol");
214  } catch (InvalidParameterName&) {
215  blkTol = defaultParams->get<MagnitudeType> ("blkTol");
216  }
217  params->set ("blkTol", blkTol);
218  }
219 
220  try {
221  singTol = params->get<MagnitudeType> ("singTol");
222  } catch (InvalidParameterName&) {
223  singTol = defaultParams->get<MagnitudeType> ("singTol");
224  params->set ("singTol", singTol);
225  }
226 
227  max_ortho_steps_ = maxNumOrthogPasses;
228  blk_tol_ = blkTol;
229  sing_tol_ = singTol;
230 
231  this->setMyParamList (params);
232  }
233 
234  Teuchos::RCP<const Teuchos::ParameterList>
235  getValidParameters () const
236  {
237  if (defaultParams_.is_null()) {
238  defaultParams_ = Belos::getIMGSDefaultParameters<ScalarType, MV, OP>();
239  }
240 
241  return defaultParams_;
242  }
243 
245 
247 
248 
250  void setBlkTol( const MagnitudeType blk_tol ) { blk_tol_ = blk_tol; }
251 
253  void setSingTol( const MagnitudeType sing_tol ) { sing_tol_ = sing_tol; }
254 
256  MagnitudeType getBlkTol() const { return blk_tol_; }
257 
259  MagnitudeType getSingTol() const { return sing_tol_; }
260 
262 
263 
265 
266 
294  void project ( MV &X, Teuchos::RCP<MV> MX,
295  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
296  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
297 
298 
301  void project ( MV &X,
302  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
303  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
304  project(X,Teuchos::null,C,Q);
305  }
306 
307 
308 
333  int normalize ( MV &X, Teuchos::RCP<MV> MX,
334  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B) const;
335 
336 
339  int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
340  return normalize(X,Teuchos::null,B);
341  }
342 
343  protected:
385  virtual int
386  projectAndNormalizeWithMxImpl (MV &X,
387  Teuchos::RCP<MV> MX,
388  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
389  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
390  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
391 
392  public:
394 
396 
400  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
401  orthonormError(const MV &X) const {
402  return orthonormError(X,Teuchos::null);
403  }
404 
409  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
410  orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const;
411 
415  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
416  orthogError(const MV &X1, const MV &X2) const {
417  return orthogError(X1,Teuchos::null,X2);
418  }
419 
424  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
425  orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const;
426 
428 
430 
431 
434  void setLabel(const std::string& label);
435 
438  const std::string& getLabel() const { return label_; }
439 
441 
443 
444 
446  static const int max_ortho_steps_default_;
448  static const MagnitudeType blk_tol_default_;
450  static const MagnitudeType sing_tol_default_;
451 
453  static const int max_ortho_steps_fast_;
455  static const MagnitudeType blk_tol_fast_;
457  static const MagnitudeType sing_tol_fast_;
458 
460 
461  private:
462 
464  int max_ortho_steps_;
466  MagnitudeType blk_tol_;
468  MagnitudeType sing_tol_;
469 
471  std::string label_;
472 #ifdef BELOS_TEUCHOS_TIME_MONITOR
473  Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, timerNorm_, timerInnerProd_;
474 #endif // BELOS_TEUCHOS_TIME_MONITOR
475 
477  mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
478 
480  int findBasis(MV &X, Teuchos::RCP<MV> MX,
481  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
482  bool completeBasis, int howMany = -1 ) const;
483 
485  bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
486  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
487  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
488 
490  bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
491  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
492  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
493 
507  int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
508  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
509  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
510  Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const;
511  };
512 
513  // Set static variables.
514  template<class ScalarType, class MV, class OP>
515  const int IMGSOrthoManager<ScalarType,MV,OP>::max_ortho_steps_default_ = 1;
516 
517  template<class ScalarType, class MV, class OP>
518  const typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
519  IMGSOrthoManager<ScalarType,MV,OP>::blk_tol_default_
520  = 10*Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
521  Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps() );
522 
523  template<class ScalarType, class MV, class OP>
524  const typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
525  IMGSOrthoManager<ScalarType,MV,OP>::sing_tol_default_
526  = 10*Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps();
527 
528  template<class ScalarType, class MV, class OP>
529  const int IMGSOrthoManager<ScalarType,MV,OP>::max_ortho_steps_fast_ = 1;
530 
531  template<class ScalarType, class MV, class OP>
532  const typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
533  IMGSOrthoManager<ScalarType,MV,OP>::blk_tol_fast_
534  = Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
535 
536  template<class ScalarType, class MV, class OP>
537  const typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
538  IMGSOrthoManager<ScalarType,MV,OP>::sing_tol_fast_
539  = Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
540 
542  // Set the label for this orthogonalization manager and create new timers if it's changed
543  template<class ScalarType, class MV, class OP>
544  void IMGSOrthoManager<ScalarType,MV,OP>::setLabel(const std::string& label)
545  {
546  if (label != label_) {
547  label_ = label;
548 #ifdef BELOS_TEUCHOS_TIME_MONITOR
549  std::stringstream ss;
550  ss << label_ + ": IMGS[" << max_ortho_steps_ << "]";
551 
552  std::string orthoLabel = ss.str() + ": Orthogonalization";
553  timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
554 
555  std::string updateLabel = ss.str() + ": Ortho (Update)";
556  timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
557 
558  std::string normLabel = ss.str() + ": Ortho (Norm)";
559  timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
560 
561  std::string ipLabel = ss.str() + ": Ortho (Inner Product)";
562  timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
563 #endif
564  }
565  }
566 
568  // Compute the distance from orthonormality
569  template<class ScalarType, class MV, class OP>
570  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
571  IMGSOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const {
572  const ScalarType ONE = SCT::one();
573  int rank = MVT::GetNumberVecs(X);
574  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
575  MatOrthoManager<ScalarType,MV,OP>::innerProd(X,X,MX,xTx);
576  for (int i=0; i<rank; i++) {
577  xTx(i,i) -= ONE;
578  }
579  return xTx.normFrobenius();
580  }
581 
583  // Compute the distance from orthogonality
584  template<class ScalarType, class MV, class OP>
585  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
586  IMGSOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const {
587  int r1 = MVT::GetNumberVecs(X1);
588  int r2 = MVT::GetNumberVecs(X2);
589  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
590  MatOrthoManager<ScalarType,MV,OP>::innerProd(X2,X1,MX1,xTx);
591  return xTx.normFrobenius();
592  }
593 
595  // Find an Op-orthonormal basis for span(X) - span(W)
596  template<class ScalarType, class MV, class OP>
597  int
598  IMGSOrthoManager<ScalarType, MV, OP>::
599  projectAndNormalizeWithMxImpl(MV &X,
600  Teuchos::RCP<MV> MX,
601  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
602  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
603  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
604  {
605  using Teuchos::Array;
606  using Teuchos::null;
607  using Teuchos::is_null;
608  using Teuchos::RCP;
609  using Teuchos::rcp;
610  using Teuchos::SerialDenseMatrix;
611  typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
612  typedef typename Array< RCP< const MV > >::size_type size_type;
613 
614 #ifdef BELOS_TEUCHOS_TIME_MONITOR
615  Teuchos::TimeMonitor orthotimer(*timerOrtho_);
616 #endif
617 
618  ScalarType ONE = SCT::one();
619  const MagnitudeType ZERO = MGT::zero();
620 
621  int nq = Q.size();
622  int xc = MVT::GetNumberVecs( X );
623  ptrdiff_t xr = MVT::GetGlobalLength( X );
624  int rank = xc;
625 
626  // If the user doesn't want to store the normalization
627  // coefficients, allocate some local memory for them. This will
628  // go away at the end of this method.
629  if (is_null (B)) {
630  B = rcp (new serial_dense_matrix_type (xc, xc));
631  }
632  // Likewise, if the user doesn't want to store the projection
633  // coefficients, allocate some local memory for them. Also make
634  // sure that all the entries of C are the right size. We're going
635  // to overwrite them anyway, so we don't have to worry about the
636  // contents (other than to resize them if they are the wrong
637  // size).
638  if (C.size() < nq)
639  C.resize (nq);
640  for (size_type k = 0; k < nq; ++k)
641  {
642  const int numRows = MVT::GetNumberVecs (*Q[k]);
643  const int numCols = xc; // Number of vectors in X
644 
645  if (is_null (C[k]))
646  C[k] = rcp (new serial_dense_matrix_type (numRows, numCols));
647  else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
648  {
649  int err = C[k]->reshape (numRows, numCols);
650  TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error,
651  "IMGS orthogonalization: failed to reshape "
652  "C[" << k << "] (the array of block "
653  "coefficients resulting from projecting X "
654  "against Q[1:" << nq << "]).");
655  }
656  }
657 
658  /****** DO NOT MODIFY *MX IF _hasOp == false ******/
659  if (this->_hasOp) {
660  if (MX == Teuchos::null) {
661  // we need to allocate space for MX
662  MX = MVT::Clone(X,MVT::GetNumberVecs(X));
663  OPT::Apply(*(this->_Op),X,*MX);
664  }
665  }
666  else {
667  // Op == I --> MX = X (ignore it if the user passed it in)
668  MX = Teuchos::rcp( &X, false );
669  }
670 
671  int mxc = MVT::GetNumberVecs( *MX );
672  ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
673 
674  // short-circuit
675  TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Belos::IMGSOrthoManager::projectAndNormalize(): X must be non-empty" );
676 
677  int numbas = 0;
678  for (int i=0; i<nq; i++) {
679  numbas += MVT::GetNumberVecs( *Q[i] );
680  }
681 
682  // check size of B
683  TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
684  "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
685  // check size of X and MX
686  TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
687  "Belos::IMGSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
688  // check size of X w.r.t. MX
689  TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
690  "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
691  // check feasibility
692  //TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument,
693  // "Belos::IMGSOrthoManager::projectAndNormalize(): Orthogonality constraints not feasible" );
694 
695  // Some flags for checking dependency returns from the internal orthogonalization methods
696  bool dep_flg = false;
697 
698  // Make a temporary copy of X and MX, just in case a block dependency is detected.
699  Teuchos::RCP<MV> tmpX, tmpMX;
700  tmpX = MVT::CloneCopy(X);
701  if (this->_hasOp) {
702  tmpMX = MVT::CloneCopy(*MX);
703  }
704 
705  if (xc == 1) {
706 
707  // Use the cheaper block orthogonalization.
708  // NOTE: Don't check for dependencies because the update has one vector.
709  dep_flg = blkOrtho1( X, MX, C, Q );
710 
711  // Normalize the new block X
712  if ( B == Teuchos::null ) {
713  B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
714  }
715  std::vector<ScalarType> diag(xc);
716  {
717 #ifdef BELOS_TEUCHOS_TIME_MONITOR
718  Teuchos::TimeMonitor normTimer( *timerNorm_ );
719 #endif
720  MVT::MvDot( X, *MX, diag );
721  }
722  (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
723 
724  if (SCT::magnitude((*B)(0,0)) > ZERO) {
725  rank = 1;
726  MVT::MvScale( X, ONE/(*B)(0,0) );
727  if (this->_hasOp) {
728  // Update MXj.
729  MVT::MvScale( *MX, ONE/(*B)(0,0) );
730  }
731  }
732  }
733  else {
734 
735  // Use the cheaper block orthogonalization.
736  dep_flg = blkOrtho( X, MX, C, Q );
737 
738  // If a dependency has been detected in this block, then perform
739  // the more expensive nonblock (single vector at a time)
740  // orthogonalization.
741  if (dep_flg) {
742  rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
743 
744  // Copy tmpX back into X.
745  MVT::Assign( *tmpX, X );
746  if (this->_hasOp) {
747  MVT::Assign( *tmpMX, *MX );
748  }
749  }
750  else {
751  // There is no dependency, so orthonormalize new block X
752  rank = findBasis( X, MX, B, false );
753  if (rank < xc) {
754  // A dependency was found during orthonormalization of X,
755  // rerun orthogonalization using more expensive single-
756  // vector orthogonalization.
757  rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
758 
759  // Copy tmpX back into X.
760  MVT::Assign( *tmpX, X );
761  if (this->_hasOp) {
762  MVT::Assign( *tmpMX, *MX );
763  }
764  }
765  }
766  } // if (xc == 1) {
767 
768  // this should not raise an std::exception; but our post-conditions oblige us to check
769  TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
770  "Belos::IMGSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
771 
772  // Return the rank of X.
773  return rank;
774  }
775 
776 
777 
779  // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
780  template<class ScalarType, class MV, class OP>
781  int IMGSOrthoManager<ScalarType, MV, OP>::normalize(
782  MV &X, Teuchos::RCP<MV> MX,
783  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
784 
785 #ifdef BELOS_TEUCHOS_TIME_MONITOR
786  Teuchos::TimeMonitor orthotimer(*timerOrtho_);
787 #endif
788 
789  // call findBasis, with the instruction to try to generate a basis of rank numvecs(X)
790  return findBasis(X, MX, B, true);
791  }
792 
793 
794 
796  template<class ScalarType, class MV, class OP>
797  void IMGSOrthoManager<ScalarType, MV, OP>::project(
798  MV &X, Teuchos::RCP<MV> MX,
799  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
800  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
801  // For the inner product defined by the operator Op or the identity (Op == 0)
802  // -> Orthogonalize X against each Q[i]
803  // Modify MX accordingly
804  //
805  // Note that when Op is 0, MX is not referenced
806  //
807  // Parameter variables
808  //
809  // X : Vectors to be transformed
810  //
811  // MX : Image of the block of vectors X by the mass matrix
812  //
813  // Q : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
814  //
815 
816 #ifdef BELOS_TEUCHOS_TIME_MONITOR
817  Teuchos::TimeMonitor orthotimer(*timerOrtho_);
818 #endif
819 
820  int xc = MVT::GetNumberVecs( X );
821  ptrdiff_t xr = MVT::GetGlobalLength( X );
822  int nq = Q.size();
823  std::vector<int> qcs(nq);
824  // short-circuit
825  if (nq == 0 || xc == 0 || xr == 0) {
826  return;
827  }
828  ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
829  // if we don't have enough C, expand it with null references
830  // if we have too many, resize to throw away the latter ones
831  // if we have exactly as many as we have Q, this call has no effect
832  C.resize(nq);
833 
834 
835  /****** DO NOT MODIFY *MX IF _hasOp == false ******/
836  if (this->_hasOp) {
837  if (MX == Teuchos::null) {
838  // we need to allocate space for MX
839  MX = MVT::Clone(X,MVT::GetNumberVecs(X));
840  OPT::Apply(*(this->_Op),X,*MX);
841  }
842  }
843  else {
844  // Op == I --> MX = X (ignore it if the user passed it in)
845  MX = Teuchos::rcp( &X, false );
846  }
847  int mxc = MVT::GetNumberVecs( *MX );
848  ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
849 
850  // check size of X and Q w.r.t. common sense
851  TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
852  "Belos::IMGSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
853  // check size of X w.r.t. MX and Q
854  TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
855  "Belos::IMGSOrthoManager::project(): Size of X not consistant with MX,Q" );
856 
857  // tally up size of all Q and check/allocate C
858  int baslen = 0;
859  for (int i=0; i<nq; i++) {
860  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
861  "Belos::IMGSOrthoManager::project(): Q lengths not mutually consistant" );
862  qcs[i] = MVT::GetNumberVecs( *Q[i] );
863  TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
864  "Belos::IMGSOrthoManager::project(): Q has less rows than columns" );
865  baslen += qcs[i];
866 
867  // check size of C[i]
868  if ( C[i] == Teuchos::null ) {
869  C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
870  }
871  else {
872  TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
873  "Belos::IMGSOrthoManager::project(): Size of Q not consistant with size of C" );
874  }
875  }
876 
877  // Use the cheaper block orthogonalization, don't check for rank deficiency.
878  blkOrtho( X, MX, C, Q );
879 
880  }
881 
883  // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that
884  // the rank is numvectors(X)
885  template<class ScalarType, class MV, class OP>
886  int IMGSOrthoManager<ScalarType, MV, OP>::findBasis(
887  MV &X, Teuchos::RCP<MV> MX,
888  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
889  bool completeBasis, int howMany ) const {
890  // For the inner product defined by the operator Op or the identity (Op == 0)
891  // -> Orthonormalize X
892  // Modify MX accordingly
893  //
894  // Note that when Op is 0, MX is not referenced
895  //
896  // Parameter variables
897  //
898  // X : Vectors to be orthonormalized
899  //
900  // MX : Image of the multivector X under the operator Op
901  //
902  // Op : Pointer to the operator for the inner product
903  //
904  //
905 
906  const ScalarType ONE = SCT::one();
907  const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
908 
909  int xc = MVT::GetNumberVecs( X );
910  ptrdiff_t xr = MVT::GetGlobalLength( X );
911 
912  if (howMany == -1) {
913  howMany = xc;
914  }
915 
916  /*******************************************************
917  * If _hasOp == false, we will not reference MX below *
918  *******************************************************/
919 
920  // if Op==null, MX == X (via pointer)
921  // Otherwise, either the user passed in MX or we will allocated and compute it
922  if (this->_hasOp) {
923  if (MX == Teuchos::null) {
924  // we need to allocate space for MX
925  MX = MVT::Clone(X,xc);
926  OPT::Apply(*(this->_Op),X,*MX);
927  }
928  }
929 
930  /* if the user doesn't want to store the coefficienets,
931  * allocate some local memory for them
932  */
933  if ( B == Teuchos::null ) {
934  B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
935  }
936 
937  int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
938  ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
939 
940  // check size of C, B
941  TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
942  "Belos::IMGSOrthoManager::findBasis(): X must be non-empty" );
943  TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
944  "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
945  TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
946  "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
947  TEUCHOS_TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument,
948  "Belos::IMGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
949  TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
950  "Belos::IMGSOrthoManager::findBasis(): Invalid howMany parameter" );
951 
952  /* xstart is which column we are starting the process with, based on howMany
953  * columns before xstart are assumed to be Op-orthonormal already
954  */
955  int xstart = xc - howMany;
956 
957  for (int j = xstart; j < xc; j++) {
958 
959  // numX is
960  // * number of currently orthonormal columns of X
961  // * the index of the current column of X
962  int numX = j;
963  bool addVec = false;
964 
965  // Get a view of the vector currently being worked on.
966  std::vector<int> index(1);
967  index[0] = numX;
968  Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
969  Teuchos::RCP<MV> MXj;
970  if ((this->_hasOp)) {
971  // MXj is a view of the current vector in MX
972  MXj = MVT::CloneViewNonConst( *MX, index );
973  }
974  else {
975  // MXj is a pointer to Xj, and MUST NOT be modified
976  MXj = Xj;
977  }
978 
979  Teuchos::RCP<MV> oldMXj;
980  if (numX > 0) {
981  oldMXj = MVT::CloneCopy( *MXj );
982  }
983 
984  // Make storage for these Gram-Schmidt iterations.
985  Teuchos::SerialDenseVector<int,ScalarType> product(numX);
986  Teuchos::SerialDenseVector<int,ScalarType> P2(1);
987  Teuchos::RCP<const MV> prevX, prevMX;
988 
989  std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
990  //
991  // Save old MXj vector and compute Op-norm
992  //
993  {
994 #ifdef BELOS_TEUCHOS_TIME_MONITOR
995  Teuchos::TimeMonitor normTimer( *timerNorm_ );
996 #endif
997  MVT::MvDot( *Xj, *MXj, oldDot );
998  }
999  // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
1000  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError,
1001  "Belos::IMGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
1002 
1003  // Perform MGS one vector at a time
1004  for (int ii=0; ii<numX; ii++) {
1005 
1006  index[0] = ii;
1007  prevX = MVT::CloneView( X, index );
1008  if (this->_hasOp) {
1009  prevMX = MVT::CloneView( *MX, index );
1010  }
1011 
1012  for (int i=0; i<max_ortho_steps_; ++i) {
1013 
1014  // product <- prevX^T MXj
1015  {
1016 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1017  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1018 #endif
1019  MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*Xj,MXj,P2);
1020  }
1021 
1022  // Xj <- Xj - prevX prevX^T MXj
1023  // = Xj - prevX product
1024  {
1025 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1026  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1027 #endif
1028  MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1029  }
1030 
1031  // Update MXj
1032  if (this->_hasOp) {
1033  // MXj <- Op*Xj_new
1034  // = Op*(Xj_old - prevX prevX^T MXj)
1035  // = MXj - prevMX product
1036 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1037  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1038 #endif
1039  MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1040  }
1041 
1042  // Set coefficients
1043  if ( i==0 )
1044  product[ii] = P2[0];
1045  else
1046  product[ii] += P2[0];
1047 
1048  } // for (int i=0; i<max_ortho_steps_; ++i)
1049 
1050  } // for (int ii=0; ii<numX; ++ii)
1051 
1052  // Compute Op-norm with old MXj
1053  if (numX > 0) {
1054 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1055  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1056 #endif
1057  MVT::MvDot( *Xj, *oldMXj, newDot );
1058  }
1059  else {
1060  newDot[0] = oldDot[0];
1061  }
1062 
1063  // Check to see if the new vector is dependent.
1064  if (completeBasis) {
1065  //
1066  // We need a complete basis, so add random vectors if necessary
1067  //
1068  if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1069 
1070  // Add a random vector and orthogonalize it against previous vectors in block.
1071  addVec = true;
1072 #ifdef ORTHO_DEBUG
1073  std::cout << "Belos::IMGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1074 #endif
1075  //
1076  Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1077  Teuchos::RCP<MV> tempMXj;
1078  MVT::MvRandom( *tempXj );
1079  if (this->_hasOp) {
1080  tempMXj = MVT::Clone( X, 1 );
1081  OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1082  }
1083  else {
1084  tempMXj = tempXj;
1085  }
1086  {
1087 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1088  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1089 #endif
1090  MVT::MvDot( *tempXj, *tempMXj, oldDot );
1091  }
1092  //
1093  // Perform MGS one vector at a time
1094  for (int ii=0; ii<numX; ii++) {
1095 
1096  index[0] = ii;
1097  prevX = MVT::CloneView( X, index );
1098  if (this->_hasOp) {
1099  prevMX = MVT::CloneView( *MX, index );
1100  }
1101 
1102  for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
1103  {
1104 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1105  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1106 #endif
1107  MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*tempXj,tempMXj,P2);
1108  }
1109  {
1110 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1111  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1112 #endif
1113  MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *tempXj );
1114  }
1115  if (this->_hasOp) {
1116 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1117  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1118 #endif
1119  MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *tempMXj );
1120  }
1121 
1122  // Set coefficients
1123  if ( num_orth==0 )
1124  product[ii] = P2[0];
1125  else
1126  product[ii] += P2[0];
1127  }
1128  }
1129 
1130  // Compute new Op-norm
1131  {
1132 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1133  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1134 #endif
1135  MVT::MvDot( *tempXj, *tempMXj, newDot );
1136  }
1137  //
1138  if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1139  // Copy vector into current column of _basisvecs
1140  MVT::Assign( *tempXj, *Xj );
1141  if (this->_hasOp) {
1142  MVT::Assign( *tempMXj, *MXj );
1143  }
1144  }
1145  else {
1146  return numX;
1147  }
1148  }
1149  }
1150  else {
1151  //
1152  // We only need to detect dependencies.
1153  //
1154  if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1155  return numX;
1156  }
1157  }
1158 
1159 
1160  // If we haven't left this method yet, then we can normalize the new vector Xj.
1161  // Normalize Xj.
1162  // Xj <- Xj / std::sqrt(newDot)
1163  ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1164  if (SCT::magnitude(diag) > ZERO) {
1165  MVT::MvScale( *Xj, ONE/diag );
1166  if (this->_hasOp) {
1167  // Update MXj.
1168  MVT::MvScale( *MXj, ONE/diag );
1169  }
1170  }
1171 
1172  // If we've added a random vector, enter a zero in the j'th diagonal element.
1173  if (addVec) {
1174  (*B)(j,j) = ZERO;
1175  }
1176  else {
1177  (*B)(j,j) = diag;
1178  }
1179 
1180  // Save the coefficients, if we are working on the original vector and not a randomly generated one
1181  if (!addVec) {
1182  for (int i=0; i<numX; i++) {
1183  (*B)(i,j) = product(i);
1184  }
1185  }
1186 
1187  } // for (j = 0; j < xc; ++j)
1188 
1189  return xc;
1190  }
1191 
1193  // Routine to compute the block orthogonalization
1194  template<class ScalarType, class MV, class OP>
1195  bool
1196  IMGSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
1197  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1198  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
1199  {
1200  int nq = Q.size();
1201  int xc = MVT::GetNumberVecs( X );
1202  const ScalarType ONE = SCT::one();
1203 
1204  std::vector<int> qcs( nq );
1205  for (int i=0; i<nq; i++) {
1206  qcs[i] = MVT::GetNumberVecs( *Q[i] );
1207  }
1208 
1209  // Perform the Gram-Schmidt transformation for a block of vectors
1210  std::vector<int> index(1);
1211  Teuchos::RCP<const MV> tempQ;
1212 
1213  Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1214  // Define the product Q^T * (Op*X)
1215  for (int i=0; i<nq; i++) {
1216 
1217  // Perform MGS one vector at a time
1218  for (int ii=0; ii<qcs[i]; ii++) {
1219 
1220  index[0] = ii;
1221  tempQ = MVT::CloneView( *Q[i], index );
1222  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1223 
1224  // Multiply Q' with MX
1225  {
1226 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1227  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1228 #endif
1229  MatOrthoManager<ScalarType,MV,OP>::innerProd(*tempQ,X,MX,tempC);
1230  }
1231  // Multiply by Q and subtract the result in X
1232  {
1233 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1234  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1235 #endif
1236  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1237  }
1238  }
1239 
1240  // Update MX, with the least number of applications of Op as possible
1241  if (this->_hasOp) {
1242  if (xc <= qcs[i]) {
1243  OPT::Apply( *(this->_Op), X, *MX);
1244  }
1245  else {
1246  // this will possibly be used again below; don't delete it
1247  MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1248  OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1249  MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1250  }
1251  }
1252  }
1253 
1254  // Do as many steps of modified Gram-Schmidt as required by max_ortho_steps_
1255  for (int j = 1; j < max_ortho_steps_; ++j) {
1256 
1257  for (int i=0; i<nq; i++) {
1258 
1259  Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],1);
1260 
1261  // Perform MGS one vector at a time
1262  for (int ii=0; ii<qcs[i]; ii++) {
1263 
1264  index[0] = ii;
1265  tempQ = MVT::CloneView( *Q[i], index );
1266  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1267  Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii, 0 );
1268 
1269  // Apply another step of modified Gram-Schmidt
1270  {
1271 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1272  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1273 #endif
1274  MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC2 );
1275  }
1276  tempC += tempC2;
1277  {
1278 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1279  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1280 #endif
1281  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1282  }
1283 
1284  }
1285 
1286  // Update MX, with the least number of applications of Op as possible
1287  if (this->_hasOp) {
1288  if (MQ[i].get()) {
1289  // MQ was allocated and computed above; use it
1290  MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1291  }
1292  else if (xc <= qcs[i]) {
1293  // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1294  OPT::Apply( *(this->_Op), X, *MX);
1295  }
1296  }
1297  } // for (int i=0; i<nq; i++)
1298  } // for (int j = 0; j < max_ortho_steps; ++j)
1299 
1300  return false;
1301  }
1302 
1304  // Routine to compute the block orthogonalization
1305  template<class ScalarType, class MV, class OP>
1306  bool
1307  IMGSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
1308  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1309  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
1310  {
1311  int nq = Q.size();
1312  int xc = MVT::GetNumberVecs( X );
1313  bool dep_flg = false;
1314  const ScalarType ONE = SCT::one();
1315 
1316  std::vector<int> qcs( nq );
1317  for (int i=0; i<nq; i++) {
1318  qcs[i] = MVT::GetNumberVecs( *Q[i] );
1319  }
1320 
1321  // Perform the Gram-Schmidt transformation for a block of vectors
1322 
1323  // Compute the initial Op-norms
1324  std::vector<ScalarType> oldDot( xc );
1325  {
1326 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1327  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1328 #endif
1329  MVT::MvDot( X, *MX, oldDot );
1330  }
1331 
1332  std::vector<int> index(1);
1333  Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1334  Teuchos::RCP<const MV> tempQ;
1335 
1336  // Define the product Q^T * (Op*X)
1337  for (int i=0; i<nq; i++) {
1338 
1339  // Perform MGS one vector at a time
1340  for (int ii=0; ii<qcs[i]; ii++) {
1341 
1342  index[0] = ii;
1343  tempQ = MVT::CloneView( *Q[i], index );
1344  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1345 
1346  // Multiply Q' with MX
1347  {
1348 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1349  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1350 #endif
1351  MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC);
1352  }
1353  // Multiply by Q and subtract the result in X
1354  {
1355 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1356  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1357 #endif
1358  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1359  }
1360  }
1361 
1362  // Update MX, with the least number of applications of Op as possible
1363  if (this->_hasOp) {
1364  if (xc <= qcs[i]) {
1365  OPT::Apply( *(this->_Op), X, *MX);
1366  }
1367  else {
1368  // this will possibly be used again below; don't delete it
1369  MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1370  OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1371  MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1372  }
1373  }
1374  }
1375 
1376  // Do as many steps of modified Gram-Schmidt as required by max_ortho_steps_
1377  for (int j = 1; j < max_ortho_steps_; ++j) {
1378 
1379  for (int i=0; i<nq; i++) {
1380  Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],xc);
1381 
1382  // Perform MGS one vector at a time
1383  for (int ii=0; ii<qcs[i]; ii++) {
1384 
1385  index[0] = ii;
1386  tempQ = MVT::CloneView( *Q[i], index );
1387  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1388  Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, xc, ii, 0 );
1389 
1390  // Apply another step of modified Gram-Schmidt
1391  {
1392 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1393  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1394 #endif
1395  MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC2 );
1396  }
1397  tempC += tempC2;
1398  {
1399 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1400  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1401 #endif
1402  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1403  }
1404  }
1405 
1406  // Update MX, with the least number of applications of Op as possible
1407  if (this->_hasOp) {
1408  if (MQ[i].get()) {
1409  // MQ was allocated and computed above; use it
1410  MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1411  }
1412  else if (xc <= qcs[i]) {
1413  // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1414  OPT::Apply( *(this->_Op), X, *MX);
1415  }
1416  }
1417  } // for (int i=0; i<nq; i++)
1418  } // for (int j = 0; j < max_ortho_steps; ++j)
1419 
1420  // Compute new Op-norms
1421  std::vector<ScalarType> newDot(xc);
1422  {
1423 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1424  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1425 #endif
1426  MVT::MvDot( X, *MX, newDot );
1427  }
1428 
1429  // Check to make sure the new block of vectors are not dependent on previous vectors
1430  for (int i=0; i<xc; i++){
1431  if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1432  dep_flg = true;
1433  break;
1434  }
1435  } // end for (i=0;...)
1436 
1437  return dep_flg;
1438  }
1439 
1440  template<class ScalarType, class MV, class OP>
1441  int
1442  IMGSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
1443  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1444  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1445  Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const
1446  {
1447  Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1448 
1449  const ScalarType ONE = SCT::one();
1450  const ScalarType ZERO = SCT::zero();
1451 
1452  int nq = Q.size();
1453  int xc = MVT::GetNumberVecs( X );
1454  std::vector<int> indX( 1 );
1455  std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1456 
1457  std::vector<int> qcs( nq );
1458  for (int i=0; i<nq; i++) {
1459  qcs[i] = MVT::GetNumberVecs( *Q[i] );
1460  }
1461 
1462  // Create pointers for the previous vectors of X that have already been orthonormalized.
1463  Teuchos::RCP<const MV> lastQ;
1464  Teuchos::RCP<MV> Xj, MXj;
1465  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1466 
1467  // Perform the Gram-Schmidt transformation for each vector in the block of vectors.
1468  for (int j=0; j<xc; j++) {
1469 
1470  bool dep_flg = false;
1471 
1472  // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
1473  if (j > 0) {
1474  std::vector<int> index( j );
1475  for (int ind=0; ind<j; ind++) {
1476  index[ind] = ind;
1477  }
1478  lastQ = MVT::CloneView( X, index );
1479 
1480  // Add these views to the Q and C arrays.
1481  Q.push_back( lastQ );
1482  C.push_back( B );
1483  qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1484  }
1485 
1486  // Get a view of the current vector in X to orthogonalize.
1487  indX[0] = j;
1488  Xj = MVT::CloneViewNonConst( X, indX );
1489  if (this->_hasOp) {
1490  MXj = MVT::CloneViewNonConst( *MX, indX );
1491  }
1492  else {
1493  MXj = Xj;
1494  }
1495 
1496  // Compute the initial Op-norms
1497  {
1498 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1499  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1500 #endif
1501  MVT::MvDot( *Xj, *MXj, oldDot );
1502  }
1503 
1504  Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1505  Teuchos::RCP<const MV> tempQ;
1506 
1507  // Define the product Q^T * (Op*X)
1508  for (int i=0; i<Q.size(); i++) {
1509 
1510  // Perform MGS one vector at a time
1511  for (int ii=0; ii<qcs[i]; ii++) {
1512 
1513  indX[0] = ii;
1514  tempQ = MVT::CloneView( *Q[i], indX );
1515  // Get a view of the current serial dense matrix
1516  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, j );
1517 
1518  // Multiply Q' with MX
1519  MatOrthoManager<ScalarType,MV,OP>::innerProd(*tempQ,*Xj,MXj,tempC);
1520 
1521  // Multiply by Q and subtract the result in Xj
1522  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *Xj );
1523  }
1524 
1525  // Update MXj, with the least number of applications of Op as possible
1526  if (this->_hasOp) {
1527  if (xc <= qcs[i]) {
1528  OPT::Apply( *(this->_Op), *Xj, *MXj);
1529  }
1530  else {
1531  // this will possibly be used again below; don't delete it
1532  MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1533  OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1534  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1535  MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1536  }
1537  }
1538  }
1539 
1540  // Do any additional steps of modified Gram-Schmidt orthogonalization
1541  for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
1542 
1543  for (int i=0; i<Q.size(); i++) {
1544  Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1545 
1546  // Perform MGS one vector at a time
1547  for (int ii=0; ii<qcs[i]; ii++) {
1548 
1549  indX[0] = ii;
1550  tempQ = MVT::CloneView( *Q[i], indX );
1551  // Get a view of the current serial dense matrix
1552  Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii );
1553 
1554  // Apply another step of modified Gram-Schmidt
1555  MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, *Xj, MXj, tempC2);
1556  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, *Xj );
1557  }
1558 
1559  // Add the coefficients into C[i]
1560  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1561  tempC += C2;
1562 
1563  // Update MXj, with the least number of applications of Op as possible
1564  if (this->_hasOp) {
1565  if (MQ[i].get()) {
1566  // MQ was allocated and computed above; use it
1567  MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1568  }
1569  else if (xc <= qcs[i]) {
1570  // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1571  OPT::Apply( *(this->_Op), *Xj, *MXj);
1572  }
1573  }
1574  } // for (int i=0; i<Q.size(); i++)
1575 
1576  } // for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps)
1577 
1578  // Compute the Op-norms after the correction step.
1579  {
1580 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1581  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1582 #endif
1583  MVT::MvDot( *Xj, *MXj, newDot );
1584  }
1585 
1586  // Check for linear dependence.
1587  if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1588  dep_flg = true;
1589  }
1590 
1591  // Normalize the new vector if it's not dependent
1592  if (!dep_flg) {
1593  ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1594 
1595  MVT::MvScale( *Xj, ONE/diag );
1596  if (this->_hasOp) {
1597  // Update MXj.
1598  MVT::MvScale( *MXj, ONE/diag );
1599  }
1600 
1601  // Enter value on diagonal of B.
1602  (*B)(j,j) = diag;
1603  }
1604  else {
1605  // Create a random vector and orthogonalize it against all previous columns of Q.
1606  Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1607  Teuchos::RCP<MV> tempMXj;
1608  MVT::MvRandom( *tempXj );
1609  if (this->_hasOp) {
1610  tempMXj = MVT::Clone( X, 1 );
1611  OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1612  }
1613  else {
1614  tempMXj = tempXj;
1615  }
1616  {
1617 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1618  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1619 #endif
1620  MVT::MvDot( *tempXj, *tempMXj, oldDot );
1621  }
1622  //
1623  for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
1624 
1625  for (int i=0; i<Q.size(); i++) {
1626  Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1627 
1628  // Perform MGS one vector at a time
1629  for (int ii=0; ii<qcs[i]; ii++) {
1630 
1631  indX[0] = ii;
1632  tempQ = MVT::CloneView( *Q[i], indX );
1633  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, product, 1, 1, ii );
1634 
1635  // Apply another step of modified Gram-Schmidt
1636  MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, *tempXj, tempMXj, tempC );
1637  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *tempXj );
1638 
1639  }
1640 
1641  // Update MXj, with the least number of applications of Op as possible
1642  if (this->_hasOp) {
1643  if (MQ[i].get()) {
1644  // MQ was allocated and computed above; use it
1645  MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1646  }
1647  else if (xc <= qcs[i]) {
1648  // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1649  OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1650  }
1651  }
1652  } // for (int i=0; i<nq; i++)
1653  } // for (int num_orth=0; num_orth<max_orth_steps_; num_orth++)
1654 
1655  // Compute the Op-norms after the correction step.
1656  {
1657 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1658  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1659 #endif
1660  MVT::MvDot( *tempXj, *tempMXj, newDot );
1661  }
1662 
1663  // Copy vector into current column of Xj
1664  if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1665  ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1666 
1667  // Enter value on diagonal of B.
1668  (*B)(j,j) = ZERO;
1669 
1670  // Copy vector into current column of _basisvecs
1671  MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1672  if (this->_hasOp) {
1673  MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1674  }
1675  }
1676  else {
1677  return j;
1678  }
1679  } // if (!dep_flg)
1680 
1681  // Remove the vectors from array
1682  if (j > 0) {
1683  Q.resize( nq );
1684  C.resize( nq );
1685  qcs.resize( nq );
1686  }
1687 
1688  } // for (int j=0; j<xc; j++)
1689 
1690  return xc;
1691  }
1692 
1693  template<class ScalarType, class MV, class OP>
1694  Teuchos::RCP<Teuchos::ParameterList> getIMGSDefaultParameters ()
1695  {
1696  using Teuchos::ParameterList;
1697  using Teuchos::parameterList;
1698  using Teuchos::RCP;
1699 
1700  RCP<ParameterList> params = parameterList ("IMGS");
1701 
1702  // Default parameter values for IMGS orthogonalization.
1703  // Documentation will be embedded in the parameter list.
1704  params->set ("maxNumOrthogPasses", IMGSOrthoManager<ScalarType, MV, OP>::max_ortho_steps_default_,
1705  "Maximum number of orthogonalization passes (includes the "
1706  "first). Default is 2, since \"twice is enough\" for Krylov "
1707  "methods.");
1708  params->set ("blkTol", IMGSOrthoManager<ScalarType, MV, OP>::blk_tol_default_,
1709  "Block reorthogonalization threshold.");
1710  params->set ("singTol", IMGSOrthoManager<ScalarType, MV, OP>::sing_tol_default_,
1711  "Singular block detection threshold.");
1712 
1713  return params;
1714  }
1715 
1716  template<class ScalarType, class MV, class OP>
1717  Teuchos::RCP<Teuchos::ParameterList> getIMGSFastParameters ()
1718  {
1719  using Teuchos::ParameterList;
1720  using Teuchos::RCP;
1721 
1722  RCP<ParameterList> params = getIMGSDefaultParameters<ScalarType, MV, OP>();
1723 
1724  params->set ("maxNumOrthogPasses",
1725  IMGSOrthoManager<ScalarType, MV, OP>::max_ortho_steps_fast_);
1726  params->set ("blkTol",
1727  IMGSOrthoManager<ScalarType, MV, OP>::blk_tol_fast_);
1728  params->set ("singTol",
1729  IMGSOrthoManager<ScalarType, MV, OP>::sing_tol_fast_);
1730 
1731  return params;
1732  }
1733 
1734 } // namespace Belos
1735 
1736 #endif // BELOS_IMGS_ORTHOMANAGER_HPP
1737 
BelosConfigDefs.hpp
Belos header file which uses auto-configuration information to include necessary C++ headers.
BelosMatOrthoManager.hpp
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
BelosMultiVecTraits.hpp
Declaration of basic traits for the multivector type.
BelosOperatorTraits.hpp
Class which defines basic traits for the operator type.
Belos::IMGSOrthoManager
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
Definition: BelosIMGSOrthoManager.hpp:84
Belos::IMGSOrthoManager::setParameterList
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
Definition: BelosIMGSOrthoManager.hpp:163
Belos::IMGSOrthoManager::orthogError
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthogError(const MV &X1, const MV &X2) const
This method computes the error in orthogonality of two multivectors, measured as the Frobenius norm o...
Definition: BelosIMGSOrthoManager.hpp:416
Belos::IMGSOrthoManager::blk_tol_fast_
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
Definition: BelosIMGSOrthoManager.hpp:455
Belos::IMGSOrthoManager::normalize
int normalize(MV &X, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B) const
This method calls normalize(X,Teuchos::null,B); see documentation for that function.
Definition: BelosIMGSOrthoManager.hpp:339
Belos::IMGSOrthoManager::IMGSOrthoManager
IMGSOrthoManager(const Teuchos::RCP< Teuchos::ParameterList > &plist, const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null)
Constructor that takes a list of parameters.
Definition: BelosIMGSOrthoManager.hpp:127
Belos::IMGSOrthoManager::normalize
int normalize(MV &X, Teuchos::RCP< MV > MX, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B) const
This method takes a multivector X and attempts to compute an orthonormal basis for ,...
Definition: BelosIMGSOrthoManager.hpp:781
Belos::IMGSOrthoManager::setLabel
void setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
Definition: BelosIMGSOrthoManager.hpp:544
Belos::IMGSOrthoManager::orthonormError
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector, measured as the Frobenius norm of ...
Definition: BelosIMGSOrthoManager.hpp:401
Belos::IMGSOrthoManager::max_ortho_steps_fast_
static const int max_ortho_steps_fast_
Max number of (re)orthogonalization steps, including the first (fast).
Definition: BelosIMGSOrthoManager.hpp:453
Belos::IMGSOrthoManager::project
void project(MV &X, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
This method calls project(X,Teuchos::null,C,Q); see documentation for that function.
Definition: BelosIMGSOrthoManager.hpp:301
Belos::IMGSOrthoManager::getLabel
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
Definition: BelosIMGSOrthoManager.hpp:438
Belos::IMGSOrthoManager::project
void project(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Given a list of (mutually and internally) orthonormal bases Q, this method takes a multivector X and ...
Definition: BelosIMGSOrthoManager.hpp:797
Belos::IMGSOrthoManager::sing_tol_default_
static const MagnitudeType sing_tol_default_
Singular block detection threshold (default).
Definition: BelosIMGSOrthoManager.hpp:450
Belos::IMGSOrthoManager::max_ortho_steps_default_
static const int max_ortho_steps_default_
Max number of (re)orthogonalization steps, including the first (default).
Definition: BelosIMGSOrthoManager.hpp:446
Belos::IMGSOrthoManager::sing_tol_fast_
static const MagnitudeType sing_tol_fast_
Singular block detection threshold (fast).
Definition: BelosIMGSOrthoManager.hpp:457
Belos::IMGSOrthoManager::blk_tol_default_
static const MagnitudeType blk_tol_default_
Block reorthogonalization threshold (default).
Definition: BelosIMGSOrthoManager.hpp:448
Belos::IMGSOrthoManager::getBlkTol
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
Definition: BelosIMGSOrthoManager.hpp:256
Belos::IMGSOrthoManager::getSingTol
MagnitudeType getSingTol() const
Return parameter for singular block detection.
Definition: BelosIMGSOrthoManager.hpp:259
Belos::IMGSOrthoManager::setSingTol
void setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
Definition: BelosIMGSOrthoManager.hpp:253
Belos::IMGSOrthoManager::IMGSOrthoManager
IMGSOrthoManager(const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null, const int max_ortho_steps=max_ortho_steps_default_, const MagnitudeType blk_tol=blk_tol_default_, const MagnitudeType sing_tol=sing_tol_default_)
Constructor specifying re-orthogonalization tolerance.
Definition: BelosIMGSOrthoManager.hpp:97
Belos::IMGSOrthoManager::setBlkTol
void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
Definition: BelosIMGSOrthoManager.hpp:250
Belos::IMGSOrthoManager::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: BelosIMGSOrthoManager.hpp:235
Belos::IMGSOrthoManager::projectAndNormalizeWithMxImpl
virtual int projectAndNormalizeWithMxImpl(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Given a set of bases Q[i] and a multivector X, this method computes an orthonormal basis for .
Definition: BelosIMGSOrthoManager.hpp:599
Belos::MatOrthoManager
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Definition: BelosMatOrthoManager.hpp:78
Belos::MatOrthoManager::innerProd
void innerProd(const MV &X, const MV &Y, Teuchos::SerialDenseMatrix< int, ScalarType > &Z) const
Provides the inner product defining the orthogonality concepts, using the provided operator.
Definition: BelosMatOrthoManager.hpp:115
Belos::MultiVecTraits
Traits class which defines basic operations on multivectors.
Definition: BelosMultiVecTraits.hpp:129
Belos::OperatorTraits
Class which defines basic traits for the operator type.
Definition: BelosOperatorTraits.hpp:110
Belos::getIMGSFastParameters
Teuchos::RCP< Teuchos::ParameterList > getIMGSFastParameters()
"Fast" but possibly unsafe or less accurate parameters.
Definition: BelosIMGSOrthoManager.hpp:1717
Belos::getIMGSDefaultParameters
Teuchos::RCP< Teuchos::ParameterList > getIMGSDefaultParameters()
"Default" parameters for robustness and accuracy.
Definition: BelosIMGSOrthoManager.hpp:1694
Generated on Wed Mar 9 2022 04:36:09 for Belos by doxygen 1.9.1