36 #ifndef ANASAZI_SADDLE_CONTAINER_HPP 37 #define ANASAZI_SADDLE_CONTAINER_HPP 40 #include "Teuchos_VerboseObject.hpp" 42 #ifdef HAVE_ANASAZI_BELOS 43 #include "BelosMultiVecTraits.hpp" 52 template <
class ScalarType,
class MV>
55 template <
class Scalar_,
class MV_,
class OP_>
friend class SaddleOperator;
59 typedef Teuchos::SerialDenseMatrix<int,ScalarType> SerialDenseMatrix;
60 const ScalarType ONE, ZERO;
62 RCP<SerialDenseMatrix> lowerRaw_;
63 std::vector<int> indices_;
65 RCP<SerialDenseMatrix> getLower()
const;
66 void setLower(
const RCP<SerialDenseMatrix> lower);
70 SaddleContainer( ) : ONE(Teuchos::ScalarTraits<ScalarType>::one()), ZERO(Teuchos::ScalarTraits<ScalarType>::zero()) { };
71 SaddleContainer(
const RCP<MV> X,
bool eye=
false );
75 SaddleContainer<ScalarType, MV> * Clone(
const int nvecs)
const;
77 SaddleContainer<ScalarType, MV> * CloneCopy()
const;
79 SaddleContainer<ScalarType, MV> * CloneCopy(
const std::vector<int> &index)
const;
80 SaddleContainer<ScalarType, MV> * CloneViewNonConst(
const std::vector<int> &index)
const;
81 const SaddleContainer<ScalarType, MV> * CloneView(
const std::vector<int> &index)
const;
82 ptrdiff_t GetGlobalLength()
const {
return MVT::GetGlobalLength(*upper_) + lowerRaw_->numRows(); };
85 void MvTimesMatAddMv(ScalarType alpha,
const SaddleContainer<ScalarType,MV> &A,
86 const Teuchos::SerialDenseMatrix<int, ScalarType> &B,
89 void MvAddMv(ScalarType alpha,
const SaddleContainer<ScalarType,MV>& A,
90 ScalarType beta,
const SaddleContainer<ScalarType,MV>& B);
92 void MvScale( ScalarType alpha );
94 void MvScale(
const std::vector<ScalarType>& alpha );
96 void MvTransMv (ScalarType alpha,
const SaddleContainer<ScalarType, MV>& A,
97 Teuchos::SerialDenseMatrix< int, ScalarType >& B)
const;
99 void MvDot (
const SaddleContainer<ScalarType, MV>& A, std::vector<ScalarType> &b)
const;
101 void MvNorm ( std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &normvec)
const;
105 void SetBlock (
const SaddleContainer<ScalarType, MV>& A,
const std::vector<int> &index);
107 void Assign (
const SaddleContainer<ScalarType, MV>&A);
111 void MvInit (ScalarType alpha);
113 void MvPrint (std::ostream &os)
const;
119 template <
class ScalarType,
class MV>
120 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > SaddleContainer<ScalarType, MV>::getLower()
const 127 int nrows = lowerRaw_->numRows();
128 int ncols = indices_.size();
130 RCP<SerialDenseMatrix> lower = rcp(
new SerialDenseMatrix(nrows,ncols,
false));
132 for(
int r=0; r<nrows; r++)
134 for(
int c=0; c<ncols; c++)
136 (*lower)(r,c) = (*lowerRaw_)(r,indices_[c]);
146 template <
class ScalarType,
class MV>
147 void SaddleContainer<ScalarType, MV>::setLower(
const RCP<SerialDenseMatrix> lower)
155 int nrows = lowerRaw_->numRows();
156 int ncols = indices_.size();
158 for(
int r=0; r<nrows; r++)
160 for(
int c=0; c<ncols; c++)
162 (*lowerRaw_)(r,indices_[c]) = (*lower)(r,c);
170 template <
class ScalarType,
class MV>
171 SaddleContainer<ScalarType, MV>::SaddleContainer(
const RCP<MV> X,
bool eye )
172 : ONE(Teuchos::ScalarTraits<ScalarType>::one()), ZERO(Teuchos::ScalarTraits<ScalarType>::zero())
183 lowerRaw_ = rcp(
new SerialDenseMatrix(nvecs,nvecs));
184 for(
int i=0; i < nvecs; i++)
185 (*lowerRaw_)(i,i) = ONE;
193 lowerRaw_ = rcp(
new SerialDenseMatrix(nvecs,nvecs));
200 template <
class ScalarType,
class MV>
201 SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::Clone(
const int nvecs)
const 203 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
206 newSC->lowerRaw_ = rcp(
new SerialDenseMatrix(lowerRaw_->numRows(),nvecs));
214 template <
class ScalarType,
class MV>
215 SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneCopy()
const 217 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
220 newSC->lowerRaw_ = getLower();
228 template <
class ScalarType,
class MV>
229 SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneCopy(
const std::vector< int > &index)
const 231 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
235 int ncols = index.size();
236 int nrows = lowerRaw_->numRows();
237 RCP<SerialDenseMatrix> lower = getLower();
238 newSC->lowerRaw_ = rcp(
new SerialDenseMatrix(nrows,ncols));
239 for(
int c=0; c < ncols; c++)
241 for(
int r=0; r < nrows; r++)
242 (*newSC->lowerRaw_)(r,c) = (*lower)(r,index[c]);
251 template <
class ScalarType,
class MV>
252 SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneViewNonConst(
const std::vector<int> &index)
const 254 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
258 newSC->lowerRaw_ = lowerRaw_;
260 if(!indices_.empty())
262 newSC->indices_.resize(index.size());
263 for(
unsigned int i=0; i<index.size(); i++)
265 newSC->indices_[i] = indices_[index[i]];
270 newSC->indices_ = index;
279 template <
class ScalarType,
class MV>
280 const SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneView(
const std::vector<int> &index)
const 282 SaddleContainer<ScalarType, MV> * newSC =
new SaddleContainer<ScalarType, MV>();
286 newSC->lowerRaw_ = lowerRaw_;
288 if(!indices_.empty())
290 newSC->indices_.resize(index.size());
291 for(
unsigned int i=0; i<index.size(); i++)
293 newSC->indices_[i] = indices_[index[i]];
298 newSC->indices_ = index;
308 template <
class ScalarType,
class MV>
309 void SaddleContainer<ScalarType, MV>::MvTimesMatAddMv(ScalarType alpha,
const SaddleContainer<ScalarType,MV> &A,
310 const Teuchos::SerialDenseMatrix<int, ScalarType> &B,
314 RCP<SerialDenseMatrix> lower = getLower();
315 RCP<SerialDenseMatrix> Alower = A.getLower();
316 lower->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,alpha,*Alower,B,beta);
323 template <
class ScalarType,
class MV>
324 void SaddleContainer<ScalarType, MV>::MvAddMv(ScalarType alpha,
const SaddleContainer<ScalarType,MV>& A,
325 ScalarType beta,
const SaddleContainer<ScalarType,MV>& B)
327 MVT::MvAddMv(alpha, *(A.upper_), beta, *(B.upper_), *upper_);
329 RCP<SerialDenseMatrix> lower = getLower();
330 RCP<SerialDenseMatrix> Alower = A.getLower();
331 RCP<SerialDenseMatrix> Blower = B.getLower();
338 lower->assign(*Alower);
344 else if(beta == -ONE)
346 else if(beta != ZERO)
348 SerialDenseMatrix scaledB(*Blower);
359 template <
class ScalarType,
class MV>
360 void SaddleContainer<ScalarType, MV>::MvScale( ScalarType alpha )
365 RCP<SerialDenseMatrix> lower = getLower();
373 template <
class ScalarType,
class MV>
374 void SaddleContainer<ScalarType, MV>::MvScale(
const std::vector<ScalarType>& alpha )
378 RCP<SerialDenseMatrix> lower = getLower();
380 int nrows = lower->numRows();
381 int ncols = lower->numCols();
383 for(
int c=0; c<ncols; c++)
385 for(
int r=0; r<nrows; r++)
386 (*lower)(r,c) *= alpha[c];
396 template <
class ScalarType,
class MV>
397 void SaddleContainer<ScalarType, MV>::MvTransMv (ScalarType alpha,
const SaddleContainer<ScalarType, MV>& A,
398 Teuchos::SerialDenseMatrix< int, ScalarType >& B)
const 401 RCP<SerialDenseMatrix> lower = getLower();
402 RCP<SerialDenseMatrix> Alower = A.getLower();
403 B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,alpha,*(Alower),*lower,ONE);
409 template <
class ScalarType,
class MV>
410 void SaddleContainer<ScalarType, MV>::MvDot (
const SaddleContainer<ScalarType, MV>& A, std::vector<ScalarType> &b)
const 414 RCP<SerialDenseMatrix> lower = getLower();
415 RCP<SerialDenseMatrix> Alower = A.getLower();
417 int nrows = lower->numRows();
418 int ncols = lower->numCols();
420 for(
int c=0; c < ncols; c++)
422 for(
int r=0; r < nrows; r++)
424 b[c] += ((*Alower)(r,c) * (*lower)(r,c));
433 template <
class ScalarType,
class MV>
434 void SaddleContainer<ScalarType, MV>::MvNorm ( std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &normvec)
const 437 MvDot(*
this,normvec);
438 for(
unsigned int i=0; i<normvec.size(); i++)
439 normvec[i] = sqrt(normvec[i]);
447 template <
class ScalarType,
class MV>
448 void SaddleContainer<ScalarType, MV>::SetBlock (
const SaddleContainer<ScalarType, MV>& A,
const std::vector<int> &index)
452 RCP<SerialDenseMatrix> lower = getLower();
453 RCP<SerialDenseMatrix> Alower = A.getLower();
455 int nrows = lower->numRows();
457 int nvecs = index.size();
458 for(
int c=0; c<nvecs; c++)
460 for(
int r=0; r<nrows; r++)
461 (*lower)(r,index[c]) = (*Alower)(r,c);
470 template <
class ScalarType,
class MV>
471 void SaddleContainer<ScalarType, MV>::Assign (
const SaddleContainer<ScalarType, MV>&A)
475 RCP<SerialDenseMatrix> lower = getLower();
476 RCP<SerialDenseMatrix> Alower = A.getLower();
487 template <
class ScalarType,
class MV>
488 void SaddleContainer<ScalarType, MV>::MvRandom ()
492 RCP<SerialDenseMatrix> lower = getLower();
500 template <
class ScalarType,
class MV>
501 void SaddleContainer<ScalarType, MV>::MvInit (ScalarType alpha)
505 RCP<SerialDenseMatrix> lower = getLower();
506 lower->putScalar(alpha);
513 template <
class ScalarType,
class MV>
514 void SaddleContainer<ScalarType, MV>::MvPrint (std::ostream &os)
const 516 RCP<SerialDenseMatrix> lower = getLower();
520 os <<
"Object SaddleContainer" << std::endl;
522 upper_->describe(*(Teuchos::VerboseObjectBase::getDefaultOStream()),Teuchos::VERB_EXTREME);
525 os <<
"Y\n" << *lower << std::endl;
530 template<
class ScalarType,
class MV >
531 class MultiVecTraits<ScalarType,
Experimental::SaddleContainer<ScalarType, MV> >
533 typedef Experimental::SaddleContainer<ScalarType,MV> SC;
536 static RCP<SC > Clone(
const SC& mv,
const int numvecs )
537 {
return rcp( const_cast<SC&>(mv).Clone(numvecs) ); }
539 static RCP<SC > CloneCopy(
const SC& mv )
540 {
return rcp( const_cast<SC&>(mv).CloneCopy() ); }
542 static RCP<SC > CloneCopy(
const SC& mv,
const std::vector<int>& index )
543 {
return rcp( const_cast<SC&>(mv).CloneCopy(index) ); }
545 static ptrdiff_t GetGlobalLength(
const SC& mv )
546 {
return mv.GetGlobalLength(); }
548 static int GetNumberVecs(
const SC& mv )
549 {
return mv.GetNumberVecs(); }
551 static void MvTimesMatAddMv( ScalarType alpha,
const SC& A,
552 const Teuchos::SerialDenseMatrix<int,ScalarType>& B,
553 ScalarType beta, SC& mv )
554 { mv.MvTimesMatAddMv(alpha, A, B, beta); }
556 static void MvAddMv( ScalarType alpha,
const SC& A, ScalarType beta,
const SC& B, SC& mv )
557 { mv.MvAddMv(alpha, A, beta, B); }
559 static void MvTransMv( ScalarType alpha,
const SC& A,
const SC& mv, Teuchos::SerialDenseMatrix<int,ScalarType>& B)
560 { mv.MvTransMv(alpha, A, B); }
562 static void MvDot(
const SC& mv,
const SC& A, std::vector<ScalarType> & b)
565 static void MvScale ( SC& mv, ScalarType alpha )
566 { mv.MvScale( alpha ); }
568 static void MvScale ( SC& mv,
const std::vector<ScalarType>& alpha )
569 { mv.MvScale( alpha ); }
571 static void MvNorm(
const SC& mv, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> & normvec)
572 { mv.MvNorm(normvec); }
574 static void SetBlock(
const SC& A,
const std::vector<int>& index, SC& mv )
575 { mv.SetBlock(A, index); }
577 static void Assign(
const SC& A, SC& mv )
580 static void MvRandom( SC& mv )
583 static void MvInit( SC& mv, ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero() )
584 { mv.MvInit(alpha); }
586 static void MvPrint(
const SC& mv, std::ostream& os )
592 #ifdef HAVE_ANASAZI_BELOS 596 template<
class ScalarType,
class MV >
597 class MultiVecTraits< ScalarType,
Anasazi::Experimental::SaddleContainer<ScalarType,MV> >
599 typedef Anasazi::Experimental::SaddleContainer<ScalarType,MV> SC;
601 static RCP<SC > Clone(
const SC& mv,
const int numvecs )
602 {
return rcp( const_cast<SC&>(mv).Clone(numvecs) ); }
604 static RCP<SC > CloneCopy(
const SC& mv )
605 {
return rcp( const_cast<SC&>(mv).CloneCopy() ); }
607 static RCP<SC > CloneCopy(
const SC& mv,
const std::vector<int>& index )
608 {
return rcp( const_cast<SC&>(mv).CloneCopy(index) ); }
610 static RCP<SC> CloneViewNonConst( SC& mv,
const std::vector<int>& index )
611 {
return rcp( mv.CloneViewNonConst(index) ); }
613 static RCP<const SC> CloneView(
const SC& mv,
const std::vector<int> & index )
614 {
return rcp( mv.CloneView(index) ); }
616 static ptrdiff_t GetGlobalLength(
const SC& mv )
617 {
return mv.GetGlobalLength(); }
619 static int GetNumberVecs(
const SC& mv )
620 {
return mv.GetNumberVecs(); }
622 static void MvTimesMatAddMv( ScalarType alpha,
const SC& A,
623 const Teuchos::SerialDenseMatrix<int,ScalarType>& B,
624 ScalarType beta, SC& mv )
625 { mv.MvTimesMatAddMv(alpha, A, B, beta); }
627 static void MvAddMv( ScalarType alpha,
const SC& A, ScalarType beta,
const SC& B, SC& mv )
628 { mv.MvAddMv(alpha, A, beta, B); }
630 static void MvTransMv( ScalarType alpha,
const SC& A,
const SC& mv, Teuchos::SerialDenseMatrix<int,ScalarType>& B)
631 { mv.MvTransMv(alpha, A, B); }
633 static void MvDot(
const SC& mv,
const SC& A, std::vector<ScalarType> & b)
636 static void MvScale ( SC& mv, ScalarType alpha )
637 { mv.MvScale( alpha ); }
639 static void MvScale ( SC& mv,
const std::vector<ScalarType>& alpha )
640 { mv.MvScale( alpha ); }
643 static void MvNorm(
const SC& mv, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> & normvec, NormType type=TwoNorm)
644 { mv.MvNorm(normvec); }
646 static void SetBlock(
const SC& A,
const std::vector<int>& index, SC& mv )
647 { mv.SetBlock(A, index); }
649 static void Assign(
const SC& A, SC& mv )
652 static void MvRandom( SC& mv )
655 static void MvInit( SC& mv, ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero() )
656 { mv.MvInit(alpha); }
658 static void MvPrint(
const SC& mv, std::ostream& os )
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.