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;
76 typedef MultiVecTraits<ScalarType,MV> MVT;
77 typedef OperatorTraits<ScalarType,MV,OP> OPT;
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.
Exception thrown to signal error in an orthogonalization manager method.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.