NewtonsMethod.h 3.59 KB
Newer Older
 Kenneth Moreland committed Aug 27, 2015 1 2 3 4 5 6 7 8 //============================================================================ // 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. // Kenneth Moreland committed Sep 20, 2017 9 // Copyright 2015 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Kenneth Moreland committed Aug 27, 2015 10 11 12 // Copyright 2015 UT-Battelle, LLC. // Copyright 2015 Los Alamos National Security. // Kenneth Moreland committed Sep 20, 2017 13 // Under the terms of Contract DE-NA0003525 with NTESS, Kenneth Moreland committed Aug 27, 2015 14 15 16 17 18 19 // 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. //============================================================================ Kenneth Moreland committed Apr 14, 2016 20 21 #ifndef vtk_m_NewtonsMethod_h #define vtk_m_NewtonsMethod_h Kenneth Moreland committed Aug 27, 2015 22 23 24 25 #include #include Kitware Robot committed May 25, 2017 26 27 namespace vtkm { Kenneth Moreland committed Aug 27, 2015 28 Sujin Philip committed Oct 10, 2017 29 30 31 32 33 34 35 36 template struct NewtonsMethodResult { bool Valid; bool Converged; vtkm::Vec Solution; }; Kenneth Moreland committed Aug 27, 2015 37 38 39 40 41 42 43 44 45 46 /// 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. /// Kenneth Moreland committed Apr 14, 2016 47 VTKM_SUPPRESS_EXEC_WARNINGS Robert Maynard committed May 26, 2017 48 49 50 template Sujin Philip committed Oct 10, 2017 52 VTKM_EXEC_CONT NewtonsMethodResult NewtonsMethod( Robert Maynard committed May 26, 2017 53 54 JacobianFunctor jacobianEvaluator, FunctionFunctor functionEvaluator, Kitware Robot committed May 25, 2017 55 56 vtkm::Vec desiredFunctionOutput, vtkm::Vec initialGuess = vtkm::Vec(ScalarType(0)), Robert Maynard committed May 26, 2017 57 58 ScalarType convergeDifference = ScalarType(1e-3), vtkm::IdComponent maxIterations = 10) Kenneth Moreland committed Aug 27, 2015 59 { Robert Maynard committed Feb 23, 2018 60 61 using VectorType = vtkm::Vec; using MatrixType = vtkm::Matrix; Kenneth Moreland committed Aug 27, 2015 62 63 64 VectorType x = initialGuess; Sujin Philip committed Oct 10, 2017 65 bool valid = false; Kenneth Moreland committed Aug 27, 2015 66 bool converged = false; Kitware Robot committed May 25, 2017 67 68 for (vtkm::IdComponent iteration = 0; !converged && (iteration < maxIterations); iteration++) { Kenneth Moreland committed Aug 27, 2015 69 70 71 72 73 74 75 76 77 78 79 80 81 // 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 = Kitware Robot committed May 25, 2017 82 vtkm::SolveLinearSystem(jacobian, currentFunctionOutput - desiredFunctionOutput, valid); Sujin Philip committed Oct 10, 2017 83 84 85 86 if (!valid) { break; } Kenneth Moreland committed Aug 27, 2015 87 88 89 90 91 x = x - deltaX; converged = true; for (vtkm::IdComponent index = 0; index < Size; index++) Kitware Robot committed May 25, 2017 92 { Kenneth Moreland committed Aug 27, 2015 93 94 converged &= (vtkm::Abs(deltaX[index]) < convergeDifference); } Kitware Robot committed May 25, 2017 95 } Kenneth Moreland committed Aug 27, 2015 96 97 // Not checking whether converged. Sujin Philip committed Oct 10, 2017 98 return { valid, converged, x }; Kenneth Moreland committed Aug 27, 2015 99 100 } Kenneth Moreland committed Apr 14, 2016 101 } // namespace vtkm Kenneth Moreland committed Aug 27, 2015 102 Kenneth Moreland committed Apr 14, 2016 103 #endif //vtk_m_NewtonsMethod_h