42 #ifndef THYRA_EXAMPLE_TRIDIAG_SERIAL_LINEAR_OP_HPP 43 #define THYRA_EXAMPLE_TRIDIAG_SERIAL_LINEAR_OP_HPP 45 #include "Thyra_LinearOpDefaultBase.hpp" 46 #include "Thyra_DefaultSpmdVectorSpace.hpp" 47 #include "Thyra_DetachedVectorView.hpp" 48 #include "Teuchos_Assert.hpp" 78 template<
class Scalar>
89 const Teuchos::ArrayView<const Scalar> &lower,
90 const Teuchos::ArrayView<const Scalar> &diag,
91 const Teuchos::ArrayView<const Scalar> &upper
118 const Teuchos::ArrayView<const Scalar> &lower,
119 const Teuchos::ArrayView<const Scalar> &diag,
120 const Teuchos::ArrayView<const Scalar> &upper
123 TEUCHOS_TEST_FOR_EXCEPT( dim < 2 );
124 space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim);
135 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
range()
const 139 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
domain()
const 157 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > space_;
158 Teuchos::Array<Scalar> lower_;
159 Teuchos::Array<Scalar> diag_;
160 Teuchos::Array<Scalar> upper_;
165 template<
class Scalar>
175 typedef Teuchos::ScalarTraits<Scalar> ST;
178 const Ordinal dim = space_->dim();
184 for (
Ordinal col_j = 0; col_j < m; ++col_j) {
189 const Teuchos::ArrayRCP<const Scalar> x = x_vec.
sv().values();
190 const Teuchos::ArrayRCP<Scalar> y = y_vec.sv().values();
194 if( beta == ST::zero() ) {
195 for(
Ordinal k = 0; k < dim; ++k ) y[k] = ST::zero();
197 else if( beta != ST::one() ) {
198 for(
Ordinal k = 0; k < dim; ++k ) y[k] *= beta;
204 y[k] += alpha * ( diag_[k]*x[k] + upper_[k]*x[k+1] );
205 for( k = 1; k < dim - 1; ++k )
206 y[k] += alpha * ( lower_[k-1]*x[k-1] + diag_[k]*x[k] + upper_[k]*x[k+1] );
207 y[k] += alpha * ( lower_[k-1]*x[k-1] + diag_[k]*x[k] );
210 y[k] += alpha * ( ST::conjugate(diag_[k])*x[k] + ST::conjugate(upper_[k])*x[k+1] );
211 for( k = 1; k < dim - 1; ++k )
212 y[k] += alpha * ( ST::conjugate(lower_[k-1])*x[k-1]
213 + ST::conjugate(diag_[k])*x[k] + ST::conjugate(upper_[k])*x[k+1] );
214 y[k] += alpha * ( ST::conjugate(lower_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] );
217 y[k] += alpha * ( diag_[k]*x[k] + lower_[k]*x[k+1] );
218 for( k = 1; k < dim - 1; ++k )
219 y[k] += alpha * ( upper_[k-1]*x[k-1] + diag_[k]*x[k] + lower_[k]*x[k+1] );
220 y[k] += alpha * ( upper_[k-1]*x[k-1] + diag_[k]*x[k] );
223 y[k] += alpha * ( ST::conjugate(diag_[k])*x[k] + ST::conjugate(lower_[k])*x[k+1] );
224 for( k = 1; k < dim - 1; ++k )
225 y[k] += alpha * ( ST::conjugate(upper_[k-1])*x[k-1]
226 + ST::conjugate(diag_[k])*x[k] + ST::conjugate(lower_[k])*x[k+1] );
227 y[k] += alpha * ( ST::conjugate(upper_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] );
230 TEUCHOS_TEST_FOR_EXCEPT(
true);
237 #endif // THYRA_EXAMPLE_TRIDIAG_SERIAL_LINEAR_OP_HPP const RTOpPack::SubVectorView< Scalar > & sv() const
Returns the explicit view as an RTOpPack::ConstSubVectorView<Scalar> object.
ExampleTridiagSerialLinearOp(const Thyra::Ordinal dim, const Teuchos::ArrayView< const Scalar > &lower, const Teuchos::ArrayView< const Scalar > &diag, const Teuchos::ArrayView< const Scalar > &upper)
initialize().
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
bool opSupportedImpl(Thyra::EOpTransp M_trans) const
void initialize(const Thyra::Ordinal dim, const Teuchos::ArrayView< const Scalar > &lower, const Teuchos::ArrayView< const Scalar > &diag, const Teuchos::ArrayView< const Scalar > &upper)
Use the non-transposed operator.
Node subclass that provides a good default implementation for the describe() function.
Create an explicit non-mutable (const) view of a VectorBase object.
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types)...
void applyImpl(const Thyra::EOpTransp M_trans, const Thyra::MultiVectorBase< Scalar > &X_in, const Teuchos::Ptr< Thyra::MultiVectorBase< Scalar > > &Y_inout, const Scalar alpha, const Scalar beta) const
Use the non-transposed operator with complex-conjugate elements (same as NOTRANS for real scalar type...
Use the transposed operator.
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
Interface for a collection of column vectors called a multi-vector.
Create an explicit mutable (non-const) view of a VectorBase object.
RCP< const VectorBase< Scalar > > col(Ordinal j) const
Calls colImpl().
ExampleTridiagSerialLinearOp()
Construct to uninitialized.
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > domain() const
Simple example subclass for serial tridiagonal matrices.
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > range() const