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