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