mirror of
https://gitlab.kitware.com/vtk/vtk-m
synced 2024-09-16 17:22:55 +00:00
c503e0ce8a
All the other math functions are in the vtkm package. This one was in vtkm::exec because it uses a callback method. This can be problematic on CUDA the the declaration of NewtonsMethod does not match the callback method. However, we now have a VTKM_SUPPRESS_EXEC_WARNINGS macro that allows a VTKM_EXEC_CONT_EXPORT function (like NewtonsMethod) to call either a VTKM_EXEC_EXPORT or VTKM_CONT_EXPORT without a warning.
98 lines
3.5 KiB
C++
98 lines
3.5 KiB
C++
//============================================================================
|
|
// 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.
|
|
//
|
|
// Copyright 2015 Sandia Corporation.
|
|
// Copyright 2015 UT-Battelle, LLC.
|
|
// Copyright 2015 Los Alamos National Security.
|
|
//
|
|
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
|
|
// the U.S. Government retains certain rights in this software.
|
|
//
|
|
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
|
|
// Laboratory (LANL), the U.S. Government retains certain rights in
|
|
// this software.
|
|
//============================================================================
|
|
#ifndef vtk_m_NewtonsMethod_h
|
|
#define vtk_m_NewtonsMethod_h
|
|
|
|
#include <vtkm/Math.h>
|
|
#include <vtkm/Matrix.h>
|
|
|
|
namespace vtkm {
|
|
|
|
/// Uses Newton's method (a.k.a. Newton-Raphson method) to solve a nonlinear
|
|
/// system of equations. This function assumes that the number of variables
|
|
/// equals the number of equations. Newton's method operates on an iterative
|
|
/// evaluate and search. Evaluations are performed using the functors passed
|
|
/// into the NewtonsMethod. The first functor returns the NxN matrix of the
|
|
/// Jacobian at a given input point. The second functor returns the N tuple
|
|
/// that is the function evaluation at the given input point. The input point
|
|
/// that evaluates to the desired output, or the closest point found, is
|
|
/// returned.
|
|
///
|
|
VTKM_SUPPRESS_EXEC_WARNINGS
|
|
template<typename ScalarType,
|
|
vtkm::IdComponent Size,
|
|
typename JacobianFunctor,
|
|
typename FunctionFunctor>
|
|
VTKM_EXEC_CONT_EXPORT
|
|
vtkm::Vec<ScalarType,Size>
|
|
NewtonsMethod(JacobianFunctor jacobianEvaluator,
|
|
FunctionFunctor functionEvaluator,
|
|
vtkm::Vec<ScalarType,Size> desiredFunctionOutput,
|
|
vtkm::Vec<ScalarType,Size> initialGuess
|
|
= vtkm::Vec<ScalarType,Size>(ScalarType(0)),
|
|
ScalarType convergeDifference = ScalarType(1e-3),
|
|
vtkm::IdComponent maxIterations = 10)
|
|
{
|
|
typedef vtkm::Vec<ScalarType,Size> VectorType;
|
|
typedef vtkm::Matrix<ScalarType,Size,Size> MatrixType;
|
|
|
|
VectorType x = initialGuess;
|
|
|
|
bool converged = false;
|
|
for (vtkm::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);
|
|
|
|
bool valid; // Ignored.
|
|
VectorType deltaX =
|
|
vtkm::SolveLinearSystem(
|
|
jacobian,
|
|
currentFunctionOutput - desiredFunctionOutput,
|
|
valid);
|
|
|
|
x = x - deltaX;
|
|
|
|
converged = true;
|
|
for (vtkm::IdComponent index = 0; index < Size; index++)
|
|
{
|
|
converged &= (vtkm::Abs(deltaX[index]) < convergeDifference);
|
|
}
|
|
}
|
|
|
|
// Not checking whether converged.
|
|
return x;
|
|
}
|
|
|
|
} // namespace vtkm
|
|
|
|
#endif //vtk_m_NewtonsMethod_h
|