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;
61 const ScalarType ONE, ZERO;
64 std::vector<int> indices_;
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;
88 void MvTimesMatAddMv(ScalarType alpha,
const SaddleContainer<ScalarType,MV> &A,
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,
102 void MvDot (
const SaddleContainer<ScalarType, MV>& A, std::vector<ScalarType> &b)
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>
130 int nrows = lowerRaw_->numRows();
131 int ncols = indices_.size();
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>
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>();
208 newSC->upper_ = MVT::Clone(*upper_,nvecs);
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>();
222 newSC->upper_ = MVT::CloneCopy(*upper_);
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>();
236 newSC->upper_ = MVT::CloneCopy(*upper_,index);
238 int ncols = index.size();
239 int nrows = lowerRaw_->numRows();
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>();
259 newSC->upper_ = MVT::CloneViewNonConst(*upper_,index);
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>();
285 newSC->upper_ = MVT::CloneViewNonConst(*upper_,index);
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>();
306 newSC->upper_ = MVT::CloneViewNonConst(*upper_,index);
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>();
332 newSC->upper_ = MVT::CloneViewNonConst(*upper_,index);
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,
354 MVT::MvTimesMatAddMv(alpha,*(A.upper_),B,beta,*upper_);
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_);
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 )
403 MVT::MvScale(*upper_, alpha);
414 template <
class ScalarType,
class MV>
415 void SaddleContainer<ScalarType, MV>::MvScale(
const std::vector<ScalarType>& alpha )
417 MVT::MvScale(*upper_, alpha);
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,
441 MVT::MvTransMv(alpha,*(A.upper_),*upper_,B);
450 template <
class ScalarType,
class MV>
451 void SaddleContainer<ScalarType, MV>::MvDot (
const SaddleContainer<ScalarType, MV>& A, std::vector<ScalarType> &b)
const 453 MVT::MvDot(*upper_, *(A.upper_), b);
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>
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)
491 MVT::SetBlock(*(A.upper_), index, *upper_);
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)
514 MVT::Assign(*(A.upper_),*(upper_));
528 template <
class ScalarType,
class MV>
529 void SaddleContainer<ScalarType, MV>::MvRandom ()
531 MVT::MvRandom(*upper_);
541 template <
class ScalarType,
class MV>
542 void SaddleContainer<ScalarType, MV>::MvInit (ScalarType alpha)
544 MVT::MvInit(*upper_,alpha);
547 lower->putScalar(alpha);
554 template <
class ScalarType,
class MV>
555 void SaddleContainer<ScalarType, MV>::MvPrint (std::ostream &os)
const 561 os <<
"Object SaddleContainer" << std::endl;
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,
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); }
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 ); }
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 )
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) ); }
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) ); }
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,
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); }
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 ); }
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 )
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.
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
OrdinalType numRows() const
Traits class which defines basic operations on multivectors.
static RCP< FancyOStream > getDefaultOStream()
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.