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