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

Generated on Thu Mar 28 2024 09:24:29 for Belos by doxygen 1.8.5