00001
00002
00003
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
00021
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
00042
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
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