BenchmarkFieldAlgorithms.cxx 29.2 KB
Newer Older
1 2 3 4 5 6 7 8
//============================================================================
//  Copyright (c) Kitware, Inc.
//  All rights reserved.
//  See LICENSE.txt for details.
//  This software is distributed WITHOUT ANY WARRANTY; without even
//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
//  PURPOSE.  See the above copyright notice for more information.
//
9
//  Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
10 11 12
//  Copyright 2014 UT-Battelle, LLC.
//  Copyright 2014 Los Alamos National Security.
//
13
//  Under the terms of Contract DE-NA0003525 with NTESS,
14 15 16 17 18 19 20 21 22 23 24 25
//  the U.S. Government retains certain rights in this software.
//
//  Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
//  Laboratory (LANL), the U.S. Government retains certain rights in
//  this software.
//============================================================================

#include <vtkm/Math.h>
#include <vtkm/VectorAnalysis.h>

#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/CellSetStructured.h>
26
#include <vtkm/cont/DynamicArrayHandle.h>
27
#include <vtkm/cont/ImplicitFunctionHandle.h>
28 29 30 31
#include <vtkm/cont/Timer.h>

#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/DispatcherMapTopology.h>
32 33
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/WorkletMapTopology.h>
34

35
#include "Benchmarker.h"
36
#include <vtkm/cont/testing/Testing.h>
37

38
#include <cctype>
39
#include <random>
40 41
#include <string>

