Anasazi  Version of the Day
AnasaziSVQBOrthoManager.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 
48 #ifndef ANASAZI_SVQB_ORTHOMANAGER_HPP
49 #define ANASAZI_SVQB_ORTHOMANAGER_HPP
50 
60 #include "AnasaziConfigDefs.hpp"
61 #include "AnasaziMultiVecTraits.hpp"
62 #include "AnasaziOperatorTraits.hpp"
63 #include "AnasaziMatOrthoManager.hpp"
64 #include "Teuchos_LAPACK.hpp"
65 
66 namespace Anasazi {
67 
68  template<class ScalarType, class MV, class OP>
69  class SVQBOrthoManager : public MatOrthoManager<ScalarType,MV,OP> {
70 
71  private:
72  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
73  typedef Teuchos::ScalarTraits<ScalarType> SCT;
74  typedef Teuchos::ScalarTraits<MagnitudeType> SCTM;
75  typedef MultiVecTraits<ScalarType,MV> MVT;
76  typedef OperatorTraits<ScalarType,MV,OP> OPT;
77  std::string dbgstr;
78 
79 
80  public:
81 
83 
84  SVQBOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null, bool debug = false );
86 
87 
89  ~SVQBOrthoManager() {};
91 
92 
93 
95 
96 
97 
136  void projectMat (
137  MV &X,
138  Teuchos::Array<Teuchos::RCP<const MV> > Q,
139  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
140  = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
141  Teuchos::RCP<MV> MX = Teuchos::null,
142  Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
143  ) const;
144 
145 
184  int normalizeMat (
185  MV &X,
186  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
187  Teuchos::RCP<MV> MX = Teuchos::null
188  ) const;
189 
190 
250  int projectAndNormalizeMat (
251  MV &X,
252  Teuchos::Array<Teuchos::RCP<const MV> > Q,
253  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
254  = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
255  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
256  Teuchos::RCP<MV> MX = Teuchos::null,
257  Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
258  ) const;
259 
261 
263 
264 
269  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
270  orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const;
271 
276  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
277  orthogErrorMat(
278  const MV &X,
279  const MV &Y,
280  Teuchos::RCP<const MV> MX = Teuchos::null,
281  Teuchos::RCP<const MV> MY = Teuchos::null
282  ) const;
283 
285 
286  private:
287 
288  MagnitudeType eps_;
289  bool debug_;
290 
291  // ! Routine to find an orthogonal/orthonormal basis
292  int findBasis(MV &X, Teuchos::RCP<MV> MX,
293  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
294  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
295  Teuchos::Array<Teuchos::RCP<const MV> > Q,
296  bool normalize_in ) const;
297  };
298 
299 
301  // Constructor
302  template<class ScalarType, class MV, class OP>
303  SVQBOrthoManager<ScalarType,MV,OP>::SVQBOrthoManager( Teuchos::RCP<const OP> Op, bool debug)
304  : MatOrthoManager<ScalarType,MV,OP>(Op), dbgstr(" *** "), debug_(debug) {
305 
306  Teuchos::LAPACK<int,MagnitudeType> lapack;
307  eps_ = lapack.LAMCH('E');
308  if (debug_) {
309  std::cout << "eps_ == " << eps_ << std::endl;
310  }
311  }
312 
313 
315  // Compute the distance from orthonormality
316  template<class ScalarType, class MV, class OP>
317  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
318  SVQBOrthoManager<ScalarType,MV,OP>::orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX) const {
319  const ScalarType ONE = SCT::one();
320  int rank = MVT::GetNumberVecs(X);
321  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
322  MatOrthoManager<ScalarType,MV,OP>::innerProdMat(X,X,xTx,MX,MX);
323  for (int i=0; i<rank; i++) {
324  xTx(i,i) -= ONE;
325  }
326  return xTx.normFrobenius();
327  }
328 
330  // Compute the distance from orthogonality
331  template<class ScalarType, class MV, class OP>
332  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
333  SVQBOrthoManager<ScalarType,MV,OP>::orthogErrorMat(
334  const MV &X,
335  const MV &Y,
336  Teuchos::RCP<const MV> MX,
337  Teuchos::RCP<const MV> MY
338  ) const {
339  int r1 = MVT::GetNumberVecs(X);
340  int r2 = MVT::GetNumberVecs(Y);
341  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
342  MatOrthoManager<ScalarType,MV,OP>::innerProdMat(X,Y,xTx,MX,MY);
343  return xTx.normFrobenius();
344  }
345 
346 
347 
349  template<class ScalarType, class MV, class OP>
350  void SVQBOrthoManager<ScalarType, MV, OP>::projectMat(
351  MV &X,
352  Teuchos::Array<Teuchos::RCP<const MV> > Q,
353  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
354  Teuchos::RCP<MV> MX,
355  Teuchos::Array<Teuchos::RCP<const MV> > MQ) const {
356  (void)MQ;
357  findBasis(X,MX,C,Teuchos::null,Q,false);
358  }
359 
360 
361 
363  // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
364  template<class ScalarType, class MV, class OP>
365  int SVQBOrthoManager<ScalarType, MV, OP>::normalizeMat(
366  MV &X,
367  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
368  Teuchos::RCP<MV> MX) const {
369  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C;
370  Teuchos::Array<Teuchos::RCP<const MV> > Q;
371  return findBasis(X,MX,C,B,Q,true);
372  }
373 
374 
375 
377  // Find an Op-orthonormal basis for span(X) - span(W)
378  template<class ScalarType, class MV, class OP>
379  int SVQBOrthoManager<ScalarType, MV, OP>::projectAndNormalizeMat(
380  MV &X,
381  Teuchos::Array<Teuchos::RCP<const MV> > Q,
382  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
383  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
384  Teuchos::RCP<MV> MX,
385  Teuchos::Array<Teuchos::RCP<const MV> > MQ) const {
386  (void)MQ;
387  return findBasis(X,MX,C,B,Q,true);
388  }
389 
390 
391 
392 
394  // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that
395  // the rank is numvectors(X)
396  //
397  // Tracking the coefficients (C[i] and B) for this code is complicated by the fact that the loop
398  // structure looks like
399  // do
400  // project
401  // do
402  // ortho
403  // end
404  // end
405  // However, the recurrence for the coefficients is not complicated:
406  // B = I
407  // C = 0
408  // do
409  // project yields newC
410  // C = C + newC*B
411  // do
412  // ortho yields newR
413  // B = newR*B
414  // end
415  // end
416  // This holds for each individual C[i] (which correspond to the list of bases we are orthogonalizing
417  // against).
418  //
419  template<class ScalarType, class MV, class OP>
420  int SVQBOrthoManager<ScalarType, MV, OP>::findBasis(
421  MV &X, Teuchos::RCP<MV> MX,
422  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
423  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
424  Teuchos::Array<Teuchos::RCP<const MV> > Q,
425  bool normalize_in) const {
426 
427  const ScalarType ONE = SCT::one();
428  const MagnitudeType MONE = SCTM::one();
429  const MagnitudeType ZERO = SCTM::zero();
430 
431  int numGS = 0,
432  numSVQB = 0,
433  numRand = 0;
434 
435  // get sizes of X,MX
436  int xc = MVT::GetNumberVecs(X);
437  ptrdiff_t xr = MVT::GetGlobalLength( X );
438 
439  // get sizes of Q[i]
440  int nq = Q.length();
441  ptrdiff_t qr = (nq == 0) ? 0 : MVT::GetGlobalLength(*Q[0]);
442  ptrdiff_t qsize = 0;
443  std::vector<int> qcs(nq);
444  for (int i=0; i<nq; i++) {
445  qcs[i] = MVT::GetNumberVecs(*Q[i]);
446  qsize += qcs[i];
447  }
448 
449  if (normalize_in == true && qsize + xc > xr) {
450  // not well-posed
451  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
452  "Anasazi::SVQBOrthoManager::findBasis(): Orthogonalization constraints not feasible" );
453  }
454 
455  // try to short-circuit as early as possible
456  if (normalize_in == false && (qsize == 0 || xc == 0)) {
457  // nothing to do
458  return 0;
459  }
460  else if (normalize_in == true && (xc == 0 || xr == 0)) {
461  // normalize requires X not empty
462  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
463  "Anasazi::SVQBOrthoManager::findBasis(): X must be non-empty" );
464  }
465 
466  // check that Q matches X
467  TEUCHOS_TEST_FOR_EXCEPTION( qsize != 0 && qr != xr , std::invalid_argument,
468  "Anasazi::SVQBOrthoManager::findBasis(): Size of X not consistant with size of Q" );
469 
470  /* If we don't have enough C, expanding it creates null references
471  * If we have too many, resizing just throws away the later ones
472  * If we have exactly as many as we have Q, this call has no effect
473  */
474  C.resize(nq);
475  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > newC(nq);
476  // check the size of the C[i] against the Q[i] and consistency between Q[i]
477  for (int i=0; i<nq; i++) {
478  // check size of Q[i]
479  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
480  "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not mutually consistant" );
481  TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
482  "Anasazi::SVQBOrthoManager::findBasis(): Q has less rows than columns" );
483  // check size of C[i]
484  if ( C[i] == Teuchos::null ) {
485  C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
486  }
487  else {
488  TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc, std::invalid_argument,
489  "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not consistant with C" );
490  }
491  // clear C[i]
492  C[i]->putScalar(ZERO);
493  newC[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(C[i]->numRows(),C[i]->numCols()) );
494  }
495 
496 
498  // Allocate necessary storage
499  // C were allocated above
500  // Allocate MX and B (if necessary)
501  // Set B = I
502  if (normalize_in == true) {
503  if ( B == Teuchos::null ) {
504  B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
505  }
506  TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
507  "Anasazi::SVQBOrthoManager::findBasis(): Size of B not consistant with X" );
508  // set B to I
509  B->putScalar(ZERO);
510  for (int i=0; i<xc; i++) {
511  (*B)(i,i) = MONE;
512  }
513  }
514  /******************************************
515  * If _hasOp == false, DO NOT MODIFY MX *
516  ******************************************
517  * if Op==null, MX == X (via pointer)
518  * Otherwise, either the user passed in MX or we will allocate and compute it
519  *
520  * workX will be a multivector of the same size as X, used to perform X*S when normalizing
521  */
522  Teuchos::RCP<MV> workX;
523  if (normalize_in) {
524  workX = MVT::Clone(X,xc);
525  }
526  if (this->_hasOp) {
527  if (MX == Teuchos::null) {
528  // we need to allocate space for MX
529  MX = MVT::Clone(X,xc);
530  OPT::Apply(*(this->_Op),X,*MX);
531  this->_OpCounter += MVT::GetNumberVecs(X);
532  }
533  }
534  else {
535  MX = Teuchos::rcpFromRef(X);
536  }
537  std::vector<ScalarType> normX(xc), invnormX(xc);
538  Teuchos::SerialDenseMatrix<int,ScalarType> XtMX(xc,xc), workU(1,1);
539  Teuchos::LAPACK<int,ScalarType> lapack;
540  /**********************************************************************
541  * allocate storage for eigenvectors,eigenvalues of X^T Op X, and for
542  * the work space needed to compute this xc-by-xc eigendecomposition
543  **********************************************************************/
544  std::vector<ScalarType> work;
545  std::vector<MagnitudeType> lambda, lambdahi, rwork;
546  if (normalize_in) {
547  // get size of work from ILAENV
548  int lwork = lapack.ILAENV(1,"hetrd","VU",xc,-1,-1,-1);
549  // lwork >= (nb+1)*n for complex
550  // lwork >= (nb+2)*n for real
551  TEUCHOS_TEST_FOR_EXCEPTION( lwork < 0, OrthoError,
552  "Anasazi::SVQBOrthoManager::findBasis(): Error code from ILAENV" );
553 
554  lwork = (lwork+2)*xc;
555  work.resize(lwork);
556  // size of rwork is max(1,3*xc-2)
557  lwork = (3*xc-2 > 1) ? 3*xc - 2 : 1;
558  rwork.resize(lwork);
559  // size of lambda is xc
560  lambda.resize(xc);
561  lambdahi.resize(xc);
562  workU.reshape(xc,xc);
563  }
564 
565  // test sizes of X,MX
566  int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
567  ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
568  TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
569  "Anasazi::SVQBOrthoManager::findBasis(): Size of X not consistant with MX" );
570 
571  // sentinel to continue the outer loop (perform another projection step)
572  bool doGramSchmidt = true;
573  // variable for testing orthonorm/orthog
574  MagnitudeType tolerance = MONE/SCTM::squareroot(eps_);
575 
576  // outer loop
577  while (doGramSchmidt) {
578 
580  // perform projection
581  if (qsize > 0) {
582 
583  numGS++;
584 
585  // Compute the norms of the vectors
586  {
587  std::vector<MagnitudeType> normX_mag(xc);
588  MatOrthoManager<ScalarType,MV,OP>::normMat(X,normX_mag,MX);
589  for (int i=0; i<xc; ++i) {
590  normX[i] = normX_mag[i];
591  invnormX[i] = (normX_mag[i] == ZERO) ? ZERO : MONE/normX_mag[i];
592  }
593  }
594  // normalize the vectors
595  MVT::MvScale(X,invnormX);
596  if (this->_hasOp) {
597  MVT::MvScale(*MX,invnormX);
598  }
599  // check that vectors are normalized now
600  if (debug_) {
601  std::vector<MagnitudeType> nrm2(xc);
602  std::cout << dbgstr << "max post-scale norm: (with/without MX) : ";
603  MagnitudeType maxpsnw = ZERO, maxpsnwo = ZERO;
604  MatOrthoManager<ScalarType,MV,OP>::normMat(X,nrm2,MX);
605  for (int i=0; i<xc; i++) {
606  maxpsnw = (nrm2[i] > maxpsnw ? nrm2[i] : maxpsnw);
607  }
608  this->norm(X,nrm2);
609  for (int i=0; i<xc; i++) {
610  maxpsnwo = (nrm2[i] > maxpsnwo ? nrm2[i] : maxpsnwo);
611  }
612  std::cout << "(" << maxpsnw << "," << maxpsnwo << ")" << std::endl;
613  }
614  // project the vectors onto the Qi
615  for (int i=0; i<nq; i++) {
616  MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*Q[i],X,*newC[i],Teuchos::null,MX);
617  }
618  // remove the components in Qi from X
619  for (int i=0; i<nq; i++) {
620  MVT::MvTimesMatAddMv(-ONE,*Q[i],*newC[i],ONE,X);
621  }
622  // un-scale the vectors
623  MVT::MvScale(X,normX);
624  // Recompute the vectors in MX
625  if (this->_hasOp) {
626  OPT::Apply(*(this->_Op),X,*MX);
627  this->_OpCounter += MVT::GetNumberVecs(X);
628  }
629 
630  //
631  // Compute largest column norm of
632  // ( C[0] )
633  // C = ( .... )
634  // ( C[nq-1] )
635  MagnitudeType maxNorm = ZERO;
636  for (int j=0; j<xc; j++) {
637  MagnitudeType sum = ZERO;
638  for (int k=0; k<nq; k++) {
639  for (int i=0; i<qcs[k]; i++) {
640  sum += SCT::magnitude((*newC[k])(i,j))*SCT::magnitude((*newC[k])(i,j));
641  }
642  }
643  maxNorm = (sum > maxNorm) ? sum : maxNorm;
644  }
645 
646  // do we perform another GS?
647  if (maxNorm < 0.36) {
648  doGramSchmidt = false;
649  }
650 
651  // unscale newC to reflect the scaling of X
652  for (int k=0; k<nq; k++) {
653  for (int j=0; j<xc; j++) {
654  for (int i=0; i<qcs[k]; i++) {
655  (*newC[k])(i,j) *= normX[j];
656  }
657  }
658  }
659  // accumulate into C
660  if (normalize_in) {
661  // we are normalizing
662  int info;
663  for (int i=0; i<nq; i++) {
664  info = C[i]->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,*newC[i],*B,ONE);
665  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::SVQBOrthoManager::findBasis(): Input error to SerialDenseMatrix::multiply.");
666  }
667  }
668  else {
669  // not normalizing
670  for (int i=0; i<nq; i++) {
671  (*C[i]) += *newC[i];
672  }
673  }
674  }
675  else { // qsize == 0... don't perform projection
676  // don't do any more outer loops; all we need is to call the normalize code below
677  doGramSchmidt = false;
678  }
679 
680 
682  // perform normalization
683  if (normalize_in) {
684 
685  MagnitudeType condT = tolerance;
686 
687  while (condT >= tolerance) {
688 
689  numSVQB++;
690 
691  // compute X^T Op X
692  MatOrthoManager<ScalarType,MV,OP>::innerProdMat(X,X,XtMX,MX,MX);
693 
694  // compute scaling matrix for XtMX: D^{.5} and D^{-.5} (D-half and D-half-inv)
695  std::vector<MagnitudeType> Dh(xc), Dhi(xc);
696  for (int i=0; i<xc; i++) {
697  Dh[i] = SCT::magnitude(SCT::squareroot(XtMX(i,i)));
698  Dhi[i] = (Dh[i] == ZERO ? ZERO : MONE/Dh[i]);
699  }
700  // scale XtMX : S = D^{-.5} * XtMX * D^{-.5}
701  for (int i=0; i<xc; i++) {
702  for (int j=0; j<xc; j++) {
703  XtMX(i,j) *= Dhi[i]*Dhi[j];
704  }
705  }
706 
707  // compute the eigenvalue decomposition of S=U*Lambda*U^T (using upper part)
708  int info;
709  lapack.HEEV('V', 'U', xc, XtMX.values(), XtMX.stride(), &lambda[0], &work[0], work.size(), &rwork[0], &info);
710  std::stringstream os;
711  os << "Anasazi::SVQBOrthoManager::findBasis(): Error code " << info << " from HEEV";
712  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OrthoError,
713  os.str() );
714  if (debug_) {
715  std::cout << dbgstr << "eigenvalues of XtMX: (";
716  for (int i=0; i<xc-1; i++) {
717  std::cout << lambda[i] << ",";
718  }
719  std::cout << lambda[xc-1] << ")" << std::endl;
720  }
721 
722  // remember, HEEV orders the eigenvalues from smallest to largest
723  // examine condition number of Lambda, compute Lambda^{-.5}
724  MagnitudeType maxLambda = lambda[xc-1],
725  minLambda = lambda[0];
726  int iZeroMax = -1;
727  for (int i=0; i<xc; i++) {
728  if (lambda[i] <= 10*eps_*maxLambda) { // finish: this was eps_*eps_*maxLambda
729  iZeroMax = i;
730  lambda[i] = ZERO;
731  lambdahi[i] = ZERO;
732  }
733  /*
734  else if (lambda[i] < eps_*maxLambda) {
735  lambda[i] = SCTM::squareroot(eps_*maxLambda);
736  lambdahi[i] = MONE/lambda[i];
737  }
738  */
739  else {
740  lambda[i] = SCTM::squareroot(lambda[i]);
741  lambdahi[i] = MONE/lambda[i];
742  }
743  }
744 
745  // compute X * D^{-.5} * U * Lambda^{-.5} and new Op*X
746  //
747  // copy X into workX
748  std::vector<int> ind(xc);
749  for (int i=0; i<xc; i++) {ind[i] = i;}
750  MVT::SetBlock(X,ind,*workX);
751  //
752  // compute D^{-.5}*U*Lambda^{-.5} into workU
753  workU.assign(XtMX);
754  for (int j=0; j<xc; j++) {
755  for (int i=0; i<xc; i++) {
756  workU(i,j) *= Dhi[i]*lambdahi[j];
757  }
758  }
759  // compute workX * workU into X
760  MVT::MvTimesMatAddMv(ONE,*workX,workU,ZERO,X);
761  //
762  // note, it seems important to apply Op exactly for large condition numbers.
763  // for small condition numbers, we can update MX "implicitly"
764  // this trick reduces the number of applications of Op
765  if (this->_hasOp) {
766  if (maxLambda >= tolerance * minLambda) {
767  // explicit update of MX
768  OPT::Apply(*(this->_Op),X,*MX);
769  this->_OpCounter += MVT::GetNumberVecs(X);
770  }
771  else {
772  // implicit update of MX
773  // copy MX into workX
774  MVT::SetBlock(*MX,ind,*workX);
775  //
776  // compute workX * workU into MX
777  MVT::MvTimesMatAddMv(ONE,*workX,workU,ZERO,*MX);
778  }
779  }
780 
781  // accumulate new B into previous B
782  // B = Lh * U^H * Dh * B
783  for (int j=0; j<xc; j++) {
784  for (int i=0; i<xc; i++) {
785  workU(i,j) = Dh[i] * (*B)(i,j);
786  }
787  }
788  info = B->multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,XtMX,workU,ZERO);
789  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::SVQBOrthoManager::findBasis(): Input error to SerialDenseMatrix::multiply.");
790  for (int j=0; j<xc ;j++) {
791  for (int i=0; i<xc; i++) {
792  (*B)(i,j) *= lambda[i];
793  }
794  }
795 
796  // check iZeroMax (rank indicator)
797  if (iZeroMax >= 0) {
798  if (debug_) {
799  std::cout << dbgstr << "augmenting multivec with " << iZeroMax+1 << " random directions" << std::endl;
800  }
801 
802  numRand++;
803  // put random info in the first iZeroMax+1 vectors of X,MX
804  std::vector<int> ind2(iZeroMax+1);
805  for (int i=0; i<iZeroMax+1; i++) {
806  ind2[i] = i;
807  }
808  Teuchos::RCP<MV> Xnull,MXnull;
809  Xnull = MVT::CloneViewNonConst(X,ind2);
810  MVT::MvRandom(*Xnull);
811  if (this->_hasOp) {
812  MXnull = MVT::CloneViewNonConst(*MX,ind2);
813  OPT::Apply(*(this->_Op),*Xnull,*MXnull);
814  this->_OpCounter += MVT::GetNumberVecs(*Xnull);
815  MXnull = Teuchos::null;
816  }
817  Xnull = Teuchos::null;
818  condT = tolerance;
819  doGramSchmidt = true;
820  break; // break from while(condT > tolerance)
821  }
822 
823  condT = SCTM::magnitude(maxLambda / minLambda);
824  if (debug_) {
825  std::cout << dbgstr << "condT: " << condT << ", tolerance = " << tolerance << std::endl;
826  }
827 
828  // if multiple passes of SVQB are necessary, then pass through outer GS loop again
829  if ((doGramSchmidt == false) && (condT > SCTM::squareroot(tolerance))) {
830  doGramSchmidt = true;
831  }
832 
833  } // end while (condT >= tolerance)
834  }
835  // end if(normalize_in)
836 
837  } // end while (doGramSchmidt)
838 
839  if (debug_) {
840  std::cout << dbgstr << "(numGS,numSVQB,numRand) : "
841  << "(" << numGS
842  << "," << numSVQB
843  << "," << numRand
844  << ")" << std::endl;
845  }
846  return xc;
847  }
848 
849 
850 } // namespace Anasazi
851 
852 #endif // ANASAZI_SVQB_ORTHOMANAGER_HPP
853 
AnasaziMultiVecTraits.hpp
Declaration of basic traits for the multivector type.
Anasazi::MultiVecTraits
Traits class which defines basic operations on multivectors.
Definition: AnasaziMultiVecTraits.hpp:127
Anasazi::SVQBOrthoManager::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: AnasaziSVQBOrthoManager.hpp:318
AnasaziMatOrthoManager.hpp
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
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:422
Anasazi::SVQBOrthoManager::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: AnasaziSVQBOrthoManager.hpp:379
Anasazi::SVQBOrthoManager
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
Definition: AnasaziSVQBOrthoManager.hpp:69
Anasazi::MatOrthoManager
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
Definition: AnasaziMatOrthoManager.hpp:76
AnasaziOperatorTraits.hpp
Virtual base class which defines basic traits for the operator type.
Anasazi::OperatorTraits
Virtual base class which defines basic traits for the operator type.
Definition: AnasaziOperatorTraits.hpp:84
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:383
Anasazi::SVQBOrthoManager::orthogErrorMat
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthogErrorMat(const MV &X, const MV &Y, Teuchos::RCP< const MV > MX=Teuchos::null, Teuchos::RCP< const MV > MY=Teuchos::null) const
This method computes the error in orthogonality of two multivectors, measured as the Frobenius norm o...
Definition: AnasaziSVQBOrthoManager.hpp:333
Anasazi::SVQBOrthoManager::SVQBOrthoManager
SVQBOrthoManager(Teuchos::RCP< const OP > Op=Teuchos::null, bool debug=false)
Constructor specifying re-orthogonalization tolerance.
Definition: AnasaziSVQBOrthoManager.hpp:303
Anasazi
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
AnasaziConfigDefs.hpp
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Anasazi::SVQBOrthoManager::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: AnasaziSVQBOrthoManager.hpp:365
Anasazi::SVQBOrthoManager::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: AnasaziSVQBOrthoManager.hpp:350
Anasazi::SVQBOrthoManager::~SVQBOrthoManager
~SVQBOrthoManager()
Destructor.
Definition: AnasaziSVQBOrthoManager.hpp:89
Generated by   doxygen 1.8.16