Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziICGSOrthoManager.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
47#ifndef ANASAZI_ICSG_ORTHOMANAGER_HPP
48#define ANASAZI_ICSG_ORTHOMANAGER_HPP
49
57#include "AnasaziConfigDefs.hpp"
61#include "Teuchos_TimeMonitor.hpp"
62#include "Teuchos_LAPACK.hpp"
63#include "Teuchos_BLAS.hpp"
64#ifdef ANASAZI_ICGS_DEBUG
65# include <Teuchos_FancyOStream.hpp>
66#endif
67
68namespace Anasazi {
69
70 template<class ScalarType, class MV, class OP>
71 class ICGSOrthoManager : public GenOrthoManager<ScalarType,MV,OP> {
72
73 private:
74 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
75 typedef Teuchos::ScalarTraits<ScalarType> SCT;
78
79 public:
80
82
83
84 ICGSOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null, int numIters = 2,
85 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps = 0.0,
86 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol = 0.20 );
87
88
92
93
95
96
178 void projectGen(
179 MV &S,
180 Teuchos::Array<Teuchos::RCP<const MV> > X,
181 Teuchos::Array<Teuchos::RCP<const MV> > Y,
182 bool isBiOrtho,
183 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
184 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
185 Teuchos::RCP<MV> MS = Teuchos::null,
186 Teuchos::Array<Teuchos::RCP<const MV> > MX = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)),
187 Teuchos::Array<Teuchos::RCP<const MV> > MY = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
188 ) const;
189
190
283 MV &S,
284 Teuchos::Array<Teuchos::RCP<const MV> > X,
285 Teuchos::Array<Teuchos::RCP<const MV> > Y,
286 bool isBiOrtho,
287 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
288 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
289 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
290 Teuchos::RCP<MV> MS = Teuchos::null,
291 Teuchos::Array<Teuchos::RCP<const MV> > MX = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)),
292 Teuchos::Array<Teuchos::RCP<const MV> > MY = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
293 ) const;
294
295
297
298
300
301
302
314 void projectMat (
315 MV &X,
316 Teuchos::Array<Teuchos::RCP<const MV> > Q,
317 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
318 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
319 Teuchos::RCP<MV> MX = Teuchos::null,
320 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
321 ) const;
322
323
332 int normalizeMat (
333 MV &X,
334 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
335 Teuchos::RCP<MV> MX = Teuchos::null
336 ) const;
337
338
348 MV &X,
349 Teuchos::Array<Teuchos::RCP<const MV> > Q,
350 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
351 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
352 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
353 Teuchos::RCP<MV> MX = Teuchos::null,
354 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
355 ) const;
356
358
360
361
366 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
367 orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const;
368
373 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
374 orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const;
375
377
379
380
382 void setNumIters( int numIters ) {
383 numIters_ = numIters;
384 TEUCHOS_TEST_FOR_EXCEPTION(numIters_ < 1,std::invalid_argument,
385 "Anasazi::ICGSOrthoManager::setNumIters(): input must be >= 1.");
386 }
387
389 int getNumIters() const { return numIters_; }
390
392
393 private:
394 MagnitudeType eps_;
395 MagnitudeType tol_;
396
398 int numIters_;
399
400 // ! Routine to find an orthonormal basis
401 int findBasis(MV &X, Teuchos::RCP<MV> MX,
402 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
403 bool completeBasis, int howMany = -1) const;
404 };
405
406
407
409 // Constructor
410 template<class ScalarType, class MV, class OP>
412 int numIters,
413 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps,
414 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol) :
415 GenOrthoManager<ScalarType,MV,OP>(Op), eps_(eps), tol_(tol)
416 {
417 setNumIters(numIters);
418 TEUCHOS_TEST_FOR_EXCEPTION(eps_ < SCT::magnitude(SCT::zero()),std::invalid_argument,
419 "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"eps\" must be non-negative.");
420 if (eps_ == 0) {
421 Teuchos::LAPACK<int,MagnitudeType> lapack;
422 eps_ = lapack.LAMCH('E');
423 eps_ = Teuchos::ScalarTraits<MagnitudeType>::pow(eps_,.50);
424 }
425 TEUCHOS_TEST_FOR_EXCEPTION(
426 tol_ < SCT::magnitude(SCT::zero()) || tol_ > SCT::magnitude(SCT::one()),
427 std::invalid_argument,
428 "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"tol\" must be in [0,1].");
429 }
430
431
432
434 // Compute the distance from orthonormality
435 template<class ScalarType, class MV, class OP>
436 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
437 ICGSOrthoManager<ScalarType,MV,OP>::orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX) const {
438 const ScalarType ONE = SCT::one();
439 int rank = MVT::GetNumberVecs(X);
440 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
442 for (int i=0; i<rank; i++) {
443 xTx(i,i) -= ONE;
444 }
445 return xTx.normFrobenius();
446 }
447
448
449
451 // Compute the distance from orthogonality
452 template<class ScalarType, class MV, class OP>
453 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
454 ICGSOrthoManager<ScalarType,MV,OP>::orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const {
455 int r1 = MVT::GetNumberVecs(X1);
456 int r2 = MVT::GetNumberVecs(X2);
457 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
459 return xTx.normFrobenius();
460 }
461
462
463
465 template<class ScalarType, class MV, class OP>
467 MV &X,
468 Teuchos::Array<Teuchos::RCP<const MV> > Q,
469 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
470 Teuchos::RCP<MV> MX,
471 Teuchos::Array<Teuchos::RCP<const MV> > MQ
472 ) const
473 {
474 projectGen(X,Q,Q,true,C,MX,MQ,MQ);
475 }
476
477
478
480 // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
481 template<class ScalarType, class MV, class OP>
483 MV &X,
484 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
485 Teuchos::RCP<MV> MX) const
486 {
487 // call findBasis(), with the instruction to try to generate a basis of rank numvecs(X)
488 // findBasis() requires MX
489
490 int xc = MVT::GetNumberVecs(X);
491 ptrdiff_t xr = MVT::GetGlobalLength(X);
492
493 // if Op==null, MX == X (via pointer)
494 // Otherwise, either the user passed in MX or we will allocated and compute it
495 if (this->_hasOp) {
496 if (MX == Teuchos::null) {
497 // we need to allocate space for MX
498 MX = MVT::Clone(X,xc);
499 OPT::Apply(*(this->_Op),X,*MX);
500 this->_OpCounter += MVT::GetNumberVecs(X);
501 }
502 }
503
504 // if the user doesn't want to store the coefficients,
505 // allocate some local memory for them
506 if ( B == Teuchos::null ) {
507 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
508 }
509
510 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
511 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
512
513 // check size of C, B
514 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
515 "Anasazi::ICGSOrthoManager::normalizeMat(): X must be non-empty" );
516 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
517 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not consistent with size of B" );
518 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
519 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not consistent with size of MX" );
520 TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument,
521 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not feasible for normalization" );
522
523 return findBasis(X, MX, *B, true );
524 }
525
526
527
529 // Find an Op-orthonormal basis for span(X) - span(W)
530 template<class ScalarType, class MV, class OP>
532 MV &X,
533 Teuchos::Array<Teuchos::RCP<const MV> > Q,
534 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
535 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
536 Teuchos::RCP<MV> MX,
537 Teuchos::Array<Teuchos::RCP<const MV> > MQ
538 ) const
539 {
540 return projectAndNormalizeGen(X,Q,Q,true,C,B,MX,MQ,MQ);
541 }
542
543
544
546 template<class ScalarType, class MV, class OP>
548 MV &S,
549 Teuchos::Array<Teuchos::RCP<const MV> > X,
550 Teuchos::Array<Teuchos::RCP<const MV> > Y,
551 bool isBiortho,
552 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
553 Teuchos::RCP<MV> MS,
554 Teuchos::Array<Teuchos::RCP<const MV> > MX,
555 Teuchos::Array<Teuchos::RCP<const MV> > MY
556 ) const
557 {
558 // For the inner product defined by the operator Op or the identity (Op == 0)
559 // -> Orthogonalize S against each Y[i], modifying it in the range of X[i]
560 // Modify MS accordingly
561 //
562 // Note that when Op is 0, MS is not referenced
563 //
564 // Parameter variables
565 //
566 // S : Multivector to be transformed
567 //
568 // MS : Image of the block vector S by the mass matrix
569 //
570 // X,Y: Bases to orthogonalize against/via.
571 //
572
573#ifdef ANASAZI_ICGS_DEBUG
574 // Get a FancyOStream from out_arg or create a new one ...
575 Teuchos::RCP<Teuchos::FancyOStream>
576 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
577 out->setShowAllFrontMatter(false).setShowProcRank(true);
578 *out << "Entering Anasazi::ICGSOrthoManager::projectGen(...)\n";
579#endif
580
581 const ScalarType ONE = SCT::one();
582 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
583 Teuchos::LAPACK<int,ScalarType> lapack;
584 Teuchos::BLAS<int,ScalarType> blas;
585
586 int sc = MVT::GetNumberVecs( S );
587 ptrdiff_t sr = MVT::GetGlobalLength( S );
588 int numxy = X.length();
589 TEUCHOS_TEST_FOR_EXCEPTION(X.length() != Y.length(),std::invalid_argument,
590 "Anasazi::ICGSOrthoManager::projectGen(): X and Y must contain the same number of multivectors.");
591 std::vector<int> xyc(numxy);
592 // short-circuit
593 if (numxy == 0 || sc == 0 || sr == 0) {
594#ifdef ANASAZI_ICGS_DEBUG
595 *out << "Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
596#endif
597 return;
598 }
599 // if we don't have enough C, expand it with null references
600 // if we have too many, resize to throw away the latter ones
601 // if we have exactly as many as we have X,Y this call has no effect
602 //
603 // same for MX, MY
604 C.resize(numxy);
605 MX.resize(numxy);
606 MY.resize(numxy);
607
608 // check size of S w.r.t. common sense
609 TEUCHOS_TEST_FOR_EXCEPTION( sc<0 || sr<0, std::invalid_argument,
610 "Anasazi::ICGSOrthoManager::projectGen(): MVT returned negative dimensions for S." );
611
612 // check size of MS
613 if (this->_hasOp == true) {
614 if (MS != Teuchos::null) {
615 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MS) != sr, std::invalid_argument,
616 "Anasazi::ICGSOrthoManager::projectGen(): MS length not consistent with S." );
617 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*MS) != sc, std::invalid_argument,
618 "Anasazi::ICGSOrthoManager::projectGen(): MS width not consistent with S." );
619 }
620 }
621
622 // tally up size of all X,Y and check/allocate C
623 ptrdiff_t sumxyc = 0;
624 int MYmissing = 0;
625 int MXmissing = 0;
626 for (int i=0; i<numxy; i++) {
627 if (X[i] != Teuchos::null && Y[i] != Teuchos::null) {
628 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*X[i]) != sr, std::invalid_argument,
629 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i << "] length not consistent with S." );
630 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*Y[i]) != sr, std::invalid_argument,
631 "Anasazi::ICGSOrthoManager::projectGen(): Y[" << i << "] length not consistent with S." );
632 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*X[i]) != MVT::GetNumberVecs(*Y[i]), std::invalid_argument,
633 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i << "] and Y[" << i << "] widths not consistent." );
634
635 xyc[i] = MVT::GetNumberVecs( *X[i] );
636 TEUCHOS_TEST_FOR_EXCEPTION( sr < static_cast<ptrdiff_t>(xyc[i]), std::invalid_argument,
637 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i << "],Y[" << i << "] have less rows than columns, and therefore cannot be full rank." );
638 sumxyc += xyc[i];
639
640 // check size of C[i]
641 if ( C[i] == Teuchos::null ) {
642 C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],sc) );
643 }
644 else {
645 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != xyc[i] || C[i]->numCols() != sc , std::invalid_argument,
646 "Anasazi::ICGSOrthoManager::projectGen(): Size of Q not consistent with size of C." );
647 }
648 // check sizes of MX[i], MY[i] if present
649 // if not present, count their absence
650 if (MX[i] != Teuchos::null) {
651 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MX[i]) != sr || MVT::GetNumberVecs(*MX[i]) != xyc[i], std::invalid_argument,
652 "Anasazi::ICGSOrthoManager::projectGen(): Size of MX[" << i << "] not consistent with size of X[" << i << "]." );
653 }
654 else {
655 MXmissing += xyc[i];
656 }
657 if (MY[i] != Teuchos::null) {
658 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MY[i]) != sr || MVT::GetNumberVecs(*MY[i]) != xyc[i], std::invalid_argument,
659 "Anasazi::ICGSOrthoManager::projectGen(): Size of MY[" << i << "] not consistent with size of Y[" << i << "]." );
660 }
661 else {
662 MYmissing += xyc[i];
663 }
664 }
665 else {
666 // if one is null and the other is not... the user may have made a mistake
667 TEUCHOS_TEST_FOR_EXCEPTION(X[i] != Teuchos::null || Y[i] != Teuchos::null, std::invalid_argument,
668 "Anasazi::ICGSOrthoManager::projectGen(): "
669 << (X[i] == Teuchos::null ? "Y[" : "X[") << i << "] was provided but "
670 << (X[i] == Teuchos::null ? "X[" : "Y[") << i << "] was not.");
671 }
672 }
673
674 // is this operation even feasible?
675 TEUCHOS_TEST_FOR_EXCEPTION(sumxyc > sr, std::invalid_argument,
676 "Anasazi::ICGSOrthoManager::projectGen(): dimension of all X[i],Y[i] is "
677 << sumxyc << ", but length of vectors is only " << sr << ". This is infeasible.");
678
679
680 /****** DO NO MODIFY *MS IF _hasOp == false
681 * if _hasOp == false, we don't need MS, MX or MY
682 *
683 * if _hasOp == true, we need MS (for S M-norms) and
684 * therefore, we must also update MS, regardless of whether
685 * it gets returned to the user (i.e., whether the user provided it)
686 * we may need to allocate and compute MX or MY
687 *
688 * let BXY denote:
689 * "X" - we have all M*X[i]
690 * "Y" - we have all M*Y[i]
691 * "B" - we have biorthogonality (for all X[i],Y[i])
692 * Encode these as values
693 * X = 1
694 * Y = 2
695 * B = 4
696 * We have 8 possibilities, 0-7
697 *
698 * We must allocate storage and computer the following, lest
699 * innerProdMat do it for us:
700 * none (0) - allocate MX, for inv(<X,Y>) and M*S
701 *
702 * for the following, we will compute M*S manually before returning
703 * B(4) BY(6) Y(2) --> updateMS = 1
704 * for the following, we will update M*S as we go, using M*X
705 * XY(3) X(1) none(0) BXY(7) BX(5) --> updateMS = 2
706 * these choices favor applications of M over allocation of memory
707 *
708 */
709 int updateMS = -1;
710 if (this->_hasOp) {
711 int whichAlloc = 0;
712 if (MXmissing == 0) {
713 whichAlloc += 1;
714 }
715 if (MYmissing == 0) {
716 whichAlloc += 2;
717 }
718 if (isBiortho) {
719 whichAlloc += 4;
720 }
721
722 switch (whichAlloc) {
723 case 2:
724 case 4:
725 case 6:
726 updateMS = 1;
727 break;
728 case 0:
729 case 1:
730 case 3:
731 case 5:
732 case 7:
733 updateMS = 2;
734 break;
735 }
736
737 // produce MS
738 if (MS == Teuchos::null) {
739#ifdef ANASAZI_ICGS_DEBUG
740 *out << "Allocating MS...\n";
741#endif
742 MS = MVT::Clone(S,MVT::GetNumberVecs(S));
743 OPT::Apply(*(this->_Op),S,*MS);
744 this->_OpCounter += MVT::GetNumberVecs(S);
745 }
746
747 // allocate the rest
748 if (whichAlloc == 0) {
749 // allocate and compute missing MX
750 for (int i=0; i<numxy; i++) {
751 if (MX[i] == Teuchos::null) {
752#ifdef ANASAZI_ICGS_DEBUG
753 *out << "Allocating MX[" << i << "]...\n";
754#endif
755 Teuchos::RCP<MV> tmpMX = MVT::Clone(*X[i],xyc[i]);
756 OPT::Apply(*(this->_Op),*X[i],*tmpMX);
757 MX[i] = tmpMX;
758 this->_OpCounter += xyc[i];
759 }
760 }
761 }
762 }
763 else {
764 // Op == I --> MS == S
765 MS = Teuchos::rcpFromRef(S);
766 updateMS = 0;
767 }
768 TEUCHOS_TEST_FOR_EXCEPTION(updateMS == -1,std::logic_error,
769 "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic.");
770
771
773 // Perform the Gram-Schmidt transformation for a block of vectors
775
776 // Compute Cholesky factorizations for the Y'*M*X
777 // YMX stores the YMX (initially) and their Cholesky factorizations (utlimately)
778 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > YMX(numxy);
779 if (isBiortho == false) {
780 for (int i=0; i<numxy; i++) {
781#ifdef ANASAZI_ICGS_DEBUG
782 *out << "Computing YMX[" << i << "] and its Cholesky factorization...\n";
783#endif
784 YMX[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],xyc[i]) );
785 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*Y[i],*X[i],*YMX[i],MY[i],MX[i]);
786#ifdef ANASAZI_ICGS_DEBUG
787 // YMX should be symmetric positive definite
788 // the cholesky will check the positive definiteness, but it looks only at the upper half
789 // we will check the symmetry by removing the upper half from the lower half, which should
790 // result in zeros
791 // also, diagonal of YMX should be real; consider the complex part as error
792 {
793 MagnitudeType err = ZERO;
794 for (int jj=0; jj<YMX[i]->numCols(); jj++) {
795 err =+ SCT::magnitude(SCT::imag((*YMX[i])(jj,jj)));
796 for (int ii=jj; ii<YMX[i]->numRows(); ii++) {
797 err += SCT::magnitude( (*YMX[i])(ii,jj) - SCT::conjugate((*YMX[i])(jj,ii)) );
798 }
799 }
800 *out << "Symmetry error in YMX[" << i << "] == " << err << "\n";
801 }
802#endif
803 // take the LU factorization
804 int info;
805 lapack.POTRF('U',YMX[i]->numRows(),YMX[i]->values(),YMX[i]->stride(),&info);
806 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
807 "Anasazi::ICGSOrthoManager::projectGen(): Error computing Cholesky factorization of Y[i]^T * M * X[i] using POTRF: returned info " << info);
808 }
809 }
810
811 // Compute the initial Op-norms
812#ifdef ANASAZI_ICGS_DEBUG
813 std::vector<MagnitudeType> oldNorms(sc);
815 *out << "oldNorms = { ";
816 std::copy(oldNorms.begin(), oldNorms.end(), std::ostream_iterator<MagnitudeType>(*out, " "));
817 *out << "}\n";
818#endif
819
820
821 // clear the C[i] and allocate Ccur
822 Teuchos::Array<Teuchos::SerialDenseMatrix<int,ScalarType> > Ccur(numxy);
823 for (int i=0; i<numxy; i++) {
824 C[i]->putScalar(ZERO);
825 Ccur[i].reshape(C[i]->numRows(),C[i]->numCols());
826 }
827
828 for (int iter=0; iter < numIters_; iter++) {
829#ifdef ANASAZI_ICGS_DEBUG
830 *out << "beginning iteration " << iter+1 << "\n";
831#endif
832
833 // Define the product Y[i]'*Op*S
834 for (int i=0; i<numxy; i++) {
835 // Compute Y[i]'*M*S
836 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*Y[i],S,Ccur[i],MY[i],MS);
837 if (isBiortho == false) {
838 // C[i] = inv(YMX[i])*(Y[i]'*M*S)
839 int info;
840 lapack.POTRS('U',YMX[i]->numCols(),Ccur[i].numCols(),
841 YMX[i]->values(),YMX[i]->stride(),
842 Ccur[i].values(),Ccur[i].stride(), &info);
843 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
844 "Anasazi::ICGSOrthoManager::projectGen(): Error code " << info << " from lapack::POTRS." );
845 }
846
847 // Multiply by X[i] and subtract the result in S
848#ifdef ANASAZI_ICGS_DEBUG
849 *out << "Applying projector P_{X[" << i << "],Y[" << i << "]}...\n";
850#endif
851 MVT::MvTimesMatAddMv( -ONE, *X[i], Ccur[i], ONE, S );
852
853 // Accumulate coeffs into previous step
854 *C[i] += Ccur[i];
855
856 // Update MS as required
857 if (updateMS == 1) {
858#ifdef ANASAZI_ICGS_DEBUG
859 *out << "Updating MS...\n";
860#endif
861 OPT::Apply( *(this->_Op), S, *MS);
862 this->_OpCounter += sc;
863 }
864 else if (updateMS == 2) {
865#ifdef ANASAZI_ICGS_DEBUG
866 *out << "Updating MS...\n";
867#endif
868 MVT::MvTimesMatAddMv( -ONE, *MX[i], Ccur[i], ONE, *MS );
869 }
870 }
871
872 // Compute new Op-norms
873#ifdef ANASAZI_ICGS_DEBUG
874 std::vector<MagnitudeType> newNorms(sc);
876 *out << "newNorms = { ";
877 std::copy(newNorms.begin(), newNorms.end(), std::ostream_iterator<MagnitudeType>(*out, " "));
878 *out << "}\n";
879#endif
880 }
881
882#ifdef ANASAZI_ICGS_DEBUG
883 *out << "Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
884#endif
885 }
886
887
888
890 // Find an Op-orthonormal basis for span(S) - span(Y)
891 template<class ScalarType, class MV, class OP>
893 MV &S,
894 Teuchos::Array<Teuchos::RCP<const MV> > X,
895 Teuchos::Array<Teuchos::RCP<const MV> > Y,
896 bool isBiortho,
897 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
898 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
899 Teuchos::RCP<MV> MS,
900 Teuchos::Array<Teuchos::RCP<const MV> > MX,
901 Teuchos::Array<Teuchos::RCP<const MV> > MY
902 ) const {
903 // For the inner product defined by the operator Op or the identity (Op == 0)
904 // -> Orthogonalize S against each Y[i], modifying it in the range of X[i]
905 // Modify MS accordingly
906 // Then construct a M-orthonormal basis for the remaining part
907 //
908 // Note that when Op is 0, MS is not referenced
909 //
910 // Parameter variables
911 //
912 // S : Multivector to be transformed
913 //
914 // MS : Image of the block vector S by the mass matrix
915 //
916 // X,Y: Bases to orthogonalize against/via.
917 //
918
919#ifdef ANASAZI_ICGS_DEBUG
920 // Get a FancyOStream from out_arg or create a new one ...
921 Teuchos::RCP<Teuchos::FancyOStream>
922 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
923 out->setShowAllFrontMatter(false).setShowProcRank(true);
924 *out << "Entering Anasazi::ICGSOrthoManager::projectAndNormalizeGen(...)\n";
925#endif
926
927 int rank;
928 int sc = MVT::GetNumberVecs( S );
929 ptrdiff_t sr = MVT::GetGlobalLength( S );
930 int numxy = X.length();
931 TEUCHOS_TEST_FOR_EXCEPTION(X.length() != Y.length(),std::invalid_argument,
932 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X and Y must contain the same number of multivectors.");
933 std::vector<int> xyc(numxy);
934 // short-circuit
935 if (sc == 0 || sr == 0) {
936#ifdef ANASAZI_ICGS_DEBUG
937 *out << "Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
938#endif
939 return 0;
940 }
941 // if we don't have enough C, expand it with null references
942 // if we have too many, resize to throw away the latter ones
943 // if we have exactly as many as we have X,Y this call has no effect
944 //
945 // same for MX, MY
946 C.resize(numxy);
947 MX.resize(numxy);
948 MY.resize(numxy);
949
950 // check size of S w.r.t. common sense
951 TEUCHOS_TEST_FOR_EXCEPTION( sc<0 || sr<0, std::invalid_argument,
952 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MVT returned negative dimensions for S." );
953
954 // check size of MS
955 if (this->_hasOp == true) {
956 if (MS != Teuchos::null) {
957 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MS) != sr, std::invalid_argument,
958 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS length not consistent with S." );
959 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*MS) != sc, std::invalid_argument,
960 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS width not consistent with S." );
961 }
962 }
963
964 // tally up size of all X,Y and check/allocate C
965 ptrdiff_t sumxyc = 0;
966 int MYmissing = 0;
967 int MXmissing = 0;
968 for (int i=0; i<numxy; i++) {
969 if (X[i] != Teuchos::null && Y[i] != Teuchos::null) {
970 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*X[i]) != sr, std::invalid_argument,
971 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i << "] length not consistent with S." );
972 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*Y[i]) != sr, std::invalid_argument,
973 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Y[" << i << "] length not consistent with S." );
974 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*X[i]) != MVT::GetNumberVecs(*Y[i]), std::invalid_argument,
975 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i << "] and Y[" << i << "] widths not consistent." );
976
977 xyc[i] = MVT::GetNumberVecs( *X[i] );
978 TEUCHOS_TEST_FOR_EXCEPTION( sr < static_cast<ptrdiff_t>(xyc[i]), std::invalid_argument,
979 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i << "],Y[" << i << "] have less rows than columns, and therefore cannot be full rank." );
980 sumxyc += xyc[i];
981
982 // check size of C[i]
983 if ( C[i] == Teuchos::null ) {
984 C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],sc) );
985 }
986 else {
987 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != xyc[i] || C[i]->numCols() != sc , std::invalid_argument,
988 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of Q not consistent with size of C." );
989 }
990 // check sizes of MX[i], MY[i] if present
991 // if not present, count their absence
992 if (MX[i] != Teuchos::null) {
993 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MX[i]) != sr || MVT::GetNumberVecs(*MX[i]) != xyc[i], std::invalid_argument,
994 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of MX[" << i << "] not consistent with size of X[" << i << "]." );
995 }
996 else {
997 MXmissing += xyc[i];
998 }
999 if (MY[i] != Teuchos::null) {
1000 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MY[i]) != sr || MVT::GetNumberVecs(*MY[i]) != xyc[i], std::invalid_argument,
1001 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of MY[" << i << "] not consistent with size of Y[" << i << "]." );
1002 }
1003 else {
1004 MYmissing += xyc[i];
1005 }
1006 }
1007 else {
1008 // if one is null and the other is not... the user may have made a mistake
1009 TEUCHOS_TEST_FOR_EXCEPTION(X[i] != Teuchos::null || Y[i] != Teuchos::null, std::invalid_argument,
1010 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): "
1011 << (X[i] == Teuchos::null ? "Y[" : "X[") << i << "] was provided but "
1012 << (X[i] == Teuchos::null ? "X[" : "Y[") << i << "] was not.");
1013 }
1014 }
1015
1016 // is this operation even feasible?
1017 TEUCHOS_TEST_FOR_EXCEPTION(sumxyc + sc > sr, std::invalid_argument,
1018 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): dimension of all X[i],Y[i] is "
1019 << sumxyc << " and requested " << sc << "-dimensional basis, but length of vectors is only "
1020 << sr << ". This is infeasible.");
1021
1022
1023 /****** DO NO MODIFY *MS IF _hasOp == false
1024 * if _hasOp == false, we don't need MS, MX or MY
1025 *
1026 * if _hasOp == true, we need MS (for S M-norms and normalization) and
1027 * therefore, we must also update MS, regardless of whether
1028 * it gets returned to the user (i.e., whether the user provided it)
1029 * we may need to allocate and compute MX or MY
1030 *
1031 * let BXY denote:
1032 * "X" - we have all M*X[i]
1033 * "Y" - we have all M*Y[i]
1034 * "B" - we have biorthogonality (for all X[i],Y[i])
1035 * Encode these as values
1036 * X = 1
1037 * Y = 2
1038 * B = 4
1039 * We have 8 possibilities, 0-7
1040 *
1041 * We must allocate storage and computer the following, lest
1042 * innerProdMat do it for us:
1043 * none (0) - allocate MX, for inv(<X,Y>) and M*S
1044 *
1045 * for the following, we will compute M*S manually before returning
1046 * B(4) BY(6) Y(2) --> updateMS = 1
1047 * for the following, we will update M*S as we go, using M*X
1048 * XY(3) X(1) none(0) BXY(7) BX(5) --> updateMS = 2
1049 * these choices favor applications of M over allocation of memory
1050 *
1051 */
1052 int updateMS = -1;
1053 if (this->_hasOp) {
1054 int whichAlloc = 0;
1055 if (MXmissing == 0) {
1056 whichAlloc += 1;
1057 }
1058 if (MYmissing == 0) {
1059 whichAlloc += 2;
1060 }
1061 if (isBiortho) {
1062 whichAlloc += 4;
1063 }
1064
1065 switch (whichAlloc) {
1066 case 2:
1067 case 4:
1068 case 6:
1069 updateMS = 1;
1070 break;
1071 case 0:
1072 case 1:
1073 case 3:
1074 case 5:
1075 case 7:
1076 updateMS = 2;
1077 break;
1078 }
1079
1080 // produce MS
1081 if (MS == Teuchos::null) {
1082#ifdef ANASAZI_ICGS_DEBUG
1083 *out << "Allocating MS...\n";
1084#endif
1085 MS = MVT::Clone(S,MVT::GetNumberVecs(S));
1086 OPT::Apply(*(this->_Op),S,*MS);
1087 this->_OpCounter += MVT::GetNumberVecs(S);
1088 }
1089
1090 // allocate the rest
1091 if (whichAlloc == 0) {
1092 // allocate and compute missing MX
1093 for (int i=0; i<numxy; i++) {
1094 if (MX[i] == Teuchos::null) {
1095#ifdef ANASAZI_ICGS_DEBUG
1096 *out << "Allocating MX[" << i << "]...\n";
1097#endif
1098 Teuchos::RCP<MV> tmpMX = MVT::Clone(*X[i],xyc[i]);
1099 OPT::Apply(*(this->_Op),*X[i],*tmpMX);
1100 MX[i] = tmpMX;
1101 this->_OpCounter += xyc[i];
1102 }
1103 }
1104 }
1105 }
1106 else {
1107 // Op == I --> MS == S
1108 MS = Teuchos::rcpFromRef(S);
1109 updateMS = 0;
1110 }
1111 TEUCHOS_TEST_FOR_EXCEPTION(updateMS == -1,std::logic_error,
1112 "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic.");
1113
1114
1115 // if the user doesn't want to store the coefficients,
1116 // allocate some local memory for them
1117 if ( B == Teuchos::null ) {
1118 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(sc,sc) );
1119 }
1120 else {
1121 // check size of B
1122 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != sc || B->numCols() != sc, std::invalid_argument,
1123 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of S must be consistent with size of B" );
1124 }
1125
1126
1127 // orthogonalize all of S against X,Y
1128 projectGen(S,X,Y,isBiortho,C,MS,MX,MY);
1129
1130 Teuchos::SerialDenseMatrix<int,ScalarType> oldCoeff(sc,1);
1131 // start working
1132 rank = 0;
1133 int numTries = 10; // each vector in X gets 10 random chances to escape degeneracy
1134 int oldrank = -1;
1135 do {
1136 int curssize = sc - rank;
1137
1138 // orthonormalize S, but quit if it is rank deficient
1139 // we can't let findBasis generated random vectors to complete the basis,
1140 // because it doesn't know about Q; we will do this ourselves below
1141#ifdef ANASAZI_ICGS_DEBUG
1142 *out << "Attempting to find orthonormal basis for X...\n";
1143#endif
1144 rank = findBasis(S,MS,*B,false,curssize);
1145
1146 if (oldrank != -1 && rank != oldrank) {
1147 // we had previously stopped before, after operating on vector oldrank
1148 // we saved its coefficients, augmented it with a random vector, and
1149 // then called findBasis() again, which proceeded to add vector oldrank
1150 // to the basis.
1151 // now, restore the saved coefficients into B
1152 for (int i=0; i<sc; i++) {
1153 (*B)(i,oldrank) = oldCoeff(i,0);
1154 }
1155 }
1156
1157 if (rank < sc) {
1158 if (rank != oldrank) {
1159 // we quit on this vector and will augment it with random below
1160 // this is the first time that we have quit on this vector
1161 // therefor, (*B)(:,rank) contains the actual coefficients of the
1162 // input vectors with respect to the previous vectors in the basis
1163 // save these values, as (*B)(:,rank) will be overwritten by our next
1164 // call to findBasis()
1165 // we will restore it after we are done working on this vector
1166 for (int i=0; i<sc; i++) {
1167 oldCoeff(i,0) = (*B)(i,rank);
1168 }
1169 }
1170 }
1171
1172 if (rank == sc) {
1173 // we are done
1174#ifdef ANASAZI_ICGS_DEBUG
1175 *out << "Finished computing basis.\n";
1176#endif
1177 break;
1178 }
1179 else {
1180 TEUCHOS_TEST_FOR_EXCEPTION( rank < oldrank, OrthoError,
1181 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): basis lost rank; this shouldn't happen");
1182
1183 if (rank != oldrank) {
1184 // we added a basis vector from random info; reset the chance counter
1185 numTries = 10;
1186 // store old rank
1187 oldrank = rank;
1188 }
1189 else {
1190 // has this vector run out of chances to escape degeneracy?
1191 if (numTries <= 0) {
1192 break;
1193 }
1194 }
1195 // use one of this vector's chances
1196 numTries--;
1197
1198 // randomize troubled direction
1199#ifdef ANASAZI_ICGS_DEBUG
1200 *out << "Inserting random vector in X[" << rank << "]. Attempt " << 10-numTries << ".\n";
1201#endif
1202 Teuchos::RCP<MV> curS, curMS;
1203 std::vector<int> ind(1);
1204 ind[0] = rank;
1205 curS = MVT::CloneViewNonConst(S,ind);
1206 MVT::MvRandom(*curS);
1207 if (this->_hasOp) {
1208#ifdef ANASAZI_ICGS_DEBUG
1209 *out << "Applying operator to random vector.\n";
1210#endif
1211 curMS = MVT::CloneViewNonConst(*MS,ind);
1212 OPT::Apply( *(this->_Op), *curS, *curMS );
1213 this->_OpCounter += MVT::GetNumberVecs(*curS);
1214 }
1215
1216 // orthogonalize against X,Y
1217 // if !this->_hasOp, the curMS will be ignored.
1218 // we don't care about these coefficients
1219 // on the contrary, we need to preserve the previous coeffs
1220 {
1221 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC(0);
1222 projectGen(*curS,X,Y,isBiortho,dummyC,curMS,MX,MY);
1223 }
1224 }
1225 } while (1);
1226
1227 // this should never raise an exception; but our post-conditions oblige us to check
1228 TEUCHOS_TEST_FOR_EXCEPTION( rank > sc || rank < 0, std::logic_error,
1229 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Debug error in rank variable." );
1230
1231#ifdef ANASAZI_ICGS_DEBUG
1232 *out << "Returning " << rank << " from Anasazi::ICGSOrthoManager::projectAndNormalizeGen(...)\n";
1233#endif
1234
1235 return rank;
1236 }
1237
1238
1239
1241 // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that
1242 // the rank is numvectors(X)
1243 template<class ScalarType, class MV, class OP>
1245 MV &X, Teuchos::RCP<MV> MX,
1246 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
1247 bool completeBasis, int howMany ) const {
1248
1249 // For the inner product defined by the operator Op or the identity (Op == 0)
1250 // -> Orthonormalize X
1251 // Modify MX accordingly
1252 //
1253 // Note that when Op is 0, MX is not referenced
1254 //
1255 // Parameter variables
1256 //
1257 // X : Vectors to be orthonormalized
1258 //
1259 // MX : Image of the multivector X under the operator Op
1260 //
1261 // Op : Pointer to the operator for the inner product
1262 //
1263
1264#ifdef ANASAZI_ICGS_DEBUG
1265 // Get a FancyOStream from out_arg or create a new one ...
1266 Teuchos::RCP<Teuchos::FancyOStream>
1267 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
1268 out->setShowAllFrontMatter(false).setShowProcRank(true);
1269 *out << "Entering Anasazi::ICGSOrthoManager::findBasis(...)\n";
1270#endif
1271
1272 const ScalarType ONE = SCT::one();
1273 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
1274
1275 int xc = MVT::GetNumberVecs( X );
1276
1277 if (howMany == -1) {
1278 howMany = xc;
1279 }
1280
1281 /*******************************************************
1282 * If _hasOp == false, we will not reference MX below *
1283 *******************************************************/
1284 TEUCHOS_TEST_FOR_EXCEPTION(this->_hasOp == true && MX == Teuchos::null, std::logic_error,
1285 "Anasazi::ICGSOrthoManager::findBasis(): calling routine did not specify MS.");
1286 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::logic_error,
1287 "Anasazi::ICGSOrthoManager::findBasis(): Invalid howMany parameter" );
1288
1289 /* xstart is which column we are starting the process with, based on howMany
1290 * columns before xstart are assumed to be Op-orthonormal already
1291 */
1292 int xstart = xc - howMany;
1293
1294 for (int j = xstart; j < xc; j++) {
1295
1296 // numX represents the number of currently orthonormal columns of X
1297 int numX = j;
1298 // j represents the index of the current column of X
1299 // these are different interpretations of the same value
1300
1301 //
1302 // set the lower triangular part of B to zero
1303 for (int i=j+1; i<xc; ++i) {
1304 B(i,j) = ZERO;
1305 }
1306
1307 // Get a view of the vector currently being worked on.
1308 std::vector<int> index(1);
1309 index[0] = j;
1310 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
1311 Teuchos::RCP<MV> MXj;
1312 if ((this->_hasOp)) {
1313 // MXj is a view of the current vector in MX
1314 MXj = MVT::CloneViewNonConst( *MX, index );
1315 }
1316 else {
1317 // MXj is a pointer to Xj, and MUST NOT be modified
1318 MXj = Xj;
1319 }
1320
1321 // Get a view of the previous vectors.
1322 std::vector<int> prev_idx( numX );
1323 Teuchos::RCP<const MV> prevX, prevMX;
1324
1325 if (numX > 0) {
1326 for (int i=0; i<numX; ++i) prev_idx[i] = i;
1327 prevX = MVT::CloneView( X, prev_idx );
1328 if (this->_hasOp) {
1329 prevMX = MVT::CloneView( *MX, prev_idx );
1330 }
1331 }
1332
1333 bool rankDef = true;
1334 /* numTrials>0 will denote that the current vector was randomized for the purpose
1335 * of finding a basis vector, and that the coefficients of that vector should
1336 * not be stored in B
1337 */
1338 for (int numTrials = 0; numTrials < 10; numTrials++) {
1339#ifdef ANASAZI_ICGS_DEBUG
1340 *out << "Trial " << numTrials << " for vector " << j << "\n";
1341#endif
1342
1343 // Make storage for these Gram-Schmidt iterations.
1344 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
1345 std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1);
1346
1347 //
1348 // Save old MXj vector and compute Op-norm
1349 //
1350 Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj );
1352#ifdef ANASAZI_ICGS_DEBUG
1353 *out << "origNorm = " << origNorm[0] << "\n";
1354#endif
1355
1356 if (numX > 0) {
1357 // Apply the first step of Gram-Schmidt
1358
1359 // product <- prevX^T MXj
1360 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*prevX,*Xj,product,Teuchos::null,MXj);
1361
1362 // Xj <- Xj - prevX prevX^T MXj
1363 // = Xj - prevX product
1364#ifdef ANASAZI_ICGS_DEBUG
1365 *out << "Orthogonalizing X[" << j << "]...\n";
1366#endif
1367 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
1368
1369 // Update MXj
1370 if (this->_hasOp) {
1371 // MXj <- Op*Xj_new
1372 // = Op*(Xj_old - prevX prevX^T MXj)
1373 // = MXj - prevMX product
1374#ifdef ANASAZI_ICGS_DEBUG
1375 *out << "Updating MX[" << j << "]...\n";
1376#endif
1377 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
1378 }
1379
1380 // Compute new Op-norm
1382 MagnitudeType product_norm = product.normOne();
1383
1384#ifdef ANASAZI_ICGS_DEBUG
1385 *out << "newNorm = " << newNorm[0] << "\n";
1386 *out << "prodoct_norm = " << product_norm << "\n";
1387#endif
1388
1389 // Check if a correction is needed.
1390 if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) {
1391#ifdef ANASAZI_ICGS_DEBUG
1392 if (product_norm/newNorm[0] >= tol_) {
1393 *out << "product_norm/newNorm == " << product_norm/newNorm[0] << "... another step of Gram-Schmidt.\n";
1394 }
1395 else {
1396 *out << "eps*origNorm == " << eps_*origNorm[0] << "... another step of Gram-Schmidt.\n";
1397 }
1398#endif
1399 // Apply the second step of Gram-Schmidt
1400 // This is the same as above
1401 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
1402 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*prevX,*Xj,P2,Teuchos::null,MXj);
1403 product += P2;
1404#ifdef ANASAZI_ICGS_DEBUG
1405 *out << "Orthogonalizing X[" << j << "]...\n";
1406#endif
1407 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1408 if ((this->_hasOp)) {
1409#ifdef ANASAZI_ICGS_DEBUG
1410 *out << "Updating MX[" << j << "]...\n";
1411#endif
1412 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1413 }
1414 // Compute new Op-norms
1416 product_norm = P2.normOne();
1417#ifdef ANASAZI_ICGS_DEBUG
1418 *out << "newNorm2 = " << newNorm2[0] << "\n";
1419 *out << "product_norm = " << product_norm << "\n";
1420#endif
1421 if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) {
1422 // we don't do another GS, we just set it to zero.
1423#ifdef ANASAZI_ICGS_DEBUG
1424 if (product_norm/newNorm2[0] >= tol_) {
1425 *out << "product_norm/newNorm2 == " << product_norm/newNorm2[0] << "... setting vector to zero.\n";
1426 }
1427 else if (newNorm[0] < newNorm2[0]) {
1428 *out << "newNorm2 > newNorm... setting vector to zero.\n";
1429 }
1430 else {
1431 *out << "eps*origNorm == " << eps_*origNorm[0] << "... setting vector to zero.\n";
1432 }
1433#endif
1434 MVT::MvInit(*Xj,ZERO);
1435 if ((this->_hasOp)) {
1436#ifdef ANASAZI_ICGS_DEBUG
1437 *out << "Setting MX[" << j << "] to zero as well.\n";
1438#endif
1439 MVT::MvInit(*MXj,ZERO);
1440 }
1441 }
1442 }
1443 } // if (numX > 0) do GS
1444
1445 // save the coefficients, if we are working on the original vector and not a randomly generated one
1446 if (numTrials == 0) {
1447 for (int i=0; i<numX; i++) {
1448 B(i,j) = product(i,0);
1449 }
1450 }
1451
1452 // Check if Xj has any directional information left after the orthogonalization.
1454 if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) {
1455#ifdef ANASAZI_ICGS_DEBUG
1456 *out << "Normalizing X[" << j << "], norm(X[" << j << "]) = " << newNorm[0] << "\n";
1457#endif
1458 // Normalize Xj.
1459 // Xj <- Xj / norm
1460 MVT::MvScale( *Xj, ONE/newNorm[0]);
1461 if (this->_hasOp) {
1462#ifdef ANASAZI_ICGS_DEBUG
1463 *out << "Normalizing M*X[" << j << "]...\n";
1464#endif
1465 // Update MXj.
1466 MVT::MvScale( *MXj, ONE/newNorm[0]);
1467 }
1468
1469 // save it, if it corresponds to the original vector and not a randomly generated one
1470 if (numTrials == 0) {
1471 B(j,j) = newNorm[0];
1472 }
1473
1474 // We are not rank deficient in this vector. Move on to the next vector in X.
1475 rankDef = false;
1476 break;
1477 }
1478 else {
1479#ifdef ANASAZI_ICGS_DEBUG
1480 *out << "Not normalizing M*X[" << j << "]...\n";
1481#endif
1482 // There was nothing left in Xj after orthogonalizing against previous columns in X.
1483 // X is rank deficient.
1484 // reflect this in the coefficients
1485 B(j,j) = ZERO;
1486
1487 if (completeBasis) {
1488 // Fill it with random information and keep going.
1489#ifdef ANASAZI_ICGS_DEBUG
1490 *out << "Inserting random vector in X[" << j << "]...\n";
1491#endif
1492 MVT::MvRandom( *Xj );
1493 if (this->_hasOp) {
1494#ifdef ANASAZI_ICGS_DEBUG
1495 *out << "Updating M*X[" << j << "]...\n";
1496#endif
1497 OPT::Apply( *(this->_Op), *Xj, *MXj );
1498 this->_OpCounter += MVT::GetNumberVecs(*Xj);
1499 }
1500 }
1501 else {
1502 rankDef = true;
1503 break;
1504 }
1505 }
1506 } // for (numTrials = 0; numTrials < 10; ++numTrials)
1507
1508 // if rankDef == true, then quit and notify user of rank obtained
1509 if (rankDef == true) {
1510 TEUCHOS_TEST_FOR_EXCEPTION( rankDef && completeBasis, OrthoError,
1511 "Anasazi::ICGSOrthoManager::findBasis(): Unable to complete basis" );
1512#ifdef ANASAZI_ICGS_DEBUG
1513 *out << "Returning early, rank " << j << " from Anasazi::ICGSOrthoManager::findBasis(...)\n";
1514#endif
1515 return j;
1516 }
1517
1518 } // for (j = 0; j < xc; ++j)
1519
1520#ifdef ANASAZI_ICGS_DEBUG
1521 *out << "Returning " << xc << " from Anasazi::ICGSOrthoManager::findBasis(...)\n";
1522#endif
1523 return xc;
1524 }
1525
1526} // namespace Anasazi
1527
1528#endif // ANASAZI_ICSG_ORTHOMANAGER_HPP
1529
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.
An implementation of the Anasazi::GenOrthoManager that performs orthogonalization using iterated clas...
void setNumIters(int numIters)
Set parameter for re-orthogonalization threshold.
ICGSOrthoManager(Teuchos::RCP< const OP > Op=Teuchos::null, int numIters=2, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType eps=0.0, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType tol=0.20)
Constructor specifying the operator defining the inner product as well as the number of orthogonaliza...
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 ,...
void projectGen(MV &S, Teuchos::Array< Teuchos::RCP< const MV > > X, Teuchos::Array< Teuchos::RCP< const MV > > Y, bool isBiOrtho, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C=Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)), Teuchos::RCP< MV > MS=Teuchos::null, Teuchos::Array< Teuchos::RCP< const MV > > MX=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null)), Teuchos::Array< Teuchos::RCP< const MV > > MY=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null))) const
Applies a series of generic projectors.
int getNumIters() const
Return parameter for re-orthogonalization threshold.
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...
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 .
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP< const MV > MX1, Teuchos::RCP< const MV > MX2) const
This method computes the error in orthogonality of two multivectors, measured as the Frobenius norm o...
int projectAndNormalizeGen(MV &S, Teuchos::Array< Teuchos::RCP< const MV > > X, Teuchos::Array< Teuchos::RCP< const MV > > Y, bool isBiOrtho, 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 > MS=Teuchos::null, Teuchos::Array< Teuchos::RCP< const MV > > MX=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null)), Teuchos::Array< Teuchos::RCP< const MV > > MY=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null))) const
Applies a series of generic projectors and returns an orthonormal basis for the residual data.
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.
Exception thrown to signal error in an orthogonalization manager method.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.