42 43 44 45
namespace vtkm
{
namespace benchmarking
{
46 47 48

#define ARRAY_SIZE (1 << 22)
#define CUBE_SIZE 256
49
static const std::string DIVIDER(40, '-');
50

51 52
enum BenchmarkName
{
53 54
  BLACK_SCHOLES = 1,
  MATH = 1 << 1,
55 56 57
  FUSED_MATH = 1 << 2,
  INTERPOLATE_FIELD = 1 << 3,
  IMPLICIT_FUNCTION = 1 << 4,
58
  ALL = BLACK_SCHOLES | MATH | FUSED_MATH | INTERPOLATE_FIELD | IMPLICIT_FUNCTION
59 60
};

61
template <typename T>
62 63 64 65
class BlackScholes : public vtkm::worklet::WorkletMapField
{
  T Riskfree;
  T Volatility;
66

67
public:
68 69 70
  using ControlSignature =
    void(FieldIn<Scalar>, FieldIn<Scalar>, FieldIn<Scalar>, FieldOut<Scalar>, FieldOut<Scalar>);
  using ExecutionSignature = void(_1, _2, _3, _4, _5);
71

72 73 74
  BlackScholes(T risk, T volatility)
    : Riskfree(risk)
    , Volatility(volatility)
75 76 77
  {
  }

78
  VTKM_EXEC
79 80
  T CumulativeNormalDistribution(T d) const
  {
81 82 83 84 85 86
    const vtkm::Float32 A1 = 0.31938153f;
    const vtkm::Float32 A2 = -0.356563782f;
    const vtkm::Float32 A3 = 1.781477937f;
    const vtkm::Float32 A4 = -1.821255978f;
    const vtkm::Float32 A5 = 1.330274429f;
    const vtkm::Float32 RSQRT2PI = 0.39894228040143267793994605993438f;
87

88
    const vtkm::Float32 df = static_cast<vtkm::Float32>(d);
89
    const vtkm::Float32 K = 1.0f / (1.0f + 0.2316419f * vtkm::Abs(df));
90

91 92
    vtkm::Float32 cnd =
      RSQRT2PI * vtkm::Exp(-0.5f * df * df) * (K * (A1 + K * (A2 + K * (A3 + K * (A4 + K * A5)))));
93

94 95
    if (df > 0.0f)
    {
96
      cnd = 1.0f - cnd;
97
    }
98

99
    return static_cast<T>(cnd);
100 101
  }

102
  template <typename U, typename V, typename W>
103 104
  VTKM_EXEC void operator()(const U& sp, const V& os, const W& oy, T& callResult, T& putResult)
    const
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
  {
    const T stockPrice = static_cast<T>(sp);
    const T optionStrike = static_cast<T>(os);
    const T optionYears = static_cast<T>(oy);

    // Black-Scholes formula for both call and put
    const T sqrtYears = vtkm::Sqrt(optionYears);
    const T volMultSqY = this->Volatility * sqrtYears;

    const T d1 = (vtkm::Log(stockPrice / optionStrike) +
                  (this->Riskfree + 0.5f * Volatility * Volatility) * optionYears) /
      (volMultSqY);
    const T d2 = d1 - volMultSqY;
    const T CNDD1 = CumulativeNormalDistribution(d1);
    const T CNDD2 = CumulativeNormalDistribution(d2);

    //Calculate Call and Put simultaneously
    T expRT = vtkm::Exp(-this->Riskfree * optionYears);
    callResult = stockPrice * CNDD1 - optionStrike * expRT * CNDD2;
    putResult = optionStrike * expRT * (1.0f - CNDD2) - stockPrice * (1.0f - CNDD1);
125 126 127 128 129 130
  }
};

class Mag : public vtkm::worklet::WorkletMapField
{
public:
131 132
  using ControlSignature = void(FieldIn<Vec3>, FieldOut<Scalar>);
  using ExecutionSignature = void(_1, _2);
133

134 135
  template <typename T, typename U>
  VTKM_EXEC void operator()(const vtkm::Vec<T, 3>& vec, U& result) const
136
  {
137
    result = static_cast<U>(vtkm::Magnitude(vec));
138 139 140 141 142 143
  }
};

class Square : public vtkm::worklet::WorkletMapField
{
public:
144 145
  using ControlSignature = void(FieldIn<Scalar>, FieldOut<Scalar>);
  using ExecutionSignature = void(_1, _2);
146

147 148
  template <typename T, typename U>
  VTKM_EXEC void operator()(T input, U& output) const
149
  {
150
    output = static_cast<U>(input * input);
151 152 153 154 155 156
  }
};

class Sin : public vtkm::worklet::WorkletMapField
{
public:
157 158
  using ControlSignature = void(FieldIn<Scalar>, FieldOut<Scalar>);
  using ExecutionSignature = void(_1, _2);
159

160 161
  template <typename T, typename U>
  VTKM_EXEC void operator()(T input, U& output) const
162
  {
163
    output = static_cast<U>(vtkm::Sin(input));
164 165 166 167 168 169
  }
};

class Cos : public vtkm::worklet::WorkletMapField
{
public:
170 171
  using ControlSignature = void(FieldIn<Scalar>, FieldOut<Scalar>);
  using ExecutionSignature = void(_1, _2);
172

173 174
  template <typename T, typename U>
  VTKM_EXEC void operator()(T input, U& output) const
175
  {
176
    output = static_cast<U>(vtkm::Cos(input));
177 178 179 180 181 182
  }
};

class FusedMath : public vtkm::worklet::WorkletMapField
{
public:
183 184
  using ControlSignature = void(FieldIn<Vec3>, FieldOut<Scalar>);
  using ExecutionSignature = void(_1, _2);
185

186 187
  template <typename T>
  VTKM_EXEC void operator()(const vtkm::Vec<T, 3>& vec, T& result) const
188 189
  {
    const T m = vtkm::Magnitude(vec);
190
    result = vtkm::Cos(vtkm::Sin(m) * vtkm::Sin(m));
191
  }
192

193 194
  template <typename T, typename U>
  VTKM_EXEC void operator()(const vtkm::Vec<T, 3>&, U&) const
195
  {
196
    this->RaiseError("Mixed types unsupported.");
197
  }
198 199 200 201 202
};

class GenerateEdges : public vtkm::worklet::WorkletMapPointToCell
{
public:
203 204
  using ControlSignature = void(CellSetIn cellset, WholeArrayOut<> edgeIds);
  using ExecutionSignature = void(PointIndices, ThreadIndices, _2);
205
  using InputDomain = _1;
206

207 208 209 210
  template <typename ConnectivityInVec, typename ThreadIndicesType, typename IdPairTableType>
  VTKM_EXEC void operator()(const ConnectivityInVec& connectivity,
                            const ThreadIndicesType threadIndices,
                            const IdPairTableType& edgeIds) const
211
  {
212
    const vtkm::Id writeOffset = (threadIndices.GetInputIndex() * 12);
213

214 215
    const vtkm::IdComponent edgeTable[24] = { 0, 1, 1, 2, 3, 2, 0, 3, 4, 5, 5, 6,
                                              7, 6, 4, 7, 0, 4, 1, 5, 2, 6, 3, 7 };
216

217
    for (vtkm::Id i = 0; i < 12; ++i)
218
    {
219 220 221
      const vtkm::Id offset = (i * 2);
      const vtkm::Id2 edge(connectivity[edgeTable[offset]], connectivity[edgeTable[offset + 1]]);
      edgeIds.Set(writeOffset + i, edge);
222 223 224 225 226 227 228
    }
  }
};

class InterpolateField : public vtkm::worklet::WorkletMapField
{
public:
229
  using ControlSignature = void(FieldIn<Id2Type> interpolation_ids,
230 231
                                FieldIn<Scalar> interpolation_weights,
                                WholeArrayIn<> inputField,
232
                                FieldOut<> output);
233
  using ExecutionSignature = void(_1, _2, _3, _4);
234
  using InputDomain = _1;
235

236
  template <typename WeightType, typename T, typename S, typename D>
237 238
  VTKM_EXEC void operator()(const vtkm::Id2& low_high,
                            const WeightType& weight,
239 240
                            const vtkm::exec::ExecutionWholeArrayConst<T, S, D>& inPortal,
                            T& result) const
241 242
  {
    //fetch the low / high values from inPortal
243
    result = vtkm::Lerp(inPortal.Get(low_high[0]), inPortal.Get(low_high[1]), weight);
244
  }
245 246

