Anasazi Version of the Day
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"
58#include "AnasaziMultiVecTraits.hpp"
59#include "AnasaziOperatorTraits.hpp"
60#include "AnasaziMatOrthoManager.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
91 ~BasicOrthoManager() {}
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
251 int projectAndNormalizeMat (
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>
315 BasicOrthoManager<ScalarType,MV,OP>::BasicOrthoManager( Teuchos::RCP<const OP> 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);
347 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(X,X,xTx,MX,MX);
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);
364 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(X1,X2,xTx,MX1,MX2);
365 return xTx.normFrobenius();
366 }
367
368
369
371 template<class ScalarType, class MV, class OP>
372 void BasicOrthoManager<ScalarType, MV, OP>::projectMat(
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 int baslen = 0;
450 for (int i=0; i<nq; i++) {
451 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
452 "Anasazi::BasicOrthoManager::projectMat(): Q lengths not mutually consistent" );
453 qcs[i] = MVT::GetNumberVecs( *Q[i] );
454 TEUCHOS_TEST_FOR_EXCEPTION( qr < static_cast<ptrdiff_t>(qcs[i]), std::invalid_argument,
455 "Anasazi::BasicOrthoManager::projectMat(): Q has less rows than columns" );
456 baslen += qcs[i];
457
458 // check size of C[i]
459 if ( C[i] == Teuchos::null ) {
460 C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
461 }
462 else {
463 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
464 "Anasazi::BasicOrthoManager::projectMat(): Size of Q not consistent with size of C" );
465 }
466 }
467
468 // Perform the Gram-Schmidt transformation for a block of vectors
469
470 // Compute the initial Op-norms
471 std::vector<ScalarType> oldDot( xc );
472 MVT::MvDot( X, *MX, oldDot );
473#ifdef ANASAZI_BASIC_ORTHO_DEBUG
474 *out << "oldDot = { ";
475 std::copy(oldDot.begin(), oldDot.end(), std::ostream_iterator<ScalarType>(*out, " "));
476 *out << "}\n";
477#endif
478
479 MQ.resize(nq);
480 // Define the product Q^T * (Op*X)
481 for (int i=0; i<nq; i++) {
482 // Multiply Q' with MX
483 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*Q[i],X,*C[i],MQ[i],MX);
484 // Multiply by Q and subtract the result in X
485#ifdef ANASAZI_BASIC_ORTHO_DEBUG
486 *out << "Applying projector P_Q[" << i << "]...\n";
487#endif
488 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
489
490 // Update MX, with the least number of applications of Op as possible
491 // Update MX. If we have MQ, use it. Otherwise, just multiply by Op
492 if (this->_hasOp) {
493 if (MQ[i] == Teuchos::null) {
494#ifdef ANASAZI_BASIC_ORTHO_DEBUG
495 *out << "Updating MX via M*X...\n";
496#endif
497 OPT::Apply( *(this->_Op), X, *MX);
498 this->_OpCounter += MVT::GetNumberVecs(X);
499 }
500 else {
501#ifdef ANASAZI_BASIC_ORTHO_DEBUG
502 *out << "Updating MX via M*Q...\n";
503#endif
504 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
505 }
506 }
507 }
508
509 // Compute new Op-norms
510 std::vector<ScalarType> newDot(xc);
511 MVT::MvDot( X, *MX, newDot );
512#ifdef ANASAZI_BASIC_ORTHO_DEBUG
513 *out << "newDot = { ";
514 std::copy(newDot.begin(), newDot.end(), std::ostream_iterator<ScalarType>(*out, " "));
515 *out << "}\n";
516#endif
517
518 // determine (individually) whether to do another step of classical Gram-Schmidt
519 for (int j = 0; j < xc; ++j) {
520
521 if ( SCT::magnitude(kappa_*newDot[j]) < SCT::magnitude(oldDot[j]) ) {
522#ifdef ANASAZI_BASIC_ORTHO_DEBUG
523 *out << "kappa_*newDot[" <<j<< "] == " << kappa_*newDot[j] << "... another step of Gram-Schmidt.\n";
524#endif
525#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
526 Teuchos::TimeMonitor lcltimer( *timerReortho_ );
527#endif
528 for (int i=0; i<nq; i++) {
529 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(), C[i]->numCols());
530
531 // Apply another step of classical Gram-Schmidt
532 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*Q[i],X,C2,MQ[i],MX);
533 *C[i] += C2;
534#ifdef ANASAZI_BASIC_ORTHO_DEBUG
535 *out << "Applying projector P_Q[" << i << "]...\n";
536#endif
537 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
538
539 // Update MX as above
540 if (this->_hasOp) {
541 if (MQ[i] == Teuchos::null) {
542#ifdef ANASAZI_BASIC_ORTHO_DEBUG
543 *out << "Updating MX via M*X...\n";
544#endif
545 OPT::Apply( *(this->_Op), X, *MX);
546 this->_OpCounter += MVT::GetNumberVecs(X);
547 }
548 else {
549#ifdef ANASAZI_BASIC_ORTHO_DEBUG
550 *out << "Updating MX via M*Q...\n";
551#endif
552 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
553 }
554 }
555 }
556 break;
557 } // if (kappa_*newDot[j] < oldDot[j])
558 } // for (int j = 0; j < xc; ++j)
559
560#ifdef ANASAZI_BASIC_ORTHO_DEBUG
561 *out << "Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
562#endif
563 }
564
565
566
568 // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
569 template<class ScalarType, class MV, class OP>
570 int BasicOrthoManager<ScalarType, MV, OP>::normalizeMat(
571 MV &X,
572 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
573 Teuchos::RCP<MV> MX) const {
574 // call findBasis(), with the instruction to try to generate a basis of rank numvecs(X)
575 // findBasis() requires MX
576
577 int xc = MVT::GetNumberVecs(X);
578 ptrdiff_t xr = MVT::GetGlobalLength(X);
579
580 // if Op==null, MX == X (via pointer)
581 // Otherwise, either the user passed in MX or we will allocated and compute it
582 if (this->_hasOp) {
583 if (MX == Teuchos::null) {
584 // we need to allocate space for MX
585 MX = MVT::Clone(X,xc);
586 OPT::Apply(*(this->_Op),X,*MX);
587 this->_OpCounter += MVT::GetNumberVecs(X);
588 }
589 }
590
591 // if the user doesn't want to store the coefficients,
592 // allocate some local memory for them
593 if ( B == Teuchos::null ) {
594 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
595 }
596
597 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
598 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
599
600 // check size of C, B
601 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
602 "Anasazi::BasicOrthoManager::normalizeMat(): X must be non-empty" );
603 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
604 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of B" );
605 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
606 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of MX" );
607 TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument,
608 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not feasible for normalization" );
609
610 return findBasis(X, MX, *B, true );
611 }
612
613
614
616 // Find an Op-orthonormal basis for span(X) - span(W)
617 template<class ScalarType, class MV, class OP>
618 int BasicOrthoManager<ScalarType, MV, OP>::projectAndNormalizeMat(
619 MV &X,
620 Teuchos::Array<Teuchos::RCP<const MV> > Q,
621 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
622 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
623 Teuchos::RCP<MV> MX,
624 Teuchos::Array<Teuchos::RCP<const MV> > MQ
625 ) const {
626
627#ifdef ANASAZI_BASIC_ORTHO_DEBUG
628 // Get a FancyOStream from out_arg or create a new one ...
629 Teuchos::RCP<Teuchos::FancyOStream>
630 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
631 out->setShowAllFrontMatter(false).setShowProcRank(true);
632 *out << "Entering Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
633#endif
634
635 int nq = Q.length();
636 int xc = MVT::GetNumberVecs( X );
637 ptrdiff_t xr = MVT::GetGlobalLength( X );
638 int rank;
639
640 /* if the user doesn't want to store the coefficients,
641 * allocate some local memory for them
642 */
643 if ( B == Teuchos::null ) {
644 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
645 }
646
647 /****** DO NO MODIFY *MX IF _hasOp == false ******/
648 if (this->_hasOp) {
649 if (MX == Teuchos::null) {
650 // we need to allocate space for MX
651#ifdef ANASAZI_BASIC_ORTHO_DEBUG
652 *out << "Allocating MX...\n";
653#endif
654 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
655 OPT::Apply(*(this->_Op),X,*MX);
656 this->_OpCounter += MVT::GetNumberVecs(X);
657 }
658 }
659 else {
660 // Op == I --> MX = X (ignore it if the user passed it in)
661 MX = Teuchos::rcpFromRef(X);
662 }
663
664 int mxc = MVT::GetNumberVecs( *MX );
665 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
666
667 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): X must be non-empty" );
668
669 ptrdiff_t numbas = 0;
670 for (int i=0; i<nq; i++) {
671 numbas += MVT::GetNumberVecs( *Q[i] );
672 }
673
674 // check size of B
675 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
676 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of B" );
677 // check size of X and MX
678 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
679 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): MVT returned negative dimensions for X,MX" );
680 // check size of X w.r.t. MX
681 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
682 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of MX" );
683 // check feasibility
684 TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument,
685 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Orthogonality constraints not feasible" );
686
687 // orthogonalize all of X against Q
688#ifdef ANASAZI_BASIC_ORTHO_DEBUG
689 *out << "Orthogonalizing X against Q...\n";
690#endif
691 projectMat(X,Q,C,MX,MQ);
692
693 Teuchos::SerialDenseMatrix<int,ScalarType> oldCoeff(xc,1);
694 // start working
695 rank = 0;
696 int numTries = 10; // each vector in X gets 10 random chances to escape degeneracy
697 int oldrank = -1;
698 do {
699 int curxsize = xc - rank;
700
701 // orthonormalize X, but quit if it is rank deficient
702 // we can't let findBasis generated random vectors to complete the basis,
703 // because it doesn't know about Q; we will do this ourselves below
704#ifdef ANASAZI_BASIC_ORTHO_DEBUG
705 *out << "Attempting to find orthonormal basis for X...\n";
706#endif
707 rank = findBasis(X,MX,*B,false,curxsize);
708
709 if (oldrank != -1 && rank != oldrank) {
710 // we had previously stopped before, after operating on vector oldrank
711 // we saved its coefficients, augmented it with a random vector, and
712 // then called findBasis() again, which proceeded to add vector oldrank
713 // to the basis.
714 // now, restore the saved coefficients into B
715 for (int i=0; i<xc; i++) {
716 (*B)(i,oldrank) = oldCoeff(i,0);
717 }
718 }
719
720 if (rank < xc) {
721 if (rank != oldrank) {
722 // we quit on this vector and will augment it with random below
723 // this is the first time that we have quit on this vector
724 // therefor, (*B)(:,rank) contains the actual coefficients of the
725 // input vectors with respect to the previous vectors in the basis
726 // save these values, as (*B)(:,rank) will be overwritten by our next
727 // call to findBasis()
728 // we will restore it after we are done working on this vector
729 for (int i=0; i<xc; i++) {
730 oldCoeff(i,0) = (*B)(i,rank);
731 }
732 }
733 }
734
735 if (rank == xc) {
736 // we are done
737#ifdef ANASAZI_BASIC_ORTHO_DEBUG
738 *out << "Finished computing basis.\n";
739#endif
740 break;
741 }
742 else {
743 TEUCHOS_TEST_FOR_EXCEPTION( rank < oldrank, OrthoError,
744 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): basis lost rank; this shouldn't happen");
745
746 if (rank != oldrank) {
747 // we added a vector to the basis; reset the chance counter
748 numTries = 10;
749 // store old rank
750 oldrank = rank;
751 }
752 else {
753 // has this vector run out of chances to escape degeneracy?
754 if (numTries <= 0) {
755 break;
756 }
757 }
758 // use one of this vector's chances
759 numTries--;
760
761 // randomize troubled direction
762#ifdef ANASAZI_BASIC_ORTHO_DEBUG
763 *out << "Randomizing X[" << rank << "]...\n";
764#endif
765 Teuchos::RCP<MV> curX, curMX;
766 std::vector<int> ind(1);
767 ind[0] = rank;
768 curX = MVT::CloneViewNonConst(X,ind);
769 MVT::MvRandom(*curX);
770 if (this->_hasOp) {
771#ifdef ANASAZI_BASIC_ORTHO_DEBUG
772 *out << "Applying operator to random vector.\n";
773#endif
774 curMX = MVT::CloneViewNonConst(*MX,ind);
775 OPT::Apply( *(this->_Op), *curX, *curMX );
776 this->_OpCounter += MVT::GetNumberVecs(*curX);
777 }
778
779 // orthogonalize against Q
780 // if !this->_hasOp, the curMX will be ignored.
781 // we don't care about these coefficients
782 // on the contrary, we need to preserve the previous coeffs
783 {
784 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC(0);
785 projectMat(*curX,Q,dummyC,curMX,MQ);
786 }
787 }
788 } while (1);
789
790 // this should never raise an exception; but our post-conditions oblige us to check
791 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
792 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Debug error in rank variable." );
793
794#ifdef ANASAZI_BASIC_ORTHO_DEBUG
795 *out << "Leaving Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
796#endif
797
798 return rank;
799 }
800
801
802
804 // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that
805 // the rank is numvectors(X)
806 template<class ScalarType, class MV, class OP>
807 int BasicOrthoManager<ScalarType, MV, OP>::findBasis(
808 MV &X, Teuchos::RCP<MV> MX,
809 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
810 bool completeBasis, int howMany ) const {
811
812 // For the inner product defined by the operator Op or the identity (Op == 0)
813 // -> Orthonormalize X
814 // Modify MX accordingly
815 //
816 // Note that when Op is 0, MX is not referenced
817 //
818 // Parameter variables
819 //
820 // X : Vectors to be orthonormalized
821 //
822 // MX : Image of the multivector X under the operator Op
823 //
824 // Op : Pointer to the operator for the inner product
825 //
826
827#ifdef ANASAZI_BASIC_ORTHO_DEBUG
828 // Get a FancyOStream from out_arg or create a new one ...
829 Teuchos::RCP<Teuchos::FancyOStream>
830 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
831 out->setShowAllFrontMatter(false).setShowProcRank(true);
832 *out << "Entering Anasazi::BasicOrthoManager::findBasis(...)\n";
833#endif
834
835 const ScalarType ONE = SCT::one();
836 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
837
838 int xc = MVT::GetNumberVecs( X );
839
840 if (howMany == -1) {
841 howMany = xc;
842 }
843
844 /*******************************************************
845 * If _hasOp == false, we will not reference MX below *
846 *******************************************************/
847 TEUCHOS_TEST_FOR_EXCEPTION(this->_hasOp == true && MX == Teuchos::null, std::logic_error,
848 "Anasazi::BasicOrthoManager::findBasis(): calling routine did not specify MS.");
849 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::logic_error,
850 "Anasazi::BasicOrthoManager::findBasis(): Invalid howMany parameter" );
851
852 /* xstart is which column we are starting the process with, based on howMany
853 * columns before xstart are assumed to be Op-orthonormal already
854 */
855 int xstart = xc - howMany;
856
857 for (int j = xstart; j < xc; j++) {
858
859 // numX represents the number of currently orthonormal columns of X
860 int numX = j;
861 // j represents the index of the current column of X
862 // these are different interpretations of the same value
863
864 //
865 // set the lower triangular part of B to zero
866 for (int i=j+1; i<xc; ++i) {
867 B(i,j) = ZERO;
868 }
869
870 // Get a view of the vector currently being worked on.
871 std::vector<int> index(1);
872 index[0] = j;
873 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
874 Teuchos::RCP<MV> MXj;
875 if ((this->_hasOp)) {
876 // MXj is a view of the current vector in MX
877 MXj = MVT::CloneViewNonConst( *MX, index );
878 }
879 else {
880 // MXj is a pointer to Xj, and MUST NOT be modified
881 MXj = Xj;
882 }
883
884 // Get a view of the previous vectors.
885 std::vector<int> prev_idx( numX );
886 Teuchos::RCP<const MV> prevX, prevMX;
887
888 if (numX > 0) {
889 for (int i=0; i<numX; ++i) prev_idx[i] = i;
890 prevX = MVT::CloneViewNonConst( X, prev_idx );
891 if (this->_hasOp) {
892 prevMX = MVT::CloneViewNonConst( *MX, prev_idx );
893 }
894 }
895
896 bool rankDef = true;
897 /* numTrials>0 will denote that the current vector was randomized for the purpose
898 * of finding a basis vector, and that the coefficients of that vector should
899 * not be stored in B
900 */
901 for (int numTrials = 0; numTrials < 10; numTrials++) {
902#ifdef ANASAZI_BASIC_ORTHO_DEBUG
903 *out << "Trial " << numTrials << " for vector " << j << "\n";
904#endif
905
906 // Make storage for these Gram-Schmidt iterations.
907 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
908 std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1);
909
910 //
911 // Save old MXj vector and compute Op-norm
912 //
913 Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj );
914 MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,origNorm,MXj);
915#ifdef ANASAZI_BASIC_ORTHO_DEBUG
916 *out << "origNorm = " << origNorm[0] << "\n";
917#endif
918
919 if (numX > 0) {
920 // Apply the first step of Gram-Schmidt
921
922 // product <- prevX^T MXj
923 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*prevX,*Xj,product,Teuchos::null,MXj);
924
925 // Xj <- Xj - prevX prevX^T MXj
926 // = Xj - prevX product
927#ifdef ANASAZI_BASIC_ORTHO_DEBUG
928 *out << "Orthogonalizing X[" << j << "]...\n";
929#endif
930 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
931
932 // Update MXj
933 if (this->_hasOp) {
934 // MXj <- Op*Xj_new
935 // = Op*(Xj_old - prevX prevX^T MXj)
936 // = MXj - prevMX product
937#ifdef ANASAZI_BASIC_ORTHO_DEBUG
938 *out << "Updating MX[" << j << "]...\n";
939#endif
940 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
941 }
942
943 // Compute new Op-norm
944 MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,newNorm,MXj);
945 MagnitudeType product_norm = product.normOne();
946
947#ifdef ANASAZI_BASIC_ORTHO_DEBUG
948 *out << "newNorm = " << newNorm[0] << "\n";
949 *out << "prodoct_norm = " << product_norm << "\n";
950#endif
951
952 // Check if a correction is needed.
953 if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) {
954#ifdef ANASAZI_BASIC_ORTHO_DEBUG
955 if (product_norm/newNorm[0] >= tol_) {
956 *out << "product_norm/newNorm == " << product_norm/newNorm[0] << "... another step of Gram-Schmidt.\n";
957 }
958 else {
959 *out << "eps*origNorm == " << eps_*origNorm[0] << "... another step of Gram-Schmidt.\n";
960 }
961#endif
962 // Apply the second step of Gram-Schmidt
963 // This is the same as above
964 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
965 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*prevX,*Xj,P2,Teuchos::null,MXj);
966 product += P2;
967#ifdef ANASAZI_BASIC_ORTHO_DEBUG
968 *out << "Orthogonalizing X[" << j << "]...\n";
969#endif
970 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
971 if ((this->_hasOp)) {
972#ifdef ANASAZI_BASIC_ORTHO_DEBUG
973 *out << "Updating MX[" << j << "]...\n";
974#endif
975 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
976 }
977 // Compute new Op-norms
978 MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,newNorm2,MXj);
979 product_norm = P2.normOne();
980#ifdef ANASAZI_BASIC_ORTHO_DEBUG
981 *out << "newNorm2 = " << newNorm2[0] << "\n";
982 *out << "product_norm = " << product_norm << "\n";
983#endif
984 if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) {
985 // we don't do another GS, we just set it to zero.
986#ifdef ANASAZI_BASIC_ORTHO_DEBUG
987 if (product_norm/newNorm2[0] >= tol_) {
988 *out << "product_norm/newNorm2 == " << product_norm/newNorm2[0] << "... setting vector to zero.\n";
989 }
990 else if (newNorm[0] < newNorm2[0]) {
991 *out << "newNorm2 > newNorm... setting vector to zero.\n";
992 }
993 else {
994 *out << "eps*origNorm == " << eps_*origNorm[0] << "... setting vector to zero.\n";
995 }
996#endif
997 MVT::MvInit(*Xj,ZERO);
998 if ((this->_hasOp)) {
999#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1000 *out << "Setting MX[" << j << "] to zero as well.\n";
1001#endif
1002 MVT::MvInit(*MXj,ZERO);
1003 }
1004 }
1005 }
1006 } // if (numX > 0) do GS
1007
1008 // save the coefficients, if we are working on the original vector and not a randomly generated one
1009 if (numTrials == 0) {
1010 for (int i=0; i<numX; i++) {
1011 B(i,j) = product(i,0);
1012 }
1013 }
1014
1015 // Check if Xj has any directional information left after the orthogonalization.
1016 MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,newNorm,MXj);
1017 if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) {
1018#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1019 *out << "Normalizing X[" << j << "], norm(X[" << j << "]) = " << newNorm[0] << "\n";
1020#endif
1021 // Normalize Xj.
1022 // Xj <- Xj / norm
1023 MVT::MvScale( *Xj, ONE/newNorm[0]);
1024 if (this->_hasOp) {
1025#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1026 *out << "Normalizing M*X[" << j << "]...\n";
1027#endif
1028 // Update MXj.
1029 MVT::MvScale( *MXj, ONE/newNorm[0]);
1030 }
1031
1032 // save it, if it corresponds to the original vector and not a randomly generated one
1033 if (numTrials == 0) {
1034 B(j,j) = newNorm[0];
1035 }
1036
1037 // We are not rank deficient in this vector. Move on to the next vector in X.
1038 rankDef = false;
1039 break;
1040 }
1041 else {
1042#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1043 *out << "Not normalizing M*X[" << j << "]...\n";
1044#endif
1045 // There was nothing left in Xj after orthogonalizing against previous columns in X.
1046 // X is rank deficient.
1047 // reflect this in the coefficients
1048 B(j,j) = ZERO;
1049
1050 if (completeBasis) {
1051 // Fill it with random information and keep going.
1052#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1053 *out << "Inserting random vector in X[" << j << "]...\n";
1054#endif
1055 MVT::MvRandom( *Xj );
1056 if (this->_hasOp) {
1057#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1058 *out << "Updating M*X[" << j << "]...\n";
1059#endif
1060 OPT::Apply( *(this->_Op), *Xj, *MXj );
1061 this->_OpCounter += MVT::GetNumberVecs(*Xj);
1062 }
1063 }
1064 else {
1065 rankDef = true;
1066 break;
1067 }
1068 }
1069 } // for (numTrials = 0; numTrials < 10; ++numTrials)
1070
1071 // if rankDef == true, then quit and notify user of rank obtained
1072 if (rankDef == true) {
1073 TEUCHOS_TEST_FOR_EXCEPTION( rankDef && completeBasis, OrthoError,
1074 "Anasazi::BasicOrthoManager::findBasis(): Unable to complete basis" );
1075#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1076 *out << "Returning early, rank " << j << " from Anasazi::BasicOrthoManager::findBasis(...)\n";
1077#endif
1078 return j;
1079 }
1080
1081 } // for (j = 0; j < xc; ++j)
1082
1083#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1084 *out << "Returning " << xc << " from Anasazi::BasicOrthoManager::findBasis(...)\n";
1085#endif
1086 return xc;
1087 }
1088
1089} // namespace Anasazi
1090
1091#endif // ANASAZI_BASIC_ORTHOMANAGER_HPP
1092
AnasaziConfigDefs.hpp
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
AnasaziMatOrthoManager.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::BasicOrthoManager
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
Definition: AnasaziBasicOrthoManager.hpp:71
Anasazi::BasicOrthoManager::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: AnasaziBasicOrthoManager.hpp:570
Anasazi::BasicOrthoManager::getKappa
Teuchos::ScalarTraits< ScalarType >::magnitudeType getKappa() const
Return parameter for re-orthogonalization threshold.
Definition: AnasaziBasicOrthoManager.hpp:289
Anasazi::BasicOrthoManager::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: AnasaziBasicOrthoManager.hpp:618
Anasazi::BasicOrthoManager::setKappa
void setKappa(typename Teuchos::ScalarTraits< ScalarType >::magnitudeType kappa)
Set parameter for re-orthogonalization threshold.
Definition: AnasaziBasicOrthoManager.hpp:286
Anasazi::BasicOrthoManager::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: AnasaziBasicOrthoManager.hpp:372
Anasazi::BasicOrthoManager::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: AnasaziBasicOrthoManager.hpp:343
Anasazi::BasicOrthoManager::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: AnasaziBasicOrthoManager.hpp:360
Anasazi::BasicOrthoManager::BasicOrthoManager
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.
Definition: AnasaziBasicOrthoManager.hpp:315
Anasazi::MatOrthoManager
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
Definition: AnasaziMatOrthoManager.hpp:76
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:422
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:383
Anasazi::MultiVecTraits
Traits class which defines basic operations on multivectors.
Definition: AnasaziMultiVecTraits.hpp:127
Anasazi::OperatorTraits
Virtual base class which defines basic traits for the operator type.
Definition: AnasaziOperatorTraits.hpp:85
Anasazi::OrthoError
Exception thrown to signal error in an orthogonalization manager method.
Definition: AnasaziOrthoManager.hpp:82
Anasazi
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
Generated by doxygen 1.9.4