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