  template <typename WeightType, typename T, typename S, typename D, typename U>
247 248 249 250
  VTKM_EXEC void operator()(const vtkm::Id2&,
                            const WeightType&,
                            const vtkm::exec::ExecutionWholeArrayConst<T, S, D>&,
                            U&) const
251 252 253
  {
    //the inPortal and result need to be the same type so this version only
    //exists to generate code when using dynamic arrays
254
    this->RaiseError("Mixed types unsupported.");
255
  }
256 257
};

258 259 260 261
template <typename ImplicitFunction>
class EvaluateImplicitFunction : public vtkm::worklet::WorkletMapField
{
public:
262 263
  using ControlSignature = void(FieldIn<Vec3>, FieldOut<Scalar>);
  using ExecutionSignature = void(_1, _2);
264

265
  EvaluateImplicitFunction(const ImplicitFunction* function)
266
    : Function(function)
267 268
  {
  }
269

270 271
  template <typename VecType, typename ScalarType>
  VTKM_EXEC void operator()(const VecType& point, ScalarType& val) const
272
  {
273
    val = this->Function->Value(point);
274 275 276
  }

private:
277
  const ImplicitFunction* Function;
278 279 280 281 282 283
};

template <typename T1, typename T2>
class Evaluate2ImplicitFunctions : public vtkm::worklet::WorkletMapField
{
public:
284 285
  using ControlSignature = void(FieldIn<Vec3>, FieldOut<Scalar>);
  using ExecutionSignature = void(_1, _2);
286

287
  Evaluate2ImplicitFunctions(const T1* f1, const T2* f2)
288 289 290 291
    : Function1(f1)
    , Function2(f2)
  {
  }
292

293 294
  template <typename VecType, typename ScalarType>
  VTKM_EXEC void operator()(const VecType& point, ScalarType& val) const
295
  {
296
    val = this->Function1->Value(point) + this->Function2->Value(point);
297 298 299
  }

private:
300 301
  const T1* Function1;
  const T2* Function2;
302 303
};

304 305 306
struct ValueTypes : vtkm::ListTagBase<vtkm::Float32, vtkm::Float64>
{
};
307

308 309 310 311
struct InterpValueTypes : vtkm::ListTagBase<vtkm::Float32,
                                            vtkm::Float64,
                                            vtkm::Vec<vtkm::Float32, 3>,
                                            vtkm::Vec<vtkm::Float64, 3>>
312 313
{
};
314
using StorageListTag = ::vtkm::cont::StorageListTagBasic;
315 316 317

/// This class runs a series of micro-benchmarks to measure
/// performance of different field operations
318 319 320
template <class DeviceAdapterTag>
class BenchmarkFieldAlgorithms
{
321
  using StorageTag = vtkm::cont::StorageTagBasic;
322

323
  using Algorithm = vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapterTag>;
324

325
  using Timer = vtkm::cont::Timer<DeviceAdapterTag>;
326

327 328
  using ValueDynamicHandle = vtkm::cont::DynamicArrayHandleBase<ValueTypes, StorageListTag>;
  using InterpDynamicHandle = vtkm::cont::DynamicArrayHandleBase<InterpValueTypes, StorageListTag>;
329
  using IdDynamicHandle =
330
    vtkm::cont::DynamicArrayHandleBase<vtkm::TypeListTagIndex, StorageListTag>;
331

332
private:
333 334 335
  template <typename Value>
  struct BenchBlackScholes
  {
336
    using ValueArrayHandle = vtkm::cont::ArrayHandle<Value, StorageTag>;
337 338 339 340 341 342 343 344 345

    ValueArrayHandle StockPrice;
    ValueArrayHandle OptionStrike;
    ValueArrayHandle OptionYears;

