// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-License-Identifier: BSD-3-Clause

//
// Note if you fix this test to fill in all the empty tests
// then remove the cppcheck suppression in VTKcppcheckSuppressions.txt
//
#include "vtkDataArrayRange.h"
#include "vtkIntArray.h"
#include "vtkMath.h"
#include "vtkMathUtilities.h"
#include "vtkSmartPointer.h"
#include "vtkType.h"
#include "vtkUnsignedCharArray.h"
#include "vtkUnsignedShortArray.h"

#include <array>
#include <numeric>
#include <vector>

#include <iostream>

static int TestPi();
static int TestDegreesFromRadians();
static int TestRound();
static int TestFloor();
static int TestCeil();
static int TestCeilLog2();
static int TestIsPowerOfTwo();
static int TestNearestPowerOfTwo();
static int TestFactorial();
static int TestBinomial();
static int TestRandom();
static int TestAddSubtract();
static int TestMultiplyScalar();
static int TestMultiplyScalar2D();
static int TestDot();
static int TestOuter();
static int TestCross();
static int TestNorm();
static int TestNormalize();
static int TestPerpendiculars();
static int TestProjectVector();
static int TestProjectVector2D();
static int TestDistance2BetweenPoints();
static int TestAngleBetweenVectors();
static int TestGaussianAmplitude();
static int TestGaussianWeight();
static int TestDot2D();
static int TestNorm2D();
static int TestNormalize2D();
static int TestDeterminant2x2();
static int TestDeterminant3x3();
static int TestLUFactor3x3();
static int TestLUSolve3x3();
static int TestLinearSolve3x3();
static int TestMultiply3x3();
static int TestMultiplyMatrix();
static int TestTranspose3x3();
static int TestInvert3x3();
static int TestInvertMatrix();
static int TestIdentity3x3();
static int TestQuaternionToMatrix3x3();
static int TestMatrix3x3ToQuaternion();
static int TestMultiplyQuaternion();
static int TestOrthogonalize3x3();
static int TestDiagonalize3x3();
static int TestSingularValueDecomposition3x3();
static int TestSolveLinearSystem();
static int TestSolveLeastSquares();
static int TestSolveHomogeneousLeastSquares();
static int TestLUSolveLinearSystemEstimateMatrixCondition();
static int TestJacobiN();
static int TestClampValue();
static int TestClampValues();
static int TestClampAndNormalizeValue();
static int TestTensorFromSymmetricTensor();
static int TestGetScalarTypeFittingRange();
static int TestGetAdjustedScalarRange();
static int TestExtentIsWithinOtherExtent();
static int TestBoundsIsWithinOtherBounds();
static int TestPointIsWithinBounds();
static int TestSolve3PointCircle();
static int TestRGBToHSV();
static int TestInf();
static int TestNegInf();
static int TestNan();
static int TestMetaMultiplyMatrix();
static int TestMetaMultiplyMatrixWithVector();
static int TestMetaDot();
static int TestMetaDeterminant();
static int TestMetaInvertMatrix();
static int TestMetaLinearSolve();

int UnitTestMath(int, char*[])
{
  int status = 0;

  status += TestPi();

  status += TestDegreesFromRadians();
  status += TestRound();
  status += TestFloor();
  status += TestCeil();
  status += TestCeilLog2();
  status += TestIsPowerOfTwo();
  status += TestNearestPowerOfTwo();
  status += TestFactorial();
  status += TestBinomial();
  status += TestRandom();
  status += TestAddSubtract();
  status += TestMultiplyScalar();
  status += TestMultiplyScalar2D();
  status += TestDot();
  status += TestOuter();
  status += TestCross();
  status += TestNorm();
  status += TestNormalize();
  status += TestPerpendiculars();
  status += TestProjectVector();
  status += TestProjectVector2D();
  status += TestDistance2BetweenPoints();
  status += TestAngleBetweenVectors();
  status += TestGaussianAmplitude();
  status += TestGaussianWeight();
  status += TestDot2D();
  status += TestNorm2D();
  status += TestNormalize2D();
  status += TestDeterminant2x2();
  status += TestDeterminant3x3();
  status += TestLUFactor3x3();
  status += TestLUSolve3x3();
  status += TestLinearSolve3x3();
  status += TestMultiply3x3();
  status += TestMultiplyMatrix();
  status += TestTranspose3x3();
  status += TestInvert3x3();
  status += TestInvertMatrix();
  status += TestIdentity3x3();
  status += TestQuaternionToMatrix3x3();
  status += TestMatrix3x3ToQuaternion();
  status += TestMultiplyQuaternion();
  status += TestOrthogonalize3x3();
  status += TestDiagonalize3x3();
  status += TestSingularValueDecomposition3x3();
  status += TestSolveLinearSystem();
  status += TestSolveLeastSquares();
  status += TestSolveHomogeneousLeastSquares();
  status += TestLUSolveLinearSystemEstimateMatrixCondition();
  status += TestJacobiN();
  status += TestClampValue();
  status += TestClampValues();
  status += TestClampAndNormalizeValue();
  status += TestTensorFromSymmetricTensor();
  status += TestGetScalarTypeFittingRange();
  status += TestGetAdjustedScalarRange();
  status += TestExtentIsWithinOtherExtent();
  status += TestBoundsIsWithinOtherBounds();
  status += TestPointIsWithinBounds();
  status += TestSolve3PointCircle();
  status += TestRGBToHSV();
  status += TestInf();
  status += TestNegInf();
  status += TestNan();
  status += TestMetaMultiplyMatrix();
  status += TestMetaMultiplyMatrixWithVector();
  status += TestMetaDot();
  status += TestMetaDeterminant();
  status += TestMetaInvertMatrix();
  status += TestMetaLinearSolve();
  if (status != 0)
  {
    return EXIT_FAILURE;
  }

  vtkSmartPointer<vtkMath> math = vtkSmartPointer<vtkMath>::New();
  math->Print(std::cout);

  return EXIT_SUCCESS;
}

