00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00034 #ifndef ANASAZI_THYRA_ADAPTER_HPP
00035 #define ANASAZI_THYRA_ADAPTER_HPP
00036
00037 #include "AnasaziMultiVecTraits.hpp"
00038 #include "AnasaziOperatorTraits.hpp"
00039 #include "AnasaziConfigDefs.hpp"
00040
00041 #include <Thyra_DetachedMultiVectorView.hpp>
00042 #include <Thyra_MultiVectorBaseDecl.hpp>
00043 #include <Thyra_MultiVectorStdOps.hpp>
00044
00045 namespace Anasazi {
00046
00048
00049
00050
00052
00061 template<class ScalarType>
00062 class MultiVecTraits< ScalarType, Thyra::MultiVectorBase<ScalarType> >
00063 {
00064 public:
00065
00068
00073 static Teuchos::RCP< Thyra::MultiVectorBase<ScalarType> > Clone( const Thyra::MultiVectorBase<ScalarType> & mv, const int numvecs )
00074 {
00075 Teuchos::RCP< Thyra::MultiVectorBase<ScalarType> > c = Thyra::createMembers( mv.range(), numvecs );
00076 return c;
00077 }
00078
00083 static Teuchos::RCP< Thyra::MultiVectorBase<ScalarType> > CloneCopy( const Thyra::MultiVectorBase< ScalarType > & mv )
00084 {
00085 int numvecs = mv.domain()->dim();
00086
00087 Teuchos::RCP< Thyra::MultiVectorBase< ScalarType > > cc = Thyra::createMembers( mv.range(), numvecs );
00088
00089 Thyra::assign(&*cc, mv);
00090 return cc;
00091 }
00092
00098 static Teuchos::RCP< Thyra::MultiVectorBase< ScalarType > > CloneCopy( const Thyra::MultiVectorBase< ScalarType > & mv, const std::vector<int>& index )
00099 {
00100 int numvecs = index.size();
00101
00102 Teuchos::RCP< Thyra::MultiVectorBase< ScalarType > > cc = Thyra::createMembers( mv.range(), numvecs );
00103
00104 Teuchos::RCP< const Thyra::MultiVectorBase< ScalarType > > view = mv.subView( numvecs, &(index[0]) );
00105
00106 Thyra::assign(&*cc, *view);
00107 return cc;
00108 }
00109
00115 static Teuchos::RCP< Thyra::MultiVectorBase< ScalarType > > CloneView( Thyra::MultiVectorBase< ScalarType > & mv, const std::vector<int>& index )
00116 {
00117 int numvecs = index.size();
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131 int lb = index[0];
00132 bool contig = true;
00133 for (int i=0; i<numvecs; i++) {
00134 if (lb+i != index[i]) contig = false;
00135 }
00136
00137 Teuchos::RCP< Thyra::MultiVectorBase< ScalarType > > cc;
00138 if (contig) {
00139 const Thyra::Range1D rng(lb,lb+numvecs-1);
00140
00141 cc = mv.subView(rng);
00142 }
00143 else {
00144
00145 cc = mv.subView( numvecs, &(index[0]) );
00146 }
00147 return cc;
00148 }
00149
00155 static Teuchos::RCP<const Thyra::MultiVectorBase< ScalarType > > CloneView( const Thyra::MultiVectorBase< ScalarType > & mv, const std::vector<int>& index )
00156 {
00157 int numvecs = index.size();
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171 int lb = index[0];
00172 bool contig = true;
00173 for (int i=0; i<numvecs; i++) {
00174 if (lb+i != index[i]) contig = false;
00175 }
00176
00177 Teuchos::RCP< const Thyra::MultiVectorBase< ScalarType > > cc;
00178 if (contig) {
00179 const Thyra::Range1D rng(lb,lb+numvecs-1);
00180
00181 cc = mv.subView(rng);
00182 }
00183 else {
00184
00185 cc = mv.subView( numvecs, &(index[0]) );
00186 }
00187 return cc;
00188 }
00189
00191
00194
00196 static int GetVecLength( const Thyra::MultiVectorBase< ScalarType > & mv )
00197 { return mv.range()->dim(); }
00198
00200 static int GetNumberVecs( const Thyra::MultiVectorBase< ScalarType > & mv )
00201 { return mv.domain()->dim(); }
00202
00204
00207
00210 static void MvTimesMatAddMv( const ScalarType alpha, const Thyra::MultiVectorBase< ScalarType > & A,
00211 const Teuchos::SerialDenseMatrix<int,ScalarType>& B,
00212 const ScalarType beta, Thyra::MultiVectorBase< ScalarType > & mv )
00213 {
00214 int m = B.numRows();
00215 int n = B.numCols();
00216
00217 Teuchos::RCP< const Thyra::MultiVectorBase< ScalarType > >
00218 B_thyra = Thyra::createMembersView(
00219 A.domain()
00220 ,RTOpPack::ConstSubMultiVectorView<ScalarType>(0,m,0,n,&B(0,0),B.stride())
00221 );
00222
00223 A.apply(Thyra::NONCONJ_ELE,*B_thyra,&mv,alpha,beta);
00224 }
00225
00228 static void MvAddMv( const ScalarType alpha, const Thyra::MultiVectorBase< ScalarType > & A,
00229 const ScalarType beta, const Thyra::MultiVectorBase< ScalarType > & B, Thyra::MultiVectorBase< ScalarType > & mv )
00230 {
00231 ScalarType coef[2], zero = Teuchos::ScalarTraits<ScalarType>::zero();
00232 const Thyra::MultiVectorBase< ScalarType > * in[2];
00233
00234 in[0] = &A;
00235 in[1] = &B;
00236 coef[0] = alpha;
00237 coef[1] = beta;
00238
00239 Thyra::linear_combination(2,coef,in,zero,&mv);
00240 }
00241
00244 static void MvTransMv( const ScalarType alpha, const Thyra::MultiVectorBase< ScalarType > & A, const Thyra::MultiVectorBase< ScalarType > & mv,
00245 Teuchos::SerialDenseMatrix<int,ScalarType>& B )
00246 {
00247
00248 int m = A.domain()->dim();
00249 int n = mv.domain()->dim();
00250
00251 Teuchos::RCP< Thyra::MultiVectorBase< ScalarType > >
00252 B_thyra = Thyra::createMembersView(
00253 A.domain()
00254 ,RTOpPack::SubMultiVectorView<ScalarType>(0,m,0,n,&B(0,0),B.stride())
00255 );
00256 A.applyTranspose(Thyra::CONJ_ELE,mv,&*B_thyra,alpha,Teuchos::ScalarTraits<ScalarType>::zero());
00257 }
00258
00262 static void MvDot( const Thyra::MultiVectorBase< ScalarType > & mv, const Thyra::MultiVectorBase< ScalarType > & A, std::vector<ScalarType> &b )
00263 { Thyra::dots(mv,A,&(b[0])); }
00264
00267 static void MvScale ( Thyra::MultiVectorBase< ScalarType > & mv, const ScalarType alpha )
00268 { Thyra::scale(alpha,&mv); }
00269
00272 static void MvScale ( Thyra::MultiVectorBase< ScalarType > & mv, const std::vector<ScalarType>& alpha )
00273 {
00274 for (unsigned int i=0; i<alpha.size(); i++) {
00275 Thyra::scale(alpha[i],mv.col(i).get());
00276 }
00277 }
00278
00280
00283
00287 static void MvNorm( const Thyra::MultiVectorBase< ScalarType > & mv, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &normvec )
00288 { Thyra::norms_2(mv,&(normvec[0])); }
00289
00291
00294
00297 static void SetBlock( const Thyra::MultiVectorBase< ScalarType > & A, const std::vector<int>& index, Thyra::MultiVectorBase< ScalarType > & mv )
00298 {
00299
00300 int numvecs = index.size();
00301 std::vector<int> indexA(numvecs);
00302 int numAcols = A.domain()->dim();
00303 for (int i=0; i<numvecs; i++) {
00304 indexA[i] = i;
00305 }
00306
00307
00308 if ( numAcols < numvecs ) {
00309
00310
00311 numvecs = numAcols;
00312 }
00313 else if ( numAcols > numvecs ) {
00314 numAcols = numvecs;
00315 indexA.resize( numAcols );
00316 }
00317
00318 Teuchos::RCP< const Thyra::MultiVectorBase< ScalarType > > relsource = A.subView( numAcols, &(indexA[0]) );
00319
00320 Teuchos::RCP< Thyra::MultiVectorBase< ScalarType > > reldest = mv.subView( numvecs, &(index[0]) );
00321
00322 Thyra::assign(&*reldest, *relsource);
00323 }
00324
00327 static void MvRandom( Thyra::MultiVectorBase< ScalarType > & mv )
00328 {
00329
00330
00331 Thyra::randomize(-Teuchos::ScalarTraits<ScalarType>::one(),
00332 Teuchos::ScalarTraits<ScalarType>::one(),
00333 &mv);
00334 }
00335
00338 static void MvInit( Thyra::MultiVectorBase< ScalarType > & mv, ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero() )
00339 { Thyra::assign(&mv,alpha); }
00340
00342
00345
00348 static void MvPrint( const Thyra::MultiVectorBase< ScalarType > & mv, std::ostream& os )
00349 {
00350 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream(Teuchos::rcp(&os,false));
00351 out->setf(std::ios_base::scientific);
00352 out->precision(16);
00353 mv.describe(*out,Teuchos::VERB_EXTREME);
00354 }
00355
00357
00358 };
00359
00361
00362
00363
00365
00375 template <class ScalarType>
00376 class OperatorTraits < ScalarType, Thyra::MultiVectorBase<ScalarType>, Thyra::LinearOpBase<ScalarType> >
00377 {
00378 public:
00379
00383 static void Apply ( const Thyra::LinearOpBase< ScalarType >& Op, const Thyra::MultiVectorBase< ScalarType > & x, Thyra::MultiVectorBase< ScalarType > & y )
00384 {
00385 Op.apply(Thyra::NONCONJ_ELE,x,&y);
00386 }
00387
00388 };
00389
00390 }
00391
00392 #endif
00393