36 #ifndef ANASAZI_SADDLE_CONTAINER_HPP 37 #define ANASAZI_SADDLE_CONTAINER_HPP 40 #include "Teuchos_VerboseObject.hpp" 42 #ifdef HAVE_ANASAZI_BELOS 43 #include "BelosConfigDefs.hpp" 44 #include "BelosMultiVecTraits.hpp" 53 template <
class ScalarType,
class MV>
56 template <
class Scalar_,
class MV_,
class OP_>
friend class SaddleOperator;
60 typedef Teuchos::SerialDenseMatrix<int,ScalarType> SerialDenseMatrix;
61 const ScalarType ONE, ZERO;
63 RCP<SerialDenseMatrix> lowerRaw_;
64 std::vector<int> indices_;
66 RCP<SerialDenseMatrix> getLower()
const;
67 void setLower(
const RCP<SerialDenseMatrix> lower);
71 SaddleContainer( ) : ONE(Teuchos::ScalarTraits<ScalarType>::one()), ZERO(Teuchos::ScalarTraits<ScalarType>::zero()) { };
72 SaddleContainer(
const RCP<MV> X,
bool eye=
false );
76 SaddleContainer<ScalarType, MV> * Clone(
const int nvecs)
const;
78 SaddleContainer<ScalarType, MV> * CloneCopy()
const;
80 SaddleContainer<ScalarType, MV> * CloneCopy(
const std::vector<int> &index)
const;
81 SaddleContainer<ScalarType, MV> * CloneViewNonConst(
const std::vector<int> &index)
const;
82 SaddleContainer<ScalarType, MV> * CloneViewNonConst(
const Teuchos::Range1D& index)
const;
83 const SaddleContainer<ScalarType, MV> * CloneView(
const std::vector<int> &index)
const;
84 const SaddleContainer<ScalarType, MV> * CloneView(
const Teuchos::Range1D& index)
const;
85 ptrdiff_t GetGlobalLength()
const {
return MVT::GetGlobalLength(*upper_) + lowerRaw_->numRows(); };
88 void MvTimesMatAddMv(ScalarType alpha,
const SaddleContainer<ScalarType,MV> &A,
89 const Teuchos::SerialDenseMatrix<int, ScalarType> &B,
92 void MvAddMv(ScalarType alpha,
const SaddleContainer<ScalarType,MV>& A,
93 ScalarType beta,
const SaddleContainer<ScalarType,MV>& B);
95 void MvScale( ScalarType alpha );
97 void MvScale(
const std::vector<ScalarType>& alpha );
99 void MvTransMv (ScalarType alpha,
const SaddleContainer<ScalarType, MV>& A,
100 Teuchos::SerialDenseMatrix< int, ScalarType >& B)
const;
102 void MvDot (
const SaddleContainer<ScalarType, MV>& A, std::vector<ScalarType> &b)
const;
104 void MvNorm ( std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &normvec)
const;
108 void SetBlock (
const SaddleContainer<ScalarType, MV>& A,
const std::vector<int> &index);
110 void Assign (
const SaddleContainer<ScalarType, MV>&A);
114 void MvInit (ScalarType alpha);
116 void MvPrint (std::ostream &os)
const;
122 template <
class ScalarType,
class MV>
123 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > SaddleContainer<ScalarType, MV>::getLower()
const 130 int nrows = lowerRaw_->numRows();
131 int ncols = indices_.size();
133 RCP<SerialDenseMatrix> lower = rcp(
new SerialDenseMatrix(nrows,ncols,
false));
135 for(
int r=0; r<nrows; r++)
137 for(
int c=0; c<ncols; c++)
139 (*lower)(r,c) = (*lowerRaw_)(r,indices_[c]);
149 template <
class ScalarType,
class MV>
150 void SaddleContainer<ScalarType, MV>::setLower(
const RCP<SerialDenseMatrix> lower)
158 int nrows = lowerRaw_->numRows();
159 int ncols = indices_.size();
161 for(
int r=0; r<nrows; r++)
163 for(
int c=0; c<ncols; c++)
165 (*lowerRaw_)(r,indices_[c]) = (*lower)(r,c);
173 template <
class ScalarType,
class MV>
174 SaddleContainer<ScalarType, MV>::SaddleContainer(
const RCP<MV> X,
bool eye )
175 : ONE(Teuchos::ScalarTraits<ScalarType>::one()), ZERO(Teuchos::ScalarTraits<ScalarType>::zero())
186 lowerRaw_ = rcp(
new SerialDenseMatrix(nvecs,nvecs));
187 for(
int i=0; i < nvecs; i++)
188 (*lowerRaw_)(i,i) = ONE;
196 lowerRaw_ = rcp(
new SerialDenseMatrix(nvecs,nvecs));
203 template <
class ScalarType,
class MV>
204 SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::Clone(
const int nvecs)
const 206 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
209 newSC->lowerRaw_ = rcp(
new SerialDenseMatrix(lowerRaw_->numRows(),nvecs));
217 template <
class ScalarType,
class MV>
218 SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneCopy()
const 220 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
223 newSC->lowerRaw_ = getLower();
231 template <
class ScalarType,
class MV>
232 SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneCopy(
const std::vector< int > &index)
const 234 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
238 int ncols = index.size();
239 int nrows = lowerRaw_->numRows();
240 RCP<SerialDenseMatrix> lower = getLower();
241 newSC->lowerRaw_ = rcp(
new SerialDenseMatrix(nrows,ncols));
242 for(
int c=0; c < ncols; c++)
244 for(
int r=0; r < nrows; r++)
245 (*newSC->lowerRaw_)(r,c) = (*lower)(r,index[c]);
254 template <
class ScalarType,
class MV>
255 SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneViewNonConst(
const std::vector<int> &index)
const 257 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
261 newSC->lowerRaw_ = lowerRaw_;
263 if(!indices_.empty())
265 newSC->indices_.resize(index.size());
266 for(
unsigned int i=0; i<index.size(); i++)
268 newSC->indices_[i] = indices_[index[i]];
273 newSC->indices_ = index;
280 template <
class ScalarType,
class MV>
281 SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneViewNonConst(
const Teuchos::Range1D& index)
const 283 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
287 newSC->lowerRaw_ = lowerRaw_;
289 newSC->indices_.resize(index.size());
290 for(
unsigned int i=0; i<index.size(); i++)
292 newSC->indices_[i] = indices_[index.lbound()+i];
301 template <
class ScalarType,
class MV>
302 const SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneView(
const std::vector<int> &index)
const 304 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
308 newSC->lowerRaw_ = lowerRaw_;
310 if(!indices_.empty())
312 newSC->indices_.resize(index.size());
313 for(
unsigned int i=0; i<index.size(); i++)
315 newSC->indices_[i] = indices_[index[i]];
320 newSC->indices_ = index;
327 template <
class ScalarType,
class MV>
328 const SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneView(
const Teuchos::Range1D& index)
const 330 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
334 newSC->lowerRaw_ = lowerRaw_;
336 newSC->indices_.resize(index.size());
337 for(
unsigned int i=0; i<index.size(); i++)
339 newSC->indices_[i] = indices_[index.lbound()+i];
349 template <
class ScalarType,
class MV>
350 void SaddleContainer<ScalarType, MV>::MvTimesMatAddMv(ScalarType alpha,
const SaddleContainer<ScalarType,MV> &A,
351 const Teuchos::SerialDenseMatrix<int, ScalarType> &B,
355 RCP<SerialDenseMatrix> lower = getLower();
356 RCP<SerialDenseMatrix> Alower = A.getLower();
357 lower->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,alpha,*Alower,B,beta);
364 template <
class ScalarType,
class MV>
365 void SaddleContainer<ScalarType, MV>::MvAddMv(ScalarType alpha,
const SaddleContainer<ScalarType,MV>& A,
366 ScalarType beta,
const SaddleContainer<ScalarType,MV>& B)
368 MVT::MvAddMv(alpha, *(A.upper_), beta, *(B.upper_), *upper_);
370 RCP<SerialDenseMatrix> lower = getLower();
371 RCP<SerialDenseMatrix> Alower = A.getLower();
372 RCP<SerialDenseMatrix> Blower = B.getLower();
379 lower->assign(*Alower);
385 else if(beta == -ONE)
387 else if(beta != ZERO)
389 SerialDenseMatrix scaledB(*Blower);
400 template <
class ScalarType,
class MV>
401 void SaddleContainer<ScalarType, MV>::MvScale( ScalarType alpha )
406 RCP<SerialDenseMatrix> lower = getLower();
414 template <
class ScalarType,
class MV>
415 void SaddleContainer<ScalarType, MV>::MvScale(
const std::vector<ScalarType>& alpha )
419 RCP<SerialDenseMatrix> lower = getLower();
421 int nrows = lower->numRows();
422 int ncols = lower->numCols();
424 for(
int c=0; c<ncols; c++)
426 for(
int r=0; r<nrows; r++)
427 (*lower)(r,c) *= alpha[c];
437 template <
class ScalarType,
class MV>
438 void SaddleContainer<ScalarType, MV>::MvTransMv (ScalarType alpha,
const SaddleContainer<ScalarType, MV>& A,
439 Teuchos::SerialDenseMatrix< int, ScalarType >& B)
const 442 RCP<SerialDenseMatrix> lower = getLower();
443 RCP<SerialDenseMatrix> Alower = A.getLower();
444 B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,alpha,*(Alower),*lower,ONE);
450 template <
class ScalarType,
class MV>
451 void SaddleContainer<ScalarType, MV>::MvDot (
const SaddleContainer<ScalarType, MV>& A, std::vector<ScalarType> &b)
const 455 RCP<SerialDenseMatrix> lower = getLower();
456 RCP<SerialDenseMatrix> Alower = A.getLower();
458 int nrows = lower->numRows();
459 int ncols = lower->numCols();
461 for(
int c=0; c < ncols; c++)
463 for(
int r=0; r < nrows; r++)
465 b[c] += ((*Alower)(r,c) * (*lower)(r,c));
474 template <
class ScalarType,
class MV>
475 void SaddleContainer<ScalarType, MV>::MvNorm ( std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &normvec)
const 478 MvDot(*
this,normvec);
479 for(
unsigned int i=0; i<normvec.size(); i++)
480 normvec[i] = sqrt(normvec[i]);
488 template <
class ScalarType,
class MV>
489 void SaddleContainer<ScalarType, MV>::SetBlock (
const SaddleContainer<ScalarType, MV>& A,
const std::vector<int> &index)
493 RCP<SerialDenseMatrix> lower = getLower();
494 RCP<SerialDenseMatrix> Alower = A.getLower();
496 int nrows = lower->numRows();
498 int nvecs = index.size();
499 for(
int c=0; c<nvecs; c++)
501 for(
int r=0; r<nrows; r++)
502 (*lower)(r,index[c]) = (*Alower)(r,c);
511 template <
class ScalarType,
class MV>
512 void SaddleContainer<ScalarType, MV>::Assign (
const SaddleContainer<ScalarType, MV>&A)
516 RCP<SerialDenseMatrix> lower = getLower();
517 RCP<SerialDenseMatrix> Alower = A.getLower();
528 template <
class ScalarType,
class MV>
529 void SaddleContainer<ScalarType, MV>::MvRandom ()
533 RCP<SerialDenseMatrix> lower = getLower();
541 template <
class ScalarType,
class MV>
542 void SaddleContainer<ScalarType, MV>::MvInit (ScalarType alpha)
546 RCP<SerialDenseMatrix> lower = getLower();
547 lower->putScalar(alpha);
554 template <
class ScalarType,
class MV>
555 void SaddleContainer<ScalarType, MV>::MvPrint (std::ostream &os)
const 557 RCP<SerialDenseMatrix> lower = getLower();
561 os <<
"Object SaddleContainer" << std::endl;
563 upper_->describe(*(Teuchos::VerboseObjectBase::getDefaultOStream()),Teuchos::VERB_EXTREME);
566 os <<
"Y\n" << *lower << std::endl;
571 template<
class ScalarType,
class MV >
572 class MultiVecTraits<ScalarType,
Experimental::SaddleContainer<ScalarType, MV> >
574 typedef Experimental::SaddleContainer<ScalarType,MV> SC;
577 static RCP<SC > Clone(
const SC& mv,
const int numvecs )
578 {
return rcp( const_cast<SC&>(mv).Clone(numvecs) ); }
580 static RCP<SC > CloneCopy(
const SC& mv )
581 {
return rcp( const_cast<SC&>(mv).CloneCopy() ); }
583 static RCP<SC > CloneCopy(
const SC& mv,
const std::vector<int>& index )
584 {
return rcp( const_cast<SC&>(mv).CloneCopy(index) ); }
586 static ptrdiff_t GetGlobalLength(
const SC& mv )
587 {
return mv.GetGlobalLength(); }
589 static int GetNumberVecs(
const SC& mv )
590 {
return mv.GetNumberVecs(); }
592 static void MvTimesMatAddMv( ScalarType alpha,
const SC& A,
593 const Teuchos::SerialDenseMatrix<int,ScalarType>& B,
594 ScalarType beta, SC& mv )
595 { mv.MvTimesMatAddMv(alpha, A, B, beta); }
597 static void MvAddMv( ScalarType alpha,
const SC& A, ScalarType beta,
const SC& B, SC& mv )
598 { mv.MvAddMv(alpha, A, beta, B); }
600 static void MvTransMv( ScalarType alpha,
const SC& A,
const SC& mv, Teuchos::SerialDenseMatrix<int,ScalarType>& B)
601 { mv.MvTransMv(alpha, A, B); }
603 static void MvDot(
const SC& mv,
const SC& A, std::vector<ScalarType> & b)
606 static void MvScale ( SC& mv, ScalarType alpha )
607 { mv.MvScale( alpha ); }
609 static void MvScale ( SC& mv,
const std::vector<ScalarType>& alpha )
610 { mv.MvScale( alpha ); }
612 static void MvNorm(
const SC& mv, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> & normvec)
613 { mv.MvNorm(normvec); }
615 static void SetBlock(
const SC& A,
const std::vector<int>& index, SC& mv )
616 { mv.SetBlock(A, index); }
618 static void Assign(
const SC& A, SC& mv )
621 static void MvRandom( SC& mv )
624 static void MvInit( SC& mv, ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero() )
625 { mv.MvInit(alpha); }
627 static void MvPrint(
const SC& mv, std::ostream& os )
633 #ifdef HAVE_ANASAZI_BELOS 637 template<
class ScalarType,
class MV >
638 class MultiVecTraits< ScalarType,
Anasazi::Experimental::SaddleContainer<ScalarType,MV> >
640 typedef Anasazi::Experimental::SaddleContainer<ScalarType,MV> SC;
642 static RCP<SC > Clone(
const SC& mv,
const int numvecs )
643 {
return rcp( const_cast<SC&>(mv).Clone(numvecs) ); }
645 static RCP<SC > CloneCopy(
const SC& mv )
646 {
return rcp( const_cast<SC&>(mv).CloneCopy() ); }
648 static RCP<SC > CloneCopy(
const SC& mv,
const std::vector<int>& index )
649 {
return rcp( const_cast<SC&>(mv).CloneCopy(index) ); }
651 static RCP<SC> CloneViewNonConst( SC& mv,
const std::vector<int>& index )
652 {
return rcp( mv.CloneViewNonConst(index) ); }
654 static RCP<SC> CloneViewNonConst( SC& mv,
const Teuchos::Range1D& index )
655 {
return rcp( mv.CloneViewNonConst(index) ); }
657 static RCP<const SC> CloneView(
const SC& mv,
const std::vector<int> & index )
658 {
return rcp( mv.CloneView(index) ); }
660 static RCP<const SC> CloneView(
const SC& mv,
const Teuchos::Range1D& index )
661 {
return rcp( mv.CloneView(index) ); }
663 static ptrdiff_t GetGlobalLength(
const SC& mv )
664 {
return mv.GetGlobalLength(); }
666 static int GetNumberVecs(
const SC& mv )
667 {
return mv.GetNumberVecs(); }
669 static void MvTimesMatAddMv( ScalarType alpha,
const SC& A,
670 const Teuchos::SerialDenseMatrix<int,ScalarType>& B,
671 ScalarType beta, SC& mv )
672 { mv.MvTimesMatAddMv(alpha, A, B, beta); }
674 static void MvAddMv( ScalarType alpha,
const SC& A, ScalarType beta,
const SC& B, SC& mv )
675 { mv.MvAddMv(alpha, A, beta, B); }
677 static void MvTransMv( ScalarType alpha,
const SC& A,
const SC& mv, Teuchos::SerialDenseMatrix<int,ScalarType>& B)
678 { mv.MvTransMv(alpha, A, B); }
680 static void MvDot(
const SC& mv,
const SC& A, std::vector<ScalarType> & b)
683 static void MvScale ( SC& mv, ScalarType alpha )
684 { mv.MvScale( alpha ); }
686 static void MvScale ( SC& mv,
const std::vector<ScalarType>& alpha )
687 { mv.MvScale( alpha ); }
690 static void MvNorm(
const SC& mv, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> & normvec, NormType type=TwoNorm)
691 { mv.MvNorm(normvec); }
693 static void SetBlock(
const SC& A,
const std::vector<int>& index, SC& mv )
694 { mv.SetBlock(A, index); }
696 static void Assign(
const SC& A, SC& mv )
699 static void MvRandom( SC& mv )
702 static void MvInit( SC& mv, ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero() )
703 { mv.MvInit(alpha); }
705 static void MvPrint(
const SC& mv, std::ostream& os )
708 #ifdef HAVE_BELOS_TSQR 709 typedef Belos::details::StubTsqrAdapter<SC> tsqr_adaptor_type;
722 #endif // HAVE_BELOS_TSQR static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
static void MvDot(const MV &mv, const MV &A, std::vector< ScalarType > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
static void Assign(const MV &A, MV &mv)
mv := A
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Traits class which defines basic operations on multivectors.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &B, Teuchos::SerialDenseMatrix< int, ScalarType > &C)
Compute C := alpha * A^H B.
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.