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" 66 template <
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;
135 template <
class ScalarType,
class MV>
136 RCP< 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]);
162 template <
class ScalarType,
class MV>
163 void 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);
186 template <
class ScalarType,
class MV>
187 SaddleContainer<ScalarType, MV>::SaddleContainer(
const RCP<MV> X,
bool eye )
188 : ONE(Teuchos::ScalarTraits<ScalarType>::one()), ZERO(Teuchos::ScalarTraits<ScalarType>::zero())
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));
216 template <
class ScalarType,
class MV>
217 SaddleContainer<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));
230 template <
class ScalarType,
class MV>
231 SaddleContainer<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();
244 template <
class ScalarType,
class MV>
245 SaddleContainer<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]);
267 template <
class ScalarType,
class MV>
268 SaddleContainer<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;
293 template <
class ScalarType,
class MV>
294 SaddleContainer<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];
314 template <
class ScalarType,
class MV>
315 const 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;
340 template <
class ScalarType,
class MV>
341 const 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];
362 template <
class ScalarType,
class MV>
363 void 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);
377 template <
class ScalarType,
class MV>
378 void 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);
413 template <
class ScalarType,
class MV>
414 void SaddleContainer<ScalarType, MV>::MvScale( ScalarType alpha )
416 MVT::MvScale(*upper_, alpha);
419 RCP<SerialDenseMatrix> lower = getLower();
427 template <
class ScalarType,
class MV>
428 void 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];
450 template <
class ScalarType,
class MV>
451 void 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);
463 template <
class ScalarType,
class MV>
464 void 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));
487 template <
class ScalarType,
class MV>
488 void 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]);
501 template <
class ScalarType,
class MV>
502 void 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);
524 template <
class ScalarType,
class MV>
525 void 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();
541 template <
class ScalarType,
class MV>
542 void SaddleContainer<ScalarType, MV>::MvRandom ()
544 MVT::MvRandom(*upper_);
546 RCP<SerialDenseMatrix> lower = getLower();
554 template <
class ScalarType,
class MV>
555 void SaddleContainer<ScalarType, MV>::MvInit (ScalarType alpha)
557 MVT::MvInit(*upper_,alpha);
559 RCP<SerialDenseMatrix> lower = getLower();
560 lower->putScalar(alpha);
567 template <
class ScalarType,
class MV>
568 void 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;
584 template<
class ScalarType,
class MV >
585 class MultiVecTraits<ScalarType,
Experimental::SaddleContainer<ScalarType, MV> >
587 typedef 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 650 template<
class ScalarType,
class MV >
651 class MultiVecTraits< ScalarType,
Anasazi::Experimental::SaddleContainer<ScalarType,MV> >
653 typedef 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 722 typedef Belos::details::StubTsqrAdapter<SC> tsqr_adaptor_type;
735 #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.