Belos  Version of the Day
BelosDGKSOrthoManager.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_DGKS_ORTHOMANAGER_HPP
48 #define BELOS_DGKS_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  DGKSOrthoManager( const std::string& label = "Belos",
90  Teuchos::RCP<const OP> Op = Teuchos::null,
91  const int max_blk_ortho = 2,
92  const MagnitudeType blk_tol = 10*MGT::squareroot( MGT::eps() ),
93  const MagnitudeType dep_tol = MGT::one()/MGT::squareroot( 2*MGT::one() ),
94  const MagnitudeType sing_tol = 10*MGT::eps() )
95  : MatOrthoManager<ScalarType,MV,OP>(Op),
96  max_blk_ortho_( max_blk_ortho ),
97  blk_tol_( blk_tol ),
98  dep_tol_( dep_tol ),
99  sing_tol_( sing_tol ),
100  label_( label )
101  {
102 #ifdef BELOS_TEUCHOS_TIME_MONITOR
103  std::string orthoLabel = label_ + ": Orthogonalization";
104  timerOrtho_ = Teuchos::TimeMonitor::getNewCounter( orthoLabel );
105 #endif
106  }
107 
109  DGKSOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& plist,
110  const std::string& label = "Belos",
111  Teuchos::RCP<const OP> Op = Teuchos::null)
112  : MatOrthoManager<ScalarType,MV,OP>(Op),
113  max_blk_ortho_ (2),
114  blk_tol_ (Teuchos::as<MagnitudeType>(10) * MGT::squareroot(MGT::eps())),
115  dep_tol_ (MGT::one() / MGT::squareroot (Teuchos::as<MagnitudeType>(2))),
116  sing_tol_ (Teuchos::as<MagnitudeType>(10) * MGT::eps()),
117  label_( label )
118  {
119  setParameterList (plist);
120 
121 #ifdef BELOS_TEUCHOS_TIME_MONITOR
122  std::string orthoLabel = label_ + ": Orthogonalization";
123  timerOrtho_ = Teuchos::TimeMonitor::getNewCounter( orthoLabel );
124 #endif
125  }
126 
130 
132 
133 
134  void
135  setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
136  {
137  using Teuchos::ParameterList;
138  using Teuchos::parameterList;
139  using Teuchos::RCP;
140 
141  RCP<const ParameterList> defaultParams = getValidParameters();
142  RCP<ParameterList> params;
143  if (plist.is_null()) {
144  // No need to validate default parameters.
145  params = parameterList (*defaultParams);
146  } else {
147  params = plist;
148  params->validateParametersAndSetDefaults (*defaultParams);
149  }
150 
151  // Using temporary variables and fetching all values before
152  // setting the output arguments ensures the strong exception
153  // guarantee for this function: if an exception is thrown, no
154  // externally visible side effects (in this case, setting the
155  // output arguments) have taken place.
156  const int maxNumOrthogPasses = params->get<int> ("maxNumOrthogPasses");
157  const MagnitudeType blkTol = params->get<MagnitudeType> ("blkTol");
158  const MagnitudeType depTol = params->get<MagnitudeType> ("depTol");
159  const MagnitudeType singTol = params->get<MagnitudeType> ("singTol");
160 
161  max_blk_ortho_ = maxNumOrthogPasses;
162  blk_tol_ = blkTol;
163  dep_tol_ = depTol;
164  sing_tol_ = singTol;
165 
166  setMyParamList (params);
167  }
168 
169  Teuchos::RCP<const Teuchos::ParameterList>
171  {
172  using Teuchos::as;
173  using Teuchos::ParameterList;
174  using Teuchos::parameterList;
175  using Teuchos::RCP;
176 
177  if (defaultParams_.is_null()) {
178  RCP<ParameterList> params = parameterList ("DGKS");
179  const MagnitudeType eps = MGT::eps ();
180 
181  // Default parameter values for DGKS orthogonalization.
182  // Documentation will be embedded in the parameter list.
183  const int defaultMaxNumOrthogPasses = 2;
184  const MagnitudeType defaultBlkTol =
185  as<MagnitudeType> (10) * MGT::squareroot (eps);
186  const MagnitudeType defaultDepTol =
187  MGT::one() / MGT::squareroot (as<MagnitudeType> (2));
188 
189  const MagnitudeType defaultSingTol = as<MagnitudeType> (10) * eps;
190 
191  params->set ("maxNumOrthogPasses", defaultMaxNumOrthogPasses,
192  "Maximum number of orthogonalization passes (includes the "
193  "first). Default is 2, since \"twice is enough\" for Krylov "
194  "methods.");
195  params->set ("blkTol", defaultBlkTol, "Block reorthogonalization "
196  "threshhold.");
197  params->set ("depTol", defaultDepTol,
198  "(Non-block) reorthogonalization threshold.");
199  params->set ("singTol", defaultSingTol, "Singular block detection "
200  "threshold.");
201  defaultParams_ = params;
202  }
203  return defaultParams_;
204  }
205 
207 
208  Teuchos::RCP<const Teuchos::ParameterList>
210  {
211  using Teuchos::ParameterList;
212  using Teuchos::RCP;
213  using Teuchos::rcp;
214 
215  RCP<const ParameterList> defaultParams = getValidParameters ();
216  // Start with a clone of the default parameters.
217  RCP<ParameterList> params = rcp (new ParameterList (*defaultParams));
218 
219  const int maxBlkOrtho = 1;
220  const MagnitudeType blkTol = MGT::zero();
221  const MagnitudeType depTol = MGT::zero();
222  const MagnitudeType singTol = MGT::zero();
223 
224  params->set ("maxNumOrthogPasses", maxBlkOrtho);
225  params->set ("blkTol", blkTol);
226  params->set ("depTol", depTol);
227  params->set ("singTol", singTol);
228 
229  return params;
230  }
231 
233 
234 
236  void setBlkTol( const MagnitudeType blk_tol ) {
237  // Update the parameter list as well.
238  Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
239  if (! params.is_null()) {
240  // If it's null, then we haven't called setParameterList()
241  // yet. It's entirely possible to construct the parameter
242  // list on demand, so we don't try to create the parameter
243  // list here.
244  params->set ("blkTol", blk_tol);
245  }
246  blk_tol_ = blk_tol;
247  }
248 
250  void setDepTol( const MagnitudeType dep_tol ) {
251  // Update the parameter list as well.
252  Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
253  if (! params.is_null()) {
254  params->set ("depTol", dep_tol);
255  }
256  dep_tol_ = dep_tol;
257  }
258 
260  void setSingTol( const MagnitudeType sing_tol ) {
261  // Update the parameter list as well.
262  Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
263  if (! params.is_null()) {
264  params->set ("singTol", sing_tol);
265  }
266  sing_tol_ = sing_tol;
267  }
268 
270  MagnitudeType getBlkTol() const { return blk_tol_; }
271 
273  MagnitudeType getDepTol() const { return dep_tol_; }
274 
276  MagnitudeType getSingTol() const { return sing_tol_; }
277 
279 
280 
282 
283 
311  void project ( MV &X, Teuchos::RCP<MV> MX,
312  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
313  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
314 
315 
318  void project ( MV &X,
319  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
320  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
321  project(X,Teuchos::null,C,Q);
322  }
323 
324 
325 
350  int normalize ( MV &X, Teuchos::RCP<MV> MX,
351  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B) const;
352 
353 
356  int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
357  return normalize(X,Teuchos::null,B);
358  }
359 
360  protected:
361 
417  virtual int
419  Teuchos::RCP<MV> MX,
420  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
421  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
422  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
423 
424  public:
426 
428 
434  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
435  orthonormError(const MV &X) const {
436  return orthonormError(X,Teuchos::null);
437  }
438 
445  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
446  orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const;
447 
451  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
452  orthogError(const MV &X1, const MV &X2) const {
453  return orthogError(X1,Teuchos::null,X2);
454  }
455 
460  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
461  orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const;
462 
464 
466 
467 
470  void setLabel(const std::string& label);
471 
474  const std::string& getLabel() const { return label_; }
475 
477 
478  private:
479 
481  int max_blk_ortho_;
483  MagnitudeType blk_tol_;
485  MagnitudeType dep_tol_;
487  MagnitudeType sing_tol_;
488 
490  std::string label_;
491 #ifdef BELOS_TEUCHOS_TIME_MONITOR
492  Teuchos::RCP<Teuchos::Time> timerOrtho_;
493 #endif // BELOS_TEUCHOS_TIME_MONITOR
494 
496  mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
497 
499  int findBasis(MV &X, Teuchos::RCP<MV> MX,
500  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
501  bool completeBasis, int howMany = -1 ) const;
502 
504  bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
505  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
506  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
507 
521  int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
522  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
523  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
524  Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const;
525  };
526 
528  // Set the label for this orthogonalization manager and create new timers if it's changed
529  template<class ScalarType, class MV, class OP>
530  void DGKSOrthoManager<ScalarType,MV,OP>::setLabel(const std::string& label)
531  {
532  if (label != label_) {
533  label_ = label;
534  std::string orthoLabel = label_ + ": Orthogonalization";
535 #ifdef BELOS_TEUCHOS_TIME_MONITOR
536  timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
537 #endif
538  }
539  }
540 
542  // Compute the distance from orthonormality
543  template<class ScalarType, class MV, class OP>
544  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
545  DGKSOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const {
546  const ScalarType ONE = SCT::one();
547  int rank = MVT::GetNumberVecs(X);
548  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
550  for (int i=0; i<rank; i++) {
551  xTx(i,i) -= ONE;
552  }
553  return xTx.normFrobenius();
554  }
555 
557  // Compute the distance from orthogonality
558  template<class ScalarType, class MV, class OP>
559  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
560  DGKSOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const {
561  int r1 = MVT::GetNumberVecs(X1);
562  int r2 = MVT::GetNumberVecs(X2);
563  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
565  return xTx.normFrobenius();
566  }
567 
569  // Find an Op-orthonormal basis for span(X) - span(W)
570  template<class ScalarType, class MV, class OP>
571  int
574  Teuchos::RCP<MV> MX,
575  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
576  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
577  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
578  {
579  using Teuchos::Array;
580  using Teuchos::null;
581  using Teuchos::is_null;
582  using Teuchos::RCP;
583  using Teuchos::rcp;
584  using Teuchos::SerialDenseMatrix;
585  typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
586  typedef typename Array< RCP< const MV > >::size_type size_type;
587 
588 #ifdef BELOS_TEUCHOS_TIME_MONITOR
589  Teuchos::TimeMonitor orthotimer(*timerOrtho_);
590 #endif
591 
592  ScalarType ONE = SCT::one();
593  ScalarType ZERO = SCT::zero();
594 
595  int nq = Q.size();
596  int xc = MVT::GetNumberVecs( X );
597  ptrdiff_t xr = MVT::GetGlobalLength( X );
598  int rank = xc;
599 
600  // If the user doesn't want to store the normalization
601  // coefficients, allocate some local memory for them. This will
602  // go away at the end of this method.
603  if (is_null (B)) {
604  B = rcp (new serial_dense_matrix_type (xc, xc));
605  }
606  // Likewise, if the user doesn't want to store the projection
607  // coefficients, allocate some local memory for them. Also make
608  // sure that all the entries of C are the right size. We're going
609  // to overwrite them anyway, so we don't have to worry about the
610  // contents (other than to resize them if they are the wrong
611  // size).
612  if (C.size() < nq)
613  C.resize (nq);
614  for (size_type k = 0; k < nq; ++k)
615  {
616  const int numRows = MVT::GetNumberVecs (*Q[k]);
617  const int numCols = xc; // Number of vectors in X
618 
619  if (is_null (C[k]))
620  C[k] = rcp (new serial_dense_matrix_type (numRows, numCols));
621  else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
622  {
623  int err = C[k]->reshape (numRows, numCols);
624  TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error,
625  "DGKS orthogonalization: failed to reshape "
626  "C[" << k << "] (the array of block "
627  "coefficients resulting from projecting X "
628  "against Q[1:" << nq << "]).");
629  }
630  }
631 
632  /****** DO NO MODIFY *MX IF _hasOp == false ******/
633  if (this->_hasOp) {
634  if (MX == Teuchos::null) {
635  // we need to allocate space for MX
636  MX = MVT::Clone(X,MVT::GetNumberVecs(X));
637  OPT::Apply(*(this->_Op),X,*MX);
638  }
639  }
640  else {
641  // Op == I --> MX = X (ignore it if the user passed it in)
642  MX = Teuchos::rcp( &X, false );
643  }
644 
645  int mxc = MVT::GetNumberVecs( *MX );
646  ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
647 
648  // short-circuit
649  TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Belos::DGKSOrthoManager::projectAndNormalize(): X must be non-empty" );
650 
651  int numbas = 0;
652  for (int i=0; i<nq; i++) {
653  numbas += MVT::GetNumberVecs( *Q[i] );
654  }
655 
656  // check size of B
657  TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
658  "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
659  // check size of X and MX
660  TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
661  "Belos::DGKSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
662  // check size of X w.r.t. MX
663  TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
664  "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
665  // check feasibility
666  //TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument,
667  // "Belos::DGKSOrthoManager::projectAndNormalize(): Orthogonality constraints not feasible" );
668 
669  // Some flags for checking dependency returns from the internal orthogonalization methods
670  bool dep_flg = false;
671 
672  // Make a temporary copy of X and MX, just in case a block dependency is detected.
673  Teuchos::RCP<MV> tmpX, tmpMX;
674  tmpX = MVT::CloneCopy(X);
675  if (this->_hasOp) {
676  tmpMX = MVT::CloneCopy(*MX);
677  }
678 
679  // Use the cheaper block orthogonalization.
680  dep_flg = blkOrtho( X, MX, C, Q );
681 
682  // If a dependency has been detected in this block, then perform
683  // the more expensive single-vector orthogonalization.
684  if (dep_flg) {
685  rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
686 
687  // Copy tmpX back into X.
688  MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
689  if (this->_hasOp) {
690  MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
691  }
692  }
693  else {
694  // There is no dependency, so orthonormalize new block X
695  rank = findBasis( X, MX, B, false );
696  if (rank < xc) {
697  // A dependency was found during orthonormalization of X,
698  // rerun orthogonalization using more expensive single-vector orthogonalization.
699  rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
700 
701  // Copy tmpX back into X.
702  MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
703  if (this->_hasOp) {
704  MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
705  }
706  }
707  }
708 
709  // this should not raise an std::exception; but our post-conditions oblige us to check
710  TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
711  "Belos::DGKSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
712 
713  // Return the rank of X.
714  return rank;
715  }
716 
717 
718 
720  // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
721  template<class ScalarType, class MV, class OP>
723  MV &X, Teuchos::RCP<MV> MX,
724  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
725 
726 #ifdef BELOS_TEUCHOS_TIME_MONITOR
727  Teuchos::TimeMonitor orthotimer(*timerOrtho_);
728 #endif
729 
730  // call findBasis, with the instruction to try to generate a basis of rank numvecs(X)
731  return findBasis(X, MX, B, true);
732  }
733 
734 
735 
737  template<class ScalarType, class MV, class OP>
739  MV &X, Teuchos::RCP<MV> MX,
740  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
741  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
742  // For the inner product defined by the operator Op or the identity (Op == 0)
743  // -> Orthogonalize X against each Q[i]
744  // Modify MX accordingly
745  //
746  // Note that when Op is 0, MX is not referenced
747  //
748  // Parameter variables
749  //
750  // X : Vectors to be transformed
751  //
752  // MX : Image of the block vector X by the mass matrix
753  //
754  // Q : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
755  //
756 
757 #ifdef BELOS_TEUCHOS_TIME_MONITOR
758  Teuchos::TimeMonitor orthotimer(*timerOrtho_);
759 #endif
760 
761  int xc = MVT::GetNumberVecs( X );
762  ptrdiff_t xr = MVT::GetGlobalLength( X );
763  int nq = Q.size();
764  std::vector<int> qcs(nq);
765  // short-circuit
766  if (nq == 0 || xc == 0 || xr == 0) {
767  return;
768  }
769  ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
770  // if we don't have enough C, expand it with null references
771  // if we have too many, resize to throw away the latter ones
772  // if we have exactly as many as we have Q, this call has no effect
773  C.resize(nq);
774 
775 
776  /****** DO NO MODIFY *MX IF _hasOp == false ******/
777  if (this->_hasOp) {
778  if (MX == Teuchos::null) {
779  // we need to allocate space for MX
780  MX = MVT::Clone(X,MVT::GetNumberVecs(X));
781  OPT::Apply(*(this->_Op),X,*MX);
782  }
783  }
784  else {
785  // Op == I --> MX = X (ignore it if the user passed it in)
786  MX = Teuchos::rcp( &X, false );
787  }
788  int mxc = MVT::GetNumberVecs( *MX );
789  ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
790 
791  // check size of X and Q w.r.t. common sense
792  TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
793  "Belos::DGKSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
794  // check size of X w.r.t. MX and Q
795  TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
796  "Belos::DGKSOrthoManager::project(): Size of X not consistant with MX,Q" );
797 
798  // tally up size of all Q and check/allocate C
799  int baslen = 0;
800  for (int i=0; i<nq; i++) {
801  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
802  "Belos::DGKSOrthoManager::project(): Q lengths not mutually consistant" );
803  qcs[i] = MVT::GetNumberVecs( *Q[i] );
804  TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
805  "Belos::DGKSOrthoManager::project(): Q has less rows than columns" );
806  baslen += qcs[i];
807 
808  // check size of C[i]
809  if ( C[i] == Teuchos::null ) {
810  C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
811  }
812  else {
813  TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
814  "Belos::DGKSOrthoManager::project(): Size of Q not consistant with size of C" );
815  }
816  }
817 
818  // Use the cheaper block orthogonalization, don't check for rank deficiency.
819  blkOrtho( X, MX, C, Q );
820 
821  }
822 
824  // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that
825  // the rank is numvectors(X)
826  template<class ScalarType, class MV, class OP>
828  MV &X, Teuchos::RCP<MV> MX,
829  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
830  bool completeBasis, int howMany ) const {
831  // For the inner product defined by the operator Op or the identity (Op == 0)
832  // -> Orthonormalize X
833  // Modify MX accordingly
834  //
835  // Note that when Op is 0, MX is not referenced
836  //
837  // Parameter variables
838  //
839  // X : Vectors to be orthonormalized
840  //
841  // MX : Image of the multivector X under the operator Op
842  //
843  // Op : Pointer to the operator for the inner product
844  //
845  //
846 
847  const ScalarType ONE = SCT::one();
848  const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
849 
850  int xc = MVT::GetNumberVecs( X );
851  ptrdiff_t xr = MVT::GetGlobalLength( X );
852 
853  if (howMany == -1) {
854  howMany = xc;
855  }
856 
857  /*******************************************************
858  * If _hasOp == false, we will not reference MX below *
859  *******************************************************/
860 
861  // if Op==null, MX == X (via pointer)
862  // Otherwise, either the user passed in MX or we will allocated and compute it
863  if (this->_hasOp) {
864  if (MX == Teuchos::null) {
865  // we need to allocate space for MX
866  MX = MVT::Clone(X,xc);
867  OPT::Apply(*(this->_Op),X,*MX);
868  }
869  }
870 
871  /* if the user doesn't want to store the coefficienets,
872  * allocate some local memory for them
873  */
874  if ( B == Teuchos::null ) {
875  B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
876  }
877 
878  int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
879  ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
880 
881  // check size of C, B
882  TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
883  "Belos::DGKSOrthoManager::findBasis(): X must be non-empty" );
884  TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
885  "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of B" );
886  TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
887  "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
888  TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument,
889  "Belos::DGKSOrthoManager::findBasis(): Size of X not feasible for normalization" );
890  TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
891  "Belos::DGKSOrthoManager::findBasis(): Invalid howMany parameter" );
892 
893  /* xstart is which column we are starting the process with, based on howMany
894  * columns before xstart are assumed to be Op-orthonormal already
895  */
896  int xstart = xc - howMany;
897 
898  for (int j = xstart; j < xc; j++) {
899 
900  // numX is
901  // * number of currently orthonormal columns of X
902  // * the index of the current column of X
903  int numX = j;
904  bool addVec = false;
905 
906  // Get a view of the vector currently being worked on.
907  std::vector<int> index(1);
908  index[0] = numX;
909  Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
910  Teuchos::RCP<MV> MXj;
911  if ((this->_hasOp)) {
912  // MXj is a view of the current vector in MX
913  MXj = MVT::CloneViewNonConst( *MX, index );
914  }
915  else {
916  // MXj is a pointer to Xj, and MUST NOT be modified
917  MXj = Xj;
918  }
919 
920  // Get a view of the previous vectors.
921  std::vector<int> prev_idx( numX );
922  Teuchos::RCP<const MV> prevX, prevMX;
923 
924  if (numX > 0) {
925  for (int i=0; i<numX; i++) {
926  prev_idx[i] = i;
927  }
928  prevX = MVT::CloneView( X, prev_idx );
929  if (this->_hasOp) {
930  prevMX = MVT::CloneView( *MX, prev_idx );
931  }
932  }
933 
934  // Make storage for these Gram-Schmidt iterations.
935  Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
936  std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
937  //
938  // Save old MXj vector and compute Op-norm
939  //
940  Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj );
941  MVT::MvDot( *Xj, *MXj, oldDot );
942  // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
943  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError,
944  "Belos::DGKSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
945 
946  if (numX > 0) {
947  // Apply the first step of Gram-Schmidt
948 
949  // product <- prevX^T MXj
950  MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*Xj,MXj,product);
951 
952  // Xj <- Xj - prevX prevX^T MXj
953  // = Xj - prevX product
954  MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
955 
956  // Update MXj
957  if (this->_hasOp) {
958  // MXj <- Op*Xj_new
959  // = Op*(Xj_old - prevX prevX^T MXj)
960  // = MXj - prevMX product
961  MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
962  }
963 
964  // Compute new Op-norm
965  MVT::MvDot( *Xj, *MXj, newDot );
966 
967  // Check if a correction is needed.
968  if ( SCT::magnitude(newDot[0]) < SCT::magnitude(dep_tol_*oldDot[0]) ) {
969  // Apply the second step of Gram-Schmidt
970  // This is the same as above
971  Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
972 
974  product += P2;
975  MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
976  if ((this->_hasOp)) {
977  MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
978  }
979  } // if (newDot[0] < dep_tol_*oldDot[0])
980 
981  } // if (numX > 0)
982 
983  // Compute Op-norm with old MXj
984  MVT::MvDot( *Xj, *oldMXj, newDot );
985 
986  // Check to see if the new vector is dependent.
987  if (completeBasis) {
988  //
989  // We need a complete basis, so add random vectors if necessary
990  //
991  if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
992 
993  // Add a random vector and orthogonalize it against previous vectors in block.
994  addVec = true;
995 #ifdef ORTHO_DEBUG
996  std::cout << "Belos::DGKSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
997 #endif
998  //
999  Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1000  Teuchos::RCP<MV> tempMXj;
1001  MVT::MvRandom( *tempXj );
1002  if (this->_hasOp) {
1003  tempMXj = MVT::Clone( X, 1 );
1004  OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1005  }
1006  else {
1007  tempMXj = tempXj;
1008  }
1009  MVT::MvDot( *tempXj, *tempMXj, oldDot );
1010  //
1011  for (int num_orth=0; num_orth<max_blk_ortho_; num_orth++){
1012  MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*tempXj,tempMXj,product);
1013  MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
1014  if (this->_hasOp) {
1015  MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
1016  }
1017  }
1018  // Compute new Op-norm
1019  MVT::MvDot( *tempXj, *tempMXj, newDot );
1020  //
1021  if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1022  // Copy vector into current column of _basisvecs
1023  MVT::MvAddMv( ONE, *tempXj, ZERO, *tempXj, *Xj );
1024  if (this->_hasOp) {
1025  MVT::MvAddMv( ONE, *tempMXj, ZERO, *tempMXj, *MXj );
1026  }
1027  }
1028  else {
1029  return numX;
1030  }
1031  }
1032  }
1033  else {
1034  //
1035  // We only need to detect dependencies.
1036  //
1037  if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1038  return numX;
1039  }
1040  }
1041 
1042  // If we haven't left this method yet, then we can normalize the new vector Xj.
1043  // Normalize Xj.
1044  // Xj <- Xj / std::sqrt(newDot)
1045  ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1046 
1047  if (SCT::magnitude(diag) > ZERO) {
1048  MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
1049  if (this->_hasOp) {
1050  // Update MXj.
1051  MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
1052  }
1053  }
1054 
1055  // If we've added a random vector, enter a zero in the j'th diagonal element.
1056  if (addVec) {
1057  (*B)(j,j) = ZERO;
1058  }
1059  else {
1060  (*B)(j,j) = diag;
1061  }
1062 
1063  // Save the coefficients, if we are working on the original vector and not a randomly generated one
1064  if (!addVec) {
1065  for (int i=0; i<numX; i++) {
1066  (*B)(i,j) = product(i,0);
1067  }
1068  }
1069 
1070  } // for (j = 0; j < xc; ++j)
1071 
1072  return xc;
1073  }
1074 
1076  // Routine to compute the block orthogonalization
1077  template<class ScalarType, class MV, class OP>
1078  bool
1079  DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
1080  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1081  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
1082  {
1083  int nq = Q.size();
1084  int xc = MVT::GetNumberVecs( X );
1085  bool dep_flg = false;
1086  const ScalarType ONE = SCT::one();
1087 
1088  std::vector<int> qcs( nq );
1089  for (int i=0; i<nq; i++) {
1090  qcs[i] = MVT::GetNumberVecs( *Q[i] );
1091  }
1092 
1093  // Perform the Gram-Schmidt transformation for a block of vectors
1094 
1095  // Compute the initial Op-norms
1096  std::vector<ScalarType> oldDot( xc );
1097  MVT::MvDot( X, *MX, oldDot );
1098 
1099  Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1100  // Define the product Q^T * (Op*X)
1101  for (int i=0; i<nq; i++) {
1102  // Multiply Q' with MX
1104  // Multiply by Q and subtract the result in X
1105  MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1106 
1107  // Update MX, with the least number of applications of Op as possible
1108  if (this->_hasOp) {
1109  if (xc <= qcs[i]) {
1110  OPT::Apply( *(this->_Op), X, *MX);
1111  }
1112  else {
1113  // this will possibly be used again below; don't delete it
1114  MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1115  OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1116  MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1117  }
1118  }
1119  }
1120 
1121  // Do as many steps of classical Gram-Schmidt as required by max_blk_ortho_
1122  for (int j = 1; j < max_blk_ortho_; ++j) {
1123 
1124  for (int i=0; i<nq; i++) {
1125  Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(), C[i]->numCols());
1126 
1127  // Apply another step of classical Gram-Schmidt
1129  *C[i] += C2;
1130  MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1131 
1132  // Update MX, with the least number of applications of Op as possible
1133  if (this->_hasOp) {
1134  if (MQ[i].get()) {
1135  // MQ was allocated and computed above; use it
1136  MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1137  }
1138  else if (xc <= qcs[i]) {
1139  // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1140  OPT::Apply( *(this->_Op), X, *MX);
1141  }
1142  }
1143  } // for (int i=0; i<nq; i++)
1144  } // for (int j = 0; j < max_blk_ortho; ++j)
1145 
1146  // Compute new Op-norms
1147  std::vector<ScalarType> newDot(xc);
1148  MVT::MvDot( X, *MX, newDot );
1149 
1150  // Check to make sure the new block of vectors are not dependent on previous vectors
1151  for (int i=0; i<xc; i++){
1152  if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1153  dep_flg = true;
1154  break;
1155  }
1156  } // end for (i=0;...)
1157 
1158  return dep_flg;
1159  }
1160 
1161 
1162  template<class ScalarType, class MV, class OP>
1163  int
1164  DGKSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
1165  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1166  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1167  Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const
1168  {
1169  Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1170 
1171  const ScalarType ONE = SCT::one();
1172  const ScalarType ZERO = SCT::zero();
1173 
1174  int nq = Q.size();
1175  int xc = MVT::GetNumberVecs( X );
1176  std::vector<int> indX( 1 );
1177  std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1178 
1179  std::vector<int> qcs( nq );
1180  for (int i=0; i<nq; i++) {
1181  qcs[i] = MVT::GetNumberVecs( *Q[i] );
1182  }
1183 
1184  // Create pointers for the previous vectors of X that have already been orthonormalized.
1185  Teuchos::RCP<const MV> lastQ;
1186  Teuchos::RCP<MV> Xj, MXj;
1187  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1188 
1189  // Perform the Gram-Schmidt transformation for each vector in the block of vectors.
1190  for (int j=0; j<xc; j++) {
1191 
1192  bool dep_flg = false;
1193 
1194  // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
1195  if (j > 0) {
1196  std::vector<int> index( j );
1197  for (int ind=0; ind<j; ind++) {
1198  index[ind] = ind;
1199  }
1200  lastQ = MVT::CloneView( X, index );
1201 
1202  // Add these views to the Q and C arrays.
1203  Q.push_back( lastQ );
1204  C.push_back( B );
1205  qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1206  }
1207 
1208  // Get a view of the current vector in X to orthogonalize.
1209  indX[0] = j;
1210  Xj = MVT::CloneViewNonConst( X, indX );
1211  if (this->_hasOp) {
1212  MXj = MVT::CloneViewNonConst( *MX, indX );
1213  }
1214  else {
1215  MXj = Xj;
1216  }
1217 
1218  // Compute the initial Op-norms
1219  MVT::MvDot( *Xj, *MXj, oldDot );
1220 
1221  Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1222  // Define the product Q^T * (Op*X)
1223  for (int i=0; i<Q.size(); i++) {
1224 
1225  // Get a view of the current serial dense matrix
1226  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1227 
1228  // Multiply Q' with MX
1229  MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*Xj,MXj,tempC);
1230  // Multiply by Q and subtract the result in Xj
1231  MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
1232 
1233  // Update MXj, with the least number of applications of Op as possible
1234  if (this->_hasOp) {
1235  if (xc <= qcs[i]) {
1236  OPT::Apply( *(this->_Op), *Xj, *MXj);
1237  }
1238  else {
1239  // this will possibly be used again below; don't delete it
1240  MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1241  OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1242  MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1243  }
1244  }
1245  }
1246 
1247  // Compute the Op-norms
1248  MVT::MvDot( *Xj, *MXj, newDot );
1249 
1250  // Do one step of classical Gram-Schmidt orthogonalization
1251  // with a second correction step if needed.
1252 
1253  if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*dep_tol_) ) {
1254 
1255  for (int i=0; i<Q.size(); i++) {
1256  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1257  Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1258 
1259  // Apply another step of classical Gram-Schmidt
1261  tempC += C2;
1262  MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
1263 
1264  // Update MXj, with the least number of applications of Op as possible
1265  if (this->_hasOp) {
1266  if (MQ[i].get()) {
1267  // MQ was allocated and computed above; use it
1268  MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1269  }
1270  else if (xc <= qcs[i]) {
1271  // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1272  OPT::Apply( *(this->_Op), *Xj, *MXj);
1273  }
1274  }
1275  } // for (int i=0; i<Q.size(); i++)
1276 
1277  // Compute the Op-norms after the correction step.
1278  MVT::MvDot( *Xj, *MXj, newDot );
1279 
1280  } // if ()
1281 
1282  // Check for linear dependence.
1283  if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1284  dep_flg = true;
1285  }
1286 
1287  // Normalize the new vector if it's not dependent
1288  if (!dep_flg) {
1289  ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1290 
1291  MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
1292  if (this->_hasOp) {
1293  // Update MXj.
1294  MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
1295  }
1296 
1297  // Enter value on diagonal of B.
1298  (*B)(j,j) = diag;
1299  }
1300  else {
1301  // Create a random vector and orthogonalize it against all previous columns of Q.
1302  Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1303  Teuchos::RCP<MV> tempMXj;
1304  MVT::MvRandom( *tempXj );
1305  if (this->_hasOp) {
1306  tempMXj = MVT::Clone( X, 1 );
1307  OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1308  }
1309  else {
1310  tempMXj = tempXj;
1311  }
1312  MVT::MvDot( *tempXj, *tempMXj, oldDot );
1313  //
1314  for (int num_orth=0; num_orth<max_blk_ortho_; num_orth++) {
1315 
1316  for (int i=0; i<Q.size(); i++) {
1317  Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1318 
1319  // Apply another step of classical Gram-Schmidt
1320  MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*tempXj,tempMXj,product);
1321  MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
1322 
1323  // Update MXj, with the least number of applications of Op as possible
1324  if (this->_hasOp) {
1325  if (MQ[i].get()) {
1326  // MQ was allocated and computed above; use it
1327  MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1328  }
1329  else if (xc <= qcs[i]) {
1330  // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1331  OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1332  }
1333  }
1334  } // for (int i=0; i<nq; i++)
1335 
1336  }
1337 
1338  // Compute the Op-norms after the correction step.
1339  MVT::MvDot( *tempXj, *tempMXj, newDot );
1340 
1341  // Copy vector into current column of Xj
1342  if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1343  ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1344 
1345  // Enter value on diagonal of B.
1346  (*B)(j,j) = ZERO;
1347 
1348  // Copy vector into current column of _basisvecs
1349  MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1350  if (this->_hasOp) {
1351  MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1352  }
1353  }
1354  else {
1355  return j;
1356  }
1357  } // if (!dep_flg)
1358 
1359  // Remove the vectors from array
1360  if (j > 0) {
1361  Q.resize( nq );
1362  C.resize( nq );
1363  qcs.resize( nq );
1364  }
1365 
1366  } // for (int j=0; j<xc; j++)
1367 
1368  return xc;
1369  }
1370 
1371 } // namespace Belos
1372 
1373 #endif // BELOS_DGKS_ORTHOMANAGER_HPP
1374 
void setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
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 ...
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 ...
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
void setDepTol(const MagnitudeType dep_tol)
Set parameter for re-orthogonalization threshhold.
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.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
DGKSOrthoManager(const Teuchos::RCP< Teuchos::ParameterList > &plist, const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null)
Constructor that takes a list of parameters.
Class which defines basic traits for the operator type.
void setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
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...
MagnitudeType getDepTol() const
Return parameter for re-orthogonalization threshhold.
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 ...
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
Exception thrown to signal error in an orthogonalization manager method.
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.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
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...
DGKSOrthoManager(const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null, const int max_blk_ortho=2, const MagnitudeType blk_tol=10 *MGT::squareroot(MGT::eps()), const MagnitudeType dep_tol=MGT::one()/MGT::squareroot(2 *MGT::one()), 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 ...
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
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.
Belos&#39;s templated virtual class for providing routines for orthogonalization and orthonormzalition of...
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
Belos header file which uses auto-configuration information to include necessary C++ headers...
MagnitudeType getSingTol() const
Return parameter for singular block detection.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters() const

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