49#ifndef ANASAZI_SADDLE_CONTAINER_HPP
50#define ANASAZI_SADDLE_CONTAINER_HPP
53#include "Teuchos_VerboseObject.hpp"
55#ifdef HAVE_ANASAZI_BELOS
56#include "BelosConfigDefs.hpp"
57#include "BelosMultiVecTraits.hpp"
66template <
class ScalarType,
class MV>
69 template <
class Scalar_,
class MV_,
class OP_>
friend class SaddleOperator;
73 typedef Teuchos::SerialDenseMatrix<int,ScalarType> SerialDenseMatrix;
74 const ScalarType ONE, ZERO;
76 RCP<SerialDenseMatrix> lowerRaw_;
77 std::vector<int> indices_;
79 RCP<SerialDenseMatrix> getLower()
const;
80 void setLower(
const RCP<SerialDenseMatrix> lower);
84 SaddleContainer( ) : ONE(Teuchos::ScalarTraits<ScalarType>::one()), ZERO(Teuchos::ScalarTraits<ScalarType>::zero()) { };
85 SaddleContainer(
const RCP<MV> X,
bool eye=
false );
89 SaddleContainer<ScalarType, MV> * Clone(
const int nvecs)
const;
91 SaddleContainer<ScalarType, MV> * CloneCopy()
const;
93 SaddleContainer<ScalarType, MV> * CloneCopy(
const std::vector<int> &index)
const;
94 SaddleContainer<ScalarType, MV> * CloneViewNonConst(
const std::vector<int> &index)
const;
95 SaddleContainer<ScalarType, MV> * CloneViewNonConst(
const Teuchos::Range1D& index)
const;
96 const SaddleContainer<ScalarType, MV> * CloneView(
const std::vector<int> &index)
const;
97 const SaddleContainer<ScalarType, MV> * CloneView(
const Teuchos::Range1D& index)
const;
98 ptrdiff_t GetGlobalLength()
const {
return MVT::GetGlobalLength(*upper_) + lowerRaw_->numRows(); };
101 void MvTimesMatAddMv(ScalarType alpha,
const SaddleContainer<ScalarType,MV> &A,
102 const Teuchos::SerialDenseMatrix<int, ScalarType> &B,
105 void MvAddMv(ScalarType alpha,
const SaddleContainer<ScalarType,MV>& A,
106 ScalarType beta,
const SaddleContainer<ScalarType,MV>& B);
108 void MvScale( ScalarType alpha );
110 void MvScale(
const std::vector<ScalarType>& alpha );
112 void MvTransMv (ScalarType alpha,
const SaddleContainer<ScalarType, MV>& A,
113 Teuchos::SerialDenseMatrix< int, ScalarType >& B)
const;
115 void MvDot (
const SaddleContainer<ScalarType, MV>& A, std::vector<ScalarType> &b)
const;
117 void MvNorm ( std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &normvec)
const;
121 void SetBlock (
const SaddleContainer<ScalarType, MV>& A,
const std::vector<int> &index);
123 void Assign (
const SaddleContainer<ScalarType, MV>&A);
127 void MvInit (ScalarType alpha);
129 void MvPrint (std::ostream &os)
const;
135template <
class ScalarType,
class MV>
136RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > SaddleContainer<ScalarType, MV>::getLower()
const
143 int nrows = lowerRaw_->numRows();
144 int ncols = indices_.size();
146 RCP<SerialDenseMatrix> lower = rcp(
new SerialDenseMatrix(nrows,ncols,
false));
148 for(
int r=0; r<nrows; r++)
150 for(
int c=0; c<ncols; c++)
152 (*lower)(r,c) = (*lowerRaw_)(r,indices_[c]);
162template <
class ScalarType,
class MV>
163void SaddleContainer<ScalarType, MV>::setLower(
const RCP<SerialDenseMatrix> lower)
171 int nrows = lowerRaw_->numRows();
172 int ncols = indices_.size();
174 for(
int r=0; r<nrows; r++)
176 for(
int c=0; c<ncols; c++)
178 (*lowerRaw_)(r,indices_[c]) = (*lower)(r,c);
186template <
class ScalarType,
class MV>
187SaddleContainer<ScalarType, MV>::SaddleContainer(
const RCP<MV> X,
bool eye )
188: ONE(Teuchos::ScalarTraits<ScalarType>::one()), ZERO(Teuchos::ScalarTraits<ScalarType>::zero())
190 int nvecs = MVT::GetNumberVecs(*X);
195 upper_ = MVT::Clone(*X, nvecs);
196 MVT::MvInit(*upper_);
199 lowerRaw_ = rcp(
new SerialDenseMatrix(nvecs,nvecs));
200 for(
int i=0; i < nvecs; i++)
201 (*lowerRaw_)(i,i) = ONE;
209 lowerRaw_ = rcp(
new SerialDenseMatrix(nvecs,nvecs));
216template <
class ScalarType,
class MV>
217SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::Clone(
const int nvecs)
const
219 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
221 newSC->upper_ = MVT::Clone(*upper_,nvecs);
222 newSC->lowerRaw_ = rcp(
new SerialDenseMatrix(lowerRaw_->numRows(),nvecs));
230template <
class ScalarType,
class MV>
231SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneCopy()
const
233 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
235 newSC->upper_ = MVT::CloneCopy(*upper_);
236 newSC->lowerRaw_ = getLower();
244template <
class ScalarType,
class MV>
245SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneCopy(
const std::vector< int > &index)
const
247 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
249 newSC->upper_ = MVT::CloneCopy(*upper_,index);
251 int ncols = index.size();
252 int nrows = lowerRaw_->numRows();
253 RCP<SerialDenseMatrix> lower = getLower();
254 newSC->lowerRaw_ = rcp(
new SerialDenseMatrix(nrows,ncols));
255 for(
int c=0; c < ncols; c++)
257 for(
int r=0; r < nrows; r++)
258 (*newSC->lowerRaw_)(r,c) = (*lower)(r,index[c]);
267template <
class ScalarType,
class MV>
268SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneViewNonConst(
const std::vector<int> &index)
const
270 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
272 newSC->upper_ = MVT::CloneViewNonConst(*upper_,index);
274 newSC->lowerRaw_ = lowerRaw_;
276 if(!indices_.empty())
278 newSC->indices_.resize(index.size());
279 for(
unsigned int i=0; i<index.size(); i++)
281 newSC->indices_[i] = indices_[index[i]];
286 newSC->indices_ = index;
293template <
class ScalarType,
class MV>
294SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneViewNonConst(
const Teuchos::Range1D& index)
const
296 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
298 newSC->upper_ = MVT::CloneViewNonConst(*upper_,index);
300 newSC->lowerRaw_ = lowerRaw_;
302 newSC->indices_.resize(index.size());
303 for(
unsigned int i=0; i<index.size(); i++)
305 newSC->indices_[i] = indices_[index.lbound()+i];
314template <
class ScalarType,
class MV>
315const SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneView(
const std::vector<int> &index)
const
317 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
319 newSC->upper_ = MVT::CloneViewNonConst(*upper_,index);
321 newSC->lowerRaw_ = lowerRaw_;
323 if(!indices_.empty())
325 newSC->indices_.resize(index.size());
326 for(
unsigned int i=0; i<index.size(); i++)
328 newSC->indices_[i] = indices_[index[i]];
333 newSC->indices_ = index;
340template <
class ScalarType,
class MV>
341const SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneView(
const Teuchos::Range1D& index)
const
343 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
345 newSC->upper_ = MVT::CloneViewNonConst(*upper_,index);
347 newSC->lowerRaw_ = lowerRaw_;
349 newSC->indices_.resize(index.size());
350 for(
unsigned int i=0; i<index.size(); i++)
352 newSC->indices_[i] = indices_[index.lbound()+i];
362template <
class ScalarType,
class MV>
363void SaddleContainer<ScalarType, MV>::MvTimesMatAddMv(ScalarType alpha,
const SaddleContainer<ScalarType,MV> &A,
364 const Teuchos::SerialDenseMatrix<int, ScalarType> &B,
367 MVT::MvTimesMatAddMv(alpha,*(A.upper_),B,beta,*upper_);
368 RCP<SerialDenseMatrix> lower = getLower();
369 RCP<SerialDenseMatrix> Alower = A.getLower();
370 lower->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,alpha,*Alower,B,beta);
377template <
class ScalarType,
class MV>
378void SaddleContainer<ScalarType, MV>::MvAddMv(ScalarType alpha,
const SaddleContainer<ScalarType,MV>& A,
379 ScalarType beta,
const SaddleContainer<ScalarType,MV>& B)
381 MVT::MvAddMv(alpha, *(A.upper_), beta, *(B.upper_), *upper_);
383 RCP<SerialDenseMatrix> lower = getLower();
384 RCP<SerialDenseMatrix> Alower = A.getLower();
385 RCP<SerialDenseMatrix> Blower = B.getLower();
392 lower->assign(*Alower);
398 else if(beta == -ONE)
400 else if(beta != ZERO)
402 SerialDenseMatrix scaledB(*Blower);
413template <
class ScalarType,
class MV>
414void SaddleContainer<ScalarType, MV>::MvScale( ScalarType alpha )
416 MVT::MvScale(*upper_, alpha);
419 RCP<SerialDenseMatrix> lower = getLower();
427template <
class ScalarType,
class MV>
428void SaddleContainer<ScalarType, MV>::MvScale(
const std::vector<ScalarType>& alpha )
430 MVT::MvScale(*upper_, alpha);
432 RCP<SerialDenseMatrix> lower = getLower();
434 int nrows = lower->numRows();
435 int ncols = lower->numCols();
437 for(
int c=0; c<ncols; c++)
439 for(
int r=0; r<nrows; r++)
440 (*lower)(r,c) *= alpha[c];
450template <
class ScalarType,
class MV>
451void SaddleContainer<ScalarType, MV>::MvTransMv (ScalarType alpha,
const SaddleContainer<ScalarType, MV>& A,
452 Teuchos::SerialDenseMatrix< int, ScalarType >& B)
const
454 MVT::MvTransMv(alpha,*(A.upper_),*upper_,B);
455 RCP<SerialDenseMatrix> lower = getLower();
456 RCP<SerialDenseMatrix> Alower = A.getLower();
457 B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,alpha,*(Alower),*lower,ONE);
463template <
class ScalarType,
class MV>
464void SaddleContainer<ScalarType, MV>::MvDot (
const SaddleContainer<ScalarType, MV>& A, std::vector<ScalarType> &b)
const
466 MVT::MvDot(*upper_, *(A.upper_), b);
468 RCP<SerialDenseMatrix> lower = getLower();
469 RCP<SerialDenseMatrix> Alower = A.getLower();
471 int nrows = lower->numRows();
472 int ncols = lower->numCols();
474 for(
int c=0; c < ncols; c++)
476 for(
int r=0; r < nrows; r++)
478 b[c] += ((*Alower)(r,c) * (*lower)(r,c));
487template <
class ScalarType,
class MV>
488void SaddleContainer<ScalarType, MV>::MvNorm ( std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &normvec)
const
491 MvDot(*
this,normvec);
492 for(
unsigned int i=0; i<normvec.size(); i++)
493 normvec[i] = sqrt(normvec[i]);
501template <
class ScalarType,
class MV>
502void SaddleContainer<ScalarType, MV>::SetBlock (
const SaddleContainer<ScalarType, MV>& A,
const std::vector<int> &index)
504 MVT::SetBlock(*(A.upper_), index, *upper_);
506 RCP<SerialDenseMatrix> lower = getLower();
507 RCP<SerialDenseMatrix> Alower = A.getLower();
509 int nrows = lower->numRows();
511 int nvecs = index.size();
512 for(
int c=0; c<nvecs; c++)
514 for(
int r=0; r<nrows; r++)
515 (*lower)(r,index[c]) = (*Alower)(r,c);
524template <
class ScalarType,
class MV>
525void SaddleContainer<ScalarType, MV>::Assign (
const SaddleContainer<ScalarType, MV>&A)
527 MVT::Assign(*(A.upper_),*(upper_));
529 RCP<SerialDenseMatrix> lower = getLower();
530 RCP<SerialDenseMatrix> Alower = A.getLower();
541template <
class ScalarType,
class MV>
542void SaddleContainer<ScalarType, MV>::MvRandom ()
544 MVT::MvRandom(*upper_);
546 RCP<SerialDenseMatrix> lower = getLower();
554template <
class ScalarType,
class MV>
555void SaddleContainer<ScalarType, MV>::MvInit (ScalarType alpha)
557 MVT::MvInit(*upper_,alpha);
559 RCP<SerialDenseMatrix> lower = getLower();
560 lower->putScalar(alpha);
567template <
class ScalarType,
class MV>
568void SaddleContainer<ScalarType, MV>::MvPrint (std::ostream &os)
const
570 RCP<SerialDenseMatrix> lower = getLower();
574 os <<
"Object SaddleContainer" << std::endl;
576 upper_->describe(*(Teuchos::VerboseObjectBase::getDefaultOStream()),Teuchos::VERB_EXTREME);
579 os <<
"Y\n" << *lower << std::endl;
584template<
class ScalarType,
class MV >
585class MultiVecTraits<ScalarType,
Experimental::SaddleContainer<ScalarType, MV> >
587typedef Experimental::SaddleContainer<ScalarType,MV> SC;
590 static RCP<SC > Clone(
const SC& mv,
const int numvecs )
591 {
return rcp(
const_cast<SC&
>(mv).Clone(numvecs) ); }
593 static RCP<SC > CloneCopy(
const SC& mv )
594 {
return rcp(
const_cast<SC&
>(mv).CloneCopy() ); }
596 static RCP<SC > CloneCopy(
const SC& mv,
const std::vector<int>& index )
597 {
return rcp(
const_cast<SC&
>(mv).CloneCopy(index) ); }
599 static ptrdiff_t GetGlobalLength(
const SC& mv )
600 {
return mv.GetGlobalLength(); }
602 static int GetNumberVecs(
const SC& mv )
603 {
return mv.GetNumberVecs(); }
605 static void MvTimesMatAddMv( ScalarType alpha,
const SC& A,
606 const Teuchos::SerialDenseMatrix<int,ScalarType>& B,
607 ScalarType beta, SC& mv )
608 { mv.MvTimesMatAddMv(alpha, A, B, beta); }
610 static void MvAddMv( ScalarType alpha,
const SC& A, ScalarType beta,
const SC& B, SC& mv )
611 { mv.MvAddMv(alpha, A, beta, B); }
613 static void MvTransMv( ScalarType alpha,
const SC& A,
const SC& mv, Teuchos::SerialDenseMatrix<int,ScalarType>& B)
614 { mv.MvTransMv(alpha, A, B); }
616 static void MvDot(
const SC& mv,
const SC& A, std::vector<ScalarType> & b)
619 static void MvScale ( SC& mv, ScalarType alpha )
620 { mv.MvScale( alpha ); }
622 static void MvScale ( SC& mv,
const std::vector<ScalarType>& alpha )
623 { mv.MvScale( alpha ); }
625 static void MvNorm(
const SC& mv, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> & normvec)
626 { mv.MvNorm(normvec); }
628 static void SetBlock(
const SC& A,
const std::vector<int>& index, SC& mv )
629 { mv.SetBlock(A, index); }
631 static void Assign(
const SC& A, SC& mv )
634 static void MvRandom( SC& mv )
637 static void MvInit( SC& mv, ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero() )
638 { mv.MvInit(alpha); }
640 static void MvPrint(
const SC& mv, std::ostream& os )
646#ifdef HAVE_ANASAZI_BELOS
650template<
class ScalarType,
class MV >
651class MultiVecTraits< ScalarType,
Anasazi::Experimental::SaddleContainer<ScalarType,MV> >
653typedef Anasazi::Experimental::SaddleContainer<ScalarType,MV> SC;
655 static RCP<SC > Clone(
const SC& mv,
const int numvecs )
656 {
return rcp(
const_cast<SC&
>(mv).Clone(numvecs) ); }
658 static RCP<SC > CloneCopy(
const SC& mv )
659 {
return rcp(
const_cast<SC&
>(mv).CloneCopy() ); }
661 static RCP<SC > CloneCopy(
const SC& mv,
const std::vector<int>& index )
662 {
return rcp(
const_cast<SC&
>(mv).CloneCopy(index) ); }
664 static RCP<SC> CloneViewNonConst( SC& mv,
const std::vector<int>& index )
665 {
return rcp( mv.CloneViewNonConst(index) ); }
667 static RCP<SC> CloneViewNonConst( SC& mv,
const Teuchos::Range1D& index )
668 {
return rcp( mv.CloneViewNonConst(index) ); }
670 static RCP<const SC> CloneView(
const SC& mv,
const std::vector<int> & index )
671 {
return rcp( mv.CloneView(index) ); }
673 static RCP<const SC> CloneView(
const SC& mv,
const Teuchos::Range1D& index )
674 {
return rcp( mv.CloneView(index) ); }
676 static ptrdiff_t GetGlobalLength(
const SC& mv )
677 {
return mv.GetGlobalLength(); }
679 static int GetNumberVecs(
const SC& mv )
680 {
return mv.GetNumberVecs(); }
682 static void MvTimesMatAddMv( ScalarType alpha,
const SC& A,
683 const Teuchos::SerialDenseMatrix<int,ScalarType>& B,
684 ScalarType beta, SC& mv )
685 { mv.MvTimesMatAddMv(alpha, A, B, beta); }
687 static void MvAddMv( ScalarType alpha,
const SC& A, ScalarType beta,
const SC& B, SC& mv )
688 { mv.MvAddMv(alpha, A, beta, B); }
690 static void MvTransMv( ScalarType alpha,
const SC& A,
const SC& mv, Teuchos::SerialDenseMatrix<int,ScalarType>& B)
691 { mv.MvTransMv(alpha, A, B); }
693 static void MvDot(
const SC& mv,
const SC& A, std::vector<ScalarType> & b)
696 static void MvScale ( SC& mv, ScalarType alpha )
697 { mv.MvScale( alpha ); }
699 static void MvScale ( SC& mv,
const std::vector<ScalarType>& alpha )
700 { mv.MvScale( alpha ); }
703 static void MvNorm(
const SC& mv, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> & normvec, NormType type=TwoNorm)
704 { mv.MvNorm(normvec); }
706 static void SetBlock(
const SC& A,
const std::vector<int>& index, SC& mv )
707 { mv.SetBlock(A, index); }
709 static void Assign(
const SC& A, SC& mv )
712 static void MvRandom( SC& mv )
715 static void MvInit( SC& mv, ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero() )
716 { mv.MvInit(alpha); }
718 static void MvPrint(
const SC& mv, std::ostream& os )
721#ifdef HAVE_BELOS_TSQR
734 typedef Belos::details::StubTsqrAdapter<SC> tsqr_adaptor_type;
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Traits class which defines basic operations on multivectors.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.