Belos Package Browser (Single Doxygen Collection)  Development
BelosIMGSOrthoManager.hpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 
47 #ifndef BELOS_IMGS_ORTHOMANAGER_HPP
48 #define BELOS_IMGS_ORTHOMANAGER_HPP
49 
57 // #define ORTHO_DEBUG
58 
59 #include "BelosConfigDefs.hpp"
60 #include "BelosMultiVecTraits.hpp"
61 #include "BelosOperatorTraits.hpp"
62 #include "BelosMatOrthoManager.hpp"
63 #include "Teuchos_SerialDenseMatrix.hpp"
64 #include "Teuchos_SerialDenseVector.hpp"
65 
66 #include "Teuchos_as.hpp"
67 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
68 #ifdef BELOS_TEUCHOS_TIME_MONITOR
69 #include "Teuchos_TimeMonitor.hpp"
70 #endif // BELOS_TEUCHOS_TIME_MONITOR
71 
72 namespace Belos {
73 
74  template<class ScalarType, class MV, class OP>
76  public MatOrthoManager<ScalarType,MV,OP>,
77  public Teuchos::ParameterListAcceptorDefaultBase
78  {
79  private:
80  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
81  typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
82  typedef Teuchos::ScalarTraits<ScalarType> SCT;
85 
86  public:
88 
89 
91  IMGSOrthoManager( const std::string& label = "Belos",
92  Teuchos::RCP<const OP> Op = Teuchos::null,
93  const int max_ortho_steps = 1,
94  const MagnitudeType blk_tol = 10*MGT::squareroot( MGT::eps() ),
95  const MagnitudeType sing_tol = 10*MGT::eps() )
96  : MatOrthoManager<ScalarType,MV,OP>(Op),
97  max_ortho_steps_( max_ortho_steps ),
98  blk_tol_( blk_tol ),
99  sing_tol_( sing_tol ),
100  label_( label )
101  {
102  std::string orthoLabel = label_ + ": Orthogonalization";
103 #ifdef BELOS_TEUCHOS_TIME_MONITOR
104  timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
105 #endif
106 
107  std::string updateLabel = label_ + ": Ortho (Update)";
108 #ifdef BELOS_TEUCHOS_TIME_MONITOR
109  timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
110 #endif
111 
112  std::string normLabel = label_ + ": Ortho (Norm)";
113 #ifdef BELOS_TEUCHOS_TIME_MONITOR
114  timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
115 #endif
116 
117  std::string ipLabel = label_ + ": Ortho (Inner Product)";
118 #ifdef BELOS_TEUCHOS_TIME_MONITOR
119  timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
120 #endif
121  }
122 
124  IMGSOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& plist,
125  const std::string& label = "Belos",
126  Teuchos::RCP<const OP> Op = Teuchos::null) :
127  MatOrthoManager<ScalarType,MV,OP>(Op),
128  max_ortho_steps_ (2),
129  blk_tol_ (10 * MGT::squareroot (MGT::eps())),
130  sing_tol_ (10 * MGT::eps()),
131  label_ (label)
132  {
133  setParameterList (plist);
134 
135  std::string orthoLabel = label_ + ": Orthogonalization";
136 #ifdef BELOS_TEUCHOS_TIME_MONITOR
137  timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
138 #endif
139  std::string updateLabel = label_ + ": Ortho (Update)";
140 #ifdef BELOS_TEUCHOS_TIME_MONITOR
141  timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
142 #endif
143  std::string normLabel = label_ + ": Ortho (Norm)";
144 #ifdef BELOS_TEUCHOS_TIME_MONITOR
145  timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
146 #endif
147  std::string ipLabel = label_ + ": Ortho (Inner Product)";
148 #ifdef BELOS_TEUCHOS_TIME_MONITOR
149  timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
150 #endif
151  }
152 
156 
158 
159  void
160  setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
161  {
162  using Teuchos::Exceptions::InvalidParameterName;
163  using Teuchos::ParameterList;
164  using Teuchos::parameterList;
165  using Teuchos::RCP;
166 
167  RCP<const ParameterList> defaultParams = getValidParameters();
168  RCP<ParameterList> params;
169  if (plist.is_null()) {
170  params = parameterList (*defaultParams);
171  } else {
172  params = plist;
173  // Some users might want to specify "blkTol" as "depTol". Due
174  // to this case, we don't invoke
175  // validateParametersAndSetDefaults on params. Instead, we go
176  // through the parameter list one parameter at a time and look
177  // for alternatives.
178  }
179 
180  // Using temporary variables and fetching all values before
181  // setting the output arguments ensures the strong exception
182  // guarantee for this function: if an exception is thrown, no
183  // externally visible side effects (in this case, setting the
184  // output arguments) have taken place.
185  int maxNumOrthogPasses;
186  MagnitudeType blkTol;
187  MagnitudeType singTol;
188 
189  try {
190  maxNumOrthogPasses = params->get<int> ("maxNumOrthogPasses");
191  } catch (InvalidParameterName&) {
192  maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
193  params->set ("maxNumOrthogPasses", maxNumOrthogPasses);
194  }
195 
196  // Handling of the "blkTol" parameter is a special case. This
197  // is because some users may prefer to call this parameter
198  // "depTol" for consistency with DGKS. However, our default
199  // parameter list calls this "blkTol", and we don't want the
200  // default list's value to override the user's value. Thus, we
201  // first check the user's parameter list for both names, and
202  // only then access the default parameter list.
203  try {
204  blkTol = params->get<MagnitudeType> ("blkTol");
205  } catch (InvalidParameterName&) {
206  try {
207  blkTol = params->get<MagnitudeType> ("depTol");
208  // "depTol" is the wrong name, so remove it and replace with
209  // "blkTol". We'll set "blkTol" below.
210  params->remove ("depTol");
211  } catch (InvalidParameterName&) {
212  blkTol = defaultParams->get<MagnitudeType> ("blkTol");
213  }
214  params->set ("blkTol", blkTol);
215  }
216 
217  try {
218  singTol = params->get<MagnitudeType> ("singTol");
219  } catch (InvalidParameterName&) {
220  singTol = defaultParams->get<MagnitudeType> ("singTol");
221  params->set ("singTol", singTol);
222  }
223 
224  max_ortho_steps_ = maxNumOrthogPasses;
225  blk_tol_ = blkTol;
226  sing_tol_ = singTol;
227 
228  setMyParamList (params);
229  }
230 
231  Teuchos::RCP<const Teuchos::ParameterList>
233  {
234  using Teuchos::as;
235  using Teuchos::ParameterList;
236  using Teuchos::parameterList;
237  using Teuchos::RCP;
238 
239  if (defaultParams_.is_null()) {
240  RCP<ParameterList> params = parameterList ("IMGS");
241 
242  // Default parameter values for IMGS orthogonalization.
243  // Documentation will be embedded in the parameter list.
244  const int defaultMaxNumOrthogPasses = 2;
245  const MagnitudeType eps = MGT::eps();
246  const MagnitudeType defaultBlkTol =
247  as<MagnitudeType> (10) * MGT::squareroot (eps);
248  const MagnitudeType defaultSingTol = as<MagnitudeType> (10) * eps;
249 
250  params->set ("maxNumOrthogPasses", defaultMaxNumOrthogPasses,
251  "Maximum number of orthogonalization passes (includes the "
252  "first). Default is 2, since \"twice is enough\" for Krylov "
253  "methods.");
254  params->set ("blkTol", defaultBlkTol, "Block reorthogonalization "
255  "threshhold.");
256  params->set ("singTol", defaultSingTol, "Singular block detection "
257  "threshold.");
258  defaultParams_ = params;
259  }
260  return defaultParams_;
261  }
263 
268  Teuchos::RCP<const Teuchos::ParameterList>
270  {
271  using Teuchos::as;
272  using Teuchos::ParameterList;
273  using Teuchos::parameterList;
274  using Teuchos::RCP;
275 
276  RCP<const ParameterList> defaultParams = getValidParameters ();
277  // Start with a clone of the default parameters.
278  RCP<ParameterList> params = parameterList (*defaultParams);
279 
280  const int maxBlkOrtho = 1; // No block reorthogonalization
281  const MagnitudeType blkTol = MGT::zero();
282  const MagnitudeType singTol = MGT::zero();
283 
284  params->set ("maxNumOrthogPasses", maxBlkOrtho);
285  params->set ("blkTol", blkTol);
286  params->set ("singTol", singTol);
287 
288  return params;
289  }
290 
292 
293 
295  void setBlkTol( const MagnitudeType blk_tol ) { blk_tol_ = blk_tol; }
296 
298  void setSingTol( const MagnitudeType sing_tol ) { sing_tol_ = sing_tol; }
299 
301  MagnitudeType getBlkTol() const { return blk_tol_; }
302 
304  MagnitudeType getSingTol() const { return sing_tol_; }
305 
307 
308 
310 
311 
339  void project ( MV &X, Teuchos::RCP<MV> MX,
340  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
341  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
342 
343 
346  void project ( MV &X,
347  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
348  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
349  project(X,Teuchos::null,C,Q);
350  }
351 
352 
353 
378  int normalize ( MV &X, Teuchos::RCP<MV> MX,
379  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B) const;
380 
381 
384  int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
385  return normalize(X,Teuchos::null,B);
386  }
387 
388  protected:
430  virtual int
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 
487  private:
494 
496  std::string label_;
497 #ifdef BELOS_TEUCHOS_TIME_MONITOR
498  Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_,
499  timerNorm_, timerScale_, timerInnerProd_;
500 #endif // BELOS_TEUCHOS_TIME_MONITOR
501 
503  mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
504 
506  int findBasis(MV &X, Teuchos::RCP<MV> MX,
507  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
508  bool completeBasis, int howMany = -1 ) const;
509 
511  bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
512  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
513  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
514 
516  bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
517  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
518  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
519 
533  int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
534  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
535  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
536  Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const;
537  };
538 
540  // Set the label for this orthogonalization manager and create new timers if it's changed
541  template<class ScalarType, class MV, class OP>
542  void IMGSOrthoManager<ScalarType,MV,OP>::setLabel(const std::string& label)
543  {
544  if (label != label_) {
545  label_ = label;
546  std::string orthoLabel = label_ + ": Orthogonalization";
547 #ifdef BELOS_TEUCHOS_TIME_MONITOR
548  timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
549 #endif
550 
551  std::string updateLabel = label_ + ": Ortho (Update)";
552 #ifdef BELOS_TEUCHOS_TIME_MONITOR
553  timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
554 #endif
555 
556  std::string normLabel = label_ + ": Ortho (Norm)";
557 #ifdef BELOS_TEUCHOS_TIME_MONITOR
558  timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
559 #endif
560 
561  std::string ipLabel = label_ + ": Ortho (Inner Product)";
562 #ifdef BELOS_TEUCHOS_TIME_MONITOR
563  timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
564 #endif
565  }
566  }
567 
569  // Compute the distance from orthonormality
570  template<class ScalarType, class MV, class OP>
571  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
572  IMGSOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const {
573  const ScalarType ONE = SCT::one();
574  int rank = MVT::GetNumberVecs(X);
575  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
577  for (int i=0; i<rank; i++) {
578  xTx(i,i) -= ONE;
579  }
580  return xTx.normFrobenius();
581  }
582 
584  // Compute the distance from orthogonality
585  template<class ScalarType, class MV, class OP>
586  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
587  IMGSOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const {
588  int r1 = MVT::GetNumberVecs(X1);
589  int r2 = MVT::GetNumberVecs(X2);
590  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
592  return xTx.normFrobenius();
593  }
594 
596  // Find an Op-orthonormal basis for span(X) - span(W)
597  template<class ScalarType, class MV, class OP>
598  int
601  Teuchos::RCP<MV> MX,
602  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
603  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
604  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
605  {
606  using Teuchos::Array;
607  using Teuchos::null;
608  using Teuchos::is_null;
609  using Teuchos::RCP;
610  using Teuchos::rcp;
611  using Teuchos::SerialDenseMatrix;
612  typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
613  typedef typename Array< RCP< const MV > >::size_type size_type;
614 
615 #ifdef BELOS_TEUCHOS_TIME_MONITOR
616  Teuchos::TimeMonitor orthotimer(*timerOrtho_);
617 #endif
618 
619  ScalarType ONE = SCT::one();
620  //ScalarType ZERO = SCT::zero();
621 
622  int nq = Q.size();
623  int xc = MVT::GetNumberVecs( X );
624  ptrdiff_t xr = MVT::GetGlobalLength( X );
625  int rank = xc;
626 
627  // If the user doesn't want to store the normalization
628  // coefficients, allocate some local memory for them. This will
629  // go away at the end of this method.
630  if (is_null (B)) {
631  B = rcp (new serial_dense_matrix_type (xc, xc));
632  }
633  // Likewise, if the user doesn't want to store the projection
634  // coefficients, allocate some local memory for them. Also make
635  // sure that all the entries of C are the right size. We're going
636  // to overwrite them anyway, so we don't have to worry about the
637  // contents (other than to resize them if they are the wrong
638  // size).
639  if (C.size() < nq)
640  C.resize (nq);
641  for (size_type k = 0; k < nq; ++k)
642  {
643  const int numRows = MVT::GetNumberVecs (*Q[k]);
644  const int numCols = xc; // Number of vectors in X
645 
646  if (is_null (C[k]))
647  C[k] = rcp (new serial_dense_matrix_type (numRows, numCols));
648  else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
649  {
650  int err = C[k]->reshape (numRows, numCols);
651  TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error,
652  "IMGS orthogonalization: failed to reshape "
653  "C[" << k << "] (the array of block "
654  "coefficients resulting from projecting X "
655  "against Q[1:" << nq << "]).");
656  }
657  }
658 
659  /****** DO NOT MODIFY *MX IF _hasOp == false ******/
660  if (this->_hasOp) {
661  if (MX == Teuchos::null) {
662  // we need to allocate space for MX
663  MX = MVT::Clone(X,MVT::GetNumberVecs(X));
664  OPT::Apply(*(this->_Op),X,*MX);
665  }
666  }
667  else {
668  // Op == I --> MX = X (ignore it if the user passed it in)
669  MX = Teuchos::rcp( &X, false );
670  }
671 
672  int mxc = MVT::GetNumberVecs( *MX );
673  ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
674 
675  // short-circuit
676  TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Belos::IMGSOrthoManager::projectAndNormalize(): X must be non-empty" );
677 
678  int numbas = 0;
679  for (int i=0; i<nq; i++) {
680  numbas += MVT::GetNumberVecs( *Q[i] );
681  }
682 
683  // check size of B
684  TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
685  "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
686  // check size of X and MX
687  TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
688  "Belos::IMGSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
689  // check size of X w.r.t. MX
690  TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
691  "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
692  // check feasibility
693  //TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument,
694  // "Belos::IMGSOrthoManager::projectAndNormalize(): Orthogonality constraints not feasible" );
695 
696  // Some flags for checking dependency returns from the internal orthogonalization methods
697  bool dep_flg = false;
698 
699  // Make a temporary copy of X and MX, just in case a block dependency is detected.
700  Teuchos::RCP<MV> tmpX, tmpMX;
701  tmpX = MVT::CloneCopy(X);
702  if (this->_hasOp) {
703  tmpMX = MVT::CloneCopy(*MX);
704  }
705 
706  if (xc == 1) {
707 
708  // Use the cheaper block orthogonalization.
709  // NOTE: Don't check for dependencies because the update has one vector.
710  dep_flg = blkOrtho1( X, MX, C, Q );
711 
712  // Normalize the new block X
713  {
714 #ifdef BELOS_TEUCHOS_TIME_MONITOR
715  Teuchos::TimeMonitor normTimer( *timerNorm_ );
716 #endif
717  if ( B == Teuchos::null ) {
718  B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
719  }
720  std::vector<ScalarType> diag(xc);
721  MVT::MvDot( X, *MX, diag );
722  (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
723  rank = 1;
724  MVT::MvScale( X, ONE/(*B)(0,0) );
725  if (this->_hasOp) {
726  // Update MXj.
727  MVT::MvScale( *MX, ONE/(*B)(0,0) );
728  }
729  }
730 
731  }
732  else {
733 
734  // Use the cheaper block orthogonalization.
735  dep_flg = blkOrtho( X, MX, C, Q );
736 
737  // If a dependency has been detected in this block, then perform
738  // the more expensive nonblock (single vector at a time)
739  // orthogonalization.
740  if (dep_flg) {
741  rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
742 
743  // Copy tmpX back into X.
744  MVT::Assign( *tmpX, X );
745  if (this->_hasOp) {
746  MVT::Assign( *tmpMX, *MX );
747  }
748  }
749  else {
750  // There is no dependency, so orthonormalize new block X
751  rank = findBasis( X, MX, B, false );
752  if (rank < xc) {
753  // A dependency was found during orthonormalization of X,
754  // rerun orthogonalization using more expensive single-
755  // vector orthogonalization.
756  rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
757 
758  // Copy tmpX back into X.
759  MVT::Assign( *tmpX, X );
760  if (this->_hasOp) {
761  MVT::Assign( *tmpMX, *MX );
762  }
763  }
764  }
765  } // if (xc == 1) {
766 
767  // this should not raise an std::exception; but our post-conditions oblige us to check
768  TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
769  "Belos::IMGSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
770 
771  // Return the rank of X.
772  return rank;
773  }
774 
775 
776 
778  // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
779  template<class ScalarType, class MV, class OP>
781  MV &X, Teuchos::RCP<MV> MX,
782  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
783 
784 #ifdef BELOS_TEUCHOS_TIME_MONITOR
785  Teuchos::TimeMonitor orthotimer(*timerOrtho_);
786 #endif
787 
788  // call findBasis, with the instruction to try to generate a basis of rank numvecs(X)
789  return findBasis(X, MX, B, true);
790  }
791 
792 
793 
795  template<class ScalarType, class MV, class OP>
797  MV &X, Teuchos::RCP<MV> MX,
798  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
799  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
800  // For the inner product defined by the operator Op or the identity (Op == 0)
801  // -> Orthogonalize X against each Q[i]
802  // Modify MX accordingly
803  //
804  // Note that when Op is 0, MX is not referenced
805  //
806  // Parameter variables
807  //
808  // X : Vectors to be transformed
809  //
810  // MX : Image of the block of vectors X by the mass matrix
811  //
812  // Q : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
813  //
814 
815 #ifdef BELOS_TEUCHOS_TIME_MONITOR
816  Teuchos::TimeMonitor orthotimer(*timerOrtho_);
817 #endif
818 
819  int xc = MVT::GetNumberVecs( X );
820  ptrdiff_t xr = MVT::GetGlobalLength( X );
821  int nq = Q.size();
822  std::vector<int> qcs(nq);
823  // short-circuit
824  if (nq == 0 || xc == 0 || xr == 0) {
825  return;
826  }
827  ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
828  // if we don't have enough C, expand it with null references
829  // if we have too many, resize to throw away the latter ones
830  // if we have exactly as many as we have Q, this call has no effect
831  C.resize(nq);
832 
833 
834  /****** DO NOT MODIFY *MX IF _hasOp == false ******/
835  if (this->_hasOp) {
836  if (MX == Teuchos::null) {
837  // we need to allocate space for MX
838  MX = MVT::Clone(X,MVT::GetNumberVecs(X));
839  OPT::Apply(*(this->_Op),X,*MX);
840  }
841  }
842  else {
843  // Op == I --> MX = X (ignore it if the user passed it in)
844  MX = Teuchos::rcp( &X, false );
845  }
846  int mxc = MVT::GetNumberVecs( *MX );
847  ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
848 
849  // check size of X and Q w.r.t. common sense
850  TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
851  "Belos::IMGSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
852  // check size of X w.r.t. MX and Q
853  TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
854  "Belos::IMGSOrthoManager::project(): Size of X not consistant with MX,Q" );
855 
856  // tally up size of all Q and check/allocate C
857  int baslen = 0;
858  for (int i=0; i<nq; i++) {
859  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
860  "Belos::IMGSOrthoManager::project(): Q lengths not mutually consistant" );
861  qcs[i] = MVT::GetNumberVecs( *Q[i] );
862  TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
863  "Belos::IMGSOrthoManager::project(): Q has less rows than columns" );
864  baslen += qcs[i];
865 
866  // check size of C[i]
867  if ( C[i] == Teuchos::null ) {
868  C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
869  }
870  else {
871  TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
872  "Belos::IMGSOrthoManager::project(): Size of Q not consistant with size of C" );
873  }
874  }
875 
876  // Use the cheaper block orthogonalization, don't check for rank deficiency.
877  blkOrtho( X, MX, C, Q );
878 
879  }
880 
882  // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that
883  // the rank is numvectors(X)
884  template<class ScalarType, class MV, class OP>
886  MV &X, Teuchos::RCP<MV> MX,
887  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
888  bool completeBasis, int howMany ) const {
889  // For the inner product defined by the operator Op or the identity (Op == 0)
890  // -> Orthonormalize X
891  // Modify MX accordingly
892  //
893  // Note that when Op is 0, MX is not referenced
894  //
895  // Parameter variables
896  //
897  // X : Vectors to be orthonormalized
898  //
899  // MX : Image of the multivector X under the operator Op
900  //
901  // Op : Pointer to the operator for the inner product
902  //
903  //
904 
905 #ifdef BELOS_TEUCHOS_TIME_MONITOR
906  Teuchos::TimeMonitor normTimer( *timerNorm_ );
907 #endif
908 
909  const ScalarType ONE = SCT::one();
910  const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
911 
912  int xc = MVT::GetNumberVecs( X );
913  ptrdiff_t xr = MVT::GetGlobalLength( X );
914 
915  if (howMany == -1) {
916  howMany = xc;
917  }
918 
919  /*******************************************************
920  * If _hasOp == false, we will not reference MX below *
921  *******************************************************/
922 
923  // if Op==null, MX == X (via pointer)
924  // Otherwise, either the user passed in MX or we will allocated and compute it
925  if (this->_hasOp) {
926  if (MX == Teuchos::null) {
927  // we need to allocate space for MX
928  MX = MVT::Clone(X,xc);
929  OPT::Apply(*(this->_Op),X,*MX);
930  }
931  }
932 
933  /* if the user doesn't want to store the coefficienets,
934  * allocate some local memory for them
935  */
936  if ( B == Teuchos::null ) {
937  B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
938  }
939 
940  int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
941  ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
942 
943  // check size of C, B
944  TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
945  "Belos::IMGSOrthoManager::findBasis(): X must be non-empty" );
946  TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
947  "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
948  TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
949  "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
950  TEUCHOS_TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument,
951  "Belos::IMGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
952  TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
953  "Belos::IMGSOrthoManager::findBasis(): Invalid howMany parameter" );
954 
955  /* xstart is which column we are starting the process with, based on howMany
956  * columns before xstart are assumed to be Op-orthonormal already
957  */
958  int xstart = xc - howMany;
959 
960  for (int j = xstart; j < xc; j++) {
961 
962  // numX is
963  // * number of currently orthonormal columns of X
964  // * the index of the current column of X
965  int numX = j;
966  bool addVec = false;
967 
968  // Get a view of the vector currently being worked on.
969  std::vector<int> index(1);
970  index[0] = numX;
971  Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
972  Teuchos::RCP<MV> MXj;
973  if ((this->_hasOp)) {
974  // MXj is a view of the current vector in MX
975  MXj = MVT::CloneViewNonConst( *MX, index );
976  }
977  else {
978  // MXj is a pointer to Xj, and MUST NOT be modified
979  MXj = Xj;
980  }
981 
982  // Make storage for these Gram-Schmidt iterations.
983  Teuchos::SerialDenseVector<int,ScalarType> product(numX);
984  Teuchos::SerialDenseVector<int,ScalarType> P2(1);
985  Teuchos::RCP<const MV> prevX, prevMX;
986 
987  std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
988  //
989  // Save old MXj vector and compute Op-norm
990  //
991  Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj );
992  MVT::MvDot( *Xj, *MXj, oldDot );
993  // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
994  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError,
995  "Belos::IMGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
996 
997  // Perform MGS one vector at a time
998  for (int ii=0; ii<numX; ii++) {
999 
1000  index[0] = ii;
1001  prevX = MVT::CloneView( X, index );
1002  if (this->_hasOp) {
1003  prevMX = MVT::CloneView( *MX, index );
1004  }
1005 
1006  for (int i=0; i<max_ortho_steps_; ++i) {
1007 
1008  // product <- prevX^T MXj
1009  {
1010 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1011  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1012 #endif
1014  }
1015 
1016  // Xj <- Xj - prevX prevX^T MXj
1017  // = Xj - prevX product
1018  {
1019 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1020  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1021 #endif
1022  MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1023  }
1024 
1025  // Update MXj
1026  if (this->_hasOp) {
1027  // MXj <- Op*Xj_new
1028  // = Op*(Xj_old - prevX prevX^T MXj)
1029  // = MXj - prevMX product
1030 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1031  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1032 #endif
1033  MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1034  }
1035 
1036  // Set coefficients
1037  if ( i==0 )
1038  product[ii] = P2[0];
1039  else
1040  product[ii] += P2[0];
1041 
1042  } // for (int i=0; i<max_ortho_steps_; ++i)
1043 
1044  } // for (int ii=0; ii<numX; ++ii)
1045 
1046  // Compute Op-norm with old MXj
1047  MVT::MvDot( *Xj, *oldMXj, newDot );
1048 
1049  // Check to see if the new vector is dependent.
1050  if (completeBasis) {
1051  //
1052  // We need a complete basis, so add random vectors if necessary
1053  //
1054  if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1055 
1056  // Add a random vector and orthogonalize it against previous vectors in block.
1057  addVec = true;
1058 #ifdef ORTHO_DEBUG
1059  std::cout << "Belos::IMGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1060 #endif
1061  //
1062  Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1063  Teuchos::RCP<MV> tempMXj;
1064  MVT::MvRandom( *tempXj );
1065  if (this->_hasOp) {
1066  tempMXj = MVT::Clone( X, 1 );
1067  OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1068  }
1069  else {
1070  tempMXj = tempXj;
1071  }
1072  MVT::MvDot( *tempXj, *tempMXj, oldDot );
1073  //
1074  // Perform MGS one vector at a time
1075  for (int ii=0; ii<numX; ii++) {
1076 
1077  index[0] = ii;
1078  prevX = MVT::CloneView( X, index );
1079  if (this->_hasOp) {
1080  prevMX = MVT::CloneView( *MX, index );
1081  }
1082 
1083  for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
1084  {
1085 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1086  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1087 #endif
1088  MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*tempXj,tempMXj,P2);
1089  }
1090  {
1091 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1092  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1093 #endif
1094  MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *tempXj );
1095  }
1096  if (this->_hasOp) {
1097 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1098  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1099 #endif
1100  MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *tempMXj );
1101  }
1102 
1103  // Set coefficients
1104  if ( num_orth==0 )
1105  product[ii] = P2[0];
1106  else
1107  product[ii] += P2[0];
1108  }
1109  }
1110 
1111  // Compute new Op-norm
1112  MVT::MvDot( *tempXj, *tempMXj, newDot );
1113  //
1114  if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1115  // Copy vector into current column of _basisvecs
1116  MVT::Assign( *tempXj, *Xj );
1117  if (this->_hasOp) {
1118  MVT::Assign( *tempMXj, *MXj );
1119  }
1120  }
1121  else {
1122  return numX;
1123  }
1124  }
1125  }
1126  else {
1127  //
1128  // We only need to detect dependencies.
1129  //
1130  if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1131  return numX;
1132  }
1133  }
1134 
1135 
1136  // If we haven't left this method yet, then we can normalize the new vector Xj.
1137  // Normalize Xj.
1138  // Xj <- Xj / std::sqrt(newDot)
1139  ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1140  {
1141  MVT::MvScale( *Xj, ONE/diag );
1142  if (this->_hasOp) {
1143  // Update MXj.
1144  MVT::MvScale( *MXj, ONE/diag );
1145  }
1146  }
1147 
1148  // If we've added a random vector, enter a zero in the j'th diagonal element.
1149  if (addVec) {
1150  (*B)(j,j) = ZERO;
1151  }
1152  else {
1153  (*B)(j,j) = diag;
1154  }
1155 
1156  // Save the coefficients, if we are working on the original vector and not a randomly generated one
1157  if (!addVec) {
1158  for (int i=0; i<numX; i++) {
1159  (*B)(i,j) = product(i);
1160  }
1161  }
1162 
1163  } // for (j = 0; j < xc; ++j)
1164 
1165  return xc;
1166  }
1167 
1169  // Routine to compute the block orthogonalization
1170  template<class ScalarType, class MV, class OP>
1171  bool
1173  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1174  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
1175  {
1176  int nq = Q.size();
1177  int xc = MVT::GetNumberVecs( X );
1178  const ScalarType ONE = SCT::one();
1179 
1180  std::vector<int> qcs( nq );
1181  for (int i=0; i<nq; i++) {
1182  qcs[i] = MVT::GetNumberVecs( *Q[i] );
1183  }
1184 
1185  // Perform the Gram-Schmidt transformation for a block of vectors
1186  std::vector<int> index(1);
1187  Teuchos::RCP<const MV> tempQ;
1188 
1189  Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1190  // Define the product Q^T * (Op*X)
1191  for (int i=0; i<nq; i++) {
1192 
1193  // Perform MGS one vector at a time
1194  for (int ii=0; ii<qcs[i]; ii++) {
1195 
1196  index[0] = ii;
1197  tempQ = MVT::CloneView( *Q[i], index );
1198  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1199 
1200  // Multiply Q' with MX
1201  {
1202 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1203  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1204 #endif
1206  }
1207  // Multiply by Q and subtract the result in X
1208  {
1209 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1210  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1211 #endif
1212  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1213  }
1214  }
1215 
1216  // Update MX, with the least number of applications of Op as possible
1217  if (this->_hasOp) {
1218  if (xc <= qcs[i]) {
1219  OPT::Apply( *(this->_Op), X, *MX);
1220  }
1221  else {
1222  // this will possibly be used again below; don't delete it
1223  MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1224  OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1225  MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1226  }
1227  }
1228  }
1229 
1230  // Do as many steps of modified Gram-Schmidt as required by max_ortho_steps_
1231  for (int j = 1; j < max_ortho_steps_; ++j) {
1232 
1233  for (int i=0; i<nq; i++) {
1234 
1235  Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],1);
1236 
1237  // Perform MGS one vector at a time
1238  for (int ii=0; ii<qcs[i]; ii++) {
1239 
1240  index[0] = ii;
1241  tempQ = MVT::CloneView( *Q[i], index );
1242  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1243  Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii, 0 );
1244 
1245  // Apply another step of modified Gram-Schmidt
1246  {
1247 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1248  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1249 #endif
1250  MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC2 );
1251  }
1252  tempC += tempC2;
1253  {
1254 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1255  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1256 #endif
1257  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1258  }
1259 
1260  }
1261 
1262  // Update MX, with the least number of applications of Op as possible
1263  if (this->_hasOp) {
1264  if (MQ[i].get()) {
1265  // MQ was allocated and computed above; use it
1266  MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1267  }
1268  else if (xc <= qcs[i]) {
1269  // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1270  OPT::Apply( *(this->_Op), X, *MX);
1271  }
1272  }
1273  } // for (int i=0; i<nq; i++)
1274  } // for (int j = 0; j < max_ortho_steps; ++j)
1275 
1276  return false;
1277  }
1278 
1280  // Routine to compute the block orthogonalization
1281  template<class ScalarType, class MV, class OP>
1282  bool
1284  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1285  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
1286  {
1287  int nq = Q.size();
1288  int xc = MVT::GetNumberVecs( X );
1289  bool dep_flg = false;
1290  const ScalarType ONE = SCT::one();
1291 
1292  std::vector<int> qcs( nq );
1293  for (int i=0; i<nq; i++) {
1294  qcs[i] = MVT::GetNumberVecs( *Q[i] );
1295  }
1296 
1297  // Perform the Gram-Schmidt transformation for a block of vectors
1298 
1299  // Compute the initial Op-norms
1300  std::vector<ScalarType> oldDot( xc );
1301  MVT::MvDot( X, *MX, oldDot );
1302 
1303  std::vector<int> index(1);
1304  Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1305  Teuchos::RCP<const MV> tempQ;
1306 
1307  // Define the product Q^T * (Op*X)
1308  for (int i=0; i<nq; i++) {
1309 
1310  // Perform MGS one vector at a time
1311  for (int ii=0; ii<qcs[i]; ii++) {
1312 
1313  index[0] = ii;
1314  tempQ = MVT::CloneView( *Q[i], index );
1315  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1316 
1317  // Multiply Q' with MX
1318  {
1319 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1320  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1321 #endif
1322  MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC);
1323  }
1324  // Multiply by Q and subtract the result in X
1325  {
1326 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1327  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1328 #endif
1329  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1330  }
1331  }
1332 
1333  // Update MX, with the least number of applications of Op as possible
1334  if (this->_hasOp) {
1335  if (xc <= qcs[i]) {
1336  OPT::Apply( *(this->_Op), X, *MX);
1337  }
1338  else {
1339  // this will possibly be used again below; don't delete it
1340  MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1341  OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1342  MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1343  }
1344  }
1345  }
1346 
1347  // Do as many steps of modified Gram-Schmidt as required by max_ortho_steps_
1348  for (int j = 1; j < max_ortho_steps_; ++j) {
1349 
1350  for (int i=0; i<nq; i++) {
1351  Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],xc);
1352 
1353  // Perform MGS one vector at a time
1354  for (int ii=0; ii<qcs[i]; ii++) {
1355 
1356  index[0] = ii;
1357  tempQ = MVT::CloneView( *Q[i], index );
1358  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1359  Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, xc, ii, 0 );
1360 
1361  // Apply another step of modified Gram-Schmidt
1362  {
1363 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1364  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1365 #endif
1366  MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC2 );
1367  }
1368  tempC += tempC2;
1369  {
1370 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1371  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1372 #endif
1373  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1374  }
1375  }
1376 
1377  // Update MX, with the least number of applications of Op as possible
1378  if (this->_hasOp) {
1379  if (MQ[i].get()) {
1380  // MQ was allocated and computed above; use it
1381  MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1382  }
1383  else if (xc <= qcs[i]) {
1384  // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1385  OPT::Apply( *(this->_Op), X, *MX);
1386  }
1387  }
1388  } // for (int i=0; i<nq; i++)
1389  } // for (int j = 0; j < max_ortho_steps; ++j)
1390 
1391  // Compute new Op-norms
1392  std::vector<ScalarType> newDot(xc);
1393  MVT::MvDot( X, *MX, newDot );
1394 
1395  // Check to make sure the new block of vectors are not dependent on previous vectors
1396  for (int i=0; i<xc; i++){
1397  if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1398  dep_flg = true;
1399  break;
1400  }
1401  } // end for (i=0;...)
1402 
1403  return dep_flg;
1404  }
1405 
1406  template<class ScalarType, class MV, class OP>
1407  int
1409  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1410  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1411  Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const
1412  {
1413  Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1414 
1415  const ScalarType ONE = SCT::one();
1416  const ScalarType ZERO = SCT::zero();
1417 
1418  int nq = Q.size();
1419  int xc = MVT::GetNumberVecs( X );
1420  std::vector<int> indX( 1 );
1421  std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1422 
1423  std::vector<int> qcs( nq );
1424  for (int i=0; i<nq; i++) {
1425  qcs[i] = MVT::GetNumberVecs( *Q[i] );
1426  }
1427 
1428  // Create pointers for the previous vectors of X that have already been orthonormalized.
1429  Teuchos::RCP<const MV> lastQ;
1430  Teuchos::RCP<MV> Xj, MXj;
1431  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1432 
1433  // Perform the Gram-Schmidt transformation for each vector in the block of vectors.
1434  for (int j=0; j<xc; j++) {
1435 
1436  bool dep_flg = false;
1437 
1438  // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
1439  if (j > 0) {
1440  std::vector<int> index( j );
1441  for (int ind=0; ind<j; ind++) {
1442  index[ind] = ind;
1443  }
1444  lastQ = MVT::CloneView( X, index );
1445 
1446  // Add these views to the Q and C arrays.
1447  Q.push_back( lastQ );
1448  C.push_back( B );
1449  qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1450  }
1451 
1452  // Get a view of the current vector in X to orthogonalize.
1453  indX[0] = j;
1454  Xj = MVT::CloneViewNonConst( X, indX );
1455  if (this->_hasOp) {
1456  MXj = MVT::CloneViewNonConst( *MX, indX );
1457  }
1458  else {
1459  MXj = Xj;
1460  }
1461 
1462  // Compute the initial Op-norms
1463  MVT::MvDot( *Xj, *MXj, oldDot );
1464 
1465  Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1466  Teuchos::RCP<const MV> tempQ;
1467 
1468  // Define the product Q^T * (Op*X)
1469  for (int i=0; i<Q.size(); i++) {
1470 
1471  // Perform MGS one vector at a time
1472  for (int ii=0; ii<qcs[i]; ii++) {
1473 
1474  indX[0] = ii;
1475  tempQ = MVT::CloneView( *Q[i], indX );
1476  // Get a view of the current serial dense matrix
1477  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, j );
1478 
1479  // Multiply Q' with MX
1480  MatOrthoManager<ScalarType,MV,OP>::innerProd(*tempQ,*Xj,MXj,tempC);
1481 
1482  // Multiply by Q and subtract the result in Xj
1483  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *Xj );
1484  }
1485 
1486  // Update MXj, with the least number of applications of Op as possible
1487  if (this->_hasOp) {
1488  if (xc <= qcs[i]) {
1489  OPT::Apply( *(this->_Op), *Xj, *MXj);
1490  }
1491  else {
1492  // this will possibly be used again below; don't delete it
1493  MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1494  OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1495  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1496  MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1497  }
1498  }
1499  }
1500 
1501  // Do any additional steps of modified Gram-Schmidt orthogonalization
1502  for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
1503 
1504  for (int i=0; i<Q.size(); i++) {
1505  Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1506 
1507  // Perform MGS one vector at a time
1508  for (int ii=0; ii<qcs[i]; ii++) {
1509 
1510  indX[0] = ii;
1511  tempQ = MVT::CloneView( *Q[i], indX );
1512  // Get a view of the current serial dense matrix
1513  Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii );
1514 
1515  // Apply another step of modified Gram-Schmidt
1516  MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, *Xj, MXj, tempC2);
1517  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, *Xj );
1518  }
1519 
1520  // Add the coefficients into C[i]
1521  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1522  tempC += C2;
1523 
1524  // Update MXj, with the least number of applications of Op as possible
1525  if (this->_hasOp) {
1526  if (MQ[i].get()) {
1527  // MQ was allocated and computed above; use it
1528  MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1529  }
1530  else if (xc <= qcs[i]) {
1531  // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1532  OPT::Apply( *(this->_Op), *Xj, *MXj);
1533  }
1534  }
1535  } // for (int i=0; i<Q.size(); i++)
1536 
1537  } // for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps)
1538 
1539  // Check for linear dependence.
1540  if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1541  dep_flg = true;
1542  }
1543 
1544  // Normalize the new vector if it's not dependent
1545  if (!dep_flg) {
1546  ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1547 
1548  MVT::MvScale( *Xj, ONE/diag );
1549  if (this->_hasOp) {
1550  // Update MXj.
1551  MVT::MvScale( *MXj, ONE/diag );
1552  }
1553 
1554  // Enter value on diagonal of B.
1555  (*B)(j,j) = diag;
1556  }
1557  else {
1558  // Create a random vector and orthogonalize it against all previous columns of Q.
1559  Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1560  Teuchos::RCP<MV> tempMXj;
1561  MVT::MvRandom( *tempXj );
1562  if (this->_hasOp) {
1563  tempMXj = MVT::Clone( X, 1 );
1564  OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1565  }
1566  else {
1567  tempMXj = tempXj;
1568  }
1569  MVT::MvDot( *tempXj, *tempMXj, oldDot );
1570  //
1571  for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
1572 
1573  for (int i=0; i<Q.size(); i++) {
1574  Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1575 
1576  // Perform MGS one vector at a time
1577  for (int ii=0; ii<qcs[i]; ii++) {
1578 
1579  indX[0] = ii;
1580  tempQ = MVT::CloneView( *Q[i], indX );
1581  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, product, 1, 1, ii );
1582 
1583  // Apply another step of modified Gram-Schmidt
1584  MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, *tempXj, tempMXj, tempC );
1585  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *tempXj );
1586 
1587  }
1588 
1589  // Update MXj, with the least number of applications of Op as possible
1590  if (this->_hasOp) {
1591  if (MQ[i].get()) {
1592  // MQ was allocated and computed above; use it
1593  MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1594  }
1595  else if (xc <= qcs[i]) {
1596  // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1597  OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1598  }
1599  }
1600  } // for (int i=0; i<nq; i++)
1601  } // for (int num_orth=0; num_orth<max_orth_steps_; num_orth++)
1602 
1603  // Compute the Op-norms after the correction step.
1604  MVT::MvDot( *tempXj, *tempMXj, newDot );
1605 
1606  // Copy vector into current column of Xj
1607  if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1608  ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1609 
1610  // Enter value on diagonal of B.
1611  (*B)(j,j) = ZERO;
1612 
1613  // Copy vector into current column of _basisvecs
1614  MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1615  if (this->_hasOp) {
1616  MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1617  }
1618  }
1619  else {
1620  return j;
1621  }
1622  } // if (!dep_flg)
1623 
1624  // Remove the vectors from array
1625  if (j > 0) {
1626  Q.resize( nq );
1627  C.resize( nq );
1628  qcs.resize( nq );
1629  }
1630 
1631  } // for (int j=0; j<xc; j++)
1632 
1633  return xc;
1634  }
1635 
1636 } // namespace Belos
1637 
1638 #endif // BELOS_IMGS_ORTHOMANAGER_HPP
1639 
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
MagnitudeType blk_tol_
Block reorthogonalization tolerance.
void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
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 ...
int blkOrthoSing(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 > > QQ) const
Project X against QQ and normalize X, one vector at a time.
Declaration of basic traits for the multivector type.
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...
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 ...
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 ...
Class which defines basic traits for the operator type.
Traits class which defines basic operations on multivectors.
int findBasis(MV &X, Teuchos::RCP< MV > MX, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > C, bool completeBasis, int howMany=-1) const
Routine to find an orthonormal basis for X.
std::string label_
Label for timers.
MagnitudeType getSingTol() const
Return parameter for singular block detection.
IMGSOrthoManager(const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null, const int max_ortho_steps=1, const MagnitudeType blk_tol=10 *MGT::squareroot(MGT::eps()), const MagnitudeType sing_tol=10 *MGT::eps())
Constructor specifying re-orthogonalization tolerance.
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::ScalarTraits< MagnitudeType > MGT
bool blkOrtho1(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Routine to compute the block orthogonalization.
Exception thrown to signal error in an orthogonalization manager method.
MagnitudeType sing_tol_
Singular block detection threshold.
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.
int max_ortho_steps_
Max number of (re)orthogonalization steps, including the first.
bool blkOrtho(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Routine to compute the block orthogonalization.
MultiVecTraits< ScalarType, MV > MVT
void setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters() const
"Fast" but possibly unsafe or less accurate parameters.
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::RCP< Teuchos::ParameterList > defaultParams_
Default parameter list.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
IMGSOrthoManager(const Teuchos::RCP< Teuchos::ParameterList > &plist, const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null)
Constructor that takes a list of parameters.
Class which defines basic traits for the operator type.
Belos&#39;s templated virtual class for providing routines for orthogonalization and orthonormzalition of...
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 ...
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...
void setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
Belos header file which uses auto-configuration information to include necessary C++ headers...
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.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.