Belos  Version of the Day
BelosIMGSOrthoManager.hpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 
47 #ifndef BELOS_IMGS_ORTHOMANAGER_HPP
48 #define BELOS_IMGS_ORTHOMANAGER_HPP
49 
57 // #define ORTHO_DEBUG
58 
59 #include "BelosConfigDefs.hpp"
60 #include "BelosMultiVecTraits.hpp"
61 #include "BelosOperatorTraits.hpp"
62 #include "BelosMatOrthoManager.hpp"
63 #include "Teuchos_SerialDenseMatrix.hpp"
64 #include "Teuchos_SerialDenseVector.hpp"
65 
66 #include "Teuchos_as.hpp"
67 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
68 #ifdef BELOS_TEUCHOS_TIME_MONITOR
69 #include "Teuchos_TimeMonitor.hpp"
70 #endif // BELOS_TEUCHOS_TIME_MONITOR
71 
72 namespace Belos {
73 
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:
489  int max_ortho_steps_;
491  MagnitudeType blk_tol_;
493  MagnitudeType sing_tol_;
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
1172  IMGSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
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
1283  IMGSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
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
1408  IMGSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
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)
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...
static void MvDot(const MV &mv, const MV &A, std::vector< ScalarType > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
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.
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 ...
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void Assign(const MV &A, MV &mv)
mv := A
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
Traits class which defines basic operations on multivectors.
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...
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
static void Apply(const OP &Op, const MV &x, MV &y, ETrans trans=NOTRANS)
Apply Op to x, putting the result into y.
Teuchos::RCP< const OP > _Op
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 ...
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
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.
Exception thrown to signal error in an orthogonalization manager method.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
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 ...
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
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 ...
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.
MagnitudeType getSingTol() const
Return parameter for singular block detection.
Class which defines basic traits for the operator type.
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.
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
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 ...
Belos&#39;s templated virtual class for providing routines for orthogonalization and orthonormzalition of...
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...
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters() const
"Fast" but possibly unsafe or less accurate parameters.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...

Generated on Thu Jul 21 2016 14:43:57 for Belos by doxygen 1.8.11