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;
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,
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,
174  Teuchos::RCP<MV> MX = Teuchos::null
175  ) const;
176 
177 
238  MV &X,
241  = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(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 
257  orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const;
258 
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,
283  bool normalize_in ) const;
284  };
285 
286 
288  // Constructor
289  template<class ScalarType, class MV, class OP>
291  : MatOrthoManager<ScalarType,MV,OP>(Op), dbgstr(" *** "), debug_(debug) {
292 
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>
306  const ScalarType ONE = SCT::one();
307  int rank = MVT::GetNumberVecs(X);
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>
321  const MV &X,
322  const MV &Y,
325  ) const {
326  int r1 = MVT::GetNumberVecs(X);
327  int r2 = MVT::GetNumberVecs(Y);
330  return xTx.normFrobenius();
331  }
332 
333 
334 
336  template<class ScalarType, class MV, class OP>
338  MV &X,
341  Teuchos::RCP<MV> MX,
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,
355  Teuchos::RCP<MV> MX) const {
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,
371  Teuchos::RCP<MV> MX,
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,
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);
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 ) {
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 ) {
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);
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 
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Declaration of basic traits for the multivector type.
void projectMat(MV &X, Teuchos::Array< Teuchos::RCP< const MV > > Q, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C=Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)), Teuchos::RCP< MV > MX=Teuchos::null, Teuchos::Array< Teuchos::RCP< const MV > > MQ=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null))) const
Given a list of mutually orthogonal and internally orthonormal bases Q, this method projects a multiv...
OrdinalType ILAENV(const OrdinalType ispec, const std::string &NAME, const std::string &OPTS, const OrdinalType N1=-1, const OrdinalType N2=-1, const OrdinalType N3=-1, const OrdinalType N4=-1) const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Virtual base class which defines basic traits for the operator type.
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...
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
void HEEV(const char JOBZ, const char UPLO, const OrdinalType n, ScalarType *A, const OrdinalType lda, MagnitudeType *W, ScalarType *WORK, const OrdinalType lwork, MagnitudeType *RWORK, OrdinalType *info) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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 ...
ScalarType LAMCH(const char CMACH) const
void normMat(const MV &X, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, Teuchos::RCP< const MV > MX=Teuchos::null) const
Provides the norm induced by the matrix-based inner product.
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
void innerProdMat(const MV &X, const MV &Y, Teuchos::SerialDenseMatrix< int, ScalarType > &Z, Teuchos::RCP< const MV > MX=Teuchos::null, Teuchos::RCP< const MV > MY=Teuchos::null) const
Provides a matrix-based inner product.
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 ...
SVQBOrthoManager(Teuchos::RCP< const OP > Op=Teuchos::null, bool debug=false)
Constructor specifying re-orthogonalization tolerance.
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 ...
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...