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;
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.
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
Exception thrown to signal error in an orthogonalization manager method.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.