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

Generated for Belos by doxygen 1.8.14