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