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