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 
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::ScalarTraits< ScalarType >::magnitudeType getKappa() const
Return parameter for re-orthogonalization threshold.
Declaration of basic traits for the multivector type.
Virtual base class which defines basic traits for the operator type.
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.
static void MvDot(const MV &mv, const MV &A, std::vector< ScalarType > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
static void Apply(const OP &Op, const MV &x, MV &y)
Application method which performs operation y = Op*x. An OperatorError exception is thrown if there i...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
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...
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...
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.
Traits class which defines basic operations on multivectors.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
Virtual base class which defines basic traits for the operator type.
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
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 ...
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
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 ...
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
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)...
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 ...
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.
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.