PlayaMultiVectorOperatorImpl.hpp

00001 /* @HEADER@ */
00002 
00003  /* @HEADER@ */
00004 
00005 #ifndef PLAYA_MULTI_VECTOR_OPERATOR_IMPL_HPP
00006 #define PLAYA_MULTI_VECTOR_OPERATOR_IMPL_HPP
00007 
00008 #include "PlayaMultiVectorOperatorDecl.hpp"
00009 #include "PlayaVectorImpl.hpp"
00010 
00011 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00012 #include "PlayaLinearOperatorImpl.hpp"
00013 #include "PlayaSimpleTransposedOpImpl.hpp"
00014 #endif
00015 
00016 namespace Playa
00017 {
00018 
00019 /*
00020  * Construct from an array of vectors and a specifier for the 
00021  * domain space. 
00022  */
00023 template <class Scalar> inline
00024 MultiVectorOperator<Scalar>
00025 ::MultiVectorOperator(const Teuchos::Array<Vector<Scalar> >& cols,
00026   const VectorSpace<Scalar>& domain)
00027   : LinearOpWithSpaces<Scalar>(domain, cols[0].space()),
00028     cols_(cols)
00029 {
00030   TEUCHOS_TEST_FOR_EXCEPTION(cols.size() == 0, std::runtime_error,
00031     "empty multivector given to MultiVectorOperator ctor");
00032   for (int i=1; i<cols.size(); i++)
00033   {
00034     TEUCHOS_TEST_FOR_EXCEPTION(cols[i].space() != cols[0].space(), std::runtime_error,
00035       "inconsistent vector spaces in  MultiVectorOperator ctor");
00036   }
00037 }
00038 
00039 
00040 /*
00041  * Apply does an element-by-element multiply between the input 
00042  * vector, x, and the diagonal values.
00043  */
00044 template <class Scalar> inline
00045 void MultiVectorOperator<Scalar>
00046 ::apply(
00047   Teuchos::ETransp transApplyType,
00048   const Vector<Scalar>& in, 
00049   Vector<Scalar> out) const
00050 {
00051   if (transApplyType == NO_TRANS)
00052   {
00053     out.zero();
00054 
00055     for (int i=0; i<cols_.size(); i++)
00056     {
00057       out.update(in[i], cols_[i]);
00058     }
00059   }
00060   else
00061   {
00062     out.zero();
00063 
00064     for (int i=0; i<cols_.size(); i++)
00065     {
00066       out[i] = in.dot(cols_[i]);
00067     }
00068   }
00069 }
00070 
00071 
00072 
00073 /* Return the kth row  */
00074 template <class Scalar> inline
00075 void MultiVectorOperator<Scalar>
00076 ::getRow(const int& k, 
00077   Teuchos::Array<int>& indices, 
00078   Teuchos::Array<Scalar>& values) const
00079 {
00080   int low = this->range()->baseGlobalNaturalIndex();
00081   indices.resize(cols_.size());
00082   values.resize(cols_.size());
00083   for (int j=0; j<cols_.size(); j++)
00084   {
00085     indices[j] = j;
00086     values[j] = cols_[j][k-low];
00087   }
00088 }
00089 
00090 
00091   
00092 template <class Scalar> inline
00093 LinearOperator<Scalar> multiVectorOperator(
00094   const Teuchos::Array<Vector<Scalar> >& cols,
00095   const VectorSpace<Scalar>& domain)
00096 {
00097   RCP<LinearOperatorBase<Scalar> > A
00098     = rcp(new MultiVectorOperator<Scalar>(cols, domain));
00099 
00100   return A;
00101 }
00102 
00103 
00104 }
00105 
00106 #endif

doxygen