41 #include "Teuchos_UnitTestHarness.hpp" 42 #include "Teuchos_UnitTestRepository.hpp" 43 #include "Teuchos_GlobalMPISession.hpp" 44 #include "Teuchos_TestingHelpers.hpp" 45 #include "Teuchos_CommHelpers.hpp" 46 #include "Teuchos_DefaultComm.hpp" 47 #include "Teuchos_Array.hpp" 48 #include "Teuchos_Comm.hpp" 53 #include "Sacado_Fad_DFad.hpp" 54 #include "Sacado_mpl_apply.hpp" 55 #include "Sacado_Random.hpp" 61 template <
typename PCEType,
typename FadType>
63 RCP<const Stokhos::CompletePolynomialBasis<int,double> >
basis;
64 RCP<Stokhos::Sparse3Tensor<int,double> >
Cijk;
72 typedef typename Sacado::mpl::apply<FadType,PCEType>::type
FadPCEType;
82 Teuchos::Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
83 for (
int i=0; i<d; i++)
91 Cijk =
basis->computeTripleProductTensor();
93 Stokhos::create_product_tensor<execution_space>(*
basis, *
Cijk);
99 rcp(
new Teuchos::ValueTypeSerializer<int,double>())));
106 template <
typename PCEType>
108 const Teuchos::Array<PCEType>& x2,
109 const std::string& tag,
110 Teuchos::FancyOStream& out) {
113 bool success = (
x.size() == x2.size());
114 out << tag <<
" PCE array size test";
119 out <<
": \n\tExpected: " <<
x.size() <<
", \n\tGot: " << x2.size()
123 for (
int i=0; i<
x.size(); i++) {
124 bool success2 = Sacado::IsEqual<PCEType>::eval(
x[i], x2[i]);
125 out << tag <<
" PCE array comparison test " << i;
130 out <<
": \n\tExpected: " <<
x[i] <<
", \n\tGot: " << x2[i] <<
"." 132 success = success && success2;
138 template<
typename Ordinal>
140 const Teuchos::Comm<Ordinal> &comm,
141 Teuchos::FancyOStream &out,
145 out <<
"\nChecking that the above test passed in all processes ...";
146 int thisResult = ( result ? 1 : 0 );
148 Teuchos::reduceAll(comm,Teuchos::REDUCE_SUM,
Ordinal(1),&thisResult,
150 const bool passed = sumResult==Teuchos::size(comm);
154 out <<
" (sumResult="<<sumResult<<
"!=numProcs="<<Teuchos::size(comm)<<
") failed\n";
164 Sacado::Random<double>
rnd;
167 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
168 comm = Teuchos::DefaultComm<Ordinal>::getComm();
172 Teuchos::Array<PCEType>
x(n), x2(n);
173 for (
int i=0; i<n; i++) {
174 x[i].reset(
setup.kokkos_cijk);
176 x[i].fastAccessCoeff(
j) =
rnd.number();
178 if (comm->getRank() == 0)
180 Teuchos::broadcast(*comm, *
setup.pce_serializer, 0, n, &x2[0]);
181 success =
checkPCEArrays(
x, x2, std::string(
"UQ::PCE")+
" Broadcast", out);
186 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
187 comm = Teuchos::DefaultComm<Ordinal>::getComm();
191 int size = comm->getSize();
192 int rank = comm->getRank();
194 Teuchos::Array<PCEType>
x(n), x2(N), x3(N);
195 for (
int i=0; i<n; i++) {
196 x[i].reset(
setup.kokkos_cijk);
198 x[i].fastAccessCoeff(
j) = (rank+1)*(i+1)*(
j+1);
200 for (
int j=0;
j<size;
j++) {
201 for (
int i=0; i<n; i++) {
202 x3[n*
j+i].reset(
setup.kokkos_cijk);
203 for (
int k=0; k<
setup.sz; k++)
204 x3[n*
j+i].fastAccessCoeff(k) = (
j+1)*(i+1)*(k+1);
207 Teuchos::gatherAll(*comm, *
setup.pce_serializer,
208 n, &
x[0], N, &x2[0]);
209 success =
checkPCEArrays(x3, x2, std::string(
"UQ::PCE")+
" Gather All", out);
214 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
215 comm = Teuchos::DefaultComm<Ordinal>::getComm();
219 int num_proc = comm->getSize();
221 Teuchos::Array<PCEType>
x(n), sums(n), sums2(n);
222 for (
int i=0; i<n; i++) {
223 x[i].reset(
setup.kokkos_cijk);
225 x[i].fastAccessCoeff(
j) = 2.0*(i+1);
227 for (
int i=0; i<n; i++) {
228 sums[i].reset(
setup.kokkos_cijk);
230 sums[i].fastAccessCoeff(
j) = 2.0*(i+1)*num_proc;
232 Teuchos::reduceAll(*comm, *
setup.pce_serializer,
233 Teuchos::REDUCE_SUM, n, &
x[0], &sums2[0]);
235 std::string(
"UQ::PCE")+
" Sum All", out);
240 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
241 comm = Teuchos::DefaultComm<Ordinal>::getComm();
245 int rank = comm->getRank();
246 int num_proc = comm->getSize();
248 Teuchos::Array<PCEType>
x(n), maxs(n), maxs2(n);
249 for (
int i=0; i<n; i++) {
250 x[i].reset(
setup.kokkos_cijk);
252 x[i].fastAccessCoeff(
j) = 2.0*(i+1)*(rank+1);
254 for (
int i=0; i<n; i++) {
255 maxs[i].reset(
setup.kokkos_cijk);
257 maxs[i].fastAccessCoeff(
j) = 2.0*(i+1)*num_proc;
259 Teuchos::reduceAll(*comm, *
setup.pce_serializer,
260 Teuchos::REDUCE_MAX, n, &
x[0], &maxs2[0]);
262 std::string(
"UQ::PCE")+
" Max All", out);
267 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
268 comm = Teuchos::DefaultComm<Ordinal>::getComm();
272 int rank = comm->getRank();
274 Teuchos::Array<PCEType>
x(n), mins(n), mins2(n);
275 for (
int i=0; i<n; i++) {
276 x[i].reset(
setup.kokkos_cijk);
278 x[i].fastAccessCoeff(
j) = 2.0*(i+1)*(rank+1);
280 for (
int i=0; i<n; i++) {
281 mins[i].reset(
setup.kokkos_cijk);
283 mins[i].fastAccessCoeff(
j) = 2.0*(i+1);
285 Teuchos::reduceAll(*comm, *
setup.pce_serializer,
286 Teuchos::REDUCE_MIN, n, &
x[0], &mins2[0]);
288 std::string(
"UQ::PCE")+
" Min All", out);
293 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
294 comm = Teuchos::DefaultComm<Ordinal>::getComm();
298 int rank = comm->getRank();
300 Teuchos::Array<PCEType>
x(n), sums(n), sums2(n);
301 for (
int i=0; i<n; i++) {
302 x[i].reset(
setup.kokkos_cijk);
304 x[i].fastAccessCoeff(
j) = 2.0*(i+1);
306 for (
int i=0; i<n; i++) {
307 sums[i].reset(
setup.kokkos_cijk);
309 sums[i].fastAccessCoeff(
j) = 2.0*(i+1)*(rank+1);
311 Teuchos::scan(*comm, *
setup.pce_serializer,
312 Teuchos::REDUCE_SUM, n, &
x[0], &sums2[0]);
314 std::string(
"UQ::PCE")+
" Scan Sum", out);
319 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
320 comm = Teuchos::DefaultComm<Ordinal>::getComm();
324 int rank = comm->getRank();
326 Teuchos::Array<PCEType>
x(n), maxs(n), maxs2(n);
327 for (
int i=0; i<n; i++) {
328 x[i].reset(
setup.kokkos_cijk);
330 x[i].fastAccessCoeff(
j) = 2.0*(i+1)*(rank+1);
332 for (
int i=0; i<n; i++) {
333 maxs[i].reset(
setup.kokkos_cijk);
335 maxs[i].fastAccessCoeff(
j) = 2.0*(i+1)*(rank+1);
337 Teuchos::scan(*comm, *
setup.pce_serializer,
338 Teuchos::REDUCE_MAX, n, &
x[0], &maxs2[0]);
340 std::string(
"UQ::PCE")+
" Scan Max", out);
345 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
346 comm = Teuchos::DefaultComm<Ordinal>::getComm();
350 int rank = comm->getRank();
352 Teuchos::Array<PCEType>
x(n), mins(n), mins2(n);
353 for (
int i=0; i<n; i++) {
354 x[i].reset(
setup.kokkos_cijk);
356 x[i].fastAccessCoeff(
j) = 2.0*(i+1)*(rank+1);
358 for (
int i=0; i<n; i++) {
359 mins[i].reset(
setup.kokkos_cijk);
361 mins[i].fastAccessCoeff(
j) = 2.0*(i+1);
363 Teuchos::scan(*comm, *
setup.pce_serializer,
364 Teuchos::REDUCE_MIN, n, &
x[0], &mins2[0]);
366 std::string(
"UQ::PCE")+
" Scan Min", out);
371 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
372 comm = Teuchos::DefaultComm<Ordinal>::getComm();
375 int num_proc = comm->getSize();
377 int rank = comm->getRank();
379 Teuchos::Array<PCEType>
x(n), x2(n);
380 for (
int i=0; i<n; i++) {
381 x[i].reset(
setup.kokkos_cijk);
383 x[i].fastAccessCoeff(
j) = 2.0*(i+1)*(
j+1);
387 if (rank == 0) Teuchos::send(*comm, *
setup.pce_serializer,
389 if (rank == 1) Teuchos::receive(*comm, *
setup.pce_serializer,
392 std::string(
"UQ::PCE")+
" Send/Receive", out);
400 typedef Sacado::mpl::apply<FadType,PCEType>::type FadPCEType;
401 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
402 comm = Teuchos::DefaultComm<Ordinal>::getComm();
407 Teuchos::Array<FadPCEType>
x(n), x2(n);
408 for (
int i=0; i<n; i++) {
410 for (
int k=0; k<
setup.sz; k++)
411 f.fastAccessCoeff(k) =
rnd.number();
412 x[i] = FadPCEType(p,
f);
413 for (
int j=0;
j<p;
j++) {
415 for (
int k=0; k<
setup.sz; k++)
416 g.fastAccessCoeff(k) =
rnd.number();
417 x[i].fastAccessDx(
j) =
g;
420 if (comm->getRank() == 0)
422 Teuchos::broadcast(*comm, *
setup.fad_pce_serializer, 0, n, &x2[0]);
424 std::string(
"DFad")+
"<"+
"UQ::PCE"+
"> Broadcast", out);
429 typedef Sacado::mpl::apply<FadType,PCEType>::type FadPCEType;
430 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
431 comm = Teuchos::DefaultComm<Ordinal>::getComm();
436 int size = comm->getSize();
437 int rank = comm->getRank();
439 Teuchos::Array<FadPCEType>
x(n), x2(N), x3(N);
440 for (
int i=0; i<n; i++) {
442 for (
int k=0; k<
setup.sz; k++)
443 f.fastAccessCoeff(k) = (rank+1)*(i+1)*(k+1);
444 x[i] = FadPCEType(p,
f);
445 for (
int j=0;
j<p;
j++) {
446 x[i].fastAccessDx(
j) =
f;
449 for (
int j=0;
j<size;
j++) {
450 for (
int i=0; i<n; i++) {
452 for (
int k=0; k<
setup.sz; k++)
453 f.fastAccessCoeff(k) = (
j+1)*(i+1)*(k+1);
454 x3[n*
j+i] = FadPCEType(p,
f);
455 for (
int k=0; k<p; k++)
459 Teuchos::gatherAll(*comm, *
setup.fad_pce_serializer,
460 n, &
x[0], N, &x2[0]);
462 std::string(
"DFad")+
"<"+
"UQ::PCE"+
"> Gather All", out);
467 typedef Sacado::mpl::apply<FadType,PCEType>::type FadPCEType;
468 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
469 comm = Teuchos::DefaultComm<Ordinal>::getComm();
474 int num_proc = comm->getSize();
476 Teuchos::Array<FadPCEType>
x(n), sums(n), sums2(n);
477 for (
int i=0; i<n; i++) {
479 for (
int k=0; k<
setup.sz; k++)
480 f.fastAccessCoeff(k) = 2.0*(i+1);
481 x[i] = FadPCEType(p,
f);
482 for (
int j=0;
j<p;
j++) {
484 for (
int k=0; k<
setup.sz; k++)
485 g.fastAccessCoeff(k) = 2.0*(i+1);
486 x[i].fastAccessDx(
j) =
g;
489 for (
int i=0; i<n; i++) {
491 for (
int k=0; k<
setup.sz; k++)
492 f.fastAccessCoeff(k) = 2.0*(i+1)*num_proc;
493 sums[i] = FadPCEType(p,
f);
494 for (
int j=0;
j<p;
j++) {
496 for (
int k=0; k<
setup.sz; k++)
497 g.fastAccessCoeff(k) = 2.0*(i+1)*num_proc;
498 sums[i].fastAccessDx(
j) =
g;
501 Teuchos::reduceAll(*comm, *
setup.fad_pce_serializer,
502 Teuchos::REDUCE_SUM, n, &
x[0], &sums2[0]);
504 std::string(
"DFad")+
"<"+
"UQ::PCE"+
"> Sum All", out);
509 typedef Sacado::mpl::apply<FadType,PCEType>::type FadPCEType;
510 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
511 comm = Teuchos::DefaultComm<Ordinal>::getComm();
516 int rank = comm->getRank();
517 int num_proc = comm->getSize();
519 Teuchos::Array<FadPCEType>
x(n), maxs(n), maxs2(n);
520 for (
int i=0; i<n; i++) {
522 for (
int k=0; k<
setup.sz; k++)
523 f.fastAccessCoeff(k) = 2.0*(i+1)*(rank+1);
524 x[i] = FadPCEType(p,
f);
525 for (
int j=0;
j<p;
j++) {
526 x[i].fastAccessDx(
j) =
f;
529 for (
int i=0; i<n; i++) {
531 for (
int k=0; k<
setup.sz; k++)
532 f.fastAccessCoeff(k) = 2.0*(i+1)*num_proc;
533 maxs[i] = FadPCEType(p,
f);
534 for (
int j=0;
j<p;
j++)
537 Teuchos::reduceAll(*comm, *
setup.fad_pce_serializer,
538 Teuchos::REDUCE_MAX, n, &
x[0], &maxs2[0]);
540 std::string(
"DFad")+
"<"+
"UQ::PCE"+
"> Max All", out);
545 typedef Sacado::mpl::apply<FadType,PCEType>::type FadPCEType;
546 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
547 comm = Teuchos::DefaultComm<Ordinal>::getComm();
552 int rank = comm->getRank();
554 Teuchos::Array<FadPCEType>
x(n), mins(n), mins2(n);
555 for (
int i=0; i<n; i++) {
557 for (
int k=0; k<
setup.sz; k++)
558 f.fastAccessCoeff(k) = 2.0*(i+1)*(rank+1);
559 x[i] = FadPCEType(p,
f);
560 for (
int j=0;
j<p;
j++) {
561 x[i].fastAccessDx(
j) =
f;
564 for (
int i=0; i<n; i++) {
566 for (
int k=0; k<
setup.sz; k++)
567 f.fastAccessCoeff(k) = 2.0*(i+1);
568 mins[i] = FadPCEType(p,
f);
569 for (
int j=0;
j<p;
j++)
572 Teuchos::reduceAll(*comm, *
setup.fad_pce_serializer,
573 Teuchos::REDUCE_MIN, n, &
x[0], &mins2[0]);
575 std::string(
"DFad")+
"<"+
"UQ::PCE"+
"> Min All", out);
580 typedef Sacado::mpl::apply<FadType,PCEType>::type FadPCEType;
581 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
582 comm = Teuchos::DefaultComm<Ordinal>::getComm();
587 int rank = comm->getRank();
589 Teuchos::Array<FadPCEType>
x(n), sums(n), sums2(n);
590 for (
int i=0; i<n; i++) {
592 for (
int k=0; k<
setup.sz; k++)
593 f.fastAccessCoeff(k) = 2.0*(i+1);
594 x[i] = FadPCEType(p,
f);
595 for (
int j=0;
j<p;
j++) {
596 x[i].fastAccessDx(
j) =
f;
599 for (
int i=0; i<n; i++) {
601 for (
int k=0; k<
setup.sz; k++)
602 f.fastAccessCoeff(k) = 2.0*(i+1)*(rank+1);
603 sums[i] = FadPCEType(p,
f);
604 for (
int j=0;
j<p;
j++)
607 Teuchos::scan(*comm, *
setup.fad_pce_serializer,
608 Teuchos::REDUCE_SUM, n, &
x[0], &sums2[0]);
610 std::string(
"DFad")+
"<"+
"UQ::PCE"+
"> Scan Sum", out);
615 typedef Sacado::mpl::apply<FadType,PCEType>::type FadPCEType;
616 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
617 comm = Teuchos::DefaultComm<Ordinal>::getComm();
622 int rank = comm->getRank();
624 Teuchos::Array<FadPCEType>
x(n), maxs(n), maxs2(n);
625 for (
int i=0; i<n; i++) {
627 for (
int k=0; k<
setup.sz; k++)
628 f.fastAccessCoeff(k) = 2.0*(i+1)*(rank+1);
629 x[i] = FadPCEType(p,
f);
630 for (
int j=0;
j<p;
j++) {
631 x[i].fastAccessDx(
j) =
f;
634 for (
int i=0; i<n; i++) {
636 for (
int k=0; k<
setup.sz; k++)
637 f.fastAccessCoeff(k) = 2.0*(i+1)*(rank+1);
638 maxs[i] = FadPCEType(p,
f);
639 for (
int j=0;
j<p;
j++)
642 Teuchos::scan(*comm, *
setup.fad_pce_serializer,
643 Teuchos::REDUCE_MAX, n, &
x[0], &maxs2[0]);
645 std::string(
"DFad")+
"<"+
"UQ::PCE"+
"> Scan Max", out);
650 typedef Sacado::mpl::apply<FadType,PCEType>::type FadPCEType;
651 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
652 comm = Teuchos::DefaultComm<Ordinal>::getComm();
657 int rank = comm->getRank();
659 Teuchos::Array<FadPCEType>
x(n), mins(n), mins2(n);
660 for (
int i=0; i<n; i++) {
662 for (
int k=0; k<
setup.sz; k++)
663 f.fastAccessCoeff(k) = 2.0*(i+1)*(rank+1);
664 x[i] = FadPCEType(p,
f);
665 for (
int j=0;
j<p;
j++) {
666 x[i].fastAccessDx(
j) =
f;
669 for (
int i=0; i<n; i++) {
671 for (
int k=0; k<
setup.sz; k++)
672 f.fastAccessCoeff(k) = 2.0*(i+1);
673 mins[i] = FadPCEType(p,
f);
674 for (
int j=0;
j<p;
j++)
677 Teuchos::scan(*comm, *
setup.fad_pce_serializer,
678 Teuchos::REDUCE_MIN, n, &
x[0], &mins2[0]);
680 std::string(
"DFad")+
"<"+
"UQ::PCE"+
"> Scan Min", out);
685 typedef Sacado::mpl::apply<FadType,PCEType>::type FadPCEType;
686 Teuchos::RCP<const Teuchos::Comm<Ordinal> >
687 comm = Teuchos::DefaultComm<Ordinal>::getComm();
690 int num_proc = comm->getSize();
692 int rank = comm->getRank();
695 Teuchos::Array<FadPCEType>
x(n), x2(n);
696 for (
int i=0; i<n; i++) {
698 for (
int k=0; k<
setup.sz; k++)
699 f.fastAccessCoeff(k) = 2.0*(i+1)*(k+1);
700 x[i] = FadPCEType(p,
f);
701 for (
int j=0;
j<p;
j++)
706 if (rank == 0) Teuchos::send(*comm, *
setup.fad_pce_serializer,
708 if (rank == 1) Teuchos::receive(*comm, *
setup.fad_pce_serializer,
711 std::string(
"DFad")+
"<"+
"UQ::PCE"+
"> Send/Receive", out);
719 Teuchos::GlobalMPISession mpiSession(&argc, &
argv);
721 Kokkos::HostSpace::execution_space::initialize();
722 if (!Kokkos::DefaultExecutionSpace::is_initialized())
723 Kokkos::DefaultExecutionSpace::initialize();
725 int res = Teuchos::UnitTestRepository::runUnitTestsFromMain(argc,
argv);
727 Kokkos::HostSpace::execution_space::finalize();
728 if (Kokkos::DefaultExecutionSpace::is_initialized())
729 Kokkos::DefaultExecutionSpace::finalize();
Kokkos::DefaultExecutionSpace execution_space
RCP< FadPCESerializerT > fad_pce_serializer
bool checkPCEArrays(const Teuchos::Array< PCEType > &x, const Teuchos::Array< PCEType > &x2, const std::string &tag, Teuchos::FancyOStream &out)
RCP< const Stokhos::CompletePolynomialBasis< int, double > > basis
RCP< PCESerializerT > pce_serializer
Teuchos::ValueTypeSerializer< int, FadPCEType > FadPCESerializerT
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
TEUCHOS_UNIT_TEST(UQ_PCE_Comm, PCE_Broadcast)
UnitTestSetup< int, double > setup
Sacado::Fad::DFad< double > FadType
expr expr expr expr fastAccessDx(i, j)) FAD_UNARYOP_MACRO(exp
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Sacado::mpl::apply< FadType, PCEType >::type FadPCEType
Legendre polynomial basis.
Sacado::UQ::PCE< storage_type > PCEType
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
int main(int argc, char *argv[])
bool checkResultOnAllProcs(const Teuchos::Comm< Ordinal > &comm, Teuchos::FancyOStream &out, const bool result)
Stokhos::DynamicStorage< int, double, execution_space > storage_type
Sacado::Random< double > rnd
PCEType::cijk_type kokkos_cijk_type
kokkos_cijk_type kokkos_cijk
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
RCP< Stokhos::Sparse3Tensor< int, double > > Cijk
Teuchos::ValueTypeSerializer< int, PCEType > PCESerializerT