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