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