NewtonsMethod.h

Go to the documentation of this file.

//============================================================================
//  The contents of this file are covered by the Viskores license. See
//  LICENSE.txt for details.
//
//  By contributing to this file, all contributors agree to the Developer
//  Certificate of Origin Version 1.1 (DCO 1.1) as stated in DCO.txt.
//============================================================================
 
//============================================================================
//  Copyright (c) Kitware, Inc.
//  All rights reserved.
//  See LICENSE.txt for details.
//
//  This software is distributed WITHOUT ANY WARRANTY; without even
//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
//  PURPOSE.  See the above copyright notice for more information.
//============================================================================
#ifndef viskores_NewtonsMethod_h
#define viskores_NewtonsMethod_h
 
#include <viskores/Math.h>
#include <viskores/Matrix.h>
 
namespace viskores
{
 
template <typename ScalarType, viskores::IdComponent Size>
struct NewtonsMethodResult
{
  bool Valid;
  bool Converged;
  viskores::Vec<ScalarType, Size> Solution;
};
 
VISKORES_SUPPRESS_EXEC_WARNINGS
template <typename ScalarType,
          viskores::IdComponent Size,
          typename JacobianFunctor,
          typename FunctionFunctor>
VISKORES_EXEC_CONT NewtonsMethodResult<ScalarType, Size> NewtonsMethod(
  JacobianFunctor jacobianEvaluator,
  FunctionFunctor functionEvaluator,
  viskores::Vec<ScalarType, Size> desiredFunctionOutput,
  viskores::Vec<ScalarType, Size> initialGuess = viskores::Vec<ScalarType, Size>(ScalarType(0)),
  ScalarType convergeDifference = ScalarType(1e-3),
  viskores::IdComponent maxIterations = 10)
{
  using VectorType = viskores::Vec<ScalarType, Size>;
  using MatrixType = viskores::Matrix<ScalarType, Size, Size>;
 
  VectorType x = initialGuess;
 
  bool valid = false;
  bool converged = false;
  for (viskores::IdComponent iteration = 0; !converged && (iteration < maxIterations); iteration++)
  {
    // For Newton's method, we solve the linear system
    //
    // Jacobian x deltaX = currentFunctionOutput - desiredFunctionOutput
    //
    // The subtraction on the right side simply makes the target of the solve
    // at zero, which is what Newton's method solves for. The deltaX tells us
    // where to move to to solve for a linear system, which we assume will be
    // closer for our nonlinear system.
 
    MatrixType jacobian = jacobianEvaluator(x);
    VectorType currentFunctionOutput = functionEvaluator(x);
 
    VectorType deltaX =
      viskores::SolveLinearSystem(jacobian, currentFunctionOutput - desiredFunctionOutput, valid);
    if (!valid)
    {
      break;
    }
 
    x = x - deltaX;
 
    converged = true;
    for (viskores::IdComponent index = 0; index < Size; index++)
    {
      converged &= (viskores::Abs(deltaX[index]) < convergeDifference);
    }
  }
 
  // Not checking whether converged.
  return { valid, converged, x };
}
 
} // namespace viskores
 
#endif //viskores_NewtonsMethod_h