fdaccc22db
Change the VTKM_CONT_EXPORT to VTKM_CONT. (Likewise for EXEC and EXEC_CONT.) Remove the inline from these macros so that they can be applied to everything, including implementations in a library. Because inline is not declared in these modifies, you have to add the keyword to functions and methods where the implementation is not inlined in the class.
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
|
|
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
|