    std::vector<Value> price;
    std::vector<Value> strike;
    std::vector<Value> years;

346
    VTKM_CONT
347 348
    BenchBlackScholes()
    {
349
      std::mt19937 rng;
350 351 352
      std::uniform_real_distribution<Value> price_range(Value(5.0f), Value(30.0f));
      std::uniform_real_distribution<Value> strike_range(Value(1.0f), Value(100.0f));
      std::uniform_real_distribution<Value> year_range(Value(0.25f), Value(10.0f));
353 354 355 356

      this->price.resize(ARRAY_SIZE);
      this->strike.resize(ARRAY_SIZE);
      this->years.resize(ARRAY_SIZE);
357
      for (std::size_t i = 0; i < ARRAY_SIZE; ++i)
358
      {
359 360 361
        this->price[i] = price_range(rng);
        this->strike[i] = strike_range(rng);
        this->years[i] = year_range(rng);
362 363 364 365 366 367 368
      }

      this->StockPrice = vtkm::cont::make_ArrayHandle(this->price);
      this->OptionStrike = vtkm::cont::make_ArrayHandle(this->strike);
      this->OptionYears = vtkm::cont::make_ArrayHandle(this->years);
    }

369
    VTKM_CONT
370 371 372 373 374 375 376 377
    vtkm::Float64 operator()()
    {
      vtkm::cont::ArrayHandle<Value> callResultHandle, putResultHandle;
      const Value RISKFREE = 0.02f;
      const Value VOLATILITY = 0.30f;

      Timer timer;
      BlackScholes<Value> worklet(RISKFREE, VOLATILITY);
378
      vtkm::worklet::DispatcherMapField<BlackScholes<Value>> dispatcher(worklet);
379

380 381
      dispatcher.Invoke(
        this->StockPrice, this->OptionStrike, this->OptionYears, callResultHandle, putResultHandle);
382 383 384 385

      return timer.GetElapsedTime();
    }

386 387
    virtual std::string Type() const { return std::string("Static"); }

388
    VTKM_CONT
389 390
    std::string Description() const
    {
391
      std::stringstream description;
392 393 394
      description << "BlackScholes "
                  << "[" << this->Type() << "] "
                  << " with a domain size of: " << ARRAY_SIZE;
395 396 397 398
      return description.str();
    }
  };

399 400 401
  template <typename Value>
  struct BenchBlackScholesDynamic : public BenchBlackScholes<Value>
  {
402

403
    VTKM_CONT
404 405 406 407 408 409 410 411 412 413 414 415
    vtkm::Float64 operator()()
    {
      ValueDynamicHandle dstocks(this->StockPrice);
      ValueDynamicHandle dstrikes(this->OptionStrike);
      ValueDynamicHandle doptions(this->OptionYears);

      vtkm::cont::ArrayHandle<Value> callResultHandle, putResultHandle;
      const Value RISKFREE = 0.02f;
      const Value VOLATILITY = 0.30f;

      Timer timer;
      BlackScholes<Value> worklet(RISKFREE, VOLATILITY);
416
      vtkm::worklet::DispatcherMapField<BlackScholes<Value>> dispatcher(worklet);
417

418
      dispatcher.Invoke(dstocks, dstrikes, doptions, callResultHandle, putResultHandle);
419 420 421 422 423 424 425

      return timer.GetElapsedTime();
    }

    virtual std::string Type() const { return std::string("Dynamic"); }
  };

426
  VTKM_MAKE_BENCHMARK(BlackScholes, BenchBlackScholes);
427 428
  VTKM_MAKE_BENCHMARK(BlackScholesDynamic, BenchBlackScholesDynamic);

429 430 431 432 433
  template <typename Value>
  struct BenchMath
  {
    std::vector<vtkm::Vec<Value, 3>> input;
    vtkm::cont::ArrayHandle<vtkm::Vec<Value, 3>, StorageTag> InputHandle;
434

435
    VTKM_CONT
436 437
    BenchMath()
    {
438 439
      std::mt19937 rng;
      std::uniform_real_distribution<Value> range;
440 441

      this->input.resize(ARRAY_SIZE);
442
      for (std::size_t i = 0; i < ARRAY_SIZE; ++i)
443
      {
444
        this->input[i] = vtkm::Vec<Value, 3>(range(rng), range(rng), range(rng));
445 446 447 448 449
      }

      this->InputHandle = vtkm::cont::make_ArrayHandle(this->input);
    }

450
    VTKM_CONT
451 452 453 454 455 456 457
    vtkm::Float64 operator()()
    {
      vtkm::cont::ArrayHandle<Value> tempHandle1;
      vtkm::cont::ArrayHandle<Value> tempHandle2;

      Timer timer;

458 459 460 461
      vtkm::worklet::DispatcherMapField<Mag>().Invoke(InputHandle, tempHandle1);
      vtkm::worklet::DispatcherMapField<Sin>().Invoke(tempHandle1, tempHandle2);
      vtkm::worklet::DispatcherMapField<Square>().Invoke(tempHandle2, tempHandle1);
      vtkm::worklet::DispatcherMapField<Cos>().Invoke(tempHandle1, tempHandle2);
462 463 464 465

      return timer.GetElapsedTime();
    }

466 467
    virtual std::string Type() const { return std::string("Static"); }

468
    VTKM_CONT
469 470
    std::string Description() const
    {
471
      std::stringstream description;
472 473 474
      description << "Magnitude -> Sine -> Square -> Cosine "
                  << "[" << this->Type() << "] "
                  << "with a domain size of: " << ARRAY_SIZE;
475 476 477 478
      return description.str();
    }
  };

479 480 481
  template <typename Value>
  struct BenchMathDynamic : public BenchMath<Value>
  {
482

483
    VTKM_CONT
484 485
    vtkm::Float64 operator()()
    {
486
      using MathTypes = vtkm::ListTagBase<vtkm::Vec<vtkm::Float32, 3>, vtkm::Vec<vtkm::Float64, 3>>;
487 488 489 490 491 492 493 494 495

      vtkm::cont::ArrayHandle<Value> temp1;
      vtkm::cont::ArrayHandle<Value> temp2;
      vtkm::cont::DynamicArrayHandleBase<MathTypes, StorageListTag> dinput(this->InputHandle);
      ValueDynamicHandle dtemp1(temp1);
      ValueDynamicHandle dtemp2(temp2);

      Timer timer;

496 497 498 499
      vtkm::worklet::DispatcherMapField<Mag>().Invoke(dinput, dtemp1);
      vtkm::worklet::DispatcherMapField<Sin>().Invoke(dtemp1, dtemp2);
      vtkm::worklet::DispatcherMapField<Square>().Invoke(dtemp2, dtemp1);
      vtkm::worklet::DispatcherMapField<Cos>().Invoke(dtemp1, dtemp2);
500 501 502 503 504 505 506

      return timer.GetElapsedTime();
    }

