Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziBasicOrthoManager.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_BASIC_ORTHOMANAGER_HPP
48#define ANASAZI_BASIC_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_BASIC_ORTHO_DEBUG
65# include <Teuchos_FancyOStream.hpp>
66#endif
67
68namespace Anasazi {
69
70 template<class ScalarType, class MV, class OP>
71 class BasicOrthoManager : public MatOrthoManager<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 BasicOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null,
85 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa = 1.41421356 /* sqrt(2) */,
86 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps = 0.0,
87 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol = 0.20 );
88
89
93
94
96
97
98
137 void projectMat (
138 MV &X,
139 Teuchos::Array<Teuchos::RCP<const MV> > Q,
140 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
141 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
142 Teuchos::RCP<MV> MX = Teuchos::null,
143 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
144 ) const;
145
146
185 int normalizeMat (
186 MV &X,
187 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
188 Teuchos::RCP<MV> MX = Teuchos::null
189 ) const;
190
191
252 MV &X,
253 Teuchos::Array<Teuchos::RCP<const MV> > Q,
254 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
255 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
256 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
257 Teuchos::RCP<MV> MX = Teuchos::null,
258 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
259 ) const;
260
262
264
265
270 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
271 orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const;
272
277 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
278 orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const;
279
281
283
284
286 void setKappa( typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa ) { kappa_ = kappa; }
287
289 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType getKappa() const { return kappa_; }
290
292
293 private:
295 MagnitudeType kappa_;
296 MagnitudeType eps_;
297 MagnitudeType tol_;
298
299 // ! Routine to find an orthonormal basis
300 int findBasis(MV &X, Teuchos::RCP<MV> MX,
301 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
302 bool completeBasis, int howMany = -1 ) const;
303
304 //
305 // Internal timers
306 //
307 Teuchos::RCP<Teuchos::Time> timerReortho_;
308
309 };
310
311
313 // Constructor
314 template<class ScalarType, class MV, class OP>
316 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa,
317 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps,
318 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol ) :
319 MatOrthoManager<ScalarType,MV,OP>(Op), kappa_(kappa), eps_(eps), tol_(tol)
320#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
321 , timerReortho_(Teuchos::TimeMonitor::getNewTimer("Anasazi::BasicOrthoManager::Re-orthogonalization"))
322#endif
323 {
324 TEUCHOS_TEST_FOR_EXCEPTION(eps_ < SCT::magnitude(SCT::zero()),std::invalid_argument,
325 "Anasazi::BasicOrthoManager::BasicOrthoManager(): argument \"eps\" must be non-negative.");
326 if (eps_ == 0) {
327 Teuchos::LAPACK<int,MagnitudeType> lapack;
328 eps_ = lapack.LAMCH('E');
329 eps_ = Teuchos::ScalarTraits<MagnitudeType>::pow(eps_,.75);
330 }
331 TEUCHOS_TEST_FOR_EXCEPTION(
332 tol_ < SCT::magnitude(SCT::zero()) || tol_ > SCT::magnitude(SCT::one()),
333 std::invalid_argument,
334 "Anasazi::BasicOrthoManager::BasicOrthoManager(): argument \"tol\" must be in [0,1].");
335 }
336
337
338
340 // Compute the distance from orthonormality
341 template<class ScalarType, class MV, class OP>
342 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
343 BasicOrthoManager<ScalarType,MV,OP>::orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX) const {
344 const ScalarType ONE = SCT::one();
345 int rank = MVT::GetNumberVecs(X);
346 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
348 for (int i=0; i<rank; i++) {
349 xTx(i,i) -= ONE;
350 }
351 return xTx.normFrobenius();
352 }
353
354
355
357 // Compute the distance from orthogonality
358 template<class ScalarType, class MV, class OP>
359 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
360 BasicOrthoManager<ScalarType,MV,OP>::orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const {
361 int r1 = MVT::GetNumberVecs(X1);
362 int r2 = MVT::GetNumberVecs(X2);
363 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
365 return xTx.normFrobenius();
366 }
367
368
369
371 template<class ScalarType, class MV, class OP>
373 MV &X,
374 Teuchos::Array<Teuchos::RCP<const MV> > Q,
375 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
376 Teuchos::RCP<MV> MX,
377 Teuchos::Array<Teuchos::RCP<const MV> > MQ
378 ) const {
379 // For the inner product defined by the operator Op or the identity (Op == 0)
380 // -> Orthogonalize X against each Q[i]
381 // Modify MX accordingly
382 //
383 // Note that when Op is 0, MX is not referenced
384 //
385 // Parameter variables
386 //
387 // X : Vectors to be transformed
388 //
389 // MX : Image of the block vector X by the mass matrix
390 //
391 // Q : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
392 //
393
394#ifdef ANASAZI_BASIC_ORTHO_DEBUG
395 // Get a FancyOStream from out_arg or create a new one ...
396 Teuchos::RCP<Teuchos::FancyOStream>
397 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
398 out->setShowAllFrontMatter(false).setShowProcRank(true);
399 *out << "Entering Anasazi::BasicOrthoManager::projectMat(...)\n";
400#endif
401
402 ScalarType ONE = SCT::one();
403
404 int xc = MVT::GetNumberVecs( X );
405 ptrdiff_t xr = MVT::GetGlobalLength( X );
406 int nq = Q.length();
407 std::vector<int> qcs(nq);
408 // short-circuit
409 if (nq == 0 || xc == 0 || xr == 0) {
410#ifdef ANASAZI_BASIC_ORTHO_DEBUG
411 *out << "Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
412#endif
413 return;
414 }
415 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
416 // if we don't have enough C, expand it with null references
417 // if we have too many, resize to throw away the latter ones
418 // if we have exactly as many as we have Q, this call has no effect
419 C.resize(nq);
420
421
422 /****** DO NO MODIFY *MX IF _hasOp == false ******/
423 if (this->_hasOp) {
424#ifdef ANASAZI_BASIC_ORTHO_DEBUG
425 *out << "Allocating MX...\n";
426#endif
427 if (MX == Teuchos::null) {
428 // we need to allocate space for MX
429 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
430 OPT::Apply(*(this->_Op),X,*MX);
431 this->_OpCounter += MVT::GetNumberVecs(X);
432 }
433 }
434 else {
435 // Op == I --> MX = X (ignore it if the user passed it in)
436 MX = Teuchos::rcpFromRef(X);
437 }
438 int mxc = MVT::GetNumberVecs( *MX );
439 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
440
441 // check size of X and Q w.r.t. common sense
442 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
443 "Anasazi::BasicOrthoManager::projectMat(): MVT returned negative dimensions for X,MX" );
444 // check size of X w.r.t. MX and Q
445 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
446 "Anasazi::BasicOrthoManager::projectMat(): Size of X not consistent with MX,Q" );
447
448 // tally up size of all Q and check/allocate C
449 for (int i=0; i<nq; i++) {
450 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
451 "Anasazi::BasicOrthoManager::projectMat(): Q lengths not mutually consistent" );
452 qcs[i] = MVT::GetNumberVecs( *Q[i] );
453 TEUCHOS_TEST_FOR_EXCEPTION( qr < static_cast<ptrdiff_t>(qcs[i]), std::invalid_argument,
454 "Anasazi::BasicOrthoManager::projectMat(): Q has less rows than columns" );
455 // check size of C[i]
456 if ( C[i] == Teuchos::null ) {
457 C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
458 }
459 else {
460 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
461 "Anasazi::BasicOrthoManager::projectMat(): Size of Q not consistent with size of C" );
462 }
463 }
464
465 // Perform the Gram-Schmidt transformation for a block of vectors
466
467 // Compute the initial Op-norms
468 std::vector<ScalarType> oldDot( xc );
469 MVT::MvDot( X, *MX, oldDot );
470#ifdef ANASAZI_BASIC_ORTHO_DEBUG
471 *out << "oldDot = { ";
472 std::copy(oldDot.begin(), oldDot.end(), std::ostream_iterator<ScalarType>(*out, " "));
473 *out << "}\n";
474#endif
475
476 MQ.resize(nq);
477 // Define the product Q^T * (Op*X)
478 for (int i=0; i<nq; i++) {
479 // Multiply Q' with MX
481 // Multiply by Q and subtract the result in X
482#ifdef ANASAZI_BASIC_ORTHO_DEBUG
483 *out << "Applying projector P_Q[" << i << "]...\n";
484#endif
485 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
486
487 // Update MX, with the least number of applications of Op as possible
488 // Update MX. If we have MQ, use it. Otherwise, just multiply by Op
489 if (this->_hasOp) {
490 if (MQ[i] == Teuchos::null) {
491#ifdef ANASAZI_BASIC_ORTHO_DEBUG
492 *out << "Updating MX via M*X...\n";
493#endif
494 OPT::Apply( *(this->_Op), X, *MX);
495 this->_OpCounter += MVT::GetNumberVecs(X);
496 }
497 else {
498#ifdef ANASAZI_BASIC_ORTHO_DEBUG
499 *out << "Updating MX via M*Q...\n";
500#endif
501 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
502 }
503 }
504 }
505
506 // Compute new Op-norms
507 std::vector<ScalarType> newDot(xc);
508 MVT::MvDot( X, *MX, newDot );
509#ifdef ANASAZI_BASIC_ORTHO_DEBUG
510 *out << "newDot = { ";
511 std::copy(newDot.begin(), newDot.end(), std::ostream_iterator<ScalarType>(*out, " "));
512 *out << "}\n";
513#endif
514
515 // determine (individually) whether to do another step of classical Gram-Schmidt
516 for (int j = 0; j < xc; ++j) {
517
518 if ( SCT::magnitude(kappa_*newDot[j]) < SCT::magnitude(oldDot[j]) ) {
519#ifdef ANASAZI_BASIC_ORTHO_DEBUG
520 *out << "kappa_*newDot[" <<j<< "] == " << kappa_*newDot[j] << "... another step of Gram-Schmidt.\n";
521#endif
522#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
523 Teuchos::TimeMonitor lcltimer( *timerReortho_ );
524#endif
525 for (int i=0; i<nq; i++) {
526 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(), C[i]->numCols());
527
528 // Apply another step of classical Gram-Schmidt
530 *C[i] += C2;
531#ifdef ANASAZI_BASIC_ORTHO_DEBUG
532 *out << "Applying projector P_Q[" << i << "]...\n";
533#endif
534 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
535
536 // Update MX as above
537 if (this->_hasOp) {
538 if (MQ[i] == Teuchos::null) {
539#ifdef ANASAZI_BASIC_ORTHO_DEBUG
540 *out << "Updating MX via M*X...\n";
541#endif
542 OPT::Apply( *(this->_Op), X, *MX);
543 this->_OpCounter += MVT::GetNumberVecs(X);
544 }
545 else {
546#ifdef ANASAZI_BASIC_ORTHO_DEBUG
547 *out << "Updating MX via M*Q...\n";
548#endif
549 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
550 }
551 }
552 }
553 break;
554 } // if (kappa_*newDot[j] < oldDot[j])
555 } // for (int j = 0; j < xc; ++j)
556
557#ifdef ANASAZI_BASIC_ORTHO_DEBUG
558 *out << "Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
559#endif
560 }
561
562
563
565 // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
566 template<class ScalarType, class MV, class OP>
568 MV &X,
569 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
570 Teuchos::RCP<MV> MX) const {
571 // call findBasis(), with the instruction to try to generate a basis of rank numvecs(X)
572 // findBasis() requires MX
573
574 int xc = MVT::GetNumberVecs(X);
575 ptrdiff_t xr = MVT::GetGlobalLength(X);
576
577 // if Op==null, MX == X (via pointer)
578 // Otherwise, either the user passed in MX or we will allocated and compute it
579 if (this->_hasOp) {
580 if (MX == Teuchos::null) {
581 // we need to allocate space for MX
582 MX = MVT::Clone(X,xc);
583 OPT::Apply(*(this->_Op),X,*MX);
584 this->_OpCounter += MVT::GetNumberVecs(X);
585 }
586 }
587
588 // if the user doesn't want to store the coefficients,
589 // allocate some local memory for them
590 if ( B == Teuchos::null ) {
591 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
592 }
593
594 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
595 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
596
597 // check size of C, B
598 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
599 "Anasazi::BasicOrthoManager::normalizeMat(): X must be non-empty" );
600 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
601 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of B" );
602 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
603 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of MX" );
604 TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument,
605 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not feasible for normalization" );
606
607 return findBasis(X, MX, *B, true );
608 }
609
610
611
613 // Find an Op-orthonormal basis for span(X) - span(W)
614 template<class ScalarType, class MV, class OP>
616 MV &X,
617 Teuchos::Array<Teuchos::RCP<const MV> > Q,
618 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
619 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
620 Teuchos::RCP<MV> MX,
621 Teuchos::Array<Teuchos::RCP<const MV> > MQ
622 ) const {
623
624#ifdef ANASAZI_BASIC_ORTHO_DEBUG
625 // Get a FancyOStream from out_arg or create a new one ...
626 Teuchos::RCP<Teuchos::FancyOStream>
627 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
628 out->setShowAllFrontMatter(false).setShowProcRank(true);
629 *out << "Entering Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
630#endif
631
632 int nq = Q.length();
633 int xc = MVT::GetNumberVecs( X );
634 ptrdiff_t xr = MVT::GetGlobalLength( X );
635 int rank;
636
637 /* if the user doesn't want to store the coefficients,
638 * allocate some local memory for them
639 */
640 if ( B == Teuchos::null ) {
641 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
642 }
643
644 /****** DO NO MODIFY *MX IF _hasOp == false ******/
645 if (this->_hasOp) {
646 if (MX == Teuchos::null) {
647 // we need to allocate space for MX
648#ifdef ANASAZI_BASIC_ORTHO_DEBUG
649 *out << "Allocating MX...\n";
650#endif
651 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
652 OPT::Apply(*(this->_Op),X,*MX);
653 this->_OpCounter += MVT::GetNumberVecs(X);
654 }
655 }
656 else {
657 // Op == I --> MX = X (ignore it if the user passed it in)
658 MX = Teuchos::rcpFromRef(X);
659 }
660
661 int mxc = MVT::GetNumberVecs( *MX );
662 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
663
664 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): X must be non-empty" );
665
666 ptrdiff_t numbas = 0;
667 for (int i=0; i<nq; i++) {
668 numbas += MVT::GetNumberVecs( *Q[i] );
669 }
670
671 // check size of B
672 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
673 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of B" );
674 // check size of X and MX
675 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
676 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): MVT returned negative dimensions for X,MX" );
677 // check size of X w.r.t. MX
678 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
679 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of MX" );
680 // check feasibility
681 TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument,
682 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Orthogonality constraints not feasible" );
683
684 // orthogonalize all of X against Q
685#ifdef ANASAZI_BASIC_ORTHO_DEBUG
686 *out << "Orthogonalizing X against Q...\n";
687#endif
688 projectMat(X,Q,C,MX,MQ);
689
690 Teuchos::SerialDenseMatrix<int,ScalarType> oldCoeff(xc,1);
691 // start working
692 rank = 0;
693 int numTries = 10; // each vector in X gets 10 random chances to escape degeneracy
694 int oldrank = -1;
695 do {
696 int curxsize = xc - rank;
697
698 // orthonormalize X, but quit if it is rank deficient
699 // we can't let findBasis generated random vectors to complete the basis,
700 // because it doesn't know about Q; we will do this ourselves below
701#ifdef ANASAZI_BASIC_ORTHO_DEBUG
702 *out << "Attempting to find orthonormal basis for X...\n";
703#endif
704 rank = findBasis(X,MX,*B,false,curxsize);
705
706 if (oldrank != -1 && rank != oldrank) {
707 // we had previously stopped before, after operating on vector oldrank
708 // we saved its coefficients, augmented it with a random vector, and
709 // then called findBasis() again, which proceeded to add vector oldrank
710 // to the basis.
711 // now, restore the saved coefficients into B
712 for (int i=0; i<xc; i++) {
713 (*B)(i,oldrank) = oldCoeff(i,0);
714 }
715 }
716
717 if (rank < xc) {
718 if (rank != oldrank) {
719 // we quit on this vector and will augment it with random below
720 // this is the first time that we have quit on this vector
721 // therefor, (*B)(:,rank) contains the actual coefficients of the
722 // input vectors with respect to the previous vectors in the basis
723 // save these values, as (*B)(:,rank) will be overwritten by our next
724 // call to findBasis()
725 // we will restore it after we are done working on this vector
726 for (int i=0; i<xc; i++) {
727 oldCoeff(i,0) = (*B)(i,rank);
728 }
729 }
730 }
731
732 if (rank == xc) {
733 // we are done
734#ifdef ANASAZI_BASIC_ORTHO_DEBUG
735 *out << "Finished computing basis.\n";
736#endif
737 break;
738 }
739 else {
740 TEUCHOS_TEST_FOR_EXCEPTION( rank < oldrank, OrthoError,
741 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): basis lost rank; this shouldn't happen");
742
743 if (rank != oldrank) {
744 // we added a vector to the basis; reset the chance counter
745 numTries = 10;
746 // store old rank
747 oldrank = rank;
748 }
749 else {
750 // has this vector run out of chances to escape degeneracy?
751 if (numTries <= 0) {
752 break;
753 }
754 }
755 // use one of this vector's chances
756 numTries--;
757
758 // randomize troubled direction
759#ifdef ANASAZI_BASIC_ORTHO_DEBUG
760 *out << "Randomizing X[" << rank << "]...\n";
761#endif
762 Teuchos::RCP<MV> curX, curMX;
763 std::vector<int> ind(1);
764 ind[0] = rank;
765 curX = MVT::CloneViewNonConst(X,ind);
766 MVT::MvRandom(*curX);
767 if (this->_hasOp) {
768#ifdef ANASAZI_BASIC_ORTHO_DEBUG
769 *out << "Applying operator to random vector.\n";
770#endif
771 curMX = MVT::CloneViewNonConst(*MX,ind);
772 OPT::Apply( *(this->_Op), *curX, *curMX );
773 this->_OpCounter += MVT::GetNumberVecs(*curX);
774 }
775
776 // orthogonalize against Q
777 // if !this->_hasOp, the curMX will be ignored.
778 // we don't care about these coefficients
779 // on the contrary, we need to preserve the previous coeffs
780 {
781 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC(0);
782 projectMat(*curX,Q,dummyC,curMX,MQ);
783 }
784 }
785 } while (1);
786
787 // this should never raise an exception; but our post-conditions oblige us to check
788 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
789 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Debug error in rank variable." );
790
791#ifdef ANASAZI_BASIC_ORTHO_DEBUG
792 *out << "Leaving Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
793#endif
794
795 return rank;
796 }
797
798
799
801 // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that
802 // the rank is numvectors(X)
803 template<class ScalarType, class MV, class OP>
805 MV &X, Teuchos::RCP<MV> MX,
806 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
807 bool completeBasis, int howMany ) const {
808
809 // For the inner product defined by the operator Op or the identity (Op == 0)
810 // -> Orthonormalize X
811 // Modify MX accordingly
812 //
813 // Note that when Op is 0, MX is not referenced
814 //
815 // Parameter variables
816 //
817 // X : Vectors to be orthonormalized
818 //
819 // MX : Image of the multivector X under the operator Op
820 //
821 // Op : Pointer to the operator for the inner product
822 //
823
824#ifdef ANASAZI_BASIC_ORTHO_DEBUG
825 // Get a FancyOStream from out_arg or create a new one ...
826 Teuchos::RCP<Teuchos::FancyOStream>
827 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
828 out->setShowAllFrontMatter(false).setShowProcRank(true);
829 *out << "Entering Anasazi::BasicOrthoManager::findBasis(...)\n";
830#endif
831
832 const ScalarType ONE = SCT::one();
833 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
834
835 int xc = MVT::GetNumberVecs( X );
836
837 if (howMany == -1) {
838 howMany = xc;
839 }
840
841 /*******************************************************
842 * If _hasOp == false, we will not reference MX below *
843 *******************************************************/
844 TEUCHOS_TEST_FOR_EXCEPTION(this->_hasOp == true && MX == Teuchos::null, std::logic_error,
845 "Anasazi::BasicOrthoManager::findBasis(): calling routine did not specify MS.");
846 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::logic_error,
847 "Anasazi::BasicOrthoManager::findBasis(): Invalid howMany parameter" );
848
849 /* xstart is which column we are starting the process with, based on howMany
850 * columns before xstart are assumed to be Op-orthonormal already
851 */
852 int xstart = xc - howMany;
853
854 for (int j = xstart; j < xc; j++) {
855
856 // numX represents the number of currently orthonormal columns of X
857 int numX = j;
858 // j represents the index of the current column of X
859 // these are different interpretations of the same value
860
861 //
862 // set the lower triangular part of B to zero
863 for (int i=j+1; i<xc; ++i) {
864 B(i,j) = ZERO;
865 }
866
867 // Get a view of the vector currently being worked on.
868 std::vector<int> index(1);
869 index[0] = j;
870 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
871 Teuchos::RCP<MV> MXj;
872 if ((this->_hasOp)) {
873 // MXj is a view of the current vector in MX
874 MXj = MVT::CloneViewNonConst( *MX, index );
875 }
876 else {
877 // MXj is a pointer to Xj, and MUST NOT be modified
878 MXj = Xj;
879 }
880
881 // Get a view of the previous vectors.
882 std::vector<int> prev_idx( numX );
883 Teuchos::RCP<const MV> prevX, prevMX;
884
885 if (numX > 0) {
886 for (int i=0; i<numX; ++i) prev_idx[i] = i;
887 prevX = MVT::CloneViewNonConst( X, prev_idx );
888 if (this->_hasOp) {
889 prevMX = MVT::CloneViewNonConst( *MX, prev_idx );
890 }
891 }
892
893 bool rankDef = true;
894 /* numTrials>0 will denote that the current vector was randomized for the purpose
895 * of finding a basis vector, and that the coefficients of that vector should
896 * not be stored in B
897 */
898 for (int numTrials = 0; numTrials < 10; numTrials++) {
899#ifdef ANASAZI_BASIC_ORTHO_DEBUG
900 *out << "Trial " << numTrials << " for vector " << j << "\n";
901#endif
902
903 // Make storage for these Gram-Schmidt iterations.
904 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
905 std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1);
906
907 //
908 // Save old MXj vector and compute Op-norm
909 //
910 Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj );
912#ifdef ANASAZI_BASIC_ORTHO_DEBUG
913 *out << "origNorm = " << origNorm[0] << "\n";
914#endif
915
916 if (numX > 0) {
917 // Apply the first step of Gram-Schmidt
918
919 // product <- prevX^T MXj
920 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*prevX,*Xj,product,Teuchos::null,MXj);
921
922 // Xj <- Xj - prevX prevX^T MXj
923 // = Xj - prevX product
924#ifdef ANASAZI_BASIC_ORTHO_DEBUG
925 *out << "Orthogonalizing X[" << j << "]...\n";
926#endif
927 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
928
929 // Update MXj
930 if (this->_hasOp) {
931 // MXj <- Op*Xj_new
932 // = Op*(Xj_old - prevX prevX^T MXj)
933 // = MXj - prevMX product
934#ifdef ANASAZI_BASIC_ORTHO_DEBUG
935 *out << "Updating MX[" << j << "]...\n";
936#endif
937 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
938 }
939
940 // Compute new Op-norm
942 MagnitudeType product_norm = product.normOne();
943
944#ifdef ANASAZI_BASIC_ORTHO_DEBUG
945 *out << "newNorm = " << newNorm[0] << "\n";
946 *out << "prodoct_norm = " << product_norm << "\n";
947#endif
948
949 // Check if a correction is needed.
950 if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) {
951#ifdef ANASAZI_BASIC_ORTHO_DEBUG
952 if (product_norm/newNorm[0] >= tol_) {
953 *out << "product_norm/newNorm == " << product_norm/newNorm[0] << "... another step of Gram-Schmidt.\n";
954 }
955 else {
956 *out << "eps*origNorm == " << eps_*origNorm[0] << "... another step of Gram-Schmidt.\n";
957 }
958#endif
959 // Apply the second step of Gram-Schmidt
960 // This is the same as above
961 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
962 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*prevX,*Xj,P2,Teuchos::null,MXj);
963 product += P2;
964#ifdef ANASAZI_BASIC_ORTHO_DEBUG
965 *out << "Orthogonalizing X[" << j << "]...\n";
966#endif
967 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
968 if ((this->_hasOp)) {
969#ifdef ANASAZI_BASIC_ORTHO_DEBUG
970 *out << "Updating MX[" << j << "]...\n";
971#endif
972 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
973 }
974 // Compute new Op-norms
976 product_norm = P2.normOne();
977#ifdef ANASAZI_BASIC_ORTHO_DEBUG
978 *out << "newNorm2 = " << newNorm2[0] << "\n";
979 *out << "product_norm = " << product_norm << "\n";
980#endif
981 if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) {
982 // we don't do another GS, we just set it to zero.
983#ifdef ANASAZI_BASIC_ORTHO_DEBUG
984 if (product_norm/newNorm2[0] >= tol_) {
985 *out << "product_norm/newNorm2 == " << product_norm/newNorm2[0] << "... setting vector to zero.\n";
986 }
987 else if (newNorm[0] < newNorm2[0]) {
988 *out << "newNorm2 > newNorm... setting vector to zero.\n";
989 }
990 else {
991 *out << "eps*origNorm == " << eps_*origNorm[0] << "... setting vector to zero.\n";
992 }
993#endif
994 MVT::MvInit(*Xj,ZERO);
995 if ((this->_hasOp)) {
996#ifdef ANASAZI_BASIC_ORTHO_DEBUG
997 *out << "Setting MX[" << j << "] to zero as well.\n";
998#endif
999 MVT::MvInit(*MXj,ZERO);
1000 }
1001 }
1002 }
1003 } // if (numX > 0) do GS
1004
1005 // save the coefficients, if we are working on the original vector and not a randomly generated one
1006 if (numTrials == 0) {
1007 for (int i=0; i<numX; i++) {
1008 B(i,j) = product(i,0);
1009 }
1010 }
1011
1012 // Check if Xj has any directional information left after the orthogonalization.
1014 if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) {
1015#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1016 *out << "Normalizing X[" << j << "], norm(X[" << j << "]) = " << newNorm[0] << "\n";
1017#endif
1018 // Normalize Xj.
1019 // Xj <- Xj / norm
1020 MVT::MvScale( *Xj, ONE/newNorm[0]);
1021 if (this->_hasOp) {
1022#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1023 *out << "Normalizing M*X[" << j << "]...\n";
1024#endif
1025 // Update MXj.
1026 MVT::MvScale( *MXj, ONE/newNorm[0]);
1027 }
1028
1029 // save it, if it corresponds to the original vector and not a randomly generated one
1030 if (numTrials == 0) {
1031 B(j,j) = newNorm[0];
1032 }
1033
1034 // We are not rank deficient in this vector. Move on to the next vector in X.
1035 rankDef = false;
1036 break;
1037 }
1038 else {
1039#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1040 *out << "Not normalizing M*X[" << j << "]...\n";
1041#endif
1042 // There was nothing left in Xj after orthogonalizing against previous columns in X.
1043 // X is rank deficient.
1044 // reflect this in the coefficients
1045 B(j,j) = ZERO;
1046
1047 if (completeBasis) {
1048 // Fill it with random information and keep going.
1049#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1050 *out << "Inserting random vector in X[" << j << "]...\n";
1051#endif
1052 MVT::MvRandom( *Xj );
1053 if (this->_hasOp) {
1054#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1055 *out << "Updating M*X[" << j << "]...\n";
1056#endif
1057 OPT::Apply( *(this->_Op), *Xj, *MXj );
1058 this->_OpCounter += MVT::GetNumberVecs(*Xj);
1059 }
1060 }
1061 else {
1062 rankDef = true;
1063 break;
1064 }
1065 }
1066 } // for (numTrials = 0; numTrials < 10; ++numTrials)
1067
1068 // if rankDef == true, then quit and notify user of rank obtained
1069 if (rankDef == true) {
1070 TEUCHOS_TEST_FOR_EXCEPTION( rankDef && completeBasis, OrthoError,
1071 "Anasazi::BasicOrthoManager::findBasis(): Unable to complete basis" );
1072#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1073 *out << "Returning early, rank " << j << " from Anasazi::BasicOrthoManager::findBasis(...)\n";
1074#endif
1075 return j;
1076 }
1077
1078 } // for (j = 0; j < xc; ++j)
1079
1080#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1081 *out << "Returning " << xc << " from Anasazi::BasicOrthoManager::findBasis(...)\n";
1082#endif
1083 return xc;
1084 }
1085
1086} // namespace Anasazi
1087
1088#endif // ANASAZI_BASIC_ORTHOMANAGER_HPP
1089
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::MatOrthoManager that performs orthogonalization using (potentially)...
int normalizeMat(MV &X, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B=Teuchos::null, Teuchos::RCP< MV > MX=Teuchos::null) const
This method takes a multivector X and attempts to compute an orthonormal basis for ,...
Teuchos::ScalarTraits< ScalarType >::magnitudeType getKappa() const
Return parameter for re-orthogonalization threshold.
int projectAndNormalizeMat(MV &X, Teuchos::Array< Teuchos::RCP< const MV > > Q, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C=Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)), Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B=Teuchos::null, Teuchos::RCP< MV > MX=Teuchos::null, Teuchos::Array< Teuchos::RCP< const MV > > MQ=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null))) const
Given a set of bases Q[i] and a multivector X, this method computes an orthonormal basis for .
void setKappa(typename Teuchos::ScalarTraits< ScalarType >::magnitudeType kappa)
Set parameter for re-orthogonalization threshold.
void projectMat(MV &X, Teuchos::Array< Teuchos::RCP< const MV > > Q, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C=Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)), Teuchos::RCP< MV > MX=Teuchos::null, Teuchos::Array< Teuchos::RCP< const MV > > MQ=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null))) const
Given a list of mutually orthogonal and internally orthonormal bases Q, this method projects a multiv...
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormErrorMat(const MV &X, Teuchos::RCP< const MV > MX=Teuchos::null) const
This method computes the error in orthonormality of a multivector, measured as the Frobenius norm of ...
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...
BasicOrthoManager(Teuchos::RCP< const OP > Op=Teuchos::null, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType kappa=1.41421356, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType eps=0.0, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType tol=0.20)
Constructor specifying re-orthogonalization tolerance.
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
void normMat(const MV &X, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, Teuchos::RCP< const MV > MX=Teuchos::null) const
Provides the norm induced by the matrix-based inner product.
void innerProdMat(const MV &X, const MV &Y, Teuchos::SerialDenseMatrix< int, ScalarType > &Z, Teuchos::RCP< const MV > MX=Teuchos::null, Teuchos::RCP< const MV > MY=Teuchos::null) const
Provides a matrix-based inner product.
Exception thrown to signal error in an orthogonalization manager method.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.