Function SolveLinearSystem

Common/Core/Testing/Cxx/UnitTestMath.cxx:2552–2619

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;
      }
    }