    virtual std::string Type() const { return std::string("Dynamic"); }
  };

507
  VTKM_MAKE_BENCHMARK(Math, BenchMath);
508
  VTKM_MAKE_BENCHMARK(MathDynamic, BenchMathDynamic);
509

510 511 512 513 514
  template <typename Value>
  struct BenchFusedMath
  {
    std::vector<vtkm::Vec<Value, 3>> input;
    vtkm::cont::ArrayHandle<vtkm::Vec<Value, 3>, StorageTag> InputHandle;
515

516
    VTKM_CONT
517 518
    BenchFusedMath()
    {
519 520
      std::mt19937 rng;
      std::uniform_real_distribution<Value> range;
521 522

      this->input.resize(ARRAY_SIZE);
523
      for (std::size_t i = 0; i < ARRAY_SIZE; ++i)
524
      {
525
        this->input[i] = vtkm::Vec<Value, 3>(range(rng), range(rng), range(rng));
526 527 528 529 530
      }

      this->InputHandle = vtkm::cont::make_ArrayHandle(this->input);
    }

531
    VTKM_CONT
532 533 534 535 536
    vtkm::Float64 operator()()
    {
      vtkm::cont::ArrayHandle<Value> result;

      Timer timer;
537
      vtkm::worklet::DispatcherMapField<FusedMath>().Invoke(this->InputHandle, result);
538 539 540
      return timer.GetElapsedTime();
    }

541 542
    virtual std::string Type() const { return std::string("Static"); }

543
    VTKM_CONT
544 545
    std::string Description() const
    {
546
      std::stringstream description;
547 548 549
      description << "Fused Magnitude -> Sine -> Square -> Cosine "
                  << "[" << this->Type() << "] "
                  << "with a domain size of: " << ARRAY_SIZE;
550 551 552 553
      return description.str();
    }
  };

554 555 556
  template <typename Value>
  struct BenchFusedMathDynamic : public BenchFusedMath<Value>
  {
557

558
    VTKM_CONT
559 560
    vtkm::Float64 operator()()
    {
561
      using MathTypes = vtkm::ListTagBase<vtkm::Vec<vtkm::Float32, 3>, vtkm::Vec<vtkm::Float64, 3>>;
562 563 564 565 566 567

      vtkm::cont::DynamicArrayHandleBase<MathTypes, StorageListTag> dinput(this->InputHandle);

      vtkm::cont::ArrayHandle<Value, StorageTag> result;

      Timer timer;
568
      vtkm::worklet::DispatcherMapField<FusedMath>().Invoke(dinput, result);
569 570 571 572 573 574
      return timer.GetElapsedTime();
    }

