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