VTK-m  2.0
NewtonsMethod.h
Go to the documentation of this file.
1 //============================================================================
2 // Copyright (c) Kitware, Inc.
4 // See LICENSE.txt for details.
5 //
6 // This software is distributed WITHOUT ANY WARRANTY; without even
7 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
9 //============================================================================
10 #ifndef vtk_m_NewtonsMethod_h
11 #define vtk_m_NewtonsMethod_h
12
13 #include <vtkm/Math.h>
14 #include <vtkm/Matrix.h>
15
16 namespace vtkm
17 {
18
19 template <typename ScalarType, vtkm::IdComponent Size>
21 {
22  bool Valid;
23  bool Converged;
25 };
26
38 template <typename ScalarType,
39  vtkm::IdComponent Size,
40  typename JacobianFunctor,
41  typename FunctionFunctor>
43  JacobianFunctor jacobianEvaluator,
44  FunctionFunctor functionEvaluator,
45  vtkm::Vec<ScalarType, Size> desiredFunctionOutput,
46  vtkm::Vec<ScalarType, Size> initialGuess = vtkm::Vec<ScalarType, Size>(ScalarType(0)),
47  ScalarType convergeDifference = ScalarType(1e-3),
48  vtkm::IdComponent maxIterations = 10)
49 {
50  using VectorType = vtkm::Vec<ScalarType, Size>;
51  using MatrixType = vtkm::Matrix<ScalarType, Size, Size>;
52
53  VectorType x = initialGuess;
54
55  bool valid = false;
56  bool converged = false;
57  for (vtkm::IdComponent iteration = 0; !converged && (iteration < maxIterations); iteration++)
58  {
59  // For Newton's method, we solve the linear system
60  //
61  // Jacobian x deltaX = currentFunctionOutput - desiredFunctionOutput
62  //
63  // The subtraction on the right side simply makes the target of the solve
64  // at zero, which is what Newton's method solves for. The deltaX tells us
65  // where to move to to solve for a linear system, which we assume will be
66  // closer for our nonlinear system.
67
68  MatrixType jacobian = jacobianEvaluator(x);
69  VectorType currentFunctionOutput = functionEvaluator(x);
70
71  VectorType deltaX =
72  vtkm::SolveLinearSystem(jacobian, currentFunctionOutput - desiredFunctionOutput, valid);
73  if (!valid)
74  {
75  break;
76  }
77
78  x = x - deltaX;
79
80  converged = true;
81  for (vtkm::IdComponent index = 0; index < Size; index++)
82  {
83  converged &= (vtkm::Abs(deltaX[index]) < convergeDifference);
84  }
85  }
86
87  // Not checking whether converged.
88  return { valid, converged, x };
89 }
90
91 } // namespace vtkm
92
93 #endif //vtk_m_NewtonsMethod_h
vtkm::NewtonsMethodResult::Valid
bool Valid
Definition: NewtonsMethod.h:22
vtkm
Groups connected points that have the same field value.
Definition: Atomic.h:19
VTKM_EXEC_CONT
#define VTKM_EXEC_CONT
Definition: ExportMacros.h:52
vtkm::IdComponent
vtkm::Int32 IdComponent
Represents a component ID (index of component in a vector).
Definition: Types.h:168
vtkm::NewtonsMethodResult
Definition: NewtonsMethod.h:20
vtkm::SolveLinearSystem
VTKM_EXEC_CONT vtkm::Vec< T, Size > SolveLinearSystem(const vtkm::Matrix< T, Size, Size > &A, const vtkm::Vec< T, Size > &b, bool &valid)
Solve the linear system Ax = b for x.
Definition: Matrix.h:446
Matrix.h
Math.h
vtkm::NewtonsMethodResult::Solution
vtkm::Vec< ScalarType, Size > Solution
Definition: NewtonsMethod.h:24
vtkm::NewtonsMethod
VTKM_SUPPRESS_EXEC_WARNINGS VTKM_EXEC_CONT NewtonsMethodResult< 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)
Uses Newton's method (a.k.a.
Definition: NewtonsMethod.h:42
vtkm::Vec< ScalarType, Size >
vtkm::Matrix
Basic Matrix type.
Definition: Matrix.h:33
vtkm::NewtonsMethodResult::Converged
bool Converged
Definition: NewtonsMethod.h:23
VTKM_SUPPRESS_EXEC_WARNINGS
#define VTKM_SUPPRESS_EXEC_WARNINGS
Definition: ExportMacros.h:53