    virtual std::string Type() const { return std::string("Dynamic"); }
  };

575
  VTKM_MAKE_BENCHMARK(FusedMath, BenchFusedMath);
576
  VTKM_MAKE_BENCHMARK(FusedMathDynamic, BenchFusedMathDynamic);
577

578 579 580
  template <typename Value>
  struct BenchEdgeInterp
  {
581 582 583
    std::vector<vtkm::Float32> weight;
    std::vector<Value> field;

584
    vtkm::cont::ArrayHandle<vtkm::Float32, StorageTag> WeightHandle;
585
    vtkm::cont::ArrayHandle<Value, StorageTag> FieldHandle;
586
    vtkm::cont::ArrayHandle<vtkm::Id2, StorageTag> EdgePairHandle;
587

588
    VTKM_CONT
589 590
    BenchEdgeInterp()
    {
591
      using CT = typename vtkm::VecTraits<Value>::ComponentType;
592
      std::mt19937 rng;
593
      std::uniform_real_distribution<vtkm::Float32> weight_range(0.0f, 1.0f);
594
      std::uniform_real_distribution<CT> field_range;
595 596 597 598 599 600 601 602

      //basically the core challenge is to generate an array whose
      //indexing pattern matches that of a edge based algorithm.
      //
      //So for this kind of problem we generate the 12 edges of each
      //cell and place them into array.
      //
      vtkm::cont::CellSetStructured<3> cellSet;
603
      cellSet.SetPointDimensions(vtkm::Id3(CUBE_SIZE, CUBE_SIZE, CUBE_SIZE));
604 605 606 607 608

      const vtkm::Id numberOfEdges = cellSet.GetNumberOfCells() * 12;
      const std::size_t esize = static_cast<std::size_t>(numberOfEdges);
      const std::size_t psize = static_cast<std::size_t>(cellSet.GetNumberOfPoints());

609 610
      this->EdgePairHandle.Allocate(numberOfEdges);
      vtkm::worklet::DispatcherMapTopology<GenerateEdges>().Invoke(cellSet, this->EdgePairHandle);
611

612 613
      this->weight.resize(esize);
      for (std::size_t i = 0; i < esize; ++i)
614
      {
615
        this->weight[i] = weight_range(rng);
616 617
      }

618 619
      this->field.resize(psize);
      for (std::size_t i = 0; i < psize; ++i)
620
      {
621
        this->field[i] = Value(field_range(rng));
622 623 624 625 626 627
      }

      this->FieldHandle = vtkm::cont::make_ArrayHandle(this->field);
      this->WeightHandle = vtkm::cont::make_ArrayHandle(this->weight);
    }

628
    VTKM_CONT
629 630 631 632 633
    vtkm::Float64 operator()()
    {
      vtkm::cont::ArrayHandle<Value> result;

      Timer timer;
634 635
      vtkm::worklet::DispatcherMapField<InterpolateField, DeviceAdapterTag>().Invoke(
        this->EdgePairHandle, this->WeightHandle, this->FieldHandle, result);
636 637 638
      return timer.GetElapsedTime();
    }

639 640
    virtual std::string Type() const { return std::string("Static"); }

641
    VTKM_CONT
642 643
    std::string Description() const
    {
644
      std::stringstream description;
645
      const std::size_t size = (CUBE_SIZE - 1) * (CUBE_SIZE - 1) * (CUBE_SIZE - 1) * 12;
646 647 648
      description << "Edge Interpolation "
                  << "[" << this->Type() << "] "
                  << "with a domain size of: " << size;
649 650 651 652
      return description.str();
    }
  };

653 654 655
  template <typename Value>
  struct BenchEdgeInterpDynamic : public BenchEdgeInterp<Value>
  {
656

657
    VTKM_CONT
658 659 660 661 662 663
    vtkm::Float64 operator()()
    {
      InterpDynamicHandle dfield(this->FieldHandle);
      InterpDynamicHandle dweight(this->WeightHandle);
      IdDynamicHandle dedges(this->EdgePairHandle);
      vtkm::cont::ArrayHandle<Value> result;
664

665
      Timer timer;
666 667
      vtkm::worklet::DispatcherMapField<InterpolateField, DeviceAdapterTag>().Invoke(
        dedges, dweight, dfield, result);
668 669
      return timer.GetElapsedTime();
    }
670

671 672
    virtual std::string Type() const { return std::string("Dynamic"); }
  };
673

674 675
  VTKM_MAKE_BENCHMARK(EdgeInterp, BenchEdgeInterp);
  VTKM_MAKE_BENCHMARK(EdgeInterpDynamic, BenchEdgeInterpDynamic);
676

677 678 679 680
  struct ImplicitFunctionBenchData
  {
    vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> Points;
    vtkm::cont::ArrayHandle<vtkm::FloatDefault> Result;
681
    vtkm::Sphere Sphere1, Sphere2;
682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703
  };

  static ImplicitFunctionBenchData MakeImplicitFunctionBenchData()
  {
    vtkm::Id count = ARRAY_SIZE;
    vtkm::FloatDefault bounds[6] = { -2.0f, 2.0f, -2.0f, 2.0f, -2.0f, 2.0f };

    ImplicitFunctionBenchData data;
    data.Points.Allocate(count);
    data.Result.Allocate(count);

    std::default_random_engine rangen;
    std::uniform_real_distribution<vtkm::FloatDefault> distx(bounds[0], bounds[1]);
    std::uniform_real_distribution<vtkm::FloatDefault> disty(bounds[2], bounds[3]);
    std::uniform_real_distribution<vtkm::FloatDefault> distz(bounds[4], bounds[5]);

    auto portal = data.Points.GetPortalControl();
    for (vtkm::Id i = 0; i < count; ++i)
    {
      portal.Set(i, vtkm::make_Vec(distx(rangen), disty(rangen), distz(rangen)));
    }

704 705
    data.Sphere1 = vtkm::Sphere({ 0.22f, 0.33f, 0.44f }, 0.55f);
    data.Sphere2 = vtkm::Sphere({ 0.22f, 0.33f, 0.11f }, 0.77f);
706 707 708 709

    return data;
  }

710
  template <typename Value>
711 712 713 714
  struct BenchImplicitFunction
  {
    BenchImplicitFunction()
      : Internal(MakeImplicitFunctionBenchData())
715 716
    {
    }
717 718 719 720

    VTKM_CONT
    vtkm::Float64 operator()()
    {
721
      using EvalWorklet = EvaluateImplicitFunction<vtkm::Sphere>;
722 723
      using EvalDispatcher = vtkm::worklet::DispatcherMapField<EvalWorklet, DeviceAdapterTag>;

724 725 726 727
      auto handle = vtkm::cont::make_ImplicitFunctionHandle(Internal.Sphere1);
      auto function =
        static_cast<const vtkm::Sphere*>(handle.PrepareForExecution(DeviceAdapterTag()));
      EvalWorklet eval(function);
728 729 730 731 732 733 734 735 736 737

      vtkm::cont::Timer<DeviceAdapterTag> timer;
      EvalDispatcher(eval).Invoke(this->Internal.Points, this->Internal.Result);
      return timer.GetElapsedTime();
    }

    VTKM_CONT
    std::string Description() const
    {
      std::stringstream description;
738 739
      description << "Implicit Function (vtkm::Sphere) on " << Internal.Points.GetNumberOfValues()
                  << " points";
740 741 742 743 744 745
      return description.str();
    }

    ImplicitFunctionBenchData Internal;
  };

746
  template <typename Value>
747
  struct BenchVirtualImplicitFunction
748
  {
749
    BenchVirtualImplicitFunction()
750
      : Internal(MakeImplicitFunctionBenchData())
751 752
    {
    }
753 754 755 756

    VTKM_CONT
    vtkm::Float64 operator()()
    {
757
      using EvalWorklet = EvaluateImplicitFunction<vtkm::ImplicitFunction>;
758 759
      using EvalDispatcher = vtkm::worklet::DispatcherMapField<EvalWorklet, DeviceAdapterTag>;

760 761
      auto sphere = vtkm::cont::make_ImplicitFunctionHandle(Internal.Sphere1);
      EvalWorklet eval(sphere.PrepareForExecution(DeviceAdapterTag()));
762 763 764 765 766 767 768 769 770 771

      vtkm::cont::Timer<DeviceAdapterTag> timer;
      EvalDispatcher(eval).Invoke(this->Internal.Points, this->Internal.Result);
      return timer.GetElapsedTime();
    }

    VTKM_CONT
    std::string Description() const
    {
      std::stringstream description;
772
      description << "Implicit Function (VirtualImplicitFunction) on "
773
                  << Internal.Points.GetNumberOfValues() << " points";
774 775 776 777 778 779
      return description.str();
    }

    ImplicitFunctionBenchData Internal;
  };

780
  template <typename Value>
781 782 783 784
  struct Bench2ImplicitFunctions
  {
    Bench2ImplicitFunctions()
      : Internal(MakeImplicitFunctionBenchData())
785 786
    {
    }
787 788 789 790

    VTKM_CONT
    vtkm::Float64 operator()()
    {
791
      using EvalWorklet = Evaluate2ImplicitFunctions<vtkm::Sphere, vtkm::Sphere>;
792 793
      using EvalDispatcher = vtkm::worklet::DispatcherMapField<EvalWorklet, DeviceAdapterTag>;

794 795 796 797 798
      auto h1 = vtkm::cont::make_ImplicitFunctionHandle(Internal.Sphere1);
      auto h2 = vtkm::cont::make_ImplicitFunctionHandle(Internal.Sphere2);
      auto f1 = static_cast<const vtkm::Sphere*>(h1.PrepareForExecution(DeviceAdapterTag()));
      auto f2 = static_cast<const vtkm::Sphere*>(h2.PrepareForExecution(DeviceAdapterTag()));
      EvalWorklet eval(f1, f2);
799 800 801 802 803 804 805 806 807 808

      vtkm::cont::Timer<DeviceAdapterTag> timer;
      EvalDispatcher(eval).Invoke(this->Internal.Points, this->Internal.Result);
      return timer.GetElapsedTime();
    }

    VTKM_CONT
    std::string Description() const
    {
      std::stringstream description;
809 810
      description << "Implicit Function 2x(vtkm::Sphere) on " << Internal.Points.GetNumberOfValues()
                  << " points";
811 812 813 814 815 816
      return description.str();
    }

    ImplicitFunctionBenchData Internal;
  };

817
  template <typename Value>
818
  struct Bench2VirtualImplicitFunctions
819
  {
820
    Bench2VirtualImplicitFunctions()
821
      : Internal(MakeImplicitFunctionBenchData())
822 823
    {
    }
824 825 826 827

    VTKM_CONT
    vtkm::Float64 operator()()
    {
828
      using EvalWorklet =
829
        Evaluate2ImplicitFunctions<vtkm::ImplicitFunction, vtkm::ImplicitFunction>;
830 831
      using EvalDispatcher = vtkm::worklet::DispatcherMapField<EvalWorklet, DeviceAdapterTag>;

832 833 834 835
      auto s1 = vtkm::cont::make_ImplicitFunctionHandle(Internal.Sphere1);
      auto s2 = vtkm::cont::make_ImplicitFunctionHandle(Internal.Sphere2);
      EvalWorklet eval(s1.PrepareForExecution(DeviceAdapterTag()),
                       s2.PrepareForExecution(DeviceAdapterTag()));
836 837 838 839 840 841 842 843 844 845

      vtkm::cont::Timer<DeviceAdapterTag> timer;
      EvalDispatcher(eval).Invoke(this->Internal.Points, this->Internal.Result);
      return timer.GetElapsedTime();
    }

    VTKM_CONT
    std::string Description() const
    {
      std::stringstream description;
846
      description << "Implicit Function 2x(VirtualImplicitFunction) on "
847
                  << Internal.Points.GetNumberOfValues() << " points";
848 849 850 851 852 853 854
      return description.str();
    }

    ImplicitFunctionBenchData Internal;
  };

  VTKM_MAKE_BENCHMARK(ImplicitFunction, BenchImplicitFunction);
855
  VTKM_MAKE_BENCHMARK(ImplicitFunctionVirtual, BenchVirtualImplicitFunction);
856
  VTKM_MAKE_BENCHMARK(ImplicitFunction2, Bench2ImplicitFunctions);
857
  VTKM_MAKE_BENCHMARK(ImplicitFunctionVirtual2, Bench2VirtualImplicitFunctions);
858

859
public:
860 861
  static VTKM_CONT int Run(int benchmarks)
  {
862 863
    std::cout << DIVIDER << "\nRunning Field Algorithm benchmarks\n";

864 865
    if (benchmarks & BLACK_SCHOLES)
    {
866 867
      std::cout << DIVIDER << "\nBenchmarking BlackScholes\n";
      VTKM_RUN_BENCHMARK(BlackScholes, ValueTypes());
868
      VTKM_RUN_BENCHMARK(BlackScholesDynamic, ValueTypes());
869 870
    }

871 872
    if (benchmarks & MATH)
    {
873 874
      std::cout << DIVIDER << "\nBenchmarking Multiple Math Worklets\n";
      VTKM_RUN_BENCHMARK(Math, ValueTypes());
875
      VTKM_RUN_BENCHMARK(MathDynamic, ValueTypes());
876 877
    }

878 879
    if (benchmarks & FUSED_MATH)
    {
880 881
      std::cout << DIVIDER << "\nBenchmarking Single Fused Math Worklet\n";
      VTKM_RUN_BENCHMARK(FusedMath, ValueTypes());
882
      VTKM_RUN_BENCHMARK(FusedMathDynamic, ValueTypes());
883 884
    }

885 886
    if (benchmarks & INTERPOLATE_FIELD)
    {
887
      std::cout << DIVIDER << "\nBenchmarking Edge Based Field InterpolationWorklet\n";
888 889
      VTKM_RUN_BENCHMARK(EdgeInterp, InterpValueTypes());
      VTKM_RUN_BENCHMARK(EdgeInterpDynamic, InterpValueTypes());
890 891
    }

892 893
    if (benchmarks & IMPLICIT_FUNCTION)
    {
894 895 896 897
      using FloatDefaultType = vtkm::ListTagBase<vtkm::FloatDefault>;

      std::cout << "\nBenchmarking Implicit Function\n";
      VTKM_RUN_BENCHMARK(ImplicitFunction, FloatDefaultType());
898
      VTKM_RUN_BENCHMARK(ImplicitFunctionVirtual, FloatDefaultType());
899
      VTKM_RUN_BENCHMARK(ImplicitFunction2, FloatDefaultType());
900
      VTKM_RUN_BENCHMARK(ImplicitFunctionVirtual2, FloatDefaultType());
901 902
    }

903 904 905 906 907 908 909 910
    return 0;
  }
};

#undef ARRAY_SIZE
}
} // namespace vtkm::benchmarking

911
int main(int argc, char* argv[])
912 913
{
  int benchmarks = 0;
914 915
  if (argc < 2)
  {
916 917
    benchmarks = vtkm::benchmarking::ALL;
  }
918 919 920 921
  else
  {
    for (int i = 1; i < argc; ++i)
    {
922
      std::string arg = argv[i];
923 924 925
      std::transform(arg.begin(), arg.end(), arg.begin(), [](char c) {
        return static_cast<char>(std::tolower(static_cast<unsigned char>(c)));
      });
926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941
      if (arg == "blackscholes")
      {
        benchmarks |= vtkm::benchmarking::BLACK_SCHOLES;
      }
      else if (arg == "math")
      {
        benchmarks |= vtkm::benchmarking::MATH;
      }
      else if (arg == "fusedmath")
      {
        benchmarks |= vtkm::benchmarking::FUSED_MATH;
      }
      else if (arg == "interpolate")
      {
        benchmarks |= vtkm::benchmarking::INTERPOLATE_FIELD;
      }
942 943 944 945
      else if (arg == "implicit_function")
      {
        benchmarks |= vtkm::benchmarking::IMPLICIT_FUNCTION;
      }
946 947 948 949 950 951 952 953 954
      else
      {
        std::cout << "Unrecognized benchmark: " << argv[i] << std::endl;
        return 1;
      }
    }
  }

  //now actually execute the benchmarks
955 956
  return vtkm::benchmarking::BenchmarkFieldAlgorithms<VTKM_DEFAULT_DEVICE_ADAPTER_TAG>::Run(
    benchmarks);
957
}