42 #ifndef THYRA_EXAMPLE_TRIDIAG_SPMD_LINEAR_OP_HPP 43 #define THYRA_EXAMPLE_TRIDIAG_SPMD_LINEAR_OP_HPP 45 #include "Thyra_LinearOpDefaultBase.hpp" 46 #include "Thyra_DefaultSpmdVectorSpace.hpp" 47 #include "Thyra_DetachedSpmdVectorView.hpp" 48 #include "Teuchos_DefaultComm.hpp" 49 #include "Teuchos_CommHelpers.hpp" 142 template<
class Scalar>
151 const Teuchos::RCP<
const Teuchos::Comm<Thyra::Ordinal> > &comm,
153 const Teuchos::ArrayView<const Scalar> &lower,
154 const Teuchos::ArrayView<const Scalar> &diag,
155 const Teuchos::ArrayView<const Scalar> &upper
157 { this->
initialize(comm, localDim, lower, diag, upper); }
181 const Teuchos::RCP<
const Teuchos::Comm<Thyra::Ordinal> > &comm,
183 const Teuchos::ArrayView<const Scalar> &lower,
184 const Teuchos::ArrayView<const Scalar> &diag,
185 const Teuchos::ArrayView<const Scalar> &upper
193 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
range()
const 197 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
domain()
const 217 Teuchos::RCP<const Teuchos::Comm<Thyra::Ordinal> > comm_;
218 Teuchos::RCP<const Thyra::SpmdVectorSpaceBase<Scalar> > space_;
219 Teuchos::Array<Scalar> lower_;
220 Teuchos::Array<Scalar> diag_;
221 Teuchos::Array<Scalar> upper_;
223 void communicate(
const bool first,
const bool last,
224 const Teuchos::ArrayView<const Scalar> &x,
225 const Teuchos::Ptr<Scalar> &x_km1,
226 const Teuchos::Ptr<Scalar> &x_kp1 )
const;
231 template<
class Scalar>
233 const Teuchos::RCP<
const Teuchos::Comm<Thyra::Ordinal> > &comm,
235 const Teuchos::ArrayView<const Scalar> &lower,
236 const Teuchos::ArrayView<const Scalar> &diag,
237 const Teuchos::ArrayView<const Scalar> &upper
240 TEUCHOS_TEST_FOR_EXCEPT( localDim < 2 );
242 space_ = Thyra::defaultSpmdVectorSpace<Scalar>(comm, localDim, -1);
249 template<
class Scalar>
259 typedef Teuchos::ScalarTraits<Scalar> ST;
261 using Teuchos::outArg;
266 const Scalar
zero = ST::zero();
267 const Ordinal localDim = space_->localSubDim();
268 const Ordinal procRank = comm_->getRank();
269 const Ordinal numProcs = comm_->getSize();
274 for (
Ordinal col_j = 0; col_j < m; ++col_j) {
279 const Teuchos::ArrayRCP<const Scalar> x = x_vec.
sv().values();
280 const Teuchos::ArrayRCP<Scalar> y = y_vec.sv().values();
283 const bool first = (procRank == 0), last = ( procRank == numProcs-1 );
287 communicate( first, last, x(), outArg(x_km1), outArg(x_kp1) );
294 y[k] = alpha * ( (first?
zero:lower_[lk]*x_km1) + diag_[k]*x[k]
295 + upper_[k]*x[k+1] );
298 for( k = 1; k < localDim - 1; ++lk, ++k )
299 y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k] + upper_[k]*x[k+1] );
301 y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k]
302 + (last?
zero:upper_[k]*x_kp1) );
306 y[k] = alpha * ( (first?
zero:lower_[lk]*x_km1) + diag_[k]*x[k]
307 + upper_[k]*x[k+1] ) + beta*y[k];
310 for( k = 1; k < localDim - 1; ++lk, ++k )
311 y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k] + upper_[k]*x[k+1] )
314 y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k]
315 + (last?
zero:upper_[k]*x_kp1) ) + beta*y[k];
323 template<
class Scalar>
325 const bool first,
const bool last,
326 const Teuchos::ArrayView<const Scalar> &x,
327 const Teuchos::Ptr<Scalar> &x_km1,
328 const Teuchos::Ptr<Scalar> &x_kp1
332 const Ordinal localDim = space_->localSubDim();
333 const Ordinal procRank = comm_->getRank();
334 const Ordinal numProcs = comm_->getSize();
335 const bool isEven = (procRank % 2 == 0);
336 const bool isOdd = !isEven;
339 if(isEven && procRank+1 < numProcs) send(*comm_, x[localDim-1], procRank+1);
340 if(isOdd && procRank-1 >= 0) receive(*comm_, procRank-1, &*x_km1);
343 if( isOdd && procRank-1 >= 0 ) send(*comm_, x[0], procRank-1);
344 if( isEven && procRank+1 < numProcs ) receive(*comm_, procRank+1, &*x_kp1);
347 if (isOdd && procRank+1 < numProcs) send(*comm_, x[localDim-1], procRank+1);
348 if (isEven && procRank-1 >= 0) receive(*comm_, procRank-1, &*x_km1);
351 if (isEven && procRank-1 >= 0) send(*comm_, x[0], procRank-1);
352 if (isOdd && procRank+1 < numProcs) receive(*comm_, procRank+1, &*x_kp1);
356 #endif // THYRA_EXAMPLE_TRIDIAG_SPMD_LINEAR_OP_HPP Create an explicit detached mutable (non-const) view of all of the local elements on this process of ...
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
ExampleTridiagSpmdLinearOp(const Teuchos::RCP< const Teuchos::Comm< Thyra::Ordinal > > &comm, const Thyra::Ordinal localDim, 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.
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
const RTOpPack::SubVectorView< Scalar > & sv() const
RCP< const LinearOpBase< Scalar > > zero(const RCP< const VectorSpaceBase< Scalar > > &range, const RCP< const VectorSpaceBase< Scalar > > &domain)
Create a zero linear operator with given range and domain spaces.
Simple example subclass for Spmd tridiagonal matrices.
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
Interface for a collection of column vectors called a multi-vector.
bool opSupported(const LinearOpBase< Scalar > &M, EOpTransp M_trans)
Determines if an operation is supported for a single scalar type.
RCP< const VectorBase< Scalar > > col(Ordinal j) const
Calls colImpl().
bool opSupportedImpl(Thyra::EOpTransp M_trans) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > domain() const
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
void initialize(const Teuchos::RCP< const Teuchos::Comm< Thyra::Ordinal > > &comm, const Thyra::Ordinal localDim, const Teuchos::ArrayView< const Scalar > &lower, const Teuchos::ArrayView< const Scalar > &diag, const Teuchos::ArrayView< const Scalar > &upper)
ExampleTridiagSpmdLinearOp()
Create an explicit detached non-mutable (const) view of all of the local elements on this process of ...
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > range() const