Thyra Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Thyra_TpetraVector_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Thyra: Interfaces and Support for Abstract Numerical Algorithms
4 //
5 // Copyright 2004 NTESS and the Thyra contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef THYRA_TPETRA_VECTOR_HPP
11 #define THYRA_TPETRA_VECTOR_HPP
12 
13 
16 #include "Kokkos_Core.hpp"
17 
18 namespace Thyra {
19 
20 
21 // Constructors/initializers/accessors
22 
23 
24 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
26 {}
27 
28 
29 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
31  const RCP<const TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraVectorSpace,
32  const RCP<Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraVector
33  )
34 {
35  initializeImpl(tpetraVectorSpace, tpetraVector);
36 }
37 
38 
39 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
41  const RCP<const TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraVectorSpace,
42  const RCP<const Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraVector
43  )
44 {
45  initializeImpl(tpetraVectorSpace, tpetraVector);
46 }
47 
48 
49 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
52 {
53  return tpetraVector_.getNonconstObj();
54 }
55 
56 
57 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
60 {
61  return tpetraVector_;
62 }
63 
64 
65 // Overridden from VectorDefaultBase
66 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
69 {
70  if (domainSpace_.is_null()) {
71  domainSpace_ = tpetraVectorSpace<Scalar>(
72  Tpetra::createLocalMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
73  1,
74  tpetraVector_.getConstObj()->getMap()->getComm()
75  )
76  );
77  }
78  return domainSpace_;
79 }
80 
81 
82 // Overridden from SpmdMultiVectorBase
83 
84 
85 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
88 {
89  return tpetraVectorSpace_;
90 }
91 
92 
93 // Overridden from SpmdVectorBase
94 
95 
96 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
98  const Ptr<ArrayRCP<Scalar> > &localValues )
99 {
100  *localValues = tpetraVector_.getNonconstObj()->get1dViewNonConst();
101 }
102 
103 
104 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
106  const Ptr<ArrayRCP<const Scalar> > &localValues ) const
107 {
108  *localValues = tpetraVector_->get1dView();
109 }
110 
111 
112 // Overridden from VectorBase
113 
114 
115 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
117  Scalar l,
118  Scalar u
119  )
120 {
121  // Tpetra randomizes with different seed for each proc, so need a global
122  // reduction to get locally-replicated random vector same on each proc.
123  if (!tpetraVector_.getNonconstObj()->isDistributed()) {
124  auto comm = tpetraVector_.getNonconstObj()->getMap()->getComm();
125  if (tpetraVector_.getConstObj()->getMap()->getComm()->getRank() == 0)
126  tpetraVector_.getNonconstObj()->randomize(l, u);
127  else
128  tpetraVector_.getNonconstObj()->putScalar(Teuchos::ScalarTraits<Scalar>::zero());
129  tpetraVector_.getNonconstObj()->reduce();
130  } else {
131  tpetraVector_.getNonconstObj()->randomize(l, u);
132  }
133 }
134 
135 
136 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
138  const VectorBase<Scalar>& x
139  )
140 {
141  auto tx = this->getConstTpetraVector(Teuchos::rcpFromRef(x));
142 
143  // If the cast succeeded, call Tpetra directly.
144  // Otherwise, fall back to the RTOp implementation.
145  if (nonnull(tx)) {
146  tpetraVector_.getNonconstObj()->abs(*tx);
147  } else {
148  VectorDefaultBase<Scalar>::absImpl(x);
149  }
150 }
151 
152 
153 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
155  const VectorBase<Scalar>& x
156  )
157 {
158  auto tx = this->getConstTpetraVector(Teuchos::rcpFromRef(x));
159 
160  // If the cast succeeded, call Tpetra directly.
161  // Otherwise, fall back to the RTOp implementation.
162  if (nonnull(tx)) {
163  tpetraVector_.getNonconstObj()->reciprocal(*tx);
164  } else {
165  VectorDefaultBase<Scalar>::reciprocalImpl(x);
166  }
167 }
168 
169 
170 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
172  const VectorBase<Scalar>& x
173  )
174 {
175  auto tx = this->getConstTpetraVector(Teuchos::rcpFromRef(x));
176 
177  // If the cast succeeded, call Tpetra directly.
178  // Otherwise, fall back to the RTOp implementation.
179  if (nonnull(tx)) {
181  tpetraVector_.getNonconstObj()->elementWiseMultiply(
182  ST::one(), *tx, *tpetraVector_.getConstObj(), ST::zero());
183  } else {
184  VectorDefaultBase<Scalar>::eleWiseScaleImpl(x);
185  }
186 }
187 
188 
189 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
192  const VectorBase<Scalar>& x
193  ) const
194 {
195  auto tx = this->getConstTpetraVector(Teuchos::rcpFromRef(x));
196 
197  // If the cast succeeded, call Tpetra directly.
198  // Otherwise, fall back to the RTOp implementation.
199  if (nonnull(tx)) {
200  // Weighted 2-norm function for Tpetra vector seems to be deprecated...
203  = Tpetra::createVector<Scalar>(tx->getMap());
204  temp->elementWiseMultiply(
205  ST::one(), *tx, *tpetraVector_.getConstObj(), ST::zero());
206  return ST::magnitude(ST::squareroot(tpetraVector_.getConstObj()->dot(*temp)));
207  } else {
208  return VectorDefaultBase<Scalar>::norm2WeightedImpl(x);
209  }
210 }
211 
212 
213 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
215  const RTOpPack::RTOpT<Scalar> &op,
216  const ArrayView<const Ptr<const VectorBase<Scalar> > > &vecs,
217  const ArrayView<const Ptr<VectorBase<Scalar> > > &targ_vecs,
218  const Ptr<RTOpPack::ReductTarget> &reduct_obj,
219  const Ordinal global_offset
220  ) const
221 {
222  SpmdVectorDefaultBase<Scalar>::applyOpImpl(op, vecs, targ_vecs, reduct_obj, global_offset);
223 }
224 
225 
226 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
229  const Range1D& rng,
231  ) const
232 {
233  SpmdVectorDefaultBase<Scalar>::acquireDetachedVectorViewImpl(rng, sub_vec);
234 }
235 
236 
237 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
240  const Range1D& rng,
242  )
243 {
244 
245  SpmdVectorDefaultBase<Scalar>::acquireNonconstDetachedVectorViewImpl(rng, sub_vec);
246 }
247 
248 
249 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
253  )
254 {
255  SpmdVectorDefaultBase<Scalar>::commitNonconstDetachedVectorViewImpl(sub_vec);
256 }
257 
258 
259 // Overridden protected functions from MultiVectorBase
260 
261 
262 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
264 {
265  tpetraVector_.getNonconstObj()->putScalar(alpha);
266 }
267 
268 
269 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
271 assignMultiVecImpl(const MultiVectorBase<Scalar>& mv)
272 {
273  auto tmv = this->getConstTpetraMultiVector(Teuchos::rcpFromRef(mv));
274 
275  // If cast succeeded, call Tpetra directly.
276  // Otherwise, fall back to the RTOp implementation.
277  if (nonnull(tmv)) {
278  tpetraVector_.getNonconstObj()->assign(*tmv);
279  } else {
280  MultiVectorDefaultBase<Scalar>::assignMultiVecImpl(mv);
281  }
282 }
283 
284 
285 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
287 {
288  tpetraVector_.getNonconstObj()->scale(alpha);
289 }
290 
291 
292 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
294  Scalar alpha,
295  const MultiVectorBase<Scalar>& mv
296  )
297 {
298  auto tmv = this->getConstTpetraMultiVector(Teuchos::rcpFromRef(mv));
299 
300  // If cast succeeded, call Tpetra directly.
301  // Otherwise, fall back to the RTOp implementation.
303  if (nonnull(tmv)) {
304  tpetraVector_.getNonconstObj()->update(alpha, *tmv, ST::one());
305  } else {
306  MultiVectorDefaultBase<Scalar>::updateImpl(alpha, mv);
307  }
308 }
309 
310 
311 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
313  const ArrayView<const Scalar>& alpha,
314  const ArrayView<const Ptr<const MultiVectorBase<Scalar> > >& mv,
315  const Scalar& beta
316  )
317 {
318 #ifdef TEUCHOS_DEBUG
319  TEUCHOS_ASSERT_EQUALITY(alpha.size(), mv.size());
320 #endif
321 
322  // Try to cast mv to Tpetra objects
325  bool allCastsSuccessful = true;
326  {
327  auto mvIter = mv.begin();
328  auto tmvIter = tmvs.begin();
329  for (; mvIter != mv.end(); ++mvIter, ++tmvIter) {
330  tmv = this->getConstTpetraMultiVector(Teuchos::rcpFromPtr(*mvIter));
331  if (nonnull(tmv)) {
332  *tmvIter = tmv;
333  } else {
334  allCastsSuccessful = false;
335  break;
336  }
337  }
338  }
339 
340  // If casts succeeded, or input arrays are size 0, call Tpetra directly.
341  // Otherwise, fall back to the RTOp implementation.
342  auto len = mv.size();
343  if (len == 0) {
344  tpetraVector_.getNonconstObj()->scale(beta);
345  } else if (len == 1 && allCastsSuccessful) {
346  tpetraVector_.getNonconstObj()->update(alpha[0], *tmvs[0], beta);
347  } else if (len == 2 && allCastsSuccessful) {
348  tpetraVector_.getNonconstObj()->update(alpha[0], *tmvs[0], alpha[1], *tmvs[1], beta);
349  } else if (allCastsSuccessful) {
351  auto tmvIter = tmvs.begin();
352  auto alphaIter = alpha.begin();
353 
354  // Check if any entry of tmvs aliases this object's wrapped vector.
355  // If so, replace that entry in the array with a copy.
356  tmv = Teuchos::null;
357  for (; tmvIter != tmvs.end(); ++tmvIter) {
358  if (tmvIter->getRawPtr() == tpetraVector_.getConstObj().getRawPtr()) {
359  if (tmv.is_null()) {
361  *tpetraVector_.getConstObj(), Teuchos::Copy));
362  }
363  *tmvIter = tmv;
364  }
365  }
366  tmvIter = tmvs.begin();
367 
368  // We add two MVs at a time, so only scale if even num MVs,
369  // and additionally do the first addition if odd num MVs.
370  if ((tmvs.size() % 2) == 0) {
371  tpetraVector_.getNonconstObj()->scale(beta);
372  } else {
373  tpetraVector_.getNonconstObj()->update(*alphaIter, *(*tmvIter), beta);
374  ++tmvIter;
375  ++alphaIter;
376  }
377  for (; tmvIter != tmvs.end(); tmvIter+=2, alphaIter+=2) {
378  tpetraVector_.getNonconstObj()->update(
379  *alphaIter, *(*tmvIter), *(alphaIter+1), *(*(tmvIter+1)), ST::one());
380  }
381  } else {
382  MultiVectorDefaultBase<Scalar>::linearCombinationImpl(alpha, mv, beta);
383  }
384 }
385 
386 
387 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
389  const MultiVectorBase<Scalar>& mv,
390  const ArrayView<Scalar>& prods
391  ) const
392 {
393  auto tmv = this->getConstTpetraMultiVector(Teuchos::rcpFromRef(mv));
394 
395  // If the cast succeeded, call Tpetra directly.
396  // Otherwise, fall back to the RTOp implementation.
397  if (nonnull(tmv)) {
398  tpetraVector_.getConstObj()->dot(*tmv, prods);
399  } else {
400  MultiVectorDefaultBase<Scalar>::dotsImpl(mv, prods);
401  }
402 }
403 
404 
405 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
407  const ArrayView<typename ScalarTraits<Scalar>::magnitudeType>& norms
408  ) const
409 {
410  tpetraVector_.getConstObj()->norm1(norms);
411 }
412 
413 
414 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
416  const ArrayView<typename ScalarTraits<Scalar>::magnitudeType>& norms
417  ) const
418 {
419  tpetraVector_.getConstObj()->norm2(norms);
420 }
421 
422 
423 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
425  const ArrayView<typename ScalarTraits<Scalar>::magnitudeType>& norms
426  ) const
427 {
428  tpetraVector_.getConstObj()->normInf(norms);
429 }
430 
431 
432 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
434  const EOpTransp M_trans,
435  const MultiVectorBase<Scalar> &X,
436  const Ptr<MultiVectorBase<Scalar> > &Y,
437  const Scalar alpha,
438  const Scalar beta
439  ) const
440 {
441  // Try to extract Tpetra objects from X and Y
442  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TMV;
443  Teuchos::RCP<const TMV> X_tpetra = this->getConstTpetraMultiVector(Teuchos::rcpFromRef(X));
444  Teuchos::RCP<TMV> Y_tpetra = this->getTpetraMultiVector(Teuchos::rcpFromPtr(Y));
445 
446  // If the cast succeeded, call Tpetra directly.
447  // Otherwise, fall back to the default implementation.
448  if (nonnull(X_tpetra) && nonnull(Y_tpetra)) {
450  TEUCHOS_TEST_FOR_EXCEPTION(ST::isComplex && (M_trans == CONJ),
451  std::logic_error,
452  "Error, conjugation without transposition is not allowed for complex scalar types!");
453 
455  switch (M_trans) {
456  case NOTRANS:
457  trans = Teuchos::NO_TRANS;
458  break;
459  case CONJ:
460  trans = Teuchos::NO_TRANS;
461  break;
462  case TRANS:
463  trans = Teuchos::TRANS;
464  break;
465  case CONJTRANS:
466  trans = Teuchos::CONJ_TRANS;
467  break;
468  }
469 
470  Y_tpetra->multiply(trans, Teuchos::NO_TRANS, alpha, *tpetraVector_.getConstObj(), *X_tpetra, beta);
471  Kokkos::fence();
472  } else {
473  VectorDefaultBase<Scalar>::applyImpl(M_trans, X, Y, alpha, beta);
474  }
475 
476 }
477 
478 
479 // private
480 
481 
482 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
483 template<class TpetraVector_t>
485  const RCP<const TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraVectorSpace,
486  const RCP<TpetraVector_t> &tpetraVector
487  )
488 {
489 #ifdef TEUCHOS_DEBUG
490  TEUCHOS_ASSERT(nonnull(tpetraVectorSpace));
491  TEUCHOS_ASSERT(nonnull(tpetraVector));
492 #endif
493  tpetraVectorSpace_ = tpetraVectorSpace;
494  tpetraVector_.initialize(tpetraVector);
495  this->updateSpmdSpace();
496 }
497 
498 
499 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
502 getTpetraMultiVector(const RCP<MultiVectorBase<Scalar> >& mv) const
503 {
504  using Teuchos::rcp_dynamic_cast;
507 
508  RCP<TMV> tmv = rcp_dynamic_cast<TMV>(mv);
509  if (nonnull(tmv)) {
510  return tmv->getTpetraMultiVector();
511  }
512 
513  RCP<TV> tv = rcp_dynamic_cast<TV>(mv);
514  if (nonnull(tv)) {
515  return tv->getTpetraVector();
516  }
517 
518  return Teuchos::null;
519 }
520 
521 
522 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
525 getConstTpetraMultiVector(const RCP<const MultiVectorBase<Scalar> >& mv) const
526 {
527  using Teuchos::rcp_dynamic_cast;
530 
531  RCP<const TMV> tmv = rcp_dynamic_cast<const TMV>(mv);
532  if (nonnull(tmv)) {
533  return tmv->getConstTpetraMultiVector();
534  }
535 
536  RCP<const TV> tv = rcp_dynamic_cast<const TV>(mv);
537  if (nonnull(tv)) {
538  return tv->getConstTpetraVector();
539  }
540 
541  return Teuchos::null;
542 }
543 
544 
545 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
548 getTpetraVector(const RCP<VectorBase<Scalar> >& v) const
549 {
551  RCP<TV> tv = Teuchos::rcp_dynamic_cast<TV>(v);
552  if (nonnull(tv)) {
553  return tv->getTpetraVector();
554  } else {
555  return Teuchos::null;
556  }
557 }
558 
559 
560 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
563 getConstTpetraVector(const RCP<const VectorBase<Scalar> >& v) const
564 {
566  RCP<const TV> tv = Teuchos::rcp_dynamic_cast<const TV>(v);
567  if (nonnull(tv)) {
568  return tv->getConstTpetraVector();
569  } else {
570  return Teuchos::null;
571  }
572 }
573 
574 
575 } // end namespace Thyra
576 
577 
578 #endif // THYRA_TPETRA_VECTOR_HPP
void constInitialize(const RCP< const TpetraVectorSpace< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &tpetraVectorSpace, const RCP< const Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &tpetraVector)
Initialize.
Concrete implementation of Thyra::MultiVector in terms of Tpetra::MultiVector.
RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetraMultiVector(const RCP< MultiVectorBase< Scalar > > &mv) const
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
Concrete implementation of an SPMD vector space for Tpetra.
void commitNonconstDetachedVectorViewImpl(RTOpPack::SubVectorView< Scalar > *sub_vec)
void initialize(const RCP< const TpetraVectorSpace< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &tpetraVectorSpace, const RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &tpetraVector)
Initialize.
RCP< const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getConstTpetraMultiVector(const RCP< const MultiVectorBase< Scalar > > &mv) const
iterator begin() const
TpetraVector()
Construct to uninitialized.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual void norms1Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
virtual void applyOpImpl(const RTOpPack::RTOpT< Scalar > &op, const ArrayView< const Ptr< const VectorBase< Scalar > > > &vecs, const ArrayView< const Ptr< VectorBase< Scalar > > > &targ_vecs, const Ptr< RTOpPack::ReductTarget > &reduct_obj, const Ordinal global_offset) const
size_type size() const
virtual void norms2Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
void getLocalVectorDataImpl(const Ptr< ArrayRCP< const Scalar > > &localValues) const
virtual void assignMultiVecImpl(const MultiVectorBase< Scalar > &mv)
void acquireNonconstDetachedVectorViewImpl(const Range1D &rng, RTOpPack::SubVectorView< Scalar > *sub_vec)
RCP< const Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getConstTpetraVector() const
Get the embedded non-const Tpetra::Vector.
virtual void reciprocalImpl(const VectorBase< Scalar > &x)
void acquireDetachedVectorViewImpl(const Range1D &rng, RTOpPack::ConstSubVectorView< Scalar > *sub_vec) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual Teuchos::ScalarTraits< Scalar >::magnitudeType norm2WeightedImpl(const VectorBase< Scalar > &x) const
virtual void linearCombinationImpl(const ArrayView< const Scalar > &alpha, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &mv, const Scalar &beta)
virtual void absImpl(const VectorBase< Scalar > &x)
void getNonconstLocalVectorDataImpl(const Ptr< ArrayRCP< Scalar > > &localValues)
Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > TpetraMultiVector_t
RCP< const SpmdVectorSpaceBase< Scalar > > spmdSpaceImpl() const
void initializeImpl(const RCP< const TpetraVectorSpace< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &tpetraVectorSpace, const RCP< TpetraVector_t > &tpetraVector)
Concrete Thyra::SpmdVectorBase using Tpetra::Vector.
virtual void dotsImpl(const MultiVectorBase< Scalar > &mv, const ArrayView< Scalar > &prods) const
RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetraVector()
Get the embedded non-const Tpetra::Vector.
bool nonnull(const boost::shared_ptr< T > &p)
virtual void assignImpl(Scalar alpha)
virtual void eleWiseScaleImpl(const VectorBase< Scalar > &x)
virtual void normsInfImpl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
RCP< const VectorSpaceBase< Scalar > > domain() const
virtual void randomizeImpl(Scalar l, Scalar u)
virtual void updateImpl(Scalar alpha, const MultiVectorBase< Scalar > &mv)
virtual void scaleImpl(Scalar alpha)
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
Teuchos_Ordinal Ordinal