Anasazi  Version of the Day
AnasaziBasicOrthoManager.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under 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 ANASAZI_BASIC_ORTHOMANAGER_HPP
48 #define ANASAZI_BASIC_ORTHOMANAGER_HPP
49 
57 #include "AnasaziConfigDefs.hpp"
61 #include "Teuchos_TimeMonitor.hpp"
62 #include "Teuchos_LAPACK.hpp"
63 #include "Teuchos_BLAS.hpp"
64 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
65 # include <Teuchos_FancyOStream.hpp>
66 #endif
67 
68 namespace Anasazi {
69 
70  template<class ScalarType, class MV, class OP>
71  class BasicOrthoManager : public MatOrthoManager<ScalarType,MV,OP> {
72 
73  private:
74  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
75  typedef Teuchos::ScalarTraits<ScalarType> SCT;
78 
79  public:
80 
82 
83  BasicOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null,
85  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa = 1.41421356 /* sqrt(2) */,
86  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps = 0.0,
87  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol = 0.20 );
88 
89 
93 
94 
96 
97 
98 
137  void projectMat (
138  MV &X,
139  Teuchos::Array<Teuchos::RCP<const MV> > Q,
140  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
141  = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
142  Teuchos::RCP<MV> MX = Teuchos::null,
143  Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
144  ) const;
145 
146 
185  int normalizeMat (
186  MV &X,
187  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
188  Teuchos::RCP<MV> MX = Teuchos::null
189  ) const;
190 
191 
252  MV &X,
253  Teuchos::Array<Teuchos::RCP<const MV> > Q,
254  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
255  = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
256  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
257  Teuchos::RCP<MV> MX = Teuchos::null,
258  Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
259  ) const;
260 
262 
264 
265 
270  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
271  orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const;
272 
277  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
278  orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const;
279 
281 
283 
284 
286  void setKappa( typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa ) { kappa_ = kappa; }
287 
289  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType getKappa() const { return kappa_; }
290 
292 
293  private:
295  MagnitudeType kappa_;
296  MagnitudeType eps_;
297  MagnitudeType tol_;
298 
299  // ! Routine to find an orthonormal basis
300  int findBasis(MV &X, Teuchos::RCP<MV> MX,
301  Teuchos::SerialDenseMatrix<int,ScalarType> &B,
302  bool completeBasis, int howMany = -1 ) const;
303 
304  //
305  // Internal timers
306  //
307  Teuchos::RCP<Teuchos::Time> timerReortho_;
308 
309  };
310 
311 
313  // Constructor
314  template<class ScalarType, class MV, class OP>
316  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa,
317  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps,
318  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol ) :
319  MatOrthoManager<ScalarType,MV,OP>(Op), kappa_(kappa), eps_(eps), tol_(tol)
320 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
321  , timerReortho_(Teuchos::TimeMonitor::getNewTimer("Anasazi::BasicOrthoManager::Re-orthogonalization"))
322 #endif
323  {
324  TEUCHOS_TEST_FOR_EXCEPTION(eps_ < SCT::magnitude(SCT::zero()),std::invalid_argument,
325  "Anasazi::BasicOrthoManager::BasicOrthoManager(): argument \"eps\" must be non-negative.");
326  if (eps_ == 0) {
327  Teuchos::LAPACK<int,MagnitudeType> lapack;
328  eps_ = lapack.LAMCH('E');
329  eps_ = Teuchos::ScalarTraits<MagnitudeType>::pow(eps_,.75);
330  }
331  TEUCHOS_TEST_FOR_EXCEPTION(
332  tol_ < SCT::magnitude(SCT::zero()) || tol_ > SCT::magnitude(SCT::one()),
333  std::invalid_argument,
334  "Anasazi::BasicOrthoManager::BasicOrthoManager(): argument \"tol\" must be in [0,1].");
335  }
336 
337 
338 
340  // Compute the distance from orthonormality
341  template<class ScalarType, class MV, class OP>
342  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
343  BasicOrthoManager<ScalarType,MV,OP>::orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX) const {
344  const ScalarType ONE = SCT::one();
345  int rank = MVT::GetNumberVecs(X);
346  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
348  for (int i=0; i<rank; i++) {
349  xTx(i,i) -= ONE;
350  }
351  return xTx.normFrobenius();
352  }
353 
354 
355 
357  // Compute the distance from orthogonality
358  template<class ScalarType, class MV, class OP>
359  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
360  BasicOrthoManager<ScalarType,MV,OP>::orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const {
361  int r1 = MVT::GetNumberVecs(X1);
362  int r2 = MVT::GetNumberVecs(X2);
363  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
365  return xTx.normFrobenius();
366  }
367 
368 
369 
371  template<class ScalarType, class MV, class OP>
373  MV &X,
374  Teuchos::Array<Teuchos::RCP<const MV> > Q,
375  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
376  Teuchos::RCP<MV> MX,
377  Teuchos::Array<Teuchos::RCP<const MV> > MQ
378  ) const {
379  // For the inner product defined by the operator Op or the identity (Op == 0)
380  // -> Orthogonalize X against each Q[i]
381  // Modify MX accordingly
382  //
383  // Note that when Op is 0, MX is not referenced
384  //
385  // Parameter variables
386  //
387  // X : Vectors to be transformed
388  //
389  // MX : Image of the block vector X by the mass matrix
390  //
391  // Q : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
392  //
393 
394 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
395  // Get a FancyOStream from out_arg or create a new one ...
396  Teuchos::RCP<Teuchos::FancyOStream>
397  out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
398  out->setShowAllFrontMatter(false).setShowProcRank(true);
399  *out << "Entering Anasazi::BasicOrthoManager::projectMat(...)\n";
400 #endif
401 
402  ScalarType ONE = SCT::one();
403 
404  int xc = MVT::GetNumberVecs( X );
405  ptrdiff_t xr = MVT::GetGlobalLength( X );
406  int nq = Q.length();
407  std::vector<int> qcs(nq);
408  // short-circuit
409  if (nq == 0 || xc == 0 || xr == 0) {
410 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
411  *out << "Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
412 #endif
413  return;
414  }
415  ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
416  // if we don't have enough C, expand it with null references
417  // if we have too many, resize to throw away the latter ones
418  // if we have exactly as many as we have Q, this call has no effect
419  C.resize(nq);
420 
421 
422  /****** DO NO MODIFY *MX IF _hasOp == false ******/
423  if (this->_hasOp) {
424 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
425  *out << "Allocating MX...\n";
426 #endif
427  if (MX == Teuchos::null) {
428  // we need to allocate space for MX
429  MX = MVT::Clone(X,MVT::GetNumberVecs(X));
430  OPT::Apply(*(this->_Op),X,*MX);
431  this->_OpCounter += MVT::GetNumberVecs(X);
432  }
433  }
434  else {
435  // Op == I --> MX = X (ignore it if the user passed it in)
436  MX = Teuchos::rcpFromRef(X);
437  }
438  int mxc = MVT::GetNumberVecs( *MX );
439  ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
440 
441  // check size of X and Q w.r.t. common sense
442  TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
443  "Anasazi::BasicOrthoManager::projectMat(): MVT returned negative dimensions for X,MX" );
444  // check size of X w.r.t. MX and Q
445  TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
446  "Anasazi::BasicOrthoManager::projectMat(): Size of X not consistent with MX,Q" );
447 
448  // tally up size of all Q and check/allocate C
449  int baslen = 0;
450  for (int i=0; i<nq; i++) {
451  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
452  "Anasazi::BasicOrthoManager::projectMat(): Q lengths not mutually consistent" );
453  qcs[i] = MVT::GetNumberVecs( *Q[i] );
454  TEUCHOS_TEST_FOR_EXCEPTION( qr < static_cast<ptrdiff_t>(qcs[i]), std::invalid_argument,
455  "Anasazi::BasicOrthoManager::projectMat(): Q has less rows than columns" );
456  baslen += qcs[i];
457 
458  // check size of C[i]
459  if ( C[i] == Teuchos::null ) {
460  C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
461  }
462  else {
463  TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
464  "Anasazi::BasicOrthoManager::projectMat(): Size of Q not consistent with size of C" );
465  }
466  }
467 
468  // Perform the Gram-Schmidt transformation for a block of vectors
469 
470  // Compute the initial Op-norms
471  std::vector<ScalarType> oldDot( xc );
472  MVT::MvDot( X, *MX, oldDot );
473 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
474  *out << "oldDot = { ";
475  std::copy(oldDot.begin(), oldDot.end(), std::ostream_iterator<ScalarType>(*out, " "));
476  *out << "}\n";
477 #endif
478 
479  MQ.resize(nq);
480  // Define the product Q^T * (Op*X)
481  for (int i=0; i<nq; i++) {
482  // Multiply Q' with MX
483  MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*Q[i],X,*C[i],MQ[i],MX);
484  // Multiply by Q and subtract the result in X
485 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
486  *out << "Applying projector P_Q[" << i << "]...\n";
487 #endif
488  MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
489 
490  // Update MX, with the least number of applications of Op as possible
491  // Update MX. If we have MQ, use it. Otherwise, just multiply by Op
492  if (this->_hasOp) {
493  if (MQ[i] == Teuchos::null) {
494 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
495  *out << "Updating MX via M*X...\n";
496 #endif
497  OPT::Apply( *(this->_Op), X, *MX);
498  this->_OpCounter += MVT::GetNumberVecs(X);
499  }
500  else {
501 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
502  *out << "Updating MX via M*Q...\n";
503 #endif
504  MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
505  }
506  }
507  }
508 
509  // Compute new Op-norms
510  std::vector<ScalarType> newDot(xc);
511  MVT::MvDot( X, *MX, newDot );
512 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
513  *out << "newDot = { ";
514  std::copy(newDot.begin(), newDot.end(), std::ostream_iterator<ScalarType>(*out, " "));
515  *out << "}\n";
516 #endif
517 
518  // determine (individually) whether to do another step of classical Gram-Schmidt
519  for (int j = 0; j < xc; ++j) {
520 
521  if ( SCT::magnitude(kappa_*newDot[j]) < SCT::magnitude(oldDot[j]) ) {
522 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
523  *out << "kappa_*newDot[" <<j<< "] == " << kappa_*newDot[j] << "... another step of Gram-Schmidt.\n";
524 #endif
525 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
526  Teuchos::TimeMonitor lcltimer( *timerReortho_ );
527 #endif
528  for (int i=0; i<nq; i++) {
529  Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(), C[i]->numCols());
530 
531  // Apply another step of classical Gram-Schmidt
533  *C[i] += C2;
534 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
535  *out << "Applying projector P_Q[" << i << "]...\n";
536 #endif
537  MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
538 
539  // Update MX as above
540  if (this->_hasOp) {
541  if (MQ[i] == Teuchos::null) {
542 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
543  *out << "Updating MX via M*X...\n";
544 #endif
545  OPT::Apply( *(this->_Op), X, *MX);
546  this->_OpCounter += MVT::GetNumberVecs(X);
547  }
548  else {
549 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
550  *out << "Updating MX via M*Q...\n";
551 #endif
552  MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
553  }
554  }
555  }
556  break;
557  } // if (kappa_*newDot[j] < oldDot[j])
558  } // for (int j = 0; j < xc; ++j)
559 
560 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
561  *out << "Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
562 #endif
563  }
564 
565 
566 
568  // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
569  template<class ScalarType, class MV, class OP>
571  MV &X,
572  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
573  Teuchos::RCP<MV> MX) const {
574  // call findBasis(), with the instruction to try to generate a basis of rank numvecs(X)
575  // findBasis() requires MX
576 
577  int xc = MVT::GetNumberVecs(X);
578  ptrdiff_t xr = MVT::GetGlobalLength(X);
579 
580  // if Op==null, MX == X (via pointer)
581  // Otherwise, either the user passed in MX or we will allocated and compute it
582  if (this->_hasOp) {
583  if (MX == Teuchos::null) {
584  // we need to allocate space for MX
585  MX = MVT::Clone(X,xc);
586  OPT::Apply(*(this->_Op),X,*MX);
587  this->_OpCounter += MVT::GetNumberVecs(X);
588  }
589  }
590 
591  // if the user doesn't want to store the coefficients,
592  // allocate some local memory for them
593  if ( B == Teuchos::null ) {
594  B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
595  }
596 
597  int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
598  ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
599 
600  // check size of C, B
601  TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
602  "Anasazi::BasicOrthoManager::normalizeMat(): X must be non-empty" );
603  TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
604  "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of B" );
605  TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
606  "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of MX" );
607  TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument,
608  "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not feasible for normalization" );
609 
610  return findBasis(X, MX, *B, true );
611  }
612 
613 
614 
616  // Find an Op-orthonormal basis for span(X) - span(W)
617  template<class ScalarType, class MV, class OP>
619  MV &X,
620  Teuchos::Array<Teuchos::RCP<const MV> > Q,
621  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
622  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
623  Teuchos::RCP<MV> MX,
624  Teuchos::Array<Teuchos::RCP<const MV> > MQ
625  ) const {
626 
627 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
628  // Get a FancyOStream from out_arg or create a new one ...
629  Teuchos::RCP<Teuchos::FancyOStream>
630  out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
631  out->setShowAllFrontMatter(false).setShowProcRank(true);
632  *out << "Entering Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
633 #endif
634 
635  int nq = Q.length();
636  int xc = MVT::GetNumberVecs( X );
637  ptrdiff_t xr = MVT::GetGlobalLength( X );
638  int rank;
639 
640  /* if the user doesn't want to store the coefficients,
641  * allocate some local memory for them
642  */
643  if ( B == Teuchos::null ) {
644  B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
645  }
646 
647  /****** DO NO MODIFY *MX IF _hasOp == false ******/
648  if (this->_hasOp) {
649  if (MX == Teuchos::null) {
650  // we need to allocate space for MX
651 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
652  *out << "Allocating MX...\n";
653 #endif
654  MX = MVT::Clone(X,MVT::GetNumberVecs(X));
655  OPT::Apply(*(this->_Op),X,*MX);
656  this->_OpCounter += MVT::GetNumberVecs(X);
657  }
658  }
659  else {
660  // Op == I --> MX = X (ignore it if the user passed it in)
661  MX = Teuchos::rcpFromRef(X);
662  }
663 
664  int mxc = MVT::GetNumberVecs( *MX );
665  ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
666 
667  TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): X must be non-empty" );
668 
669  ptrdiff_t numbas = 0;
670  for (int i=0; i<nq; i++) {
671  numbas += MVT::GetNumberVecs( *Q[i] );
672  }
673 
674  // check size of B
675  TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
676  "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of B" );
677  // check size of X and MX
678  TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
679  "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): MVT returned negative dimensions for X,MX" );
680  // check size of X w.r.t. MX
681  TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
682  "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of MX" );
683  // check feasibility
684  TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument,
685  "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Orthogonality constraints not feasible" );
686 
687  // orthogonalize all of X against Q
688 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
689  *out << "Orthogonalizing X against Q...\n";
690 #endif
691  projectMat(X,Q,C,MX,MQ);
692 
693  Teuchos::SerialDenseMatrix<int,ScalarType> oldCoeff(xc,1);
694  // start working
695  rank = 0;
696  int numTries = 10; // each vector in X gets 10 random chances to escape degeneracy
697  int oldrank = -1;
698  do {
699  int curxsize = xc - rank;
700 
701  // orthonormalize X, but quit if it is rank deficient
702  // we can't let findBasis generated random vectors to complete the basis,
703  // because it doesn't know about Q; we will do this ourselves below
704 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
705  *out << "Attempting to find orthonormal basis for X...\n";
706 #endif
707  rank = findBasis(X,MX,*B,false,curxsize);
708 
709  if (oldrank != -1 && rank != oldrank) {
710  // we had previously stopped before, after operating on vector oldrank
711  // we saved its coefficients, augmented it with a random vector, and
712  // then called findBasis() again, which proceeded to add vector oldrank
713  // to the basis.
714  // now, restore the saved coefficients into B
715  for (int i=0; i<xc; i++) {
716  (*B)(i,oldrank) = oldCoeff(i,0);
717  }
718  }
719 
720  if (rank < xc) {
721  if (rank != oldrank) {
722  // we quit on this vector and will augment it with random below
723  // this is the first time that we have quit on this vector
724  // therefor, (*B)(:,rank) contains the actual coefficients of the
725  // input vectors with respect to the previous vectors in the basis
726  // save these values, as (*B)(:,rank) will be overwritten by our next
727  // call to findBasis()
728  // we will restore it after we are done working on this vector
729  for (int i=0; i<xc; i++) {
730  oldCoeff(i,0) = (*B)(i,rank);
731  }
732  }
733  }
734 
735  if (rank == xc) {
736  // we are done
737 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
738  *out << "Finished computing basis.\n";
739 #endif
740  break;
741  }
742  else {
743  TEUCHOS_TEST_FOR_EXCEPTION( rank < oldrank, OrthoError,
744  "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): basis lost rank; this shouldn't happen");
745 
746  if (rank != oldrank) {
747  // we added a vector to the basis; reset the chance counter
748  numTries = 10;
749  // store old rank
750  oldrank = rank;
751  }
752  else {
753  // has this vector run out of chances to escape degeneracy?
754  if (numTries <= 0) {
755  break;
756  }
757  }
758  // use one of this vector's chances
759  numTries--;
760 
761  // randomize troubled direction
762 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
763  *out << "Randomizing X[" << rank << "]...\n";
764 #endif
765  Teuchos::RCP<MV> curX, curMX;
766  std::vector<int> ind(1);
767  ind[0] = rank;
768  curX = MVT::CloneViewNonConst(X,ind);
769  MVT::MvRandom(*curX);
770  if (this->_hasOp) {
771 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
772  *out << "Applying operator to random vector.\n";
773 #endif
774  curMX = MVT::CloneViewNonConst(*MX,ind);
775  OPT::Apply( *(this->_Op), *curX, *curMX );
776  this->_OpCounter += MVT::GetNumberVecs(*curX);
777  }
778 
779  // orthogonalize against Q
780  // if !this->_hasOp, the curMX will be ignored.
781  // we don't care about these coefficients
782  // on the contrary, we need to preserve the previous coeffs
783  {
784  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC(0);
785  projectMat(*curX,Q,dummyC,curMX,MQ);
786  }
787  }
788  } while (1);
789 
790  // this should never raise an exception; but our post-conditions oblige us to check
791  TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
792  "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Debug error in rank variable." );
793 
794 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
795  *out << "Leaving Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
796 #endif
797 
798  return rank;
799  }
800 
801 
802 
804  // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that
805  // the rank is numvectors(X)
806  template<class ScalarType, class MV, class OP>
808  MV &X, Teuchos::RCP<MV> MX,
809  Teuchos::SerialDenseMatrix<int,ScalarType> &B,
810  bool completeBasis, int howMany ) const {
811 
812  // For the inner product defined by the operator Op or the identity (Op == 0)
813  // -> Orthonormalize X
814  // Modify MX accordingly
815  //
816  // Note that when Op is 0, MX is not referenced
817  //
818  // Parameter variables
819  //
820  // X : Vectors to be orthonormalized
821  //
822  // MX : Image of the multivector X under the operator Op
823  //
824  // Op : Pointer to the operator for the inner product
825  //
826 
827 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
828  // Get a FancyOStream from out_arg or create a new one ...
829  Teuchos::RCP<Teuchos::FancyOStream>
830  out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
831  out->setShowAllFrontMatter(false).setShowProcRank(true);
832  *out << "Entering Anasazi::BasicOrthoManager::findBasis(...)\n";
833 #endif
834 
835  const ScalarType ONE = SCT::one();
836  const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
837 
838  int xc = MVT::GetNumberVecs( X );
839 
840  if (howMany == -1) {
841  howMany = xc;
842  }
843 
844  /*******************************************************
845  * If _hasOp == false, we will not reference MX below *
846  *******************************************************/
847  TEUCHOS_TEST_FOR_EXCEPTION(this->_hasOp == true && MX == Teuchos::null, std::logic_error,
848  "Anasazi::BasicOrthoManager::findBasis(): calling routine did not specify MS.");
849  TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::logic_error,
850  "Anasazi::BasicOrthoManager::findBasis(): Invalid howMany parameter" );
851 
852  /* xstart is which column we are starting the process with, based on howMany
853  * columns before xstart are assumed to be Op-orthonormal already
854  */
855  int xstart = xc - howMany;
856 
857  for (int j = xstart; j < xc; j++) {
858 
859  // numX represents the number of currently orthonormal columns of X
860  int numX = j;
861  // j represents the index of the current column of X
862  // these are different interpretations of the same value
863 
864  //
865  // set the lower triangular part of B to zero
866  for (int i=j+1; i<xc; ++i) {
867  B(i,j) = ZERO;
868  }
869 
870  // Get a view of the vector currently being worked on.
871  std::vector<int> index(1);
872  index[0] = j;
873  Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
874  Teuchos::RCP<MV> MXj;
875  if ((this->_hasOp)) {
876  // MXj is a view of the current vector in MX
877  MXj = MVT::CloneViewNonConst( *MX, index );
878  }
879  else {
880  // MXj is a pointer to Xj, and MUST NOT be modified
881  MXj = Xj;
882  }
883 
884  // Get a view of the previous vectors.
885  std::vector<int> prev_idx( numX );
886  Teuchos::RCP<const MV> prevX, prevMX;
887 
888  if (numX > 0) {
889  for (int i=0; i<numX; ++i) prev_idx[i] = i;
890  prevX = MVT::CloneViewNonConst( X, prev_idx );
891  if (this->_hasOp) {
892  prevMX = MVT::CloneViewNonConst( *MX, prev_idx );
893  }
894  }
895 
896  bool rankDef = true;
897  /* numTrials>0 will denote that the current vector was randomized for the purpose
898  * of finding a basis vector, and that the coefficients of that vector should
899  * not be stored in B
900  */
901  for (int numTrials = 0; numTrials < 10; numTrials++) {
902 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
903  *out << "Trial " << numTrials << " for vector " << j << "\n";
904 #endif
905 
906  // Make storage for these Gram-Schmidt iterations.
907  Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
908  std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1);
909 
910  //
911  // Save old MXj vector and compute Op-norm
912  //
913  Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj );
915 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
916  *out << "origNorm = " << origNorm[0] << "\n";
917 #endif
918 
919  if (numX > 0) {
920  // Apply the first step of Gram-Schmidt
921 
922  // product <- prevX^T MXj
923  MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*prevX,*Xj,product,Teuchos::null,MXj);
924 
925  // Xj <- Xj - prevX prevX^T MXj
926  // = Xj - prevX product
927 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
928  *out << "Orthogonalizing X[" << j << "]...\n";
929 #endif
930  MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
931 
932  // Update MXj
933  if (this->_hasOp) {
934  // MXj <- Op*Xj_new
935  // = Op*(Xj_old - prevX prevX^T MXj)
936  // = MXj - prevMX product
937 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
938  *out << "Updating MX[" << j << "]...\n";
939 #endif
940  MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
941  }
942 
943  // Compute new Op-norm
945  MagnitudeType product_norm = product.normOne();
946 
947 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
948  *out << "newNorm = " << newNorm[0] << "\n";
949  *out << "prodoct_norm = " << product_norm << "\n";
950 #endif
951 
952  // Check if a correction is needed.
953  if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) {
954 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
955  if (product_norm/newNorm[0] >= tol_) {
956  *out << "product_norm/newNorm == " << product_norm/newNorm[0] << "... another step of Gram-Schmidt.\n";
957  }
958  else {
959  *out << "eps*origNorm == " << eps_*origNorm[0] << "... another step of Gram-Schmidt.\n";
960  }
961 #endif
962  // Apply the second step of Gram-Schmidt
963  // This is the same as above
964  Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
965  MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*prevX,*Xj,P2,Teuchos::null,MXj);
966  product += P2;
967 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
968  *out << "Orthogonalizing X[" << j << "]...\n";
969 #endif
970  MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
971  if ((this->_hasOp)) {
972 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
973  *out << "Updating MX[" << j << "]...\n";
974 #endif
975  MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
976  }
977  // Compute new Op-norms
979  product_norm = P2.normOne();
980 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
981  *out << "newNorm2 = " << newNorm2[0] << "\n";
982  *out << "product_norm = " << product_norm << "\n";
983 #endif
984  if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) {
985  // we don't do another GS, we just set it to zero.
986 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
987  if (product_norm/newNorm2[0] >= tol_) {
988  *out << "product_norm/newNorm2 == " << product_norm/newNorm2[0] << "... setting vector to zero.\n";
989  }
990  else if (newNorm[0] < newNorm2[0]) {
991  *out << "newNorm2 > newNorm... setting vector to zero.\n";
992  }
993  else {
994  *out << "eps*origNorm == " << eps_*origNorm[0] << "... setting vector to zero.\n";
995  }
996 #endif
997  MVT::MvInit(*Xj,ZERO);
998  if ((this->_hasOp)) {
999 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
1000  *out << "Setting MX[" << j << "] to zero as well.\n";
1001 #endif
1002  MVT::MvInit(*MXj,ZERO);
1003  }
1004  }
1005  }
1006  } // if (numX > 0) do GS
1007 
1008  // save the coefficients, if we are working on the original vector and not a randomly generated one
1009  if (numTrials == 0) {
1010  for (int i=0; i<numX; i++) {
1011  B(i,j) = product(i,0);
1012  }
1013  }
1014 
1015  // Check if Xj has any directional information left after the orthogonalization.
1017  if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) {
1018 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
1019  *out << "Normalizing X[" << j << "], norm(X[" << j << "]) = " << newNorm[0] << "\n";
1020 #endif
1021  // Normalize Xj.
1022  // Xj <- Xj / norm
1023  MVT::MvScale( *Xj, ONE/newNorm[0]);
1024  if (this->_hasOp) {
1025 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
1026  *out << "Normalizing M*X[" << j << "]...\n";
1027 #endif
1028  // Update MXj.
1029  MVT::MvScale( *MXj, ONE/newNorm[0]);
1030  }
1031 
1032  // save it, if it corresponds to the original vector and not a randomly generated one
1033  if (numTrials == 0) {
1034  B(j,j) = newNorm[0];
1035  }
1036 
1037  // We are not rank deficient in this vector. Move on to the next vector in X.
1038  rankDef = false;
1039  break;
1040  }
1041  else {
1042 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
1043  *out << "Not normalizing M*X[" << j << "]...\n";
1044 #endif
1045  // There was nothing left in Xj after orthogonalizing against previous columns in X.
1046  // X is rank deficient.
1047  // reflect this in the coefficients
1048  B(j,j) = ZERO;
1049 
1050  if (completeBasis) {
1051  // Fill it with random information and keep going.
1052 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
1053  *out << "Inserting random vector in X[" << j << "]...\n";
1054 #endif
1055  MVT::MvRandom( *Xj );
1056  if (this->_hasOp) {
1057 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
1058  *out << "Updating M*X[" << j << "]...\n";
1059 #endif
1060  OPT::Apply( *(this->_Op), *Xj, *MXj );
1061  this->_OpCounter += MVT::GetNumberVecs(*Xj);
1062  }
1063  }
1064  else {
1065  rankDef = true;
1066  break;
1067  }
1068  }
1069  } // for (numTrials = 0; numTrials < 10; ++numTrials)
1070 
1071  // if rankDef == true, then quit and notify user of rank obtained
1072  if (rankDef == true) {
1073  TEUCHOS_TEST_FOR_EXCEPTION( rankDef && completeBasis, OrthoError,
1074  "Anasazi::BasicOrthoManager::findBasis(): Unable to complete basis" );
1075 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
1076  *out << "Returning early, rank " << j << " from Anasazi::BasicOrthoManager::findBasis(...)\n";
1077 #endif
1078  return j;
1079  }
1080 
1081  } // for (j = 0; j < xc; ++j)
1082 
1083 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
1084  *out << "Returning " << xc << " from Anasazi::BasicOrthoManager::findBasis(...)\n";
1085 #endif
1086  return xc;
1087  }
1088 
1089 } // namespace Anasazi
1090 
1091 #endif // ANASAZI_BASIC_ORTHOMANAGER_HPP
1092 
Teuchos::ScalarTraits< ScalarType >::magnitudeType getKappa() const
Return parameter for re-orthogonalization threshold.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
int projectAndNormalizeMat(MV &X, Teuchos::Array< Teuchos::RCP< const MV > > Q, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C=Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)), Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B=Teuchos::null, Teuchos::RCP< MV > MX=Teuchos::null, Teuchos::Array< Teuchos::RCP< const MV > > MQ=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null))) const
Given a set of bases Q[i] and a multivector X, this method computes an orthonormal basis for ...
Declaration of basic traits for the multivector type.
Virtual base class which defines basic traits for the operator type.
void projectMat(MV &X, Teuchos::Array< Teuchos::RCP< const MV > > Q, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C=Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)), Teuchos::RCP< MV > MX=Teuchos::null, Teuchos::Array< Teuchos::RCP< const MV > > MQ=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null))) const
Given a list of mutually orthogonal and internally orthonormal bases Q, this method projects a multiv...
int normalizeMat(MV &X, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B=Teuchos::null, Teuchos::RCP< MV > MX=Teuchos::null) const
This method takes a multivector X and attempts to compute an orthonormal basis for ...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
void setKappa(typename Teuchos::ScalarTraits< ScalarType >::magnitudeType kappa)
Set parameter for re-orthogonalization threshold.
void normMat(const MV &X, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, Teuchos::RCP< const MV > MX=Teuchos::null) const
Provides the norm induced by the matrix-based inner product.
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
void innerProdMat(const MV &X, const MV &Y, Teuchos::SerialDenseMatrix< int, ScalarType > &Z, Teuchos::RCP< const MV > MX=Teuchos::null, Teuchos::RCP< const MV > MY=Teuchos::null) const
Provides a matrix-based inner product.
Exception thrown to signal error in an orthogonalization manager method.
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormErrorMat(const MV &X, Teuchos::RCP< const MV > MX=Teuchos::null) const
This method computes the error in orthonormality of a multivector, measured as the Frobenius norm of ...
BasicOrthoManager(Teuchos::RCP< const OP > Op=Teuchos::null, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType kappa=1.41421356, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType eps=0.0, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType tol=0.20)
Constructor specifying re-orthogonalization tolerance.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP< const MV > MX1, Teuchos::RCP< const MV > MX2) const
This method computes the error in orthogonality of two multivectors, measured as the Frobenius norm o...