Thyra Version of the Day
Loading...
Searching...
No Matches
ExampleTridiagSpmdLinearOp.hpp
1// @HEADER
2// ***********************************************************************
3//
4// Thyra: Interfaces and Support for Abstract Numerical Algorithms
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef THYRA_EXAMPLE_TRIDIAG_SPMD_LINEAR_OP_HPP
43#define THYRA_EXAMPLE_TRIDIAG_SPMD_LINEAR_OP_HPP
44
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"
50
51
142template<class Scalar>
144public:
145
148
152 const Thyra::Ordinal localDim,
156 )
157 { this->initialize(comm, localDim, lower, diag, upper); }
158
180 void initialize(
182 const Thyra::Ordinal localDim,
186 );
187
188protected:
189
190 // Overridden from LinearOpBase
191
195
199
202 {
203 return (M_trans == Thyra::NOTRANS);
204 }
205
207 void applyImpl(
208 const Thyra::EOpTransp M_trans,
211 const Scalar alpha,
212 const Scalar beta
213 ) const;
214
215private:
216
219 Teuchos::Array<Scalar> lower_; // size = ( procRank == 0 ? localDim - 1 : localDim )
220 Teuchos::Array<Scalar> diag_; // size = localDim
221 Teuchos::Array<Scalar> upper_; // size = ( procRank == numProc-1 ? localDim - 1 : localDim )
222
223 void communicate( const bool first, const bool last,
225 const Teuchos::Ptr<Scalar> &x_km1,
226 const Teuchos::Ptr<Scalar> &x_kp1 ) const;
227
228}; // end class ExampleTridiagSpmdLinearOp
229
230
231template<class Scalar>
234 const Thyra::Ordinal localDim,
238 )
239{
240 TEUCHOS_TEST_FOR_EXCEPT( localDim < 2 );
241 comm_ = comm;
242 space_ = Thyra::defaultSpmdVectorSpace<Scalar>(comm, localDim, -1);
243 lower_ = lower;
244 diag_ = diag;
245 upper_ = upper;
246}
247
248
249template<class Scalar>
251 const Thyra::EOpTransp M_trans,
254 const Scalar alpha,
255 const Scalar beta
256 ) const
257{
258
260 typedef Thyra::Ordinal Ordinal;
261 using Teuchos::outArg;
262
263 TEUCHOS_ASSERT(this->opSupported(M_trans));
264
265 // Get constants
266 const Scalar zero = ST::zero();
267 const Ordinal localDim = space_->localSubDim();
268 const Ordinal procRank = comm_->getRank();
269 const Ordinal numProcs = comm_->getSize();
270 const Ordinal m = X_in.domain()->dim();
271
272 // Loop over the input columns
273
274 for (Ordinal col_j = 0; col_j < m; ++col_j) {
275
276 // Get access the the elements of column j
278 Thyra::DetachedSpmdVectorView<Scalar> y_vec(Y_inout->col(col_j));
279 const Teuchos::ArrayRCP<const Scalar> x = x_vec.sv().values();
280 const Teuchos::ArrayRCP<Scalar> y = y_vec.sv().values();
281
282 // Determine what process we are
283 const bool first = (procRank == 0), last = ( procRank == numProcs-1 );
284
285 // Communicate ghost elements
286 Scalar x_km1, x_kp1;
287 communicate( first, last, x(), outArg(x_km1), outArg(x_kp1) );
288
289 // Perform operation (if beta==0 then we must be careful since y could
290 // be uninitialized on input!)
291 Thyra::Ordinal k = 0, lk = 0;
292 if( beta == zero ) {
293 // First local row
294 y[k] = alpha * ( (first?zero:lower_[lk]*x_km1) + diag_[k]*x[k]
295 + upper_[k]*x[k+1] );
296 if(!first) ++lk;
297 // Middle local rows
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] );
300 // Last local row
301 y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k]
302 + (last?zero:upper_[k]*x_kp1) );
303 }
304 else {
305 // First local row
306 y[k] = alpha * ( (first?zero:lower_[lk]*x_km1) + diag_[k]*x[k]
307 + upper_[k]*x[k+1] ) + beta*y[k];
308 if(!first) ++lk;
309 // Middle local rows
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] )
312 + beta*y[k];
313 // Last local row
314 y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k]
315 + (last?zero:upper_[k]*x_kp1) ) + beta*y[k];
316 }
317
318 }
319
320}
321
322
323template<class Scalar>
325 const bool first, const bool last,
327 const Teuchos::Ptr<Scalar> &x_km1,
328 const Teuchos::Ptr<Scalar> &x_kp1
329 ) const
330{
331 typedef Thyra::Ordinal Ordinal;
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;
337 // 1) Send x[localDim-1] from each even process forward to the next odd
338 // process where it is received in x_km1.
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);
341 // 2) Send x[0] from each odd process backward to the previous even
342 // process where it is received in x_kp1.
343 if( isOdd && procRank-1 >= 0 ) send(*comm_, x[0], procRank-1);
344 if( isEven && procRank+1 < numProcs ) receive(*comm_, procRank+1, &*x_kp1);
345 // 3) Send x[localDim-1] from each odd process forward to the next even
346 // process where it is received in x_km1.
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);
349 // 4) Send x[0] from each even process backward to the previous odd
350 // process where it is received in x_kp1.
351 if (isEven && procRank-1 >= 0) send(*comm_, x[0], procRank-1);
352 if (isOdd && procRank+1 < numProcs) receive(*comm_, procRank+1, &*x_kp1);
353}
354
355
356#endif // THYRA_EXAMPLE_TRIDIAG_SPMD_LINEAR_OP_HPP
Simple example subclass for Spmd tridiagonal matrices.
bool opSupportedImpl(Thyra::EOpTransp M_trans) const
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)
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
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > domain() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > range() const
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)
const ArrayRCP< const Scalar > values() const
const ArrayRCP< Scalar > values() const
Create an explicit detached non-mutable (const) view of all of the local elements on this process of ...
const RTOpPack::ConstSubVectorView< Scalar > & sv() const
Create an explicit detached mutable (non-const) view of all of the local elements on this process of ...
const RTOpPack::SubVectorView< Scalar > & sv() const
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
Node subclass that provides a good default implementation for the describe() function.
Interface for a collection of column vectors called a multi-vector.
RCP< const VectorBase< Scalar > > col(Ordinal j) const
Calls colImpl().
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
@ NOTRANS
Use the non-transposed operator.