PlayaAnasaziAdapter.hpp
00001
00002
00003
00004
00005 #ifndef ANASAZI_PLAYA_ADAPTER_HPP
00006 #define ANASAZI_PLAYA_ADAPTER_HPP
00007
00008 #include "AnasaziMultiVecTraits.hpp"
00009 #include "AnasaziOperatorTraits.hpp"
00010 #include "AnasaziConfigDefs.hpp"
00011
00012
00013 #include "PlayaVectorImpl.hpp"
00014 #include "PlayaLinearOperatorImpl.hpp"
00015 #include "PlayaSimpleTransposedOpImpl.hpp"
00016 #include "PlayaLinearCombinationImpl.hpp"
00017
00018 #include "Teuchos_Array.hpp"
00019 #include "Teuchos_RCP.hpp"
00020
00021 #include "PlayaTabs.hpp"
00022
00023 namespace Anasazi
00024 {
00025 using Playa::Vector;
00026 using Teuchos::RCP;
00027 using Teuchos::Array;
00028
00029 class SimpleMV
00030 {
00031 public:
00032 SimpleMV() : data_() {}
00033
00034 SimpleMV(int n) : data_(rcp(new Array<Vector<double> >(n))) {}
00035
00036 SimpleMV(const Array<Vector<double> >& data)
00037 : data_(rcp(new Array<Vector<double> >(data.size())))
00038 {
00039 for (int i=0; i<data.size(); i++)
00040 {
00041 (*data_)[i] = data[i].copy();
00042 }
00043 }
00044
00045 SimpleMV clone() const
00046 {
00047 return SimpleMV(*data_);
00048 }
00049
00050 Vector<double>& operator[](int i) {return (*data_)[i];}
00051
00052 const Vector<double>& operator[](int i) const {return (*data_)[i];}
00053
00054 int size() const {return data_->size();}
00055
00056 void resize(int n)
00057 {
00058 data_->resize(n);
00059 }
00060
00061 void randomize()
00062 {
00063 for (int i=0; i<size(); i++) (*data_)[i].randomize();
00064 }
00065
00066 private:
00067 RCP<Array<Vector<double> > > data_;
00068 };
00069
00070
00071
00072 inline std::ostream& operator<<(std::ostream& os, const SimpleMV& mv)
00073 {
00074 os << "MV (size=" << mv.size() << ")" << std::endl;
00075 for (int i=0; i<mv.size(); i++)
00076 {
00077 os << "ptr=" << mv[i].ptr().get() << std::endl;
00078 os << mv[i] << std::endl;
00079 }
00080
00081 return os;
00082 }
00083
00085 template<>
00086 class MultiVecTraits<double, SimpleMV >
00087 {
00088 public:
00089 typedef SimpleMV _MV;
00090 typedef Teuchos::ScalarTraits<double> SCT;
00091
00092 static double one() {static double rtn = SCT::one(); return rtn;}
00093 static double zero() {static double rtn = SCT::zero(); return rtn;}
00094
00097
00100 static RCP<_MV> Clone( const _MV & mv, const int numvecs )
00101 {
00102
00103 TEUCHOS_TEST_FOR_EXCEPT(mv.size() <= 0);
00104 TEUCHOS_TEST_FOR_EXCEPT(numvecs <= 0);
00105
00106 RCP<_MV> rtn = rcp(new _MV(numvecs));
00107 for (int i=0; i<numvecs; i++)
00108 {
00109 (*rtn)[i] = mv[0].copy();
00110 (*rtn)[i].setToConstant(zero());
00111 }
00112 return rtn;
00113 }
00114
00118 static RCP< _MV > CloneCopy( const _MV & mv )
00119 {
00120
00121 int numvecs = mv.size();
00122 TEUCHOS_TEST_FOR_EXCEPT(numvecs <= 0);
00123
00124
00125 RCP<_MV> rtn = rcp(new _MV(numvecs));
00126 for (int i=0; i<numvecs; i++)
00127 {
00128 (*rtn)[i] = mv[i].copy();
00129 }
00130 return rtn;
00131 }
00132
00136 static RCP< _MV > CloneCopy( const _MV & mv, const std::vector<int>& index )
00137 {
00138
00139 int numvecs = index.size();
00140 TEUCHOS_TEST_FOR_EXCEPT(numvecs <= 0);
00141 TEUCHOS_TEST_FOR_EXCEPT((int) index.size() > mv.size());
00142
00143
00144
00145
00146 RCP< _MV > rtn = rcp(new _MV(numvecs));
00147
00148 for (int i=0; i<numvecs; i++)
00149 {
00150 (*rtn)[i] = mv[index[i]].copy();
00151 }
00152 return rtn;
00153 }
00154
00158 static RCP< _MV > CloneViewNonConst( _MV & mv, const std::vector<int>& index )
00159 {
00160 int numvecs = index.size();
00161
00162
00163 TEUCHOS_TEST_FOR_EXCEPT(numvecs <= 0);
00164 TEUCHOS_TEST_FOR_EXCEPT((int) index.size() > mv.size());
00165
00166
00167
00168
00169 RCP< _MV > rtn = rcp(new _MV(numvecs));
00170
00171 for (int i=0; i<numvecs; i++)
00172 {
00173 (*rtn)[i] = mv[index[i]];
00174 }
00175
00176 return rtn;
00177 }
00178
00182 static RCP<const _MV > CloneView( const _MV & mv, const std::vector<int>& index )
00183 {
00184
00185 int numvecs = index.size();
00186
00187
00188 TEUCHOS_TEST_FOR_EXCEPT(numvecs <= 0);
00189 TEUCHOS_TEST_FOR_EXCEPT((int) index.size() > mv.size());
00190
00191
00192
00193
00194 RCP< const _MV > rtn = rcp(new _MV(numvecs));
00195
00196 for (int i=0; i<numvecs; i++)
00197 {
00198 (*(rcp_const_cast<_MV>(rtn)))[i] = mv[index[i]];
00199 }
00200 return rtn;
00201 }
00202
00204
00207
00209 static int GetVecLength( const _MV & mv )
00210 {
00211 TEUCHOS_TEST_FOR_EXCEPT(mv.size() <= 0);
00212 return mv[0].space().dim();
00213 }
00214
00216 static int GetNumberVecs( const _MV & mv )
00217 {
00218
00219 return mv.size();
00220 }
00221
00223
00226
00229 static void MvTimesMatAddMv( const double alpha, const _MV & A,
00230 const Teuchos::SerialDenseMatrix<int,double>& B,
00231 const double beta, _MV & mv )
00232 {
00233
00234 int n = B.numCols();
00235
00236
00237 TEUCHOS_TEST_FOR_EXCEPT(mv.size() != n);
00238
00239 for (int j=0; j<mv.size(); j++)
00240 {
00241 Vector<double> tmp;
00242 if (beta==one())
00243 {
00244 tmp = mv[j].copy();
00245 }
00246 else if (beta==zero())
00247 {
00248 tmp = mv[j].copy();
00249 tmp.setToConstant(zero());
00250 }
00251 else
00252 {
00253 tmp = beta * mv[j];
00254 }
00255 if (alpha != zero())
00256 {
00257 for (int i=0; i<A.size(); i++)
00258 {
00259 tmp = tmp + alpha*B(i,j)*A[i];
00260 }
00261 }
00262 mv[j].acceptCopyOf(tmp);
00263 }
00264 }
00265
00268 static void MvAddMv( const double alpha, const _MV & A,
00269 const double beta, const _MV & B, _MV & mv )
00270 {
00271 TEUCHOS_TEST_FOR_EXCEPT(A.size() != B.size());
00272 mv.resize(A.size());
00273 for (int i=0; i<A.size(); i++)
00274 {
00275 if (alpha==zero() && beta != zero()) mv[i]=beta*B[i];
00276 else if (beta==zero() && alpha != zero()) mv[i]=alpha*A[i];
00277 else if (alpha!=zero() && beta!=zero())
00278 mv[i]= alpha*A[i] + beta*B[i] ;
00279 else
00280 {
00281 mv[i].acceptCopyOf(A[i]);
00282 mv[i].setToConstant(zero());
00283 }
00284 }
00285 }
00286
00289 static void MvTransMv( const double alpha, const _MV & A, const _MV & mv,
00290 Teuchos::SerialDenseMatrix<int,double>& B )
00291 {
00292
00293 int m = A.size();
00294 int n = mv.size();
00295
00296
00297 for (int i=0; i<m; i++)
00298 {
00299 for (int j=0; j<n; j++)
00300 {
00301 B(i,j) = alpha * (A[i] * mv[j]);
00302 }
00303 }
00304
00305 }
00306
00310 static void MvDot( const _MV & mv, const _MV & A, std::vector<double> &b )
00311 {
00312
00313 TEUCHOS_TEST_FOR_EXCEPT(mv.size() != A.size());
00314 b.resize(A.size());
00315 for (int i=0; i<mv.size(); i++)
00316 b[i] = mv[i] * A[i];
00317 }
00318
00321 static void MvScale ( _MV & mv, const double alpha )
00322 {
00323
00324 for (int i=0; i<mv.size(); i++) mv[i].scale(alpha);
00325 }
00326
00329 static void MvScale ( _MV & mv, const std::vector<double>& alpha )
00330 {
00331
00332 for (int i=0; i<mv.size(); i++) mv[i].scale(alpha[i]);
00333 }
00334
00336
00339
00341 static void MvNorm( const _MV & mv,
00342 std::vector<Teuchos::ScalarTraits<double>::magnitudeType> &normvec )
00343 {
00344
00345 normvec.resize(mv.size());
00346 for (int i=0; i<mv.size(); i++)
00347 {
00348 normvec[i] = mv[i].norm2();
00349
00350 }
00351
00352 }
00353
00355
00358
00361 static void SetBlock( const _MV & A, const std::vector<int>& index, _MV & mv )
00362 {
00363
00364 TEUCHOS_TEST_FOR_EXCEPT(A.size() < (int) index.size());
00365
00366
00367 for (unsigned int i=0; i<index.size(); i++)
00368 {
00369 mv[index[i]].acceptCopyOf(A[i]);
00370 }
00371 }
00372
00375 static void MvRandom( _MV & mv )
00376 {
00377 for (int i=0; i<mv.size(); i++) mv[i].randomize();
00378 }
00379
00382 static void MvInit( _MV & mv, double alpha = Teuchos::ScalarTraits<double>::zero() )
00383 {
00384
00385 for (int i=0; i<mv.size(); i++) mv[i].setToConstant(alpha);
00386 }
00387
00389
00392
00394 static void MvPrint( const _MV & mv, std::ostream& os )
00395 {
00396 os << mv << std::endl;
00397 }
00398
00400
00402 static bool detectRepeatedIndex(const std::vector<int>& index)
00403 {
00404 std::set<int> s;
00405 bool rtn = false;
00406
00407 for (unsigned int i=0; i<index.size(); i++)
00408 {
00409 if (s.find(index[i]) != s.end())
00410 {
00411
00412 rtn = true;
00413 }
00414 s.insert(index[i]);
00415 }
00416
00417 return rtn;
00418 }
00419
00420 };
00421
00422
00426 template <>
00427 class OperatorTraits < double, SimpleMV, LinearOperator<double> >
00428 {
00429 public:
00430 typedef SimpleMV _MV;
00433 static void Apply (
00434 const LinearOperator< double >& Op,
00435 const _MV & x,
00436 _MV & y )
00437 {
00438
00439 y.resize(x.size());
00440 for (int i=0; i<x.size(); i++)
00441 {
00442
00443 y[i].acceptCopyOf(Op * x[i]);
00444
00445
00446
00447
00448
00449 }
00450 }
00451
00452 };
00453
00454
00455
00456 }
00457
00458 #endif