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>();
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();
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>();
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,
351 const Teuchos::SerialDenseMatrix<int, ScalarType> &B,
354 MVT::MvTimesMatAddMv(alpha,*(A.upper_),B,beta,*upper_);
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 )
403 MVT::MvScale(*upper_, alpha);
406 RCP<SerialDenseMatrix> lower = getLower();
414 template <
class ScalarType,
class MV>
415 void SaddleContainer<ScalarType, MV>::MvScale(
const std::vector<ScalarType>& alpha )
417 MVT::MvScale(*upper_, 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 441 MVT::MvTransMv(alpha,*(A.upper_),*upper_,B);
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 453 MVT::MvDot(*upper_, *(A.upper_), b);
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)
491 MVT::SetBlock(*(A.upper_), index, *upper_);
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)
514 MVT::Assign(*(A.upper_),*(upper_));
516 RCP<SerialDenseMatrix> lower = getLower();
517 RCP<SerialDenseMatrix> Alower = A.getLower();
528 template <
class ScalarType,
class MV>
529 void SaddleContainer<ScalarType, MV>::MvRandom ()
531 MVT::MvRandom(*upper_);
533 RCP<SerialDenseMatrix> lower = getLower();
541 template <
class ScalarType,
class MV>
542 void SaddleContainer<ScalarType, MV>::MvInit (ScalarType alpha)
544 MVT::MvInit(*upper_,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.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Traits class which defines basic operations on multivectors.
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.