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