Anasazi Version of the Day
Loading...
Searching...
No Matches
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"
64#include "Teuchos_LAPACK.hpp"
65
66namespace 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
85 SVQBOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null, bool debug = false );
86
87
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
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
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);
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
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);
343 return xTx.normFrobenius();
344 }
345
346
347
349 template<class ScalarType, class MV, class OP>
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>
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>
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>
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);
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;
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
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
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.
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.