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