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, 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 
35 #ifndef ANASAZI_SVQB_ORTHOMANAGER_HPP
36 #define ANASAZI_SVQB_ORTHOMANAGER_HPP
37 
47 #include "AnasaziConfigDefs.hpp"
51 #include "Teuchos_LAPACK.hpp"
52 
53 namespace Anasazi {
54 
55  template<class ScalarType, class MV, class OP>
56  class SVQBOrthoManager : public MatOrthoManager<ScalarType,MV,OP> {
57 
58  private:
59  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
60  typedef Teuchos::ScalarTraits<ScalarType> SCT;
61  typedef Teuchos::ScalarTraits<MagnitudeType> SCTM;
64  std::string dbgstr;
65 
66 
67  public:
68 
70 
71  SVQBOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null, bool debug = false );
73 
74 
78 
79 
80 
82 
83 
84 
123  void projectMat (
124  MV &X,
125  Teuchos::Array<Teuchos::RCP<const MV> > Q,
126  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
127  = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
128  Teuchos::RCP<MV> MX = Teuchos::null,
129  Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
130  ) const;
131 
132 
171  int normalizeMat (
172  MV &X,
173  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
174  Teuchos::RCP<MV> MX = Teuchos::null
175  ) const;
176 
177 
238  MV &X,
239  Teuchos::Array<Teuchos::RCP<const MV> > Q,
240  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
241  = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
242  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
243  Teuchos::RCP<MV> MX = Teuchos::null,
244  Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
245  ) const;
246 
248 
250 
251 
256  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
257  orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const;
258 
263  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
265  const MV &X,
266  const MV &Y,
267  Teuchos::RCP<const MV> MX = Teuchos::null,
268  Teuchos::RCP<const MV> MY = Teuchos::null
269  ) const;
270 
272 
273  private:
274 
275  MagnitudeType eps_;
276  bool debug_;
277 
278  // ! Routine to find an orthogonal/orthonormal basis
279  int findBasis(MV &X, Teuchos::RCP<MV> MX,
280  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
281  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
282  Teuchos::Array<Teuchos::RCP<const MV> > Q,
283  bool normalize_in ) const;
284  };
285 
286 
288  // Constructor
289  template<class ScalarType, class MV, class OP>
290  SVQBOrthoManager<ScalarType,MV,OP>::SVQBOrthoManager( Teuchos::RCP<const OP> Op, bool debug)
291  : MatOrthoManager<ScalarType,MV,OP>(Op), dbgstr(" *** "), debug_(debug) {
292 
293  Teuchos::LAPACK<int,MagnitudeType> lapack;
294  eps_ = lapack.LAMCH('E');
295  if (debug_) {
296  std::cout << "eps_ == " << eps_ << std::endl;
297  }
298  }
299 
300 
302  // Compute the distance from orthonormality
303  template<class ScalarType, class MV, class OP>
304  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
305  SVQBOrthoManager<ScalarType,MV,OP>::orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX) const {
306  const ScalarType ONE = SCT::one();
307  int rank = MVT::GetNumberVecs(X);
308  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
310  for (int i=0; i<rank; i++) {
311  xTx(i,i) -= ONE;
312  }
313  return xTx.normFrobenius();
314  }
315 
317  // Compute the distance from orthogonality
318  template<class ScalarType, class MV, class OP>
319  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
321  const MV &X,
322  const MV &Y,
323  Teuchos::RCP<const MV> MX,
324  Teuchos::RCP<const MV> MY
325  ) const {
326  int r1 = MVT::GetNumberVecs(X);
327  int r2 = MVT::GetNumberVecs(Y);
328  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
330  return xTx.normFrobenius();
331  }
332 
333 
334 
336  template<class ScalarType, class MV, class OP>
338  MV &X,
339  Teuchos::Array<Teuchos::RCP<const MV> > Q,
340  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
341  Teuchos::RCP<MV> MX,
342  Teuchos::Array<Teuchos::RCP<const MV> > MQ) const {
343  (void)MQ;
344  findBasis(X,MX,C,Teuchos::null,Q,false);
345  }
346 
347 
348 
350  // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
351  template<class ScalarType, class MV, class OP>
353  MV &X,
354  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
355  Teuchos::RCP<MV> MX) const {
356  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C;
357  Teuchos::Array<Teuchos::RCP<const MV> > Q;
358  return findBasis(X,MX,C,B,Q,true);
359  }
360 
361 
362 
364  // Find an Op-orthonormal basis for span(X) - span(W)
365  template<class ScalarType, class MV, class OP>
367  MV &X,
368  Teuchos::Array<Teuchos::RCP<const MV> > Q,
369  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
370  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
371  Teuchos::RCP<MV> MX,
372  Teuchos::Array<Teuchos::RCP<const MV> > MQ) const {
373  (void)MQ;
374  return findBasis(X,MX,C,B,Q,true);
375  }
376 
377 
378 
379 
381  // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that
382  // the rank is numvectors(X)
383  //
384  // Tracking the coefficients (C[i] and B) for this code is complicated by the fact that the loop
385  // structure looks like
386  // do
387  // project
388  // do
389  // ortho
390  // end
391  // end
392  // However, the recurrence for the coefficients is not complicated:
393  // B = I
394  // C = 0
395  // do
396  // project yields newC
397  // C = C + newC*B
398  // do
399  // ortho yields newR
400  // B = newR*B
401  // end
402  // end
403  // This holds for each individual C[i] (which correspond to the list of bases we are orthogonalizing
404  // against).
405  //
406  template<class ScalarType, class MV, class OP>
408  MV &X, Teuchos::RCP<MV> MX,
409  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
410  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
411  Teuchos::Array<Teuchos::RCP<const MV> > Q,
412  bool normalize_in) const {
413 
414  const ScalarType ONE = SCT::one();
415  const MagnitudeType MONE = SCTM::one();
416  const MagnitudeType ZERO = SCTM::zero();
417 
418  int numGS = 0,
419  numSVQB = 0,
420  numRand = 0;
421 
422  // get sizes of X,MX
423  int xc = MVT::GetNumberVecs(X);
424  ptrdiff_t xr = MVT::GetGlobalLength( X );
425 
426  // get sizes of Q[i]
427  int nq = Q.length();
428  ptrdiff_t qr = (nq == 0) ? 0 : MVT::GetGlobalLength(*Q[0]);
429  ptrdiff_t qsize = 0;
430  std::vector<int> qcs(nq);
431  for (int i=0; i<nq; i++) {
432  qcs[i] = MVT::GetNumberVecs(*Q[i]);
433  qsize += qcs[i];
434  }
435 
436  if (normalize_in == true && qsize + xc > xr) {
437  // not well-posed
438  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
439  "Anasazi::SVQBOrthoManager::findBasis(): Orthogonalization constraints not feasible" );
440  }
441 
442  // try to short-circuit as early as possible
443  if (normalize_in == false && (qsize == 0 || xc == 0)) {
444  // nothing to do
445  return 0;
446  }
447  else if (normalize_in == true && (xc == 0 || xr == 0)) {
448  // normalize requires X not empty
449  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
450  "Anasazi::SVQBOrthoManager::findBasis(): X must be non-empty" );
451  }
452 
453  // check that Q matches X
454  TEUCHOS_TEST_FOR_EXCEPTION( qsize != 0 && qr != xr , std::invalid_argument,
455  "Anasazi::SVQBOrthoManager::findBasis(): Size of X not consistant with size of Q" );
456 
457  /* If we don't have enough C, expanding it creates null references
458  * If we have too many, resizing just throws away the later ones
459  * If we have exactly as many as we have Q, this call has no effect
460  */
461  C.resize(nq);
462  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > newC(nq);
463  // check the size of the C[i] against the Q[i] and consistency between Q[i]
464  for (int i=0; i<nq; i++) {
465  // check size of Q[i]
466  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
467  "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not mutually consistant" );
468  TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
469  "Anasazi::SVQBOrthoManager::findBasis(): Q has less rows than columns" );
470  // check size of C[i]
471  if ( C[i] == Teuchos::null ) {
472  C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
473  }
474  else {
475  TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc, std::invalid_argument,
476  "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not consistant with C" );
477  }
478  // clear C[i]
479  C[i]->putScalar(ZERO);
480  newC[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(C[i]->numRows(),C[i]->numCols()) );
481  }
482 
483 
485  // Allocate necessary storage
486  // C were allocated above
487  // Allocate MX and B (if necessary)
488  // Set B = I
489  if (normalize_in == true) {
490  if ( B == Teuchos::null ) {
491  B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
492  }
493  TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
494  "Anasazi::SVQBOrthoManager::findBasis(): Size of B not consistant with X" );
495  // set B to I
496  B->putScalar(ZERO);
497  for (int i=0; i<xc; i++) {
498  (*B)(i,i) = MONE;
499  }
500  }
501  /******************************************
502  * If _hasOp == false, DO NOT MODIFY MX *
503  ******************************************
504  * if Op==null, MX == X (via pointer)
505  * Otherwise, either the user passed in MX or we will allocate and compute it
506  *
507  * workX will be a multivector of the same size as X, used to perform X*S when normalizing
508  */
509  Teuchos::RCP<MV> workX;
510  if (normalize_in) {
511  workX = MVT::Clone(X,xc);
512  }
513  if (this->_hasOp) {
514  if (MX == Teuchos::null) {
515  // we need to allocate space for MX
516  MX = MVT::Clone(X,xc);
517  OPT::Apply(*(this->_Op),X,*MX);
518  this->_OpCounter += MVT::GetNumberVecs(X);
519  }
520  }
521  else {
522  MX = Teuchos::rcpFromRef(X);
523  }
524  std::vector<ScalarType> normX(xc), invnormX(xc);
525  Teuchos::SerialDenseMatrix<int,ScalarType> XtMX(xc,xc), workU(1,1);
526  Teuchos::LAPACK<int,ScalarType> lapack;
527  /**********************************************************************
528  * allocate storage for eigenvectors,eigenvalues of X^T Op X, and for
529  * the work space needed to compute this xc-by-xc eigendecomposition
530  **********************************************************************/
531  std::vector<ScalarType> work;
532  std::vector<MagnitudeType> lambda, lambdahi, rwork;
533  if (normalize_in) {
534  // get size of work from ILAENV
535  int lwork = lapack.ILAENV(1,"hetrd","VU",xc,-1,-1,-1);
536  // lwork >= (nb+1)*n for complex
537  // lwork >= (nb+2)*n for real
538  TEUCHOS_TEST_FOR_EXCEPTION( lwork < 0, OrthoError,
539  "Anasazi::SVQBOrthoManager::findBasis(): Error code from ILAENV" );
540 
541  lwork = (lwork+2)*xc;
542  work.resize(lwork);
543  // size of rwork is max(1,3*xc-2)
544  lwork = (3*xc-2 > 1) ? 3*xc - 2 : 1;
545  rwork.resize(lwork);
546  // size of lambda is xc
547  lambda.resize(xc);
548  lambdahi.resize(xc);
549  workU.reshape(xc,xc);
550  }
551 
552  // test sizes of X,MX
553  int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
554  ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
555  TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
556  "Anasazi::SVQBOrthoManager::findBasis(): Size of X not consistant with MX" );
557 
558  // sentinel to continue the outer loop (perform another projection step)
559  bool doGramSchmidt = true;
560  // variable for testing orthonorm/orthog
561  MagnitudeType tolerance = MONE/SCTM::squareroot(eps_);
562 
563  // outer loop
564  while (doGramSchmidt) {
565 
567  // perform projection
568  if (qsize > 0) {
569 
570  numGS++;
571 
572  // Compute the norms of the vectors
573  {
574  std::vector<MagnitudeType> normX_mag(xc);
576  for (int i=0; i<xc; ++i) {
577  normX[i] = normX_mag[i];
578  invnormX[i] = (normX_mag[i] == ZERO) ? ZERO : MONE/normX_mag[i];
579  }
580  }
581  // normalize the vectors
582  MVT::MvScale(X,invnormX);
583  if (this->_hasOp) {
584  MVT::MvScale(*MX,invnormX);
585  }
586  // check that vectors are normalized now
587  if (debug_) {
588  std::vector<MagnitudeType> nrm2(xc);
589  std::cout << dbgstr << "max post-scale norm: (with/without MX) : ";
590  MagnitudeType maxpsnw = ZERO, maxpsnwo = ZERO;
592  for (int i=0; i<xc; i++) {
593  maxpsnw = (nrm2[i] > maxpsnw ? nrm2[i] : maxpsnw);
594  }
595  this->norm(X,nrm2);
596  for (int i=0; i<xc; i++) {
597  maxpsnwo = (nrm2[i] > maxpsnwo ? nrm2[i] : maxpsnwo);
598  }
599  std::cout << "(" << maxpsnw << "," << maxpsnwo << ")" << std::endl;
600  }
601  // project the vectors onto the Qi
602  for (int i=0; i<nq; i++) {
603  MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*Q[i],X,*newC[i],Teuchos::null,MX);
604  }
605  // remove the components in Qi from X
606  for (int i=0; i<nq; i++) {
607  MVT::MvTimesMatAddMv(-ONE,*Q[i],*newC[i],ONE,X);
608  }
609  // un-scale the vectors
610  MVT::MvScale(X,normX);
611  // Recompute the vectors in MX
612  if (this->_hasOp) {
613  OPT::Apply(*(this->_Op),X,*MX);
614  this->_OpCounter += MVT::GetNumberVecs(X);
615  }
616 
617  //
618  // Compute largest column norm of
619  // ( C[0] )
620  // C = ( .... )
621  // ( C[nq-1] )
622  MagnitudeType maxNorm = ZERO;
623  for (int j=0; j<xc; j++) {
624  MagnitudeType sum = ZERO;
625  for (int k=0; k<nq; k++) {
626  for (int i=0; i<qcs[k]; i++) {
627  sum += SCT::magnitude((*newC[k])(i,j))*SCT::magnitude((*newC[k])(i,j));
628  }
629  }
630  maxNorm = (sum > maxNorm) ? sum : maxNorm;
631  }
632 
633  // do we perform another GS?
634  if (maxNorm < 0.36) {
635  doGramSchmidt = false;
636  }
637 
638  // unscale newC to reflect the scaling of X
639  for (int k=0; k<nq; k++) {
640  for (int j=0; j<xc; j++) {
641  for (int i=0; i<qcs[k]; i++) {
642  (*newC[k])(i,j) *= normX[j];
643  }
644  }
645  }
646  // accumulate into C
647  if (normalize_in) {
648  // we are normalizing
649  int info;
650  for (int i=0; i<nq; i++) {
651  info = C[i]->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,*newC[i],*B,ONE);
652  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::SVQBOrthoManager::findBasis(): Input error to SerialDenseMatrix::multiply.");
653  }
654  }
655  else {
656  // not normalizing
657  for (int i=0; i<nq; i++) {
658  (*C[i]) += *newC[i];
659  }
660  }
661  }
662  else { // qsize == 0... don't perform projection
663  // don't do any more outer loops; all we need is to call the normalize code below
664  doGramSchmidt = false;
665  }
666 
667 
669  // perform normalization
670  if (normalize_in) {
671 
672  MagnitudeType condT = tolerance;
673 
674  while (condT >= tolerance) {
675 
676  numSVQB++;
677 
678  // compute X^T Op X
680 
681  // compute scaling matrix for XtMX: D^{.5} and D^{-.5} (D-half and D-half-inv)
682  std::vector<MagnitudeType> Dh(xc), Dhi(xc);
683  for (int i=0; i<xc; i++) {
684  Dh[i] = SCT::magnitude(SCT::squareroot(XtMX(i,i)));
685  Dhi[i] = (Dh[i] == ZERO ? ZERO : MONE/Dh[i]);
686  }
687  // scale XtMX : S = D^{-.5} * XtMX * D^{-.5}
688  for (int i=0; i<xc; i++) {
689  for (int j=0; j<xc; j++) {
690  XtMX(i,j) *= Dhi[i]*Dhi[j];
691  }
692  }
693 
694  // compute the eigenvalue decomposition of S=U*Lambda*U^T (using upper part)
695  int info;
696  lapack.HEEV('V', 'U', xc, XtMX.values(), XtMX.stride(), &lambda[0], &work[0], work.size(), &rwork[0], &info);
697  std::stringstream os;
698  os << "Anasazi::SVQBOrthoManager::findBasis(): Error code " << info << " from HEEV";
699  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OrthoError,
700  os.str() );
701  if (debug_) {
702  std::cout << dbgstr << "eigenvalues of XtMX: (";
703  for (int i=0; i<xc-1; i++) {
704  std::cout << lambda[i] << ",";
705  }
706  std::cout << lambda[xc-1] << ")" << std::endl;
707  }
708 
709  // remember, HEEV orders the eigenvalues from smallest to largest
710  // examine condition number of Lambda, compute Lambda^{-.5}
711  MagnitudeType maxLambda = lambda[xc-1],
712  minLambda = lambda[0];
713  int iZeroMax = -1;
714  for (int i=0; i<xc; i++) {
715  if (lambda[i] <= 10*eps_*maxLambda) { // finish: this was eps_*eps_*maxLambda
716  iZeroMax = i;
717  lambda[i] = ZERO;
718  lambdahi[i] = ZERO;
719  }
720  /*
721  else if (lambda[i] < eps_*maxLambda) {
722  lambda[i] = SCTM::squareroot(eps_*maxLambda);
723  lambdahi[i] = MONE/lambda[i];
724  }
725  */
726  else {
727  lambda[i] = SCTM::squareroot(lambda[i]);
728  lambdahi[i] = MONE/lambda[i];
729  }
730  }
731 
732  // compute X * D^{-.5} * U * Lambda^{-.5} and new Op*X
733  //
734  // copy X into workX
735  std::vector<int> ind(xc);
736  for (int i=0; i<xc; i++) {ind[i] = i;}
737  MVT::SetBlock(X,ind,*workX);
738  //
739  // compute D^{-.5}*U*Lambda^{-.5} into workU
740  workU.assign(XtMX);
741  for (int j=0; j<xc; j++) {
742  for (int i=0; i<xc; i++) {
743  workU(i,j) *= Dhi[i]*lambdahi[j];
744  }
745  }
746  // compute workX * workU into X
747  MVT::MvTimesMatAddMv(ONE,*workX,workU,ZERO,X);
748  //
749  // note, it seems important to apply Op exactly for large condition numbers.
750  // for small condition numbers, we can update MX "implicitly"
751  // this trick reduces the number of applications of Op
752  if (this->_hasOp) {
753  if (maxLambda >= tolerance * minLambda) {
754  // explicit update of MX
755  OPT::Apply(*(this->_Op),X,*MX);
756  this->_OpCounter += MVT::GetNumberVecs(X);
757  }
758  else {
759  // implicit update of MX
760  // copy MX into workX
761  MVT::SetBlock(*MX,ind,*workX);
762  //
763  // compute workX * workU into MX
764  MVT::MvTimesMatAddMv(ONE,*workX,workU,ZERO,*MX);
765  }
766  }
767 
768  // accumulate new B into previous B
769  // B = Lh * U^H * Dh * B
770  for (int j=0; j<xc; j++) {
771  for (int i=0; i<xc; i++) {
772  workU(i,j) = Dh[i] * (*B)(i,j);
773  }
774  }
775  info = B->multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,XtMX,workU,ZERO);
776  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::SVQBOrthoManager::findBasis(): Input error to SerialDenseMatrix::multiply.");
777  for (int j=0; j<xc ;j++) {
778  for (int i=0; i<xc; i++) {
779  (*B)(i,j) *= lambda[i];
780  }
781  }
782 
783  // check iZeroMax (rank indicator)
784  if (iZeroMax >= 0) {
785  if (debug_) {
786  std::cout << dbgstr << "augmenting multivec with " << iZeroMax+1 << " random directions" << std::endl;
787  }
788 
789  numRand++;
790  // put random info in the first iZeroMax+1 vectors of X,MX
791  std::vector<int> ind2(iZeroMax+1);
792  for (int i=0; i<iZeroMax+1; i++) {
793  ind2[i] = i;
794  }
795  Teuchos::RCP<MV> Xnull,MXnull;
796  Xnull = MVT::CloneViewNonConst(X,ind2);
797  MVT::MvRandom(*Xnull);
798  if (this->_hasOp) {
799  MXnull = MVT::CloneViewNonConst(*MX,ind2);
800  OPT::Apply(*(this->_Op),*Xnull,*MXnull);
801  this->_OpCounter += MVT::GetNumberVecs(*Xnull);
802  MXnull = Teuchos::null;
803  }
804  Xnull = Teuchos::null;
805  condT = tolerance;
806  doGramSchmidt = true;
807  break; // break from while(condT > tolerance)
808  }
809 
810  condT = SCTM::magnitude(maxLambda / minLambda);
811  if (debug_) {
812  std::cout << dbgstr << "condT: " << condT << ", tolerance = " << tolerance << std::endl;
813  }
814 
815  // if multiple passes of SVQB are necessary, then pass through outer GS loop again
816  if ((doGramSchmidt == false) && (condT > SCTM::squareroot(tolerance))) {
817  doGramSchmidt = true;
818  }
819 
820  } // end while (condT >= tolerance)
821  }
822  // end if(normalize_in)
823 
824  } // end while (doGramSchmidt)
825 
826  if (debug_) {
827  std::cout << dbgstr << "(numGS,numSVQB,numRand) : "
828  << "(" << numGS
829  << "," << numSVQB
830  << "," << numRand
831  << ")" << std::endl;
832  }
833  return xc;
834  }
835 
836 
837 } // namespace Anasazi
838 
839 #endif // ANASAZI_SVQB_ORTHOMANAGER_HPP
840 
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
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 ...
void norm(const MV &X, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec) const
Implements the interface OrthoManager::norm().
Declaration of basic traits for the multivector type.
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 ...
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 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...
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
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 ...
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
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...
Traits class which defines basic operations on multivectors.
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...
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.
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
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 ...
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...
SVQBOrthoManager(Teuchos::RCP< const OP > Op=Teuchos::null, bool debug=false)
Constructor specifying re-orthogonalization tolerance.
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.