// Validate by comparing to atan/4
int TestPi()
{
  int status = 0;
  std::cout << "Pi..";
  if (vtkMath::Pi() != std::atan(1.0) * 4.0)
  {
    std::cout << "Expected " << vtkMath::Pi() << " but got " << std::atan(1.0) * 4.0;
    ++status;
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate against RadiansFromDegrees
int TestDegreesFromRadians()
{
  int status = 0;
  std::cout << "DegreesFromRadians..";

  unsigned int numSamples = 1000;
  for (unsigned int i = 0; i < numSamples; ++i)
  {
    float floatDegrees = vtkMath::Random(-180.0, 180.0);
    float floatRadians = vtkMath::RadiansFromDegrees(floatDegrees);
    float result = vtkMath::DegreesFromRadians(floatRadians);
    if (!vtkMathUtilities::FuzzyCompare(
          result, floatDegrees, std::numeric_limits<float>::epsilon() * 128.0f))
    {
      std::cout << "Float Expected " << floatDegrees << " but got " << result << " difference is "
                << result - floatDegrees << " ";
      std::cout << "eps ratio is: "
                << (result - floatDegrees) / std::numeric_limits<float>::epsilon() << std::endl;
      ++status;
    }
  }
  for (unsigned int i = 0; i < numSamples; ++i)
  {
    double doubleDegrees = vtkMath::Random(-180.0, 180.0);
    double doubleRadians = vtkMath::RadiansFromDegrees(doubleDegrees);
    double result = vtkMath::DegreesFromRadians(doubleRadians);
    if (!vtkMathUtilities::FuzzyCompare(
          result, doubleDegrees, std::numeric_limits<double>::epsilon() * 256.0))
    {
      std::cout << " Double Expected " << doubleDegrees << " but got " << result
                << " difference is " << result - doubleDegrees;
      std::cout << " eps ratio is: "
                << (result - doubleDegrees) / std::numeric_limits<double>::epsilon() << std::endl;
      ++status;
    }
  }
  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with http://en.wikipedia.org/wiki/Rounding#Rounding_to_integer
int TestRound()
{
  int status = 0;
  std::cout << "Round..";
  int result;
  {
    std::vector<float> values;
    std::vector<int> expecteds;

    values.push_back(23.67f);
    expecteds.push_back(24);
    values.push_back(23.50f);
    expecteds.push_back(24);
    values.push_back(23.35f);
    expecteds.push_back(23);
    values.push_back(23.00f);
    expecteds.push_back(23);
    values.push_back(0.00f);
    expecteds.push_back(0);
    values.push_back(-23.00f);
    expecteds.push_back(-23);
    values.push_back(-23.35f);
    expecteds.push_back(-23);
    values.push_back(-23.50f);
    expecteds.push_back(-24);
    values.push_back(-23.67f);
    expecteds.push_back(-24);
    for (size_t i = 0; i < values.size(); ++i)
    {
      result = vtkMath::Round(values[i]);
      if (result != expecteds[i])
      {
        std::cout << " Float Round(" << values[i] << ") got " << result << " but expected "
                  << expecteds[i];
        ++status;
      }
    }
  }
  {
    std::vector<double> values;
    std::vector<int> expecteds;

    values.push_back(23.67);
    expecteds.push_back(24);
    values.push_back(23.50);
    expecteds.push_back(24);
    values.push_back(23.35);
    expecteds.push_back(23);
    values.push_back(23.00);
    expecteds.push_back(23);
    values.push_back(0.00);
    expecteds.push_back(0);
    values.push_back(-23.00);
    expecteds.push_back(-23);
    values.push_back(-23.35);
    expecteds.push_back(-23);
    values.push_back(-23.50);
    expecteds.push_back(-24);
    values.push_back(-23.67);
    expecteds.push_back(-24);
    for (size_t i = 0; i < values.size(); ++i)
    {
      result = vtkMath::Round(values[i]);
      if (result != expecteds[i])
      {
        std::cout << " Double Round(" << values[i] << ") got " << result << " but expected "
                  << expecteds[i];
        ++status;
      }
    }
  }
  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with http://en.wikipedia.org/wiki/Floor_and_ceiling_functions
int TestFloor()
{
  int status = 0;
  std::cout << "Floor..";

  int result;
  std::vector<double> values;
  std::vector<int> expecteds;

  values.push_back(2.4);
  expecteds.push_back(2);
  values.push_back(2.7);
  expecteds.push_back(2);
  values.push_back(-2.7);
  expecteds.push_back(-3);
  values.push_back(-2.0);
  expecteds.push_back(-2);
  for (size_t i = 0; i < values.size(); ++i)
  {
    result = vtkMath::Floor(values[i]);
    if (result != expecteds[i])
    {
      std::cout << " Floor(" << values[i] << ") got " << result << " but expected " << expecteds[i];
      ++status;
    }
  }
  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with http://en.wikipedia.org/wiki/Floor_and_ceiling_functions
int TestCeil()
{
  int status = 0;
  std::cout << "Ceil..";

  int result;
  std::vector<double> values;
  std::vector<int> expecteds;

  values.push_back(2.4);
  expecteds.push_back(3);
  values.push_back(2.7);
  expecteds.push_back(3);
  values.push_back(-2.7);
  expecteds.push_back(-2);
  values.push_back(-2.0);
  expecteds.push_back(-2);
  for (size_t i = 0; i < values.size(); ++i)
  {
    result = vtkMath::Ceil(values[i]);
    if (result != expecteds[i])
    {
      std::cout << " Ceil(" << values[i] << ") got " << result << " but expected " << expecteds[i];
      ++status;
    }
  }
  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate by powers of 2 perturbations
int TestCeilLog2()
{
  int status = 0;
  std::cout << "CeilLog2..";

  int result;
  std::vector<vtkTypeUInt64> values;
  std::vector<int> expecteds;

  for (unsigned int p = 0; p < 30; ++p)
  {
    vtkTypeUInt64 shifted = (2 << p) + 1;
    values.push_back(shifted);
    expecteds.push_back(p + 2);
    shifted = (2 << p);
    values.push_back(shifted);
    expecteds.push_back(p + 1);
  }
  for (size_t i = 0; i < values.size(); ++i)
  {
    result = vtkMath::CeilLog2(values[i]);
    if (result != expecteds[i])
    {
      std::cout << " CeilLog2(" << values[i] << ") got " << result << " but expected "
                << expecteds[i];
      ++status;
    }
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate by powers of 2 perturbations
int TestIsPowerOfTwo()
{
  int status = 0;
  std::cout << "IsPowerOfTwo..";
  bool result;

  std::vector<vtkTypeUInt64> values;
  std::vector<bool> expecteds;
  int largestPower = std::numeric_limits<vtkTypeUInt64>::digits;
  vtkTypeUInt64 shifted = 1;
  for (int p = 1; p < largestPower - 1; ++p)
  {
    shifted *= 2;
    values.push_back(shifted);
    expecteds.push_back(true);
    if (shifted != 2)
    {
      values.push_back(shifted - 1);
      expecteds.push_back(false);
    }
    if (shifted < std::numeric_limits<vtkTypeUInt64>::max() - 1)
    {
      values.push_back(shifted + 1);
      expecteds.push_back(false);
    }
  }
  for (size_t i = 0; i < values.size(); ++i)
  {
    result = vtkMath::IsPowerOfTwo(values[i]);
    if (result != expecteds[i])
    {
      std::cout << " IsPowerOfTwo(" << values[i] << ") got " << result << " but expected "
                << expecteds[i];
      ++status;
    }
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate by powers of 2 perturbations
int TestNearestPowerOfTwo()
{
  int status = 0;
  std::cout << "NearestPowerOfTwo..";

  std::vector<vtkTypeUInt64> values;
  std::vector<int> expecteds;

  values.push_back(0);
  expecteds.push_back(1);

  int numDigits = std::numeric_limits<int>::digits;
  vtkTypeUInt64 shifted = 1;
  for (int p = 0; p < numDigits; ++p)
  {
    values.push_back(shifted);
    expecteds.push_back(shifted);
    if (shifted <= INT_MAX / 2)
    {
      values.push_back(shifted + 1);
      expecteds.push_back(shifted * 2);
    }
    if (shifted != 2)
    {
      values.push_back(shifted - 1);
      expecteds.push_back(shifted);
    }

    shifted *= 2;
  }

  values.push_back(INT_MAX);
  expecteds.push_back(INT_MIN);

  for (size_t i = 0; i < values.size(); ++i)
  {
    int result = vtkMath::NearestPowerOfTwo(values[i]);
    if (result != expecteds[i])
    {
      std::cout << " NearestPowerOfTwo(" << values[i] << ") got " << result << " but expected "
                << expecteds[i];
      ++status;
    }
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate by alternate computation
int TestFactorial()
{
  int status = 0;
  std::cout << "Factorial..";

  std::vector<int> values;
  std::vector<vtkTypeInt64> expecteds;
  vtkTypeInt64 expected = 1;
  for (int f = 2; f < 10; ++f)
  {
    expected *= f;
    values.push_back(f);
    expecteds.push_back(expected);
  }
  for (size_t i = 0; i < values.size(); ++i)
  {
    int result = vtkMath::Factorial(values[i]);
    if (result != expecteds[i])
    {
      std::cout << " Factorial(" << values[i] << ") got " << result << " but expected "
                << expecteds[i];
      ++status;
    }
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate by alternative computations
int TestBinomial()
{
  int status = 0;
  int m;
  int n;
  std::cout << "Binomial..";

  std::vector<int> mvalues;
  std::vector<int> nvalues;

  std::vector<vtkTypeInt64> expecteds;
  double expected;
  for (m = 1; m < 31; ++m)
  {
    for (n = 1; n <= m; ++n)
    {
      mvalues.push_back(m);
      nvalues.push_back(n);
      expected = 1;
      for (int i = 1; i <= n; ++i)
      {
        expected *= static_cast<double>(m - i + 1) / i;
      }
      expecteds.push_back(static_cast<vtkTypeInt64>(expected));
    }
  }

  for (size_t i = 0; i < mvalues.size(); ++i)
  {
    int result = vtkMath::Binomial(mvalues[i], nvalues[i]);
    if (result != expecteds[i])
    {
      std::cout << " Binomial(" << mvalues[i] << ", " << nvalues[i] << ") got " << result
                << " but expected " << expecteds[i];
      ++status;
    }
  }

  // Now test the combination iterator
  m = 6;
  n = 3;
  int more = 1;
  int count = 0;
  int* comb;
  // First, m < n should produce 0
  comb = vtkMath::BeginCombination(n, m);
  if (comb != nullptr)
  {
    ++status;
    std::cout << " Combinations(" << n << ", " << m << ") should return 0 "
              << " but got " << comb;
  }
  comb = vtkMath::BeginCombination(m, n);
  while (more)
  {
    ++count;
    more = vtkMath::NextCombination(m, n, comb);
  }
  vtkMath::FreeCombination(comb);
  if (count != vtkMath::Binomial(m, n))
  {
    ++status;
    std::cout << " Combinations(" << m << ", " << n << ") got " << count << " but expected "
              << vtkMath::Binomial(m, n);
  }
  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// No validation
int TestRandom()
{
  int status = 0;
  std::cout << "Random..";
  // Not really a test of randomness, just covering code
  int n = 1000;
  vtkMath::RandomSeed(8775070);
  vtkMath::GetSeed(); // just for coverage
  double accum = 0.0;
  for (int i = 0; i < n; ++i)
  {
    float random = vtkMath::Random();
    accum += random;
    if (random < 0.0 || random > 1.0)
    {
      std::cout << "Random(): " << random << " out of range" << std::endl;
      ++status;
    }
    random = vtkMath::Gaussian();
    accum += random;

    random = vtkMath::Gaussian(0.0, 1.0);
    accum += random;

    random = vtkMath::Random(-1000.0, 1000.0);
    accum += random;
    if (random < -1000.0 || random > 1000.0)
    {
      std::cout << "Random (-1000.0, 1000.0): " << random << " out of range" << std::endl;
      ++status;
    }
  }
  if (accum == 0.0)
  {
    ++status;
  }
  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

template <typename T>
int AddSubtract()
{
  int status = 0;
  T da[3], db[3], dc[3], dd[3];
  for (int n = 0; n < 100000; ++n)
  {
    for (int i = 0; i < 3; ++i)
    {
      da[i] = vtkMath::Random(-10.0, 10.0);
      db[i] = vtkMath::Random(-10.0, 10.0);
    }
    vtkMath::Add(da, db, dc);
    vtkMath::Subtract(dc, db, dd);
    for (int i = 0; i < 3; ++i)
    {
      if (!vtkMathUtilities::FuzzyCompare(
            da[i], dd[i], std::numeric_limits<T>::epsilon() * (T)256.0))
      {
        std::cout << " Add/Subtract got " << dd[i] << " but expected " << da[i];
        ++status;
      }
    }
  }
  return status;
}

// Validate by a + b - b = a
int TestAddSubtract()
{
  int status = 0;
  std::cout << "AddSubtract..";

  status += AddSubtract<double>();
  status += AddSubtract<float>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

template <typename T>
int MultiplyScalar()
{
  int status = 0;
  // first T
  T da[3], db[3];
  for (int n = 0; n < 100000; ++n)
  {
    for (int i = 0; i < 3; ++i)
    {
      da[i] = vtkMath::Random(-10.0, 10.0);
      db[i] = da[i];
    }
    T scale = vtkMath::Random();
    vtkMath::MultiplyScalar(da, scale);

    for (int i = 0; i < 3; ++i)
    {
      if (!vtkMathUtilities::FuzzyCompare(
            da[i], db[i] * scale, std::numeric_limits<T>::epsilon() * (T)256.0))
      {
        std::cout << " MultiplyScalar got " << da[i] << " but expected " << db[i] * scale;
        ++status;
      }
    }
  }

  return status;
}

int TestMultiplyScalar()
{
  int status = 0;
  std::cout << "MultiplyScalar..";

  status += MultiplyScalar<double>();
  status += MultiplyScalar<float>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

int TestMultiplyScalar2D()
{
  int status = 0;
  std::cout << "MultiplyScalar2D..";

  // now 2D
  // first double
  double da[2], db[2];
  for (int n = 0; n < 100000; ++n)
  {
    for (int i = 0; i < 2; ++i)
    {
      da[i] = vtkMath::Random(-10.0, 10.0);
      db[i] = da[i];
    }
    double scale = vtkMath::Random();
    vtkMath::MultiplyScalar2D(da, scale);

    for (int i = 0; i < 2; ++i)
    {
      if (!vtkMathUtilities::FuzzyCompare(
            da[i], db[i] * scale, std::numeric_limits<double>::epsilon() * 256.0))
      {
        std::cout << " MultiplyScalar2D got " << da[i] << " but expected " << db[i] * scale;
        ++status;
      }
    }
  }

  // then float
  float fa[2], fb[2];
  for (int n = 0; n < 100000; ++n)
  {
    for (int i = 0; i < 2; ++i)
    {
      fa[i] = vtkMath::Random(-10.0, 10.0);
      fb[i] = fa[i];
    }
    float scale = vtkMath::Random();
    vtkMath::MultiplyScalar2D(fa, scale);

    for (int i = 0; i < 2; ++i)
    {
      if (!vtkMathUtilities::FuzzyCompare(
            fa[i], fb[i] * scale, std::numeric_limits<float>::epsilon() * 256.0f))
      {
        std::cout << " MultiplyScalar2D got " << fa[i] << " but expected " << fb[i] * scale;
        ++status;
      }
    }
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

class valueDouble3D
{
public:
  valueDouble3D() = default;
  valueDouble3D(double aa[3], double bb[3])
  {
    for (int i = 0; i < 3; ++i)
    {
      a[i] = aa[i];
      b[i] = bb[i];
    }
  }
  double a[3];
  double b[3];
};

class valueFloat3D
{
public:
  valueFloat3D() = default;
  valueFloat3D(float aa[3], float bb[3])
  {
    for (int i = 0; i < 3; ++i)
    {
      a[i] = aa[i];
      b[i] = bb[i];
    }
  }
  float a[3];
  float b[3];
};
int TestDot()
{
  int status = 0;
  std::cout << "Dot..";

  {
    std::vector<valueDouble3D> values;
    std::vector<double> expecteds;
    for (int n = 0; n < 100; ++n)
    {
      valueDouble3D v;
      double dot = 0.0;
      for (int i = 0; i < 3; ++i)
      {
        v.a[i] = vtkMath::Random();
        v.b[i] = vtkMath::Random();
        dot += (v.a[i] * v.b[i]);
      }
      values.push_back(v);
      expecteds.push_back(dot);
    }
    valueDouble3D test;
    test.a[0] = 0.0;
    test.a[1] = 0.0;
    test.a[2] = 1.0;
    test.b[0] = 1.0;
    test.b[1] = 0.0;
    test.b[2] = 0.0;
    values.push_back(test);
    expecteds.push_back(0.0);
    test.a[0] = 0.0;
    test.a[1] = 0.0;
    test.a[2] = 1.0;
    test.b[0] = 0.0;
    test.b[1] = 1.0;
    test.b[2] = 0.0;
    values.push_back(test);
    expecteds.push_back(0.0);
    test.a[0] = 1.0;
    test.a[1] = 0.0;
    test.a[2] = 0.0;
    test.b[0] = 0.0;
    test.b[1] = 1.0;
    test.b[2] = 0.0;
    values.push_back(test);
    expecteds.push_back(0.0);

    for (size_t i = 0; i < values.size(); ++i)
    {
      double result = vtkMath::Dot(values[i].a, values[i].b);
      if (!vtkMathUtilities::FuzzyCompare(
            result, expecteds[i], std::numeric_limits<double>::epsilon() * 128.0))
      {
        std::cout << " Dot got " << result << " but expected " << expecteds[i];
        ++status;
      }
    }
  }

  // now float
  {
    std::vector<valueFloat3D> values;
    std::vector<float> expecteds;
    for (int n = 0; n < 100; ++n)
    {
      valueFloat3D v;
      float dot = 0.0;
      for (int i = 0; i < 3; ++i)
      {
        v.a[i] = vtkMath::Random();
        v.b[i] = vtkMath::Random();
        dot += (v.a[i] * v.b[i]);
      }
      values.push_back(v);
      expecteds.push_back(dot);
    }
    valueFloat3D test;
    test.a[0] = 0.0;
    test.a[1] = 0.0;
    test.a[2] = 1.0;
    test.b[0] = 1.0;
    test.b[1] = 0.0;
    test.b[2] = 0.0;
    values.push_back(test);
    expecteds.push_back(0.0);
    test.a[0] = 0.0;
    test.a[1] = 0.0;
    test.a[2] = 1.0;
    test.b[0] = 0.0;
    test.b[1] = 1.0;
    test.b[2] = 0.0;
    values.push_back(test);
    expecteds.push_back(0.0);
    test.a[0] = 1.0;
    test.a[1] = 0.0;
    test.a[2] = 0.0;
    test.b[0] = 0.0;
    test.b[1] = 1.0;
    test.b[2] = 0.0;
    values.push_back(test);
    expecteds.push_back(0.0);

    for (size_t i = 0; i < values.size(); ++i)
    {
      float result = vtkMath::Dot(values[i].a, values[i].b);
      if (!vtkMathUtilities::FuzzyCompare(
            result, expecteds[i], std::numeric_limits<float>::epsilon() * 128.0f))
      {
        std::cout << " Dot got " << result << " but expected " << expecteds[i];
        ++status;
      }
    }
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

int TestOuter()
{
  int status = 0;
  std::cout << "Outer..";

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Verify by anticommutative property
template <typename T>
int Cross()
{
  int status = 0;
  T a[3];
  T b[3];
  T c[3];
  T d[3];

  for (int n = 0; n < 1000; ++n)
  {
    for (int i = 0; i < 3; ++i)
    {
      a[i] = vtkMath::Random(-1.0, 1.0);
      b[i] = vtkMath::Random(-1.0, 1.0);
    }
    vtkMath::Cross(a, b, c);
    vtkMath::MultiplyScalar(b, (T)-1.0);
    vtkMath::Cross(b, a, d);
    // a x b = -b x a
    for (int i = 0; i < 3; ++i)
    {
      if (!vtkMathUtilities::FuzzyCompare(c[i], d[i], std::numeric_limits<T>::epsilon() * (T)128.0))
      {
        std::cout << " Cross expected " << c[i] << " but got " << d[i];
        std::cout << "eps ratio is: " << (c[i] - d[i]) / std::numeric_limits<T>::epsilon()
                  << std::endl;
        ++status;
      }
    }
  }
  return status;
}

int TestCross()
{
  int status = 0;
  std::cout << "Cross..";

  status += Cross<double>();
  status += Cross<float>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

template <typename T, int NDimension>
int Norm()
{
  int status = 0;
  T x[NDimension];

  for (int n = 0; n < 1000; ++n)
  {
    for (int i = 0; i < NDimension; ++i)
    {
      x[i] = (T)vtkMath::Random(-10.0, 10.0);
    }

    T norm = vtkMath::Norm(x, NDimension);

    for (int i = 0; i < NDimension; ++i)
    {
      x[i] /= norm;
    }

    T unitNorm = vtkMath::Norm(x, NDimension);
    if (!vtkMathUtilities::FuzzyCompare(
          unitNorm, (T)1.0, std::numeric_limits<T>::epsilon() * (T)128.0))
    {
      std::cout << "Norm Expected " << 1.0 << " but got " << unitNorm;
      std::cout << " eps ratio is: " << ((T)1.0 - unitNorm) / std::numeric_limits<T>::epsilon()
                << std::endl;
      ++status;
    }
  }

  return status;
}

int TestNorm()
{
  int status = 0;
  std::cout << "Norm..";

  status += Norm<double, 1>();
  status += Norm<double, 3>();
  status += Norm<double, 1000>();
  status += Norm<float, 1>();
  status += Norm<float, 3>();
  status += Norm<float, 1000>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate compute Norm, should be 1.0

template <typename T>
int Normalize()
{
  int status = 0;
  for (int n = 0; n < 1000; ++n)
  {
    T a[3];
    for (int i = 0; i < 3; ++i)
    {
      a[i] = vtkMath::Random(-10000.0, 10000.0);
    }
    vtkMath::Normalize(a);
    T value = vtkMath::Norm(a);
    T expected = 1.0;
    if (!vtkMathUtilities::FuzzyCompare(
          value, expected, std::numeric_limits<T>::epsilon() * (T)128.0))
    {
      std::cout << " Normalize expected " << expected << " but got " << value;
      std::cout << "eps ratio is: " << value - expected / std::numeric_limits<T>::epsilon()
                << std::endl;
      ++status;
    }
  }
  return status;
}

int TestNormalize()
{
  int status = 0;
  std::cout << "Normalize..";

  status += Normalize<double>();
  status += Normalize<float>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

int TestPerpendiculars()
{
  int status = 0;
  std::cout << "Perpendiculars..";
  {
    // first double
    double x[3], y[3], z[3];
    std::vector<valueDouble3D> values;
    std::vector<double> expecteds;
    for (int n = 0; n < 100; ++n)
    {
      for (int i = 0; i < 3; ++i)
      {
        x[i] = vtkMath::Random(-10.0, 10.0);
      }
      vtkMath::Perpendiculars(x, y, z, vtkMath::Random(-vtkMath::Pi(), vtkMath::Pi()));
      {
        valueDouble3D value(x, y);
        values.push_back(value);
        expecteds.push_back(0.0);
      }
      {
        valueDouble3D value(x, z);
        values.push_back(value);
        expecteds.push_back(0.0);
      }
      {
        valueDouble3D value(y, z);
        values.push_back(value);
        expecteds.push_back(0.0);
      }
      vtkMath::Perpendiculars(x, y, z, 0.0);
      {
        valueDouble3D value(x, y);
        values.push_back(value);
        expecteds.push_back(0.0);
      }
    }
    for (size_t i = 0; i < values.size(); ++i)
    {
      double test = vtkMath::Dot(values[i].a, values[i].b);
      if (!vtkMathUtilities::FuzzyCompare(
            expecteds[i], test, std::numeric_limits<double>::epsilon() * 256.0))
      {
        std::cout << " Perpendiculars got " << test << " but expected " << expecteds[i];
        ++status;
      }
    }
  }
  {
    // then floats
    float x[3], y[3], z[3];
    std::vector<valueFloat3D> values;
    std::vector<float> expecteds;
    for (int n = 0; n < 100; ++n)
    {
      for (int i = 0; i < 3; ++i)
      {
        x[i] = vtkMath::Random(-10.0, 10.0);
      }
      vtkMath::Perpendiculars(x, y, z, vtkMath::Random(-vtkMath::Pi(), vtkMath::Pi()));
      {
        valueFloat3D value(x, y);
        values.push_back(value);
        expecteds.push_back(0.0);
      }
      {
        valueFloat3D value(x, z);
        values.push_back(value);
        expecteds.push_back(0.0);
      }
      {
        valueFloat3D value(y, z);
        values.push_back(value);
        expecteds.push_back(0.0);
      }
      vtkMath::Perpendiculars(x, y, z, 0.0);
      {
        valueFloat3D value(x, y);
        values.push_back(value);
        expecteds.push_back(0.0);
      }
    }
    for (size_t i = 0; i < values.size(); ++i)
    {
      float test = vtkMath::Dot(values[i].a, values[i].b);
      if (!vtkMathUtilities::FuzzyCompare(
            expecteds[i], test, std::numeric_limits<float>::epsilon() * 256.0f))
      {
        std::cout << " Perpendiculars got " << test << " but expected " << expecteds[i];
        ++status;
      }
    }
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

template <typename T>
int ProjectVector()
{
  int status = 0;
  T a[3], b[3], c[3];
  for (int i = 0; i < 3; ++i)
  {
    a[i] = 0.0;
    b[i] = 0.0;
  }
  if (vtkMath::ProjectVector(a, b, c))
  {
    std::cout << "ProjectVector of a 0 vector should return false ";
    ++status;
  }
  return status;
}

// Validate 0 vector case. TestMath does the rest
int TestProjectVector()
{
  int status = 0;
  std::cout << "ProjectVector..";

  status += ProjectVector<double>();
  status += ProjectVector<float>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

template <typename T>
int ProjectVector2D()
{
  int status = 0;
  T a[2], b[2], c[2];
  for (int i = 0; i < 2; ++i)
  {
    a[i] = 0.0;
    b[i] = 0.0;
  }
  if (vtkMath::ProjectVector2D(a, b, c))
  {
    std::cout << "ProjectVector2D of a 0 vector should return false ";
    ++status;
  }
  return status;
}

// Validate 0 vector case. TestMath does the rest
int TestProjectVector2D()
{
  int status = 0;
  std::cout << "ProjectVector2D..";

  status += ProjectVector2D<double>();
  status += ProjectVector2D<float>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

int TestDistance2BetweenPoints()
{
  int status = 0;
  std::cout << "Distance2BetweenPoints..";

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

int TestAngleBetweenVectors()
{
  int status = 0;
  std::cout << "AngleBetweenVectors..";

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

int TestGaussianAmplitude()
{
  int status = 0;
  std::cout << "GaussianAmplitude..";

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

int TestGaussianWeight()
{
  int status = 0;
  std::cout << "GaussianWeight..";

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

class valueDouble2D
{
public:
  double a[2];
  double b[2];
};
class valueFloat2D
{
public:
  float a[2];
  float b[2];
};
int TestDot2D()
{
  int status = 0;
  std::cout << "Dot2D..";

  {
    std::vector<valueDouble2D> values;
    std::vector<double> expecteds;
    for (int n = 0; n < 100; ++n)
    {
      valueDouble2D v;
      double dot = 0.0;
      for (int i = 0; i < 2; ++i)
      {
        v.a[i] = vtkMath::Random();
        v.b[i] = vtkMath::Random();
        dot += (v.a[i] * v.b[i]);
      }
      values.push_back(v);
      expecteds.push_back(dot);
    }
    valueDouble2D test;
    test.a[0] = 1.0;
    test.a[1] = 0.0;
    test.b[0] = 0.0;
    test.b[1] = 1.0;
    values.push_back(test);
    expecteds.push_back(0.0);

    for (size_t i = 0; i < values.size(); ++i)
    {
      double result = vtkMath::Dot2D(values[i].a, values[i].b);
      if (!vtkMathUtilities::FuzzyCompare(
            result, expecteds[i], std::numeric_limits<double>::epsilon() * 128.0))
      {
        std::cout << " Dot got " << result << " but expected " << expecteds[i];
        ++status;
      }
    }
  }

  // now float
  {
    std::vector<valueFloat2D> values;
    std::vector<float> expecteds;
    for (int n = 0; n < 100; ++n)
    {
      valueFloat2D v;
      float dot = 0.0;
      for (int i = 0; i < 2; ++i)
      {
        v.a[i] = vtkMath::Random();
        v.b[i] = vtkMath::Random();
        dot += (v.a[i] * v.b[i]);
      }
      values.push_back(v);
      expecteds.push_back(dot);
    }
    valueFloat2D test;
    test.a[0] = 0.0;
    test.a[1] = 1.0;
    test.b[0] = 1.0;
    test.b[1] = 0.0;
    values.push_back(test);
    expecteds.push_back(0.0);

    for (size_t i = 0; i < values.size(); ++i)
    {
      float result = vtkMath::Dot2D(values[i].a, values[i].b);
      if (!vtkMathUtilities::FuzzyCompare(
            result, expecteds[i], std::numeric_limits<float>::epsilon() * 128.0f))
      {
        std::cout << " Dot got " << result << " but expected " << expecteds[i];
        ++status;
      }
    }
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

int TestNorm2D()
{
  int status = 0;
  std::cout << "Norm2D..";

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

int TestNormalize2D()
{
  int status = 0;
  std::cout << "Normalize2D..";

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

int TestDeterminant2x2()
{
  int status = 0;
  std::cout << "Determinant2x2..";
  // Frank Matrix
  {
    double a[2][2];
    for (int i = 1; i <= 2; ++i)
    {
      for (int j = 1; j <= 2; ++j)
      {
        if (j < i - 2)
        {
          a[i - 1][j - 1] = 0.0;
        }
        else if (j == (i - 1))
        {
          a[i - 1][j - 1] = 2 + 1 - i;
        }
        else
        {
          a[i - 1][j - 1] = 2 + 1 - j;
        }
      }
    }
    if (vtkMath::Determinant2x2(a[0], a[1]) != 1.0)
    {
      std::cout << "Determinant2x2 expected " << 1.0 << " but got "
                << vtkMath::Determinant2x2(a[0], a[1]) << std::endl;
      ++status;
    }
  }
  {
    float a[2][2];
    for (int i = 1; i <= 2; ++i)
    {
      for (int j = 1; j <= 2; ++j)
      {
        if (j < i - 2)
        {
          a[i - 1][j - 1] = 0.0;
        }
        else if (j == (i - 1))
        {
          a[i - 1][j - 1] = 2 + 1 - i;
        }
        else
        {
          a[i - 1][j - 1] = 2 + 1 - j;
        }
      }
    }
    if (vtkMath::Determinant2x2(a[0], a[1]) != 1.0)
    {
      std::cout << "Determinant2x2 expected " << 1.0 << " but got "
                << vtkMath::Determinant2x2(a[0], a[1]) << std::endl;
      ++status;
    }
  }
  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

int TestDeterminant3x3()
{
  int status = 0;
  std::cout << "Determinant3x3..";

  // Frank Matrix
  {
    double a[3][3];
    for (int i = 1; i <= 3; ++i)
    {
      for (int j = 1; j <= 3; ++j)
      {
        if (j < i - 3)
        {
          a[i - 1][j - 1] = 0.0;
        }
        else if (j == (i - 1))
        {
          a[i - 1][j - 1] = 3 + 1 - i;
        }
        else
        {
          a[i - 1][j - 1] = 3 + 1 - j;
        }
      }
    }
    if (vtkMath::Determinant3x3(a[0], a[1], a[2]) != 1.0)
    {
      std::cout << "Determinant3x3 expected " << 1.0 << " but got "
                << vtkMath::Determinant3x3(a[0], a[1], a[2]) << std::endl;
      ++status;
    }
  }
  {
    float a[3][3];
    for (int i = 1; i <= 3; ++i)
    {
      for (int j = 1; j <= 3; ++j)
      {
        if (j < i - 3)
        {
          a[i - 1][j - 1] = 0.0;
        }
        else if (j == (i - 1))
        {
          a[i - 1][j - 1] = 3 + 1 - i;
        }
        else
        {
          a[i - 1][j - 1] = 3 + 1 - j;
        }
      }
    }
    if (vtkMath::Determinant3x3(a[0], a[1], a[2]) != 1.0)
    {
      std::cout << "Determinant3x3 expected " << 1.0 << " but got "
                << vtkMath::Determinant3x3(a[0], a[1], a[2]) << std::endl;
      ++status;
    }
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

template <typename T>
int LUFactor3x3()
{
  int status = 0;

  T A[3][3];
  int index[3];

  for (int n = 0; n < 1000; ++n)
  {
    for (int i = 0; i < 3; ++i)
    {
      for (int j = 0; j < 3; ++j)
      {
        A[i][j] = vtkMath::Random(-10.0, 10.0);
      }
    }
    vtkMath::LUFactor3x3(A, index);
  }
  return status;
}

// Just for coverage, validated as part of TestLUSolve3x3
int TestLUFactor3x3()
{
  int status = 0;
  std::cout << "LUFactor3x3..";

  status += LUFactor3x3<double>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}
template <typename T>
int LUSolve3x3()
{
  int status = 0;

  // Generate a Hilbert Matrix
  T mat[3][3];
  int index[3];
  T lhs[3];
  T rhs[3];

  for (int n = 0; n < 1000; ++n)
  {
    for (int i = 0; i < 3; ++i)
    {
      lhs[i] = vtkMath::Random(-1.0, 1.0);
    }

    for (int i = 1; i <= 3; ++i)
    {
      rhs[i - 1] = 0.0;
      for (int j = 1; j <= 3; ++j)
      {
        mat[i - 1][j - 1] = 1.0 / (i + j - 1);
        rhs[i - 1] += (mat[i - 1][j - 1] * lhs[j - 1]);
      }
    }
    vtkMath::LUFactor3x3(mat, index);
    vtkMath::LUSolve3x3(mat, index, rhs);
    for (int i = 0; i < 3; ++i)
    {
      if (!vtkMathUtilities::FuzzyCompare(
            lhs[i], rhs[i], std::numeric_limits<T>::epsilon() * (T)256.0))
      {
        std::cout << " LUSolve3x3(T) expected " << lhs[i] << " but got " << rhs[i];
        ++status;
      }
    }
  }
  return status;
}

int TestLUSolve3x3()
{
  int status = 0;
  std::cout << "LUSolve3x3..";

  status += LUSolve3x3<double>();
  status += LUSolve3x3<float>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

template <typename T>
int LinearSolve3x3()
{
  int status = 0;

  // Generate a Hilbert Matrix
  T mat[3][3];
  T lhs[3];
  T rhs[3];
  T solution[3];

  for (int n = 0; n < 2; ++n)
  {
    for (int i = 0; i < 3; ++i)
    {
      lhs[i] = vtkMath::Random(-1.0, 1.0);
    }

    for (int i = 1; i <= 3; ++i)
    {
      rhs[i - 1] = 0.0;
      for (int j = 1; j <= 3; ++j)
      {
        mat[i - 1][j - 1] = 1.0 / (i + j - 1);
        rhs[i - 1] += (mat[i - 1][j - 1] * lhs[j - 1]);
      }
    }
    vtkMath::LinearSolve3x3(mat, rhs, solution);

    for (int i = 0; i < 3; ++i)
    {
      if (!vtkMathUtilities::FuzzyCompare(
            lhs[i], solution[i], std::numeric_limits<T>::epsilon() * (T)512.0))
      {
        std::cout << " LinearSolve3x3(T) expected " << lhs[i] << " but got " << solution[i];
        ++status;
      }
    }
  }
  return status;
}

int TestLinearSolve3x3()
{
  int status = 0;
  std::cout << "LinearSolve3x3..";

  status += LinearSolve3x3<double>();
  status += LinearSolve3x3<float>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

template <typename T>
int Multiply3x3()
{
  int status = 0;
  T A[3][3];
  T V[3];
  T U[3];

  for (int i = 0; i < 3; ++i)
  {
    for (int j = 0; j < 3; ++j)
    {
      A[i][j] = vtkMath::Random(-10.0, 10.0);
    }
    V[i] = vtkMath::Random(-10.0, 10.0);
  }

  vtkMath::Multiply3x3(A, V, U);

  return status;
}

int TestMultiply3x3()
{
  int status = 0;
  std::cout << "Multiply3x3..";

  status += Multiply3x3<double>();
  status += Multiply3x3<float>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// For coverage only. Validated as part of TestInvertMatrix
int TestMultiplyMatrix()
{
  int status = 0;
  std::cout << "MultiplyMatrix..";

  double a[3][3] = { { 1.0, 2.0, 3.0 }, { 4.0, 5.0, 6.0 }, { 7.0, 8.0, 9.0 } };
  double b[3][3] = { { 1.0, 1.0, 1.0 }, { 1.0, 1.0, 1.0 }, { 1.0, 1.0, 1.0 } };
  double c[3][3];

  double* aa[3];
  double* bb[3];
  double* cc[3];
  for (int i = 0; i < 3; ++i)
  {
    aa[i] = a[i];
    bb[i] = b[i];
    cc[i] = c[i];
  }
  vtkMath::MultiplyMatrix(aa, bb, 3, 3, 3, 3, cc);

  // WARNING: Number of columns of A must match number of rows of B.
  vtkMath::MultiplyMatrix(aa, bb, 3, 2, 3, 3, cc);

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

int TestTranspose3x3()
{
  int status = 0;
  std::cout << "Transpose3x3..";

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

int TestInvert3x3()
{
  int status = 0;
  std::cout << "Invert3x3..";
  {
    // Generate a Hilbert Matrix
    double mat[3][3];
    double matI[3][3];
    double expected[3][3] = { { 9.0, -36.0, 30.0 }, { -36.0, 192.0, -180.0 },
      { 30.0, -180.0, 180.0 } };

    for (int i = 1; i <= 3; ++i)
    {
      for (int j = 1; j <= 3; ++j)
      {
        mat[i - 1][j - 1] = 1.0 / (i + j - 1);
      }
    }
    vtkMath::Invert3x3(mat, matI);
    for (int i = 0; i < 3; ++i)
    {
      for (int j = 0; j < 3; ++j)
      {
        if (!vtkMathUtilities::FuzzyCompare(
              matI[i][j], expected[i][j], std::numeric_limits<double>::epsilon() * 16384.0))
        {
          std::cout << " Invert3x3(double) expected " << expected[i][j] << " but got "
                    << matI[i][j];
          ++status;
        }
      }
    }
  }
  {
    // Generate a Hilbert Matrix
    float mat[3][3];
    float matI[3][3];
    float expected[3][3] = { { 9.0f, -36.0f, 30.0f }, { -36.0f, 192.0f, -180.0f },
      { 30.0f, -180.0f, 180.0f } };

    for (int i = 1; i <= 3; ++i)
    {
      for (int j = 1; j <= 3; ++j)
      {
        mat[i - 1][j - 1] = 1.0 / (i + j - 1);
      }
    }
    vtkMath::Invert3x3(mat, matI);
    for (int i = 0; i < 3; ++i)
    {
      for (int j = 0; j < 3; ++j)
      {
        if (!vtkMathUtilities::FuzzyCompare(
              matI[i][j], expected[i][j], std::numeric_limits<float>::epsilon() * 8192.0f))
        {
          std::cout << " Invert3x3(single) expected " << expected[i][j] << " but got "
                    << matI[i][j];
          ++status;
        }
      }
    }
  }
  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

template <typename T, int NDimension>
int InvertMatrix()
{
  int status = 0;

  // Generate a Hilbert Matrix
  T** mat = new T*[NDimension];
  T** orig = new T*[NDimension];
  T** matI = new T*[NDimension];
  T** ident = new T*[NDimension];
  int* tmp1 = new int[NDimension];
  T* tmp2 = new T[NDimension];
  for (int i = 1; i <= NDimension; ++i)
  {
    mat[i - 1] = new T[NDimension];
    matI[i - 1] = new T[NDimension];
    orig[i - 1] = new T[NDimension];
    ident[i - 1] = new T[NDimension];
    for (int j = 1; j <= NDimension; ++j)
    {
      orig[i - 1][j - 1] = mat[i - 1][j - 1] = 1.0 / (i + j - 1);
    }
  }
  if (vtkMath::InvertMatrix(mat, matI, NDimension, tmp1, tmp2) == 0)
  {
    delete[] mat;
    delete[] orig;
    delete[] matI;
    delete[] ident;
    delete[] tmp1;
    delete[] tmp2;
    return status;
  }
  vtkMath::MultiplyMatrix(orig, matI, NDimension, NDimension, NDimension, NDimension, ident);

  T expected;
  for (int i = 0; i < NDimension; ++i)
  {
    for (int j = 0; j < NDimension; ++j)
    {
      if (i == j)
      {
        expected = 1.0;
      }
      else
      {
        expected = 0.0;
      }
      if (!vtkMathUtilities::FuzzyCompare(
            *(ident[i] + j), expected, std::numeric_limits<T>::epsilon() * (T)100000.0))
      {
        std::cout << " InvertMatrix(T) expected " << expected << " but got " << *(ident[i] + j);
        std::cout << "eps ratio is: "
                  << (*(ident[i] + j) - expected) / std::numeric_limits<T>::epsilon() << std::endl;
        ++status;
      }
    }
  }
  for (int i = 0; i < NDimension; i++)
  {
    delete[] mat[i];
    delete[] orig[i];
    delete[] matI[i];
    delete[] ident[i];
  }
  delete[] mat;
  delete[] orig;
  delete[] matI;
  delete[] ident;
  delete[] tmp1;
  delete[] tmp2;
  return status;
}
int TestInvertMatrix()
{
  int status = 0;
  status += InvertMatrix<double, 3>();
  status += InvertMatrix<double, 4>();
  status += InvertMatrix<double, 5>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

int TestIdentity3x3()
{
  int status = 0;
  std::cout << "Identity3x3..";

  double m[3][3];
  vtkMath::Identity3x3(m);

  double expected;
  for (int i = 0; i < 3; i++)
  {
    for (int j = 0; j < 3; j++)
    {
      if (i == j)
      {
        expected = 1.0;
      }
      else
      {
        expected = 0.0;
      }
      if (expected != m[i][j])
      {
        std::cout << " Identity expected " << expected << " but got " << m[i][j] << std::endl;
        ++status;
      }
    }
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

template <typename T>
int Matrix3x3ToQuaternion()
{
  int status = 0;

  T A[3][3];
  T quat[4];

  for (int n = 0; n < 1000; ++n)
  {
    for (int i = 0; i < 3; ++i)
    {
      for (int j = 0; j < 3; ++j)
      {
        A[i][j] = vtkMath::Random(-1.0, 1.0);
      }
    }
    vtkMath::Matrix3x3ToQuaternion(A, quat);
  }
  return status;
}
int TestMatrix3x3ToQuaternion()
{
  int status = 0;
  std::cout << "Matrix3x3ToQuaternion..";

  status += Matrix3x3ToQuaternion<double>();
  status += Matrix3x3ToQuaternion<float>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

template <typename T>
int QuaternionToMatrix3x3()
{
  int status = 0;

  T A[3][3];
  T quat[4];

  for (int n = 0; n < 1000; ++n)
  {
    quat[0] = vtkMath::Random(-vtkMath::Pi(), vtkMath::Pi());
    for (int i = 1; i < 4; ++i)
    {
      quat[i] = vtkMath::Random(-10.0, 10.0);
    }
    vtkMath::QuaternionToMatrix3x3(quat, A);
  }
  return status;
}

int TestQuaternionToMatrix3x3()
{
  int status = 0;
  std::cout << "QuaternionToMatrix3x3..";

  status += QuaternionToMatrix3x3<double>();
  status += QuaternionToMatrix3x3<float>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

template <typename T>
int MultiplyQuaternion()
{
  int status = 0;

  T q1[4];
  T q2[4];
  T q3[4];
  for (int n = 0; n < 1000; ++n)
  {
    q1[0] = vtkMath::Random(-vtkMath::Pi(), vtkMath::Pi());
    q2[0] = vtkMath::Random(-vtkMath::Pi(), vtkMath::Pi());
    vtkMath::MultiplyQuaternion(q1, q2, q3);
  }

  return status;
}
int TestMultiplyQuaternion()
{
  int status = 0;
  std::cout << "MultiplyQuaternion..";

  status += MultiplyQuaternion<double>();
  status += MultiplyQuaternion<float>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

int TestOrthogonalize3x3()
{
  int status = 0;
  std::cout << "Orthogonalize3x3..";

  {
    // Generate a random matrix
    double mat[3][3];
    double matO[3][3];
    double matI[3][3];

    for (int n = 0; n < 1000; ++n)
    {
      for (int i = 0; i < 3; ++i)
      {
        for (int j = 0; j < 3; ++j)
        {
          mat[i][j] = vtkMath::Random();
        }
      }
      vtkMath::Orthogonalize3x3(mat, matO);
      vtkMath::Transpose3x3(matO, mat);
      vtkMath::Multiply3x3(mat, matO, matI);

      double identity[3][3];
      vtkMath::Identity3x3(identity);
      for (int i = 0; i < 3; ++i)
      {
        for (int j = 0; j < 3; ++j)
        {
          if (!vtkMathUtilities::FuzzyCompare(
                matI[i][j], identity[i][j], std::numeric_limits<double>::epsilon() * 128.0))
          {
            std::cout << " Orthogonalize3x3 expected " << identity[i][j] << " but got "
                      << matI[i][j];
            ++status;
          }
        }
      }
    }
  }

  {
    // Generate a random matrix
    float mat[3][3];
    float matO[3][3];
    float matI[3][3];

    for (int n = 0; n < 1000; ++n)
    {
      for (int i = 0; i < 3; ++i)
      {
        for (int j = 0; j < 3; ++j)
        {
          mat[i][j] = vtkMath::Random();
        }
      }
      vtkMath::Orthogonalize3x3(mat, matO);
      vtkMath::Transpose3x3(matO, mat);
      vtkMath::Multiply3x3(mat, matO, matI);

      float identity[3][3];
      vtkMath::Identity3x3(identity);
      for (int i = 0; i < 3; ++i)
      {
        for (int j = 0; j < 3; ++j)
        {
          if (!vtkMathUtilities::FuzzyCompare(
                matI[i][j], identity[i][j], std::numeric_limits<float>::epsilon() * 128.0f))
          {
            std::cout << " Orthogonalize3x3 expected " << identity[i][j] << " but got "
                      << matI[i][j];
            ++status;
          }
        }
      }
    }
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

template <typename T>
int Diagonalize3x3()
{
  int status = 0;

  T mat[3][3];
  T eigenVector[3][3], eigenVectorT[3][3];
  T temp[3][3];
  T result[3][3];
  T eigen[3];

  for (int n = 0; n < 0; ++n)
  {
    for (int i = 0; i < 3; ++i)
    {
      for (int j = i; j < 3; ++j)
      {
        mat[i][j] = mat[j][i] = vtkMath::Random(-1.0, 1.0);
      }
    }

    vtkMath::Diagonalize3x3(mat, eigen, eigenVector);

    // Pt * A * P = diagonal matrix with eigenvalues on diagonal
    vtkMath::Multiply3x3(mat, eigenVector, temp);
    vtkMath::Invert3x3(eigenVector, eigenVectorT);
    vtkMath::Multiply3x3(eigenVectorT, temp, result);
    T expected;
    for (int i = 0; i < 3; ++i)
    {
      for (int j = 0; j < 3; ++j)
      {
        if (i == j)
        {
          expected = eigen[i];
        }
        else
        {
          expected = 0.0;
        }
        if (!vtkMathUtilities::FuzzyCompare(
              result[i][j], expected, std::numeric_limits<T>::epsilon() * (T)128.0))
        {
          std::cout << " Diagonalize3x3 expected " << expected << " but got " << result[i][j];
          ++status;
        }
      }
    }
  }

  // Now test for 2 and 3 equal eigenvalues
  vtkMath::Identity3x3(mat);
  mat[0][0] = 5.0;
  mat[1][1] = 5.0;
  mat[2][2] = 1.0;

  vtkMath::Diagonalize3x3(mat, eigen, eigenVector);
  std::cout << "eigen: " << eigen[0] << "," << eigen[1] << "," << eigen[2] << std::endl;

  vtkMath::Identity3x3(mat);
  mat[0][0] = 2.0;
  mat[1][1] = 2.0;
  mat[2][2] = 2.0;

  vtkMath::Diagonalize3x3(mat, eigen, eigenVector);
  std::cout << "eigen: " << eigen[0] << "," << eigen[1] << "," << eigen[2] << std::endl;
  return status;
}

// Validate Pt * A * P = diagonal matrix with eigenvalues on diagonal
int TestDiagonalize3x3()
{
  int status = 0;
  std::cout << "Diagonalize3x3..";

  status += Diagonalize3x3<double>();
  status += Diagonalize3x3<float>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

template <typename T>
int SingularValueDecomposition3x3()
{
  int status = 0;

  T a[3][3];
  T orig[3][3];
  T u[3][3];
  T w[3];
  T vt[3][3];

  for (int n = 0; n < 1000; ++n)
  {
    for (int i = 0; i < 3; ++i)
    {
      for (int j = 0; j < 3; ++j)
      {
        orig[i][j] = a[i][j] = vtkMath::Random(-10.0, 10.0);
      }
    }
    vtkMath::SingularValueDecomposition3x3(a, u, w, vt);

    T m[3][3];
    T W[3][3];
    vtkMath::Identity3x3(W);
    W[0][0] = w[0];
    W[1][1] = w[1];
    W[2][2] = w[2];
    vtkMath::Multiply3x3(u, W, m);
    vtkMath::Multiply3x3(m, vt, m);

    for (int i = 0; i < 3; ++i)
    {
      for (int j = 0; j < 3; ++j)
      {
        if (!vtkMathUtilities::FuzzyCompare(
              m[i][j], orig[i][j], std::numeric_limits<T>::epsilon() * (T)128.0))
        {
          std::cout << " SingularValueDecomposition3x3 expected " << orig[i][j] << " but got "
                    << m[i][j];
          std::cout << " eps ratio is: "
                    << (m[i][j] - orig[i][j]) / std::numeric_limits<T>::epsilon() << std::endl;
          ++status;
        }
      }
    }
  }
  return status;
}
// Validate u * w * vt = m
int TestSingularValueDecomposition3x3()
{
  int status = 0;
  std::cout << "SingularValueDecomposition3x3..";

  status += SingularValueDecomposition3x3<double>();
  status += SingularValueDecomposition3x3<float>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

template <typename T, int NDimension>
int SolveLinearSystem()
{
  int status = 0;

  for (int n = 0; n < 100; ++n)
  {
    // Generate a Random Matrix
    T** mat = new T*[NDimension];
    T* lhs = new T[NDimension];
    T* rhs = new T[NDimension];

    for (int i = 0; i < NDimension; ++i)
    {
      mat[i] = new T[NDimension];
      lhs[i] = vtkMath::Random(-1.0, 1.0);
      for (int j = 0; j < NDimension; ++j)
      {
        *(mat[i] + j) = vtkMath::Random(-1.0, 1.0);
      }
    }

    for (int i = 0; i < NDimension; ++i)
    {
      rhs[i] = 0.0;
      for (int j = 0; j < NDimension; ++j)
      {
        rhs[i] += (*(mat[i] + j) * lhs[j]);
      }
    }
    vtkMath::SolveLinearSystem(mat, rhs, NDimension);

    for (int i = 0; i < NDimension; ++i)
    {
      if (!vtkMathUtilities::FuzzyCompare(
            lhs[i], rhs[i], std::numeric_limits<double>::epsilon() * 32768.0))
      {
        std::cout << " SolveLinearSystem(double) expected " << lhs[i] << " but got " << rhs[i];
        std::cout << " eps ratio is: " << (lhs[i] - rhs[i]) / std::numeric_limits<T>::epsilon()
                  << std::endl;
        ++status;
      }
    }

    if (NDimension == 1 || NDimension == 2)
    {
      for (int i = 0; i < NDimension; ++i)
      {
        for (int j = 0; j < NDimension; ++j)
        {
          *(mat[i] + j) = 0.0;
        }
      }
      if (vtkMath::SolveLinearSystem(mat, rhs, NDimension) != 0.0)
      {
        std::cout << " SolveLinearSystem for a zero matrix expected " << 0 << " but got 1";
        ++status;
      }
    }
    for (int i = 0; i < NDimension; i++)
    {
      delete[] mat[i];
    }
    delete[] mat;
    delete[] rhs;
    delete[] lhs;
  }
  return status;
}

// Validate with a known left hand side
int TestSolveLinearSystem()
{
  int status = 0;
  std::cout << "SolveLinearSystem..";

  status += SolveLinearSystem<double, 1>();
  status += SolveLinearSystem<double, 2>();
  status += SolveLinearSystem<double, 3>();
  status += SolveLinearSystem<double, 50>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with a known solution
int TestSolveLeastSquares()
{
  int status = 0;
  std::cout << "SolveLeastSquares..";

  double** m = new double*[2];
  double** x = new double*[3];
  double** y = new double*[3];

  for (int i = 0; i < 2; ++i)
  {
    m[i] = new double[1];
  }
  for (int i = 0; i < 3; ++i)
  {
    x[i] = new double[2];
  }
  x[0][0] = 1;
  x[0][1] = 4;
  x[1][0] = 1;
  x[1][1] = 2;
  x[2][0] = 2;
  x[2][1] = 3;

  for (int i = 0; i < 3; ++i)
  {
    y[i] = new double[1];
  }
  y[0][0] = -2;
  y[1][0] = 6;
  y[2][0] = 1;

  vtkMath::SolveLeastSquares(3, x, 2, y, 1, m);

  std::vector<double> results;
  std::vector<double> expecteds;
  expecteds.push_back(3.0);
  results.push_back(m[0][0]);
  expecteds.push_back(-1.0);
  results.push_back(m[1][0]);

  for (size_t i = 0; i < results.size(); ++i)
  {
    if (!vtkMathUtilities::FuzzyCompare(
          results[i], expecteds[i], std::numeric_limits<double>::epsilon() * 128.0))
    {
      std::cout << " Solve Least Squares got " << results[i] << " but expected " << expecteds[i];
      ++status;
    }
  }

  // Now make one solution homogeous
  y[0][0] = 0.0;
  vtkMath::SolveLeastSquares(3, x, 2, y, 1, m);

  // Now make all homogeous
  y[0][0] = 0.0;
  y[1][0] = 0.0;
  y[2][0] = 0.0;
  vtkMath::SolveLeastSquares(3, x, 2, y, 1, m);

  // Insufficient number of samples. Underdetermined.
  if (vtkMath::SolveLeastSquares(1, x, 2, y, 1, m) != 0)
  {
    std::cout << " Solve Least Squares got " << 1 << " but expected " << 0;
    ++status;
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }

  for (int i = 0; i < 3; i++)
  {
    delete[] x[i];
    delete[] y[i];
  }
  for (int i = 0; i < 2; i++)
  {
    delete[] m[i];
  }
  delete[] x;
  delete[] y;
  delete[] m;

  return status;
}

// Only warning cases validate
// No validation, just coverage
int TestSolveHomogeneousLeastSquares()
{
  int status = 0;
  std::cout << "SolveHomogeneousLeastSquares..";

  double** m = new double*[2];
  double** x = new double*[3];
  double** y = new double*[3];

  for (int i = 0; i < 2; ++i)
  {
    m[i] = new double[1];
  }
  for (int i = 0; i < 3; ++i)
  {
    x[i] = new double[2];
    y[i] = new double[1];
  }
  x[0][0] = 1;
  x[0][1] = 2;
  x[1][0] = 2;
  x[1][1] = 4;
  x[2][0] = 3;
  x[2][1] = 6;

  vtkMath::SolveHomogeneousLeastSquares(3, x, 1, m);

  vtkMath::MultiplyMatrix(x, m, 3, 2, 2, 1, y);

  std::vector<double> results;
  std::vector<double> expecteds;

  for (size_t i = 0; i < results.size(); ++i)
  {
    if (!vtkMathUtilities::FuzzyCompare(
          results[i], expecteds[i], std::numeric_limits<double>::epsilon() * 128.0))
    {
      std::cout << " SolveHomogeneousLeastSquares got " << results[i] << " but expected "
                << expecteds[i];
      ++status;
    }
  }

  // Insufficient number of samples. Underdetermined.
  if (vtkMath::SolveHomogeneousLeastSquares(3, x, 5, m) != 0)
  {
    std::cout << " SolveHomogeneousLeastSquares got " << 1 << " but expected " << 0;
    ++status;
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }

  for (int i = 0; i < 3; i++)
  {
    delete[] x[i];
    delete[] y[i];
  }
  for (int i = 0; i < 2; i++)
  {
    delete[] m[i];
  }
  delete[] x;
  delete[] y;
  delete[] m;

  return status;
}

template <typename T, int NDimension>
int LUSolveLinearSystemEstimateMatrixCondition()
{
  int status = 0;

  // Generate a Hilbert Matrix
  T** mat = new T*[NDimension];
  int index[NDimension];

  for (int i = 1; i <= NDimension; ++i)
  {
    mat[i - 1] = new T[NDimension];
    for (int j = 1; j <= NDimension; ++j)
    {
      mat[i - 1][j - 1] = 1.0 / (i + j - 1);
    }
  }
  vtkMath::LUFactorLinearSystem(mat, index, NDimension);
  T condition = vtkMath::EstimateMatrixCondition(mat, NDimension);
  std::cout << "Condition is: " << condition << std::endl;

  T expected = condition;
  if (!vtkMathUtilities::FuzzyCompare(
        condition, expected, std::numeric_limits<T>::epsilon() * (T)128.0))
  {
    std::cout << " EstimateMatrixCondition(T) expected " << expected << " but got " << condition;
    std::cout << "eps ratio is: " << condition - expected / std::numeric_limits<T>::epsilon()
              << std::endl;
    ++status;
  }
  for (int i = 0; i < NDimension; i++)
  {
    delete[] mat[i];
  }
  delete[] mat;
  return status;
}

// Validate by observation that the condition of a hilbert matrix
// increases with dimension
int TestLUSolveLinearSystemEstimateMatrixCondition()
{
  int status = 0;
  std::cout << "LUSolveLinearSystemEstimateMatrixCondition..";
  status += LUSolveLinearSystemEstimateMatrixCondition<double, 10>();
  status += LUSolveLinearSystemEstimateMatrixCondition<double, 8>();
  status += LUSolveLinearSystemEstimateMatrixCondition<double, 6>();
  status += LUSolveLinearSystemEstimateMatrixCondition<double, 4>();
  status += LUSolveLinearSystemEstimateMatrixCondition<double, 3>();
  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

template <typename T, int NDimension>
int JacobiN()
{
  int status = 0;

  T mat[NDimension][NDimension];
  T orig[NDimension][NDimension];
  T eigenVector[NDimension][NDimension], eigenVectorT[NDimension][NDimension];
  T temp[NDimension][NDimension];
  T result[NDimension][NDimension];
  T eigen[NDimension];

  for (int n = 0; n < 10; ++n)
  {
    for (int i = 0; i < NDimension; ++i)
    {
      for (int j = i; j < NDimension; ++j)
      {
        mat[i][j] = mat[j][i] = vtkMath::Random(0.0, 1.0);
        orig[i][j] = orig[j][i] = mat[i][j];
      }
    }

    // convert to jacobiN format
    T* origJ[NDimension];
    T* matJ[NDimension];
    T* eigenVectorJ[NDimension];
    T* eigenVectorTJ[NDimension];
    T* resultJ[NDimension];
    T* tempJ[NDimension];
    for (int i = 0; i < NDimension; ++i)
    {
      matJ[i] = mat[i];
      origJ[i] = orig[i];
      eigenVectorJ[i] = eigenVector[i];
      eigenVectorTJ[i] = eigenVectorT[i];
      tempJ[i] = temp[i];
      resultJ[i] = result[i];
    }

    if (NDimension == 3)
    {
      vtkMath::Jacobi(matJ, eigen, eigenVectorJ);
    }
    else
    {
      vtkMath::JacobiN(matJ, NDimension, eigen, eigenVectorJ);
    }

    // P^-1 * A * P = diagonal matrix with eigenvalues on diagonal
    vtkMath::MultiplyMatrix(
      origJ, eigenVectorJ, NDimension, NDimension, NDimension, NDimension, tempJ);
    vtkMath::InvertMatrix(eigenVectorJ, eigenVectorTJ, NDimension);
    vtkMath::MultiplyMatrix(
      eigenVectorTJ, tempJ, NDimension, NDimension, NDimension, NDimension, resultJ);
    T expected;
    for (int i = 0; i < NDimension; ++i)
    {
      for (int j = 0; j < NDimension; ++j)
      {
        if (i == j)
        {
          expected = eigen[i];
        }
        else
        {
          expected = 0.0;
        }
        if (!vtkMathUtilities::FuzzyCompare(
              result[i][j], expected, std::numeric_limits<T>::epsilon() * (T)256.0))
        {
          std::cout << " JacobiN expected " << expected << " but got " << result[i][j];
          std::cout << "eps ratio is: "
                    << (result[i][j] - expected) / std::numeric_limits<T>::epsilon() << std::endl;
          ++status;
        }
      }
    }
  }
  return status;
}

// Validate P^-1 * A * P = diagonal matrix with eigenvalues on diagonal
int TestJacobiN()
{
  int status = 0;
  std::cout << "JacobiN..";
  status += JacobiN<double, 3>();
  status += JacobiN<double, 10>();
  status += JacobiN<double, 50>();

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

template <typename T>
int RGBToHSV()
{
  int status = 0;
  T R, G, B;
  T H, S, V;
  T CR, CG, CB;
  for (int n = 0; n < 1000; ++n)
  {
    std::vector<T> values;
    std::vector<T> expecteds;
    R = (T)vtkMath::Random(0.0, 1.0);
    G = (T)vtkMath::Random(0.0, 1.0);
    B = (T)vtkMath::Random(0.0, 1.0);

    vtkMath::RGBToHSV(R, G, B, &H, &S, &V);
    vtkMath::HSVToRGB(H, S, V, &CR, &CG, &CB);
    values.push_back(CR);
    values.push_back(CG);
    values.push_back(CB);
    expecteds.push_back(R);
    expecteds.push_back(G);
    expecteds.push_back(B);

    for (size_t i = 0; i < values.size(); ++i)
    {
      if (!vtkMathUtilities::FuzzyCompare(
            values[i], expecteds[i], std::numeric_limits<T>::epsilon() * (T)128.0))
      {
        std::cout << " RGBToHSV got " << values[i] << " but expected " << expecteds[i];
        std::cout << " eps ratio is: "
                  << (values[i] - expecteds[i]) / std::numeric_limits<T>::epsilon() << std::endl;
        ++status;
      }
    }
  }
  return status;
}

// Validate by rgb->hsv->rgb
int TestRGBToHSV()
{
  int status = 0;
  std::cout << "RGBToHSV..";

  status += RGBToHSV<double>();
  status += RGBToHSV<float>();
  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with known solutions
int TestClampValue()
{
  int status = 0;
  std::cout << "ClampValue..";

  double value;
  double clampedValue;
  double range[2] = { -1.0, 1.0 };

  value = -800.0;
  clampedValue = vtkMath::ClampValue(value, range[0], range[1]);
  if (clampedValue != range[0])
  {
    std::cout << " ClampValue expected " << range[0] << " but got " << value;
    ++status;
  }

  value = 900.0;
  clampedValue = vtkMath::ClampValue(value, range[0], range[1]);
  if (clampedValue != range[1])
  {
    std::cout << " ClampValue expected " << range[1] << " but got " << value;
    ++status;
  }

  value = 0.0;
  clampedValue = vtkMath::ClampValue(value, range[0], range[1]);
  if (clampedValue != 0.0)
  {
    std::cout << " ClampValue expected " << 0.0 << " but got " << value;
    ++status;
  }

  value = -100.0;
  vtkMath::ClampValue(&value, range);
  if (value != range[0])
  {
    std::cout << " ClampValue expected " << range[0] << " but got " << value;
    ++status;
  }
  value = 100.0;
  vtkMath::ClampValue(&value, range);
  if (value != range[1])
  {
    std::cout << " ClampValue expected " << range[1] << " but got " << value;
    ++status;
  }
  value = -100.0;
  vtkMath::ClampValue(value, range, &clampedValue);
  if (clampedValue != range[0])
  {
    std::cout << " ClampValue expected " << range[0] << " but got " << clampedValue;
    ++status;
  }

  value = 100.0;
  vtkMath::ClampValue(value, range, &clampedValue);
  if (clampedValue != range[1])
  {
    std::cout << " ClampValue expected " << range[1] << " but got " << clampedValue;
    ++status;
  }

  value = 0.0;
  vtkMath::ClampValue(value, range, &clampedValue);
  if (clampedValue != value)
  {
    std::cout << " ClampValue expected " << value << " but got " << clampedValue;
    ++status;
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with known solutions
int TestClampValues()
{
  int status = 0;
  std::cout << "ClampValues..";

  double values[1000];
  double clampedValues[1000];
  for (int n = 0; n < 1000; ++n)
  {
    values[n] = vtkMath::Random(-2.0, 2.0);
  }
  double range[2] = { -1.0, 1.0 };
  vtkMath::ClampValues(values, 1000, range, clampedValues);
  vtkMath::ClampValues(values, 1000, range);

  for (int n = 0; n < 1000; ++n)
  {
    if (values[n] != clampedValues[n])
    {
      ++status;
    }
  }

  vtkMath::ClampValues(nullptr, 1000, nullptr);
  vtkMath::ClampValues(nullptr, 1000, nullptr, nullptr);

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with known solutions
int TestClampAndNormalizeValue()
{
  int status = 0;
  std::cout << "ClampAndNormalizeValue..";

  double value;
  double result;
  double range[3] = { -1.0, 1.0 };

  value = -100.0;
  result = vtkMath::ClampAndNormalizeValue(value, range);
  if (result != 0.0)
  {
    std::cout << " ClampAndNormalizeValue expected " << 0.0 << " but got " << result;
    ++status;
  }
  value = 100.0;
  result = vtkMath::ClampAndNormalizeValue(value, range);
  if (result != 1.0)
  {
    std::cout << " ClampAndNormalizeValue expected " << 1.0 << " but got " << result;
    ++status;
  }

  range[0] = 0.0;
  range[1] = 1.0;
  value = 0.5;
  result = vtkMath::ClampAndNormalizeValue(value, range);
  if (result != 0.5)
  {
    std::cout << " ClampValue expected " << 0.5 << " but got " << result;
    ++status;
  }

  range[0] = 1.0;
  range[1] = 1.0;
  value = 1.0;
  result = vtkMath::ClampAndNormalizeValue(value, range);
  if (result != 0.0)
  {
    std::cout << " ClampValue expected " << 0.0 << " but got " << result;
    ++status;
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate by checking symmetric tensor values
// are in correct places
int TestTensorFromSymmetricTensor()
{
  int status = 0;
  std::cout << "TensorFromSymmetricTensor..";
  double symmTensor[9];
  for (int i = 0; i < 6; i++)
  {
    symmTensor[i] = vtkMath::Random();
  }
  double tensor[9];
  vtkMath::TensorFromSymmetricTensor(symmTensor, tensor);
  if (tensor[0] != symmTensor[0] || tensor[1] != symmTensor[3] || tensor[2] != symmTensor[5] ||
    tensor[3] != symmTensor[3] || tensor[4] != symmTensor[1] || tensor[5] != symmTensor[4] ||
    tensor[6] != symmTensor[5] || tensor[7] != symmTensor[4] || tensor[8] != symmTensor[2])
  {
    std::cout << " Unexpected results from TensorFromSymmetricTensor " << std::endl;
    ++status;
  }

  vtkMath::TensorFromSymmetricTensor(symmTensor);
  for (int i = 0; i < 9; i++)
  {
    if (symmTensor[i] != tensor[i])
    {
      std::cout << " Unexpected results from in place TensorFromSymmetricTensor " << std::endl;
      ++status;
      break;
    }
  }
  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate by checking ranges with numeric_limits
int TestGetScalarTypeFittingRange()
{
  int status = 0;
  std::cout << "GetScalarTypeFittingRange..";

  double rangeMin;
  double rangeMax;

  rangeMin = (double)std::numeric_limits<char>::min();
  rangeMax = (double)std::numeric_limits<char>::max();
  if (vtkMath::GetScalarTypeFittingRange(rangeMin, rangeMax, 1.0, 0.0) != VTK_CHAR)
  {
    std::cout << " Bad fitting range for VTK_CHAR" << std::endl;
    ++status;
  }

#if CHAR_MIN < 0
  auto expected_schar = VTK_CHAR;
  auto expected_uchar = VTK_UNSIGNED_CHAR;
#else
  auto expected_schar = VTK_SIGNED_CHAR;
  auto expected_uchar = VTK_CHAR;
#endif

  rangeMin = (double)std::numeric_limits<signed char>::min();
  rangeMax = (double)std::numeric_limits<signed char>::max();
  if (vtkMath::GetScalarTypeFittingRange(rangeMin, rangeMax, 1.0, 0.0) != expected_schar)
  {
    std::cout << " Bad fitting range for VTK_CHAR" << std::endl;
    ++status;
  }

  rangeMin = (double)std::numeric_limits<unsigned char>::min();
  rangeMax = (double)std::numeric_limits<unsigned char>::max();
  if (vtkMath::GetScalarTypeFittingRange(rangeMin, rangeMax, 1.0, 0.0) != expected_uchar)
  {
    std::cout << " Bad fitting range for VTK_UNSIGNED_CHAR " << std::endl;
    ++status;
  }

  rangeMin = (double)std::numeric_limits<short>::min();
  rangeMax = (double)std::numeric_limits<short>::max();
  if (vtkMath::GetScalarTypeFittingRange(rangeMin, rangeMax, 1.0, 0.0) != VTK_SHORT)
  {
    std::cout << " Bad fitting range for VTK_SHORT" << std::endl;
    ++status;
  }

  rangeMin = (double)std::numeric_limits<unsigned short>::min();
  rangeMax = (double)std::numeric_limits<unsigned short>::max();
  if (vtkMath::GetScalarTypeFittingRange(rangeMin, rangeMax, 1.0, 0.0) != VTK_UNSIGNED_SHORT)
  {
    std::cout << " Bad fitting range for VTK_UNSIGNED_SHORT" << std::endl;
    ++status;
  }

  rangeMin = (double)std::numeric_limits<int>::min();
  rangeMax = (double)std::numeric_limits<int>::max();
  if (vtkMath::GetScalarTypeFittingRange(rangeMin, rangeMax, 1.0, 0.0) != VTK_INT)
  {
    std::cout << " Bad fitting range for VTK_INT" << std::endl;
    ++status;
  }

  rangeMin = (double)std::numeric_limits<unsigned int>::min();
  rangeMax = (double)std::numeric_limits<unsigned int>::max();
  if (vtkMath::GetScalarTypeFittingRange(rangeMin, rangeMax, 1.0, 0.0) != VTK_UNSIGNED_INT)
  {
    std::cout << " Bad fitting range for VTK_UNSIGNED_INT" << std::endl;
    ++status;
  }

  rangeMin = (double)std::numeric_limits<long>::min();
  rangeMax = (double)std::numeric_limits<long>::max();
  int scalarType = vtkMath::GetScalarTypeFittingRange(rangeMin, rangeMax, 1.0, 0.0);
  if (sizeof(long) == sizeof(int))
  {
    if (scalarType != VTK_INT)
    {
      std::cout << " Bad fitting range for VTK_LONG" << std::endl;
      std::cout << " Expected " << VTK_INT << " but got " << scalarType;
      ++status;
    }
  }
  else
  {
    if (scalarType != VTK_LONG)
    {
      std::cout << " Bad fitting range for VTK_LONG" << std::endl;
      std::cout << " Expected " << VTK_LONG << " but got " << scalarType;
      ++status;
    }
  }

  rangeMin = (double)std::numeric_limits<unsigned long>::min();
  rangeMax = (double)std::numeric_limits<unsigned long>::max();
  scalarType = vtkMath::GetScalarTypeFittingRange(rangeMin, rangeMax, 1.0, 0.0);
  if (sizeof(unsigned long) == sizeof(unsigned int))
  {
    if (scalarType != VTK_UNSIGNED_INT)
    {
      std::cout << " Bad fitting range for VTK_UNSIGNED_LONG" << std::endl;
      std::cout << " Expected " << VTK_UNSIGNED_INT << " but got " << scalarType;
      ++status;
    }
  }
  else
  {
    if (scalarType != VTK_UNSIGNED_LONG)
    {
      std::cout << " Bad fitting range for VTK_UNSIGNED_LONG" << std::endl;
      std::cout << " Expected " << VTK_UNSIGNED_LONG << " but got " << scalarType;
      ++status;
    }
  }

  rangeMin = (double)std::numeric_limits<short>::min();
  rangeMax = (double)std::numeric_limits<short>::max();
  if (vtkMath::GetScalarTypeFittingRange(rangeMin, rangeMax, 1.1, 0.0) != VTK_FLOAT)
  {
    std::cout << " Bad fitting range for VTK_FLOAT" << std::endl;
    ++status;
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with known solutions
int TestGetAdjustedScalarRange()
{
  int status = 0;
  std::cout << "GetAdjustedScalarRange..";

  vtkSmartPointer<vtkUnsignedCharArray> uc = vtkSmartPointer<vtkUnsignedCharArray>::New();
  uc->SetNumberOfComponents(3);
  uc->SetNumberOfTuples(100);
  for (int i = 0; i < 100; ++i)
  {
    for (int j = 0; j < 3; ++j)
    {
      uc->SetComponent(i, j,
        vtkMath::Random(
          std::numeric_limits<unsigned char>::min(), std::numeric_limits<unsigned char>::max()));
    }
  }
  double range[2];
  vtkMath::GetAdjustedScalarRange(uc, 1, range);
  if (range[0] != uc->GetDataTypeMin() || range[1] != uc->GetDataTypeMax())
  {
    std::cout << " GetAdjustedScalarRange(unsigned char) expected " << uc->GetDataTypeMin() << ", "
              << uc->GetDataTypeMax() << " but got " << range[0] << ", " << range[1] << std::endl;
    ++status;
  }

  vtkSmartPointer<vtkUnsignedShortArray> us = vtkSmartPointer<vtkUnsignedShortArray>::New();
  us->SetNumberOfComponents(3);
  us->SetNumberOfTuples(10000);
  for (int i = 0; i < 10000; ++i)
  {
    us->SetComponent(i, 0,
      vtkMath::Random(
        std::numeric_limits<unsigned short>::min(), std::numeric_limits<unsigned short>::max()));
    us->SetComponent(i, 1,
      vtkMath::Random(std::numeric_limits<unsigned short>::min(),
        std::numeric_limits<unsigned char>::max() + 100));
    us->SetComponent(i, 2,
      vtkMath::Random(
        std::numeric_limits<unsigned short>::min(), std::numeric_limits<unsigned char>::max()));
  }
  vtkMath::GetAdjustedScalarRange(us, 0, range);
  if (range[0] != us->GetDataTypeMin() || range[1] != us->GetDataTypeMax())
  {
    std::cout << " GetAdjustedScalarRange(unsigned short) expected " << us->GetDataTypeMin() << ", "
              << us->GetDataTypeMax() << " but got " << range[0] << ", " << range[1] << std::endl;
    ++status;
  }

  vtkMath::GetAdjustedScalarRange(us, 1, range);
  if (range[0] != us->GetDataTypeMin() || range[1] != 4095.0)
  {
    std::cout << " GetAdjustedScalarRange(unsigned short) expected " << us->GetDataTypeMin() << ", "
              << 4095.0 << " but got " << range[0] << ", " << range[1] << std::endl;
    ++status;
  }

  vtkMath::GetAdjustedScalarRange(us, 2, range);
  if (range[0] != us->GetDataTypeMin() || range[1] >= uc->GetDataTypeMax())
  {
    std::cout << " GetAdjustedScalarRange(unsigned short) expected " << us->GetDataTypeMin() << ", "
              << ">= " << uc->GetDataTypeMax() << " but got " << range[0] << ", " << range[1]
              << std::endl;
    ++status;
  }

  // Test nullptr array
  if (vtkMath::GetAdjustedScalarRange(nullptr, 1000, nullptr))
  {
    std::cout << " GetAdjustedScalarRange with a nullptr array expected " << 0 << " but got " << 1
              << std::endl;
    ++status;
  }
  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with known solutions
int TestExtentIsWithinOtherExtent()
{
  int status = 0;
  std::cout << "ExtentIsWithinOtherExtent..";

  if (vtkMath::ExtentIsWithinOtherExtent(nullptr, nullptr))
  {
    std::cout << " ExtentIsWithinOtherExtent expected 0 but got 1" << std::endl;
    ++status;
  }

  int extent1[6];
  int extent2[6];
  extent1[0] = 100;
  extent1[1] = 101;
  extent1[2] = 100;
  extent1[3] = 101;
  extent1[4] = 100;
  extent1[5] = 101;

  extent2[0] = 100;
  extent2[1] = 101;
  extent2[2] = 100;
  extent2[3] = 101;
  extent2[4] = 100;
  extent2[5] = 101;

  if (!vtkMath::ExtentIsWithinOtherExtent(extent1, extent2))
  {
    std::cout << " ExtentIsWithinOtherExtent expected 1 but got 0" << std::endl;
    ++status;
  }

  extent1[0] = 99;
  extent1[1] = 101;
  if (vtkMath::ExtentIsWithinOtherExtent(extent1, extent2))
  {
    std::cout << " ExtentIsWithinOtherExtent expected 0 but got 1" << std::endl;
    ++status;
  }

  extent1[0] = 98;
  extent1[1] = 99;
  if (vtkMath::ExtentIsWithinOtherExtent(extent1, extent2))
  {
    std::cout << " ExtentIsWithinOtherExtent expected 0 but got 1" << std::endl;
    ++status;
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with known solutions
int TestBoundsIsWithinOtherBounds()
{
  int status = 0;
  std::cout << "BoundsIsWithinOtherBounds..";

  if (vtkMath::BoundsIsWithinOtherBounds(nullptr, nullptr, nullptr))
  {
    std::cout << " BoundsIsWithinOtherBounds expected 0 but got 1" << std::endl;
    ++status;
  }

  double delta[3];
  delta[0] = delta[1] = delta[2] = std::numeric_limits<double>::epsilon();

  double bounds1[6];
  double bounds2[6];
  bounds1[0] = 1 - delta[0];
  bounds1[1] = 2 + delta[0];
  bounds1[2] = 1;
  bounds1[3] = 2;
  bounds1[4] = 1;
  bounds1[5] = 2;

  bounds2[0] = 1;
  bounds2[1] = 2;
  bounds2[2] = 1;
  bounds2[3] = 2;
  bounds2[4] = 1;
  bounds2[5] = 2;

  if (!vtkMath::BoundsIsWithinOtherBounds(bounds1, bounds2, delta))
  {
    std::cout << " BoundsIsWithinOtherBounds expected 1 but got 0" << std::endl;
    ++status;
  }

  bounds1[0] = 1 - 2.0 * delta[0];
  bounds1[1] = 2 + 2.0 * delta[0];
  if (vtkMath::BoundsIsWithinOtherBounds(bounds1, bounds2, delta))
  {
    std::cout << " BoundsIsWithinOtherBounds expected 0 but got 1" << std::endl;
    ++status;
  }

  bounds1[0] = 1 - 4.0 * delta[0];
  bounds1[1] = 1 - 2.0 * delta[0];
  if (vtkMath::BoundsIsWithinOtherBounds(bounds1, bounds2, delta))
  {
    std::cout << " BoundsIsWithinOtherBounds expected 0 but got 1" << std::endl;
    ++status;
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with known solutions
int TestPointIsWithinBounds()
{
  int status = 0;
  std::cout << "PointIsWithinBounds..";

  if (vtkMath::PointIsWithinBounds(nullptr, nullptr, nullptr))
  {
    std::cout << " PointIsWithinBounds expected 0 but got 1" << std::endl;
    ++status;
  }

  double delta[3];
  delta[0] = std::numeric_limits<double>::epsilon();
  delta[1] = std::numeric_limits<double>::epsilon() * 2.0;
  delta[2] = std::numeric_limits<double>::epsilon() * 256.0;

  double bounds1[6];
  bounds1[0] = 1.0;
  bounds1[1] = 2.0;
  bounds1[2] = 1.0;
  bounds1[3] = 2.0;
  bounds1[4] = 1.0;
  bounds1[5] = 2.0;

  double point[3];
  point[0] = bounds1[0] - delta[0];
  point[1] = bounds1[2] - delta[1];
  point[2] = bounds1[4];

  if (!vtkMath::PointIsWithinBounds(point, bounds1, delta))
  {
    std::cout << " PointIsWithinBounds expected 1 but got 0" << std::endl;
    ++status;
  }

  point[0] = bounds1[0] - delta[0];
  point[1] = bounds1[2] - delta[1];
  point[2] = bounds1[4] - 2.0 * delta[2];

  if (vtkMath::PointIsWithinBounds(point, bounds1, delta))
  {
    std::cout << " PointIsWithinBounds expected 0 but got 1" << std::endl;
    ++status;
  }

  point[0] = bounds1[1] + delta[0];
  point[1] = bounds1[3] + delta[1];
  point[2] = bounds1[5] + 2.0 * delta[2];

  if (vtkMath::PointIsWithinBounds(point, bounds1, delta))
  {
    std::cout << " PointIsOtherBounds expected 0 but got 1" << std::endl;
    ++status;
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with with alternative solution
int TestSolve3PointCircle()
{
  int status = 0;
  std::cout << "Solve3PointCircle..";

  for (int n = 0; n < 1000; ++n)
  {
    double A[3], B[3], C[3];
    double center[3];
    double a[3], b[3], aMinusb[3], aCrossb[3];

    for (int i = 0; i < 3; ++i)
    {
      A[i] = vtkMath::Random(-1.0, 1.0);
      B[i] = vtkMath::Random(-1.0, 1.0);
      C[i] = vtkMath::Random(-1.0, 1.0);
    }

    vtkMath::Subtract(A, C, a);
    vtkMath::Subtract(B, C, b);
    vtkMath::Subtract(a, b, aMinusb);
    vtkMath::Cross(a, b, aCrossb);

    double expectedRadius;
    expectedRadius = (vtkMath::Norm(a) * vtkMath::Norm(b) * vtkMath::Norm(aMinusb)) /
      (2.0 * vtkMath::Norm(aCrossb));

    double radius;
    radius = vtkMath::Solve3PointCircle(A, B, C, center);
    if (!vtkMathUtilities::FuzzyCompare(
          radius, expectedRadius, std::numeric_limits<double>::epsilon() * 1024.0))
    {
      std::cout << " Solve3PointCircle radius expected " << expectedRadius << " but got " << radius;
      std::cout << "eps ratio is: "
                << (expectedRadius - radius) / std::numeric_limits<double>::epsilon() << std::endl;
      ++status;
    }

    double ab[3], ba[3];
    double abMinusba[3];
    double abMinusbaCrossaCrossb[3];

    vtkMath::Subtract(B, C, ab);
    vtkMath::Subtract(A, C, ba);
    vtkMath::MultiplyScalar(ab, vtkMath::Norm(a) * vtkMath::Norm(a));
    vtkMath::MultiplyScalar(ba, vtkMath::Norm(b) * vtkMath::Norm(b));
    vtkMath::Subtract(ab, ba, abMinusba);
    vtkMath::Cross(abMinusba, aCrossb, abMinusbaCrossaCrossb);

    double expectedCenter[3];
    vtkMath::MultiplyScalar(
      abMinusbaCrossaCrossb, 1.0 / (2.0 * vtkMath::Norm(aCrossb) * vtkMath::Norm(aCrossb)));
    vtkMath::Add(abMinusbaCrossaCrossb, C, expectedCenter);

    for (int i = 0; i < 3; ++i)
    {
      if (!vtkMathUtilities::FuzzyCompare(
            center[i], expectedCenter[i], std::numeric_limits<double>::epsilon() * 1024.0))
      {
        std::cout << " Solve3PointCircle center expected " << expectedCenter[i] << " but got "
                  << center[i];
        std::cout << "eps ratio is: "
                  << (expectedCenter[i] - center[i]) / std::numeric_limits<double>::epsilon()
                  << std::endl;
        ++status;
      }
    }
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Tested by TestMath
int TestInf()
{
  int status = 0;
  std::cout << "Inf..";

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Tested by TestMath
int TestNegInf()
{
  int status = 0;
  std::cout << "NegInf..";

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Tested by TestMath
int TestNan()
{
  int status = 0;
  std::cout << "Nan..";

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with known solutions
int TestMetaMultiplyMatrix()
{
  int status = 0;
  std::cout << "MetaMuliplyMatrix..";

  vtkNew<vtkIntArray> arrayM1, arrayM2, arrayM3;

  arrayM1->SetNumberOfComponents(25);
  arrayM1->SetNumberOfTuples(1);
  arrayM2->SetNumberOfComponents(10);
  arrayM2->SetNumberOfTuples(1);
  arrayM3->SetNumberOfComponents(25);
  arrayM3->SetNumberOfTuples(1);

  auto rangeM1 = vtk::DataArrayTupleRange<25>(arrayM1);
  auto rangeM2 = vtk::DataArrayTupleRange<10>(arrayM2);
  auto rangeM3 = vtk::DataArrayTupleRange<25>(arrayM3);

  auto M1 = rangeM1[0];
  auto M2 = rangeM2[0];
  auto M3 = rangeM3[0];

  // We fill M1 and M2 with the beginning of the sequence of natural integers.
  std::iota(M1.begin(), M1.end(), 0);
  std::iota(M2.begin(), M2.end(), 0);

  //  0  1  2  3  4      0  1     60  70
  //  5  6  7  8  9      2  3     160 195
  //  10 11 12 13 14  x  4  5  =  260 320
  //  15 16 17 18 19     6  7     360 445
  //  20 21 22 23 24     8  9     460 570
  vtkMath::MultiplyMatrix<5, 5, 2>(M1, M2, M3);

  if (M3[0] != 60 || M3[1] != 70 || M3[2] != 160 || M3[3] != 195 || M3[4] != 260 || M3[5] != 320 ||
    M3[6] != 360 || M3[7] != 445 || M3[8] != 460 || M3[9] != 570)
  {
    status = 1;
  }

  //                         0  5  10 15 20
  //    0  1  2  3  4        1  6  11 16 21       60 160 260 360 460
  //    5  6  7  8  9   x    2  7  12 17 22    =  70 195 320 445 570
  //                         3  8  13 18 23
  //                         4  9  14 19 24
  vtkMath::MultiplyMatrix<2, 5, 5, vtkMatrixUtilities::Layout::Transpose,
    vtkMatrixUtilities::Layout::Transpose>(M2, M1, M3);

  if (M3[0] != 60 || M3[5] != 70 || M3[1] != 160 || M3[6] != 195 || M3[2] != 260 || M3[7] != 320 ||
    M3[3] != 360 || M3[8] != 445 || M3[4] != 460 || M3[9] != 570)
  {
    status = 1;
  }

  // 0  4  8  12 16      0  1     240 280
  // 1  5  9  13 17      2  3     260 305
  // 2  6  10 14 18   x  4  5  =  280 330
  // 3  7  11 15 19      6  7     300 355
  //                     8  9
  vtkMath::MultiplyMatrix<4, 5, 2, vtkMatrixUtilities::Layout::Transpose,
    vtkMatrixUtilities::Layout::Identity>(M1, M2, M3);

  if (M3[0] != 240 || M3[1] != 280 || M3[2] != 260 || M3[3] != 305 || M3[4] != 280 ||
    M3[5] != 330 || M3[6] != 300 || M3[7] != 355)
  {
    status = 1;
  }

  // 0  1  2  3    0     14
  // 4  5  6  7  x 1  =  38
  // 8  9  10 11   2     62
  //               3
  vtkMath::MultiplyMatrix<3, 4, 1>(M1, M2, M3);

  if (M3[0] != 14 || M3[1] != 38 || M3[2] != 62)
  {
    status = 1;
  }

  // 0                   0  0  0  0  0
  // 1                   0  1  2  3  4
  // 2  *  0 1 2 3 4  =  0  2  4  6  8
  // 3                   0  3  6  9  12
  // 4                   0  4  8  12 16
  vtkMath::MultiplyMatrix<5, 1, 5, vtkMatrixUtilities::Layout::Transpose,
    vtkMatrixUtilities::Layout::Transpose>(M1, M2, M3);

  if (M3[0] != 0 || M3[1] != 0 || M3[2] != 0 || M3[3] != 0 || M3[4] != 0 || M3[5] != 0 ||
    M3[6] != 1 || M3[7] != 2 || M3[8] != 3 || M3[9] != 4 || M3[10] != 0 || M3[11] != 2 ||
    M3[12] != 4 || M3[13] != 6 || M3[14] != 8 || M3[15] != 0 || M3[16] != 3 || M3[17] != 6 ||
    M3[18] != 9 || M3[19] != 12 || M3[20] != 0 || M3[21] != 4 || M3[22] != 8 || M3[23] != 12 ||
    M3[24] != 16)
  {
    status = 1;
  }

  //  0  1  2  3  4      0             0  1  4  9  16
  //  5  6  7  8  9        1           0  6  14 24 36
  //  10 11 12 13 14  x      2      =  0  11 24 39 56
  //  15 16 17 18 19           3       0  16 34 54 76
  //  20 21 22 23 24             4     0  21 44 69 96
  vtkMath::MultiplyMatrix<5, 5, 5, vtkMatrixUtilities::Layout::Identity,
    vtkMatrixUtilities::Layout::Diag>(M1, M2, M3);

  if (M3[0] != 0 || M3[1] != 1 || M3[2] != 4 || M3[3] != 9 || M3[4] != 16 || M3[5] != 0 ||
    M3[6] != 6 || M3[7] != 14 || M3[8] != 24 || M3[9] != 36 || M3[10] != 0 || M3[11] != 11 ||
    M3[12] != 24 || M3[13] != 39 || M3[14] != 56 || M3[15] != 0 || M3[16] != 16 || M3[17] != 34 ||
    M3[18] != 54 || M3[19] != 76 || M3[20] != 0 || M3[21] != 21 || M3[22] != 44 || M3[23] != 69 ||
    M3[24] != 96)
  {
    status = 1;
  }

  //  0            0  5  10 15 20     0  0  0  0  0
  //    1          1  6  11 16 21     1  6  11 16 21
  //      2     x  2  7  12 17 22  =  4  14 24 34 56
  //        3      3  8  13 18 23     9  24 39 54 76
  //          4    4  9  14 19 24     16 36 56 76 96
  vtkMath::MultiplyMatrix<5, 5, 5, vtkMatrixUtilities::Layout::Diag,
    vtkMatrixUtilities::Layout::Transpose>(M2, M1, M3);

  if (M3[0] != 0 || M3[1] != 0 || M3[2] != 0 || M3[3] != 0 || M3[4] != 0 || M3[5] != 1 ||
    M3[6] != 6 || M3[7] != 11 || M3[8] != 16 || M3[9] != 21 || M3[10] != 4 || M3[11] != 14 ||
    M3[12] != 24 || M3[13] != 34 || M3[14] != 44 || M3[15] != 9 || M3[16] != 24 || M3[17] != 39 ||
    M3[18] != 54 || M3[19] != 69 || M3[20] != 16 || M3[21] != 36 || M3[22] != 56 || M3[23] != 76 ||
    M3[24] != 96)
  {
    status = 1;
  }

  //  0            0             0
  //    1            1             1
  //      2     x      2      =      4
  //        3            3             9
  //          4            4             16
  vtkMath::MultiplyMatrix<5, 5, 5, vtkMatrixUtilities::Layout::Diag,
    vtkMatrixUtilities::Layout::Diag>(M1, M2, M3);

  if (M3[0] != 0 || M3[1] != 1 || M3[2] != 4 || M3[3] != 9 || M3[4] != 16)
  {
    status = 1;
  }

  //  0 0   0 1 2   0 0 0
  //  0 1 x 3 4 5 = 3 4 5
  vtkMath::MultiplyMatrix<2, 2, 3, vtkMatrixUtilities::Layout::Diag>(M1, M2, M3);

  if (M3[0] != 0 || M3[1] != 0 || M3[2] != 0 || M3[3] != 3 || M3[4] != 4 || M3[5] != 5)
  {
    status = 1;
  }

  //  0 0   0 1   0 0
  //  0 1 x 2 3 = 2 3
  //  0 0         0 0
  vtkMath::MultiplyMatrix<3, 2, 2, vtkMatrixUtilities::Layout::Diag>(M1, M2, M3);

  if (M3[0] != 0 || M3[1] != 0 || M3[2] != 2 || M3[3] != 3 || M3[4] != 0 || M3[5] != 0)
  {
    status = 1;
  }

  // 0 1 2   0 0 0 0   0  1  4  0
  // 3 4 5 x 0 1 0 0 = 0  4  10 0
  // 6 7 8   0 0 2 0   0  7  16 0
  vtkMath::MultiplyMatrix<3, 3, 4, vtkMatrixUtilities::Layout::Identity,
    vtkMatrixUtilities::Layout::Diag>(M1, M2, M3);

  if (M3[0] != 0 || M3[1] != 1 || M3[2] != 4 || M3[3] != 0 || M3[4] != 0 || M3[5] != 4 ||
    M3[6] != 10 || M3[7] != 0 || M3[8] != 0 || M3[9] != 7 || M3[10] != 16 || M3[11] != 0)
  {
    status = 1;
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with known solutions
int TestMetaMultiplyMatrixWithVector()
{
  int status = 0;
  std::cout << "MetaMuliplyMatrixWithVector..";

  std::array<int, 9> M = { 0, 1, 2, 3, 4, 5, 6, 7, 8 };
  int x[3] = { 10, 100, 1000 }, y[3];

  // Testing with regular array output
  vtkMath::MultiplyMatrixWithVector<3, 3>(M, x, y);

  if (y[0] != 2100 || y[1] != 5430 || y[2] != 8760)
  {
    status = 1;
  }

  // Testing with the transpose
  vtkMath::MultiplyMatrixWithVector<3, 3, vtkMatrixUtilities::Layout::Transpose>(M, x, y);

  if (y[0] != 6300 || y[1] != 7410 || y[2] != 8520)
  {
    status = 1;
  }

  int* tmp = &(y[0]);
  // Testing with pointer output, only for compiling purpose.
  vtkMath::MultiplyMatrixWithVector<3, 3>(M, x, tmp);

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with known solutions
int TestMetaDot()
{
  int status = 0;
  std::cout << "MetaDot..";

  std::array<int, 5> x = { 0, 1, 2, 3, 4 };
  std::array<int, 5> y = { 0, 1, 2, 3, 4 };
  int M[5] = { 2, 4, 6, 8, 9 };

  if (vtkMath::Dot<int, 5>(x, y) != 30)
  {
    status = 1;
  }

  if (vtkMath::Dot<int, 5, vtkMatrixUtilities::Layout::Diag>(x, M, y) != 244)
  {
    status = 1;
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with known solutions
int TestMetaDeterminant()
{
  int status = 0;
  std::cout << "MetaDeterminant..";

  int M[9] = { 1, 2, 3, 3, 2, 1, 2, 3, 1 };

  // | 1 2 3 |
  // | 3 2 1 | = 12
  // | 2 3 1 |
  if (vtkMath::Determinant<3>(M) != 12)
  {
    status = 1;
  }

  // | 1 0 0 |
  // | 0 2 0 | = 6
  // | 0 0 3 |
  if (vtkMath::Determinant<3, vtkMatrixUtilities::Layout::Diag>(M) != 6)
  {
    status = 1;
  }

  // | 1 2 |
  // | 3 3 | = -3
  if (vtkMath::Determinant<2>(M) != -3)
  {
    status = 1;
  }

  // | 1 0 |
  // | 0 2 | = 2
  if (vtkMath::Determinant<2, vtkMatrixUtilities::Layout::Diag>(M) != 2)
  {
    status = 1;
  }

  // | 1 | = 1
  if (vtkMath::Determinant<1>(M) != 1)
  {
    status = 1;
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with known solutions
int TestMetaInvertMatrix()
{
  int status = 0;
  std::cout << "MetaInvertMatrix..";

  int M[9] = { 1, 1, 0, 0, 0, 1, 1, 0, 1 };
  int invM[9];
  double diag[3] = { 1.0, 2.0, 3.0 };
  double invDiag[3];
  int N[4] = { 1, -1, 0, 1 };
  int invN[4];

  //     1 1 0     0 -1  1
  // inv 0 0 1  =  1  1 -1
  //     1 0 1     0  1  0
  vtkMath::InvertMatrix<3>(M, invM);
  if (invM[0] != 0 || invM[1] != -1 || invM[2] != 1 || invM[3] != 1 || invM[4] != 1 ||
    invM[5] != -1 || invM[6] != 0 || invM[7] != 1 || invM[8] != 0)
  {
    std::cout << "a" << std::endl;
    status = 1;
  }

  //     1 0 1     0  1  0
  // inv 1 0 0  = -1  1  1
  //     0 1 1     1 -1  0
  vtkMath::InvertMatrix<3, vtkMatrixUtilities::Layout::Transpose>(M, invM);
  if (invM[0] != 0 || invM[1] != 1 || invM[2] != 0 || invM[3] != -1 || invM[4] != 1 ||
    invM[5] != 1 || invM[6] != 1 || invM[7] != -1 || invM[8] != 0)
  {
    std::cout << "b" << std::endl;
    status = 1;
  }

  //     1         1
  // inv   2    =    1/2
  //         3           1/3
  vtkMath::InvertMatrix<3, vtkMatrixUtilities::Layout::Diag>(diag, invDiag);
  if (invDiag[0] != 1.0 || invDiag[1] != 1.0 / 2.0 || invDiag[2] != 1.0 / 3.0)
  {
    std::cout << "c" << std::endl;
    status = 1;
  }

  //     1 -1     1 1
  // inv 0  1  =  0 1
  vtkMath::InvertMatrix<2>(N, invN);
  if (invN[0] != 1 || invN[1] != 1 || invN[2] != 0 || invN[3] != 1)
  {
    std::cout << "d" << std::endl;
    status = 1;
  }

  //      1 0     1 0
  // inv -1 1  =  1 1
  vtkMath::InvertMatrix<2, vtkMatrixUtilities::Layout::Transpose>(N, invN);
  if (invN[0] != 1 || invN[1] != 0 || invN[2] != 1 || invN[3] != 1)
  {
    std::cout << "e" << std::endl;
    status = 1;
  }

  //     1      1
  // inv   -1 =   -1
  vtkMath::InvertMatrix<2, vtkMatrixUtilities::Layout::Diag>(N, invN);
  if (invN[0] != 1 || invN[1] != -1)
  {
    std::cout << "f" << std::endl;
    status = 1;
  }

  // inv 1 = 1
  vtkMath::InvertMatrix<1>(N, invN);
  if (invN[0] != 1)
  {
    std::cout << "g" << std::endl;
    status = 1;
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

// Validate with known solutions
int TestMetaLinearSolve()
{
  int status = 0;
  std::cout << "MetaLinearSolve..";

  int M[9] = { 1, 1, 0, 1, 0, 1, 0, 0, 1 };
  int X[3] = { 1, 2, 3 };
  int Y[3];

  // 1 1 0   -1     1
  // 0 1 0 x  2  =  2
  // 0 0 1    3     3
  vtkMath::LinearSolve<3, 3>(M, X, Y);
  if (Y[0] != -1 || Y[1] != 2 || Y[2] != 3)
  {
    status = 1;
  }

  // 1 0 0    2     1
  // 1 1 0 x -1  =  2
  // 0 0 1    4     3
  vtkMath::LinearSolve<3, 3, vtkMatrixUtilities::Layout::Transpose>(M, X, Y);
  if (Y[0] != 2 || Y[1] != -1 || Y[2] != 4)
  {
    status = 1;
  }

  // 1        1     1
  //   2   x  1  =  2
  //     3    1     3
  vtkMath::LinearSolve<3, 3, vtkMatrixUtilities::Layout::Diag>(X, X, Y);
  if (Y[0] != 1 || Y[1] != 1 || Y[2] != 1)
  {
    status = 1;
  }

  int N[4] = { -1, 1, 0, 1 };

  // -1  1   1   1
  //  0  1 x 2 = 2
  vtkMath::LinearSolve<2, 2>(N, X, Y);
  if (Y[0] != 1 || Y[1] != 2)
  {
    status = 1;
  }

  // -1  0   -1   1
  //  1  1 x  3 = 2
  vtkMath::LinearSolve<2, 2, vtkMatrixUtilities::Layout::Transpose>(N, X, Y);
  if (Y[0] != -1 || Y[1] != 3)
  {
    std::cout << Y[0] << ", " << Y[1] << std::endl;
    status = 1;
  }

  // -1      -1   1
  //     1 x  2 = 2
  vtkMath::LinearSolve<2, 2, vtkMatrixUtilities::Layout::Diag>(N, X, Y);
  if (Y[0] != -1 || Y[1] != 2)
  {
    status = 1;
  }

  if (status)
  {
    std::cout << "..FAILED" << std::endl;
  }
  else
  {
    std::cout << ".PASSED" << std::endl;
  }
  return status;
}

static_assert(vtkMatrixUtilities::MatrixIs2DArray<double**>(), "Failed");
static_assert(vtkMatrixUtilities::MatrixIs2DArray<double[][3]>(), "Failed");
static_assert(vtkMatrixUtilities::MatrixIs2DArray<double[3][3]>(), "Failed");
static_assert(vtkMatrixUtilities::MatrixIs2DArray<double*[]>(), "Failed");
static_assert(!vtkMatrixUtilities::MatrixIs2DArray<double*>(), "Failed");
static_assert(!vtkMatrixUtilities::MatrixIs2DArray<double[]>(), "Failed");
static_assert(!vtkMatrixUtilities::MatrixIs2DArray<double>(), "Failed");
static_assert(!vtkMatrixUtilities::MatrixIs2DArray<double***>(), "Failed");