HybridNonLinearSolver.h revision c981c48f5bc9aefeffc0bcb0cc3934c2fae179dd
1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// -*- coding: utf-8 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// vim: set fileencoding=utf-8 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org> 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_HYBRIDNONLINEARSOLVER_H 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_HYBRIDNONLINEARSOLVER_H 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace HybridNonLinearSolverSpace { 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum Status { 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Running = -1, 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ImproperInputParameters = 0, 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RelativeErrorTooSmall = 1, 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath TooManyFunctionEvaluation = 2, 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath TolTooSmall = 3, 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath NotMakingProgressJacobian = 4, 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath NotMakingProgressIterations = 5, 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath UserAsked = 6 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \ingroup NonLinearOptimization_Module 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Finds a zero of a system of n 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * nonlinear functions in n variables by a modification of the Powell 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * hybrid method ("dogleg"). 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The user must provide a subroutine which calculates the 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * functions. The Jacobian is either provided by the user, or approximated 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * using a forward-difference method. 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename FunctorType, typename Scalar=double> 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass HybridNonLinearSolver 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef DenseIndex Index; 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolver(FunctorType &_functor) 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : functor(_functor) { nfev=njev=iter = 0; fnorm= 0.; useExternalScaling=false;} 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath struct Parameters { 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Parameters() 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : factor(Scalar(100.)) 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath , maxfev(1000) 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath , xtol(internal::sqrt(NumTraits<Scalar>::epsilon())) 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath , nb_of_subdiagonals(-1) 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath , nb_of_superdiagonals(-1) 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath , epsfcn(Scalar(0.)) {} 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar factor; 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index maxfev; // maximum number of function evaluation 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar xtol; 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index nb_of_subdiagonals; 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index nb_of_superdiagonals; 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar epsfcn; 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix< Scalar, Dynamic, 1 > FVectorType; 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix< Scalar, Dynamic, Dynamic > JacobianType; 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* TODO: if eigen provides a triangular storage, use it here */ 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix< Scalar, Dynamic, Dynamic > UpperTriangularType; 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolverSpace::Status hybrj1( 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath FVectorType &x, 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Scalar tol = internal::sqrt(NumTraits<Scalar>::epsilon()) 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ); 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolverSpace::Status solveInit(FVectorType &x); 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolverSpace::Status solveOneStep(FVectorType &x); 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolverSpace::Status solve(FVectorType &x); 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolverSpace::Status hybrd1( 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath FVectorType &x, 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Scalar tol = internal::sqrt(NumTraits<Scalar>::epsilon()) 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ); 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolverSpace::Status solveNumericalDiffInit(FVectorType &x); 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolverSpace::Status solveNumericalDiffOneStep(FVectorType &x); 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolverSpace::Status solveNumericalDiff(FVectorType &x); 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void resetParameters(void) { parameters = Parameters(); } 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Parameters parameters; 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath FVectorType fvec, qtf, diag; 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobianType fjac; 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath UpperTriangularType R; 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index nfev; 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index njev; 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index iter; 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar fnorm; 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool useExternalScaling; 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathprivate: 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath FunctorType &functor; 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index n; 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar sum; 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool sing; 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar temp; 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar delta; 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool jeval; 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index ncsuc; 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar ratio; 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar pnorm, xnorm, fnorm1; 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index nslow1, nslow2; 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index ncfail; 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar actred, prered; 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath FVectorType wa1, wa2, wa3, wa4; 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolver& operator=(const HybridNonLinearSolver&); 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename FunctorType, typename Scalar> 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathHybridNonLinearSolverSpace::Status 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathHybridNonLinearSolver<FunctorType,Scalar>::hybrj1( 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath FVectorType &x, 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Scalar tol 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ) 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath n = x.size(); 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* check the input parameters for errors. */ 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (n <= 0 || tol < 0.) 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::ImproperInputParameters; 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath resetParameters(); 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath parameters.maxfev = 100*(n+1); 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath parameters.xtol = tol; 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath diag.setConstant(n, 1.); 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath useExternalScaling = true; 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return solve(x); 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename FunctorType, typename Scalar> 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathHybridNonLinearSolverSpace::Status 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathHybridNonLinearSolver<FunctorType,Scalar>::solveInit(FVectorType &x) 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath n = x.size(); 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa1.resize(n); wa2.resize(n); wa3.resize(n); wa4.resize(n); 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec.resize(n); 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath qtf.resize(n); 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac.resize(n, n); 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (!useExternalScaling) 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath diag.resize(n); 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'"); 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* Function Body */ 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath nfev = 0; 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath njev = 0; 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* check the input parameters for errors. */ 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0. ) 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::ImproperInputParameters; 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (useExternalScaling) 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j = 0; j < n; ++j) 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (diag[j] <= 0.) 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::ImproperInputParameters; 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* evaluate the function at the starting point */ 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* and calculate its norm. */ 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath nfev = 1; 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ( functor(x, fvec) < 0) 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::UserAsked; 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fnorm = fvec.stableNorm(); 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* initialize iteration counter and monitors. */ 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath iter = 1; 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ncsuc = 0; 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ncfail = 0; 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath nslow1 = 0; 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath nslow2 = 0; 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::Running; 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename FunctorType, typename Scalar> 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathHybridNonLinearSolverSpace::Status 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathHybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(FVectorType &x) 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(x.size()==n); // check the caller is not cheating us 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index j; 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::vector<JacobiRotation<Scalar> > v_givens(n), w_givens(n); 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jeval = true; 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* calculate the jacobian matrix. */ 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ( functor.df(x, fjac) < 0) 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::UserAsked; 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++njev; 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa2 = fjac.colwise().blueNorm(); 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* on the first iteration and if external scaling is not used, scale according */ 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* to the norms of the columns of the initial jacobian. */ 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (iter == 1) { 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (!useExternalScaling) 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (j = 0; j < n; ++j) 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath diag[j] = (wa2[j]==0.) ? 1. : wa2[j]; 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* on the first iteration, calculate the norm of the scaled x */ 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* and initialize the step bound delta. */ 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath xnorm = diag.cwiseProduct(x).stableNorm(); 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath delta = parameters.factor * xnorm; 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (delta == 0.) 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath delta = parameters.factor; 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* compute the qr factorization of the jacobian. */ 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HouseholderQR<JacobianType> qrfac(fjac); // no pivoting: 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* copy the triangular factor of the qr factorization into r. */ 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath R = qrfac.matrixQR(); 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* accumulate the orthogonal factor in fjac. */ 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac = qrfac.householderQ(); 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* form (q transpose)*fvec and store in qtf. */ 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath qtf = fjac.transpose() * fvec; 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* rescale if necessary. */ 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (!useExternalScaling) 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath diag = diag.cwiseMax(wa2); 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while (true) { 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* determine the direction p. */ 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::dogleg<Scalar>(R, diag, qtf, delta, wa1); 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* store the direction p and x + p. calculate the norm of p. */ 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa1 = -wa1; 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa2 = x + wa1; 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pnorm = diag.cwiseProduct(wa1).stableNorm(); 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* on the first iteration, adjust the initial step bound. */ 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (iter == 1) 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath delta = (std::min)(delta,pnorm); 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* evaluate the function at x + p and calculate its norm. */ 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ( functor(wa2, wa4) < 0) 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::UserAsked; 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++nfev; 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fnorm1 = wa4.stableNorm(); 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* compute the scaled actual reduction. */ 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath actred = -1.; 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (fnorm1 < fnorm) /* Computing 2nd power */ 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath actred = 1. - internal::abs2(fnorm1 / fnorm); 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* compute the scaled predicted reduction. */ 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa3 = R.template triangularView<Upper>()*wa1 + qtf; 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp = wa3.stableNorm(); 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath prered = 0.; 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (temp < fnorm) /* Computing 2nd power */ 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath prered = 1. - internal::abs2(temp / fnorm); 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* compute the ratio of the actual to the predicted reduction. */ 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ratio = 0.; 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (prered > 0.) 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ratio = actred / prered; 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* update the step bound. */ 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ratio < Scalar(.1)) { 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ncsuc = 0; 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++ncfail; 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath delta = Scalar(.5) * delta; 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } else { 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ncfail = 0; 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++ncsuc; 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ratio >= Scalar(.5) || ncsuc > 1) 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath delta = (std::max)(delta, pnorm / Scalar(.5)); 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (internal::abs(ratio - 1.) <= Scalar(.1)) { 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath delta = pnorm / Scalar(.5); 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* test for successful iteration. */ 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ratio >= Scalar(1e-4)) { 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* successful iteration. update x, fvec, and their norms. */ 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x = wa2; 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa2 = diag.cwiseProduct(x); 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec = wa4; 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath xnorm = wa2.stableNorm(); 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fnorm = fnorm1; 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++iter; 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* determine the progress of the iteration. */ 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++nslow1; 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (actred >= Scalar(.001)) 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath nslow1 = 0; 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (jeval) 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++nslow2; 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (actred >= Scalar(.1)) 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath nslow2 = 0; 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* test for convergence. */ 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (delta <= parameters.xtol * xnorm || fnorm == 0.) 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::RelativeErrorTooSmall; 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* tests for termination and stringent tolerances. */ 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (nfev >= parameters.maxfev) 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::TooManyFunctionEvaluation; 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (Scalar(.1) * (std::max)(Scalar(.1) * delta, pnorm) <= NumTraits<Scalar>::epsilon() * xnorm) 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::TolTooSmall; 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (nslow2 == 5) 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::NotMakingProgressJacobian; 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (nslow1 == 10) 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::NotMakingProgressIterations; 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* criterion for recalculating jacobian. */ 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ncfail == 2) 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; // leave inner loop and go for the next outer loop iteration 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* calculate the rank one modification to the jacobian */ 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* and update qtf if necessary. */ 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa1 = diag.cwiseProduct( diag.cwiseProduct(wa1)/pnorm ); 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa2 = fjac.transpose() * wa4; 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ratio >= Scalar(1e-4)) 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath qtf = wa2; 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa2 = (wa2-wa3)/pnorm; 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* compute the qr factorization of the updated jacobian. */ 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::r1updt<Scalar>(R, wa1, v_givens, w_givens, wa2, wa3, &sing); 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::r1mpyq<Scalar>(n, n, fjac.data(), v_givens, w_givens); 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::r1mpyq<Scalar>(1, n, qtf.data(), v_givens, w_givens); 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jeval = false; 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::Running; 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename FunctorType, typename Scalar> 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathHybridNonLinearSolverSpace::Status 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathHybridNonLinearSolver<FunctorType,Scalar>::solve(FVectorType &x) 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolverSpace::Status status = solveInit(x); 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (status==HybridNonLinearSolverSpace::ImproperInputParameters) 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return status; 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while (status==HybridNonLinearSolverSpace::Running) 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath status = solveOneStep(x); 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return status; 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename FunctorType, typename Scalar> 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathHybridNonLinearSolverSpace::Status 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathHybridNonLinearSolver<FunctorType,Scalar>::hybrd1( 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath FVectorType &x, 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Scalar tol 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ) 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath n = x.size(); 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* check the input parameters for errors. */ 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (n <= 0 || tol < 0.) 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::ImproperInputParameters; 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath resetParameters(); 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath parameters.maxfev = 200*(n+1); 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath parameters.xtol = tol; 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath diag.setConstant(n, 1.); 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath useExternalScaling = true; 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return solveNumericalDiff(x); 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename FunctorType, typename Scalar> 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathHybridNonLinearSolverSpace::Status 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathHybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(FVectorType &x) 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath n = x.size(); 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (parameters.nb_of_subdiagonals<0) parameters.nb_of_subdiagonals= n-1; 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (parameters.nb_of_superdiagonals<0) parameters.nb_of_superdiagonals= n-1; 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa1.resize(n); wa2.resize(n); wa3.resize(n); wa4.resize(n); 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath qtf.resize(n); 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac.resize(n, n); 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec.resize(n); 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (!useExternalScaling) 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath diag.resize(n); 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'"); 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* Function Body */ 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath nfev = 0; 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath njev = 0; 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* check the input parameters for errors. */ 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.nb_of_subdiagonals< 0 || parameters.nb_of_superdiagonals< 0 || parameters.factor <= 0. ) 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::ImproperInputParameters; 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (useExternalScaling) 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j = 0; j < n; ++j) 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (diag[j] <= 0.) 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::ImproperInputParameters; 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* evaluate the function at the starting point */ 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* and calculate its norm. */ 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath nfev = 1; 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ( functor(x, fvec) < 0) 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::UserAsked; 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fnorm = fvec.stableNorm(); 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* initialize iteration counter and monitors. */ 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath iter = 1; 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ncsuc = 0; 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ncfail = 0; 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath nslow1 = 0; 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath nslow2 = 0; 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::Running; 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename FunctorType, typename Scalar> 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathHybridNonLinearSolverSpace::Status 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathHybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(FVectorType &x) 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(x.size()==n); // check the caller is not cheating us 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index j; 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::vector<JacobiRotation<Scalar> > v_givens(n), w_givens(n); 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jeval = true; 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (parameters.nb_of_subdiagonals<0) parameters.nb_of_subdiagonals= n-1; 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (parameters.nb_of_superdiagonals<0) parameters.nb_of_superdiagonals= n-1; 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* calculate the jacobian matrix. */ 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (internal::fdjac1(functor, x, fvec, fjac, parameters.nb_of_subdiagonals, parameters.nb_of_superdiagonals, parameters.epsfcn) <0) 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::UserAsked; 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath nfev += (std::min)(parameters.nb_of_subdiagonals+parameters.nb_of_superdiagonals+ 1, n); 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa2 = fjac.colwise().blueNorm(); 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* on the first iteration and if external scaling is not used, scale according */ 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* to the norms of the columns of the initial jacobian. */ 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (iter == 1) { 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (!useExternalScaling) 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (j = 0; j < n; ++j) 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath diag[j] = (wa2[j]==0.) ? 1. : wa2[j]; 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* on the first iteration, calculate the norm of the scaled x */ 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* and initialize the step bound delta. */ 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath xnorm = diag.cwiseProduct(x).stableNorm(); 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath delta = parameters.factor * xnorm; 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (delta == 0.) 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath delta = parameters.factor; 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* compute the qr factorization of the jacobian. */ 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HouseholderQR<JacobianType> qrfac(fjac); // no pivoting: 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* copy the triangular factor of the qr factorization into r. */ 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath R = qrfac.matrixQR(); 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* accumulate the orthogonal factor in fjac. */ 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac = qrfac.householderQ(); 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* form (q transpose)*fvec and store in qtf. */ 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath qtf = fjac.transpose() * fvec; 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* rescale if necessary. */ 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (!useExternalScaling) 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath diag = diag.cwiseMax(wa2); 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while (true) { 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* determine the direction p. */ 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::dogleg<Scalar>(R, diag, qtf, delta, wa1); 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* store the direction p and x + p. calculate the norm of p. */ 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa1 = -wa1; 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa2 = x + wa1; 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pnorm = diag.cwiseProduct(wa1).stableNorm(); 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* on the first iteration, adjust the initial step bound. */ 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (iter == 1) 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath delta = (std::min)(delta,pnorm); 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* evaluate the function at x + p and calculate its norm. */ 487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ( functor(wa2, wa4) < 0) 488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::UserAsked; 489c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++nfev; 490c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fnorm1 = wa4.stableNorm(); 491c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* compute the scaled actual reduction. */ 493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath actred = -1.; 494c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (fnorm1 < fnorm) /* Computing 2nd power */ 495c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath actred = 1. - internal::abs2(fnorm1 / fnorm); 496c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 497c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* compute the scaled predicted reduction. */ 498c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa3 = R.template triangularView<Upper>()*wa1 + qtf; 499c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp = wa3.stableNorm(); 500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath prered = 0.; 501c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (temp < fnorm) /* Computing 2nd power */ 502c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath prered = 1. - internal::abs2(temp / fnorm); 503c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 504c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* compute the ratio of the actual to the predicted reduction. */ 505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ratio = 0.; 506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (prered > 0.) 507c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ratio = actred / prered; 508c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 509c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* update the step bound. */ 510c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ratio < Scalar(.1)) { 511c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ncsuc = 0; 512c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++ncfail; 513c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath delta = Scalar(.5) * delta; 514c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } else { 515c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ncfail = 0; 516c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++ncsuc; 517c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ratio >= Scalar(.5) || ncsuc > 1) 518c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath delta = (std::max)(delta, pnorm / Scalar(.5)); 519c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (internal::abs(ratio - 1.) <= Scalar(.1)) { 520c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath delta = pnorm / Scalar(.5); 521c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 522c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 523c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 524c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* test for successful iteration. */ 525c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ratio >= Scalar(1e-4)) { 526c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* successful iteration. update x, fvec, and their norms. */ 527c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x = wa2; 528c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa2 = diag.cwiseProduct(x); 529c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec = wa4; 530c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath xnorm = wa2.stableNorm(); 531c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fnorm = fnorm1; 532c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++iter; 533c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 534c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 535c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* determine the progress of the iteration. */ 536c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++nslow1; 537c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (actred >= Scalar(.001)) 538c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath nslow1 = 0; 539c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (jeval) 540c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++nslow2; 541c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (actred >= Scalar(.1)) 542c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath nslow2 = 0; 543c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 544c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* test for convergence. */ 545c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (delta <= parameters.xtol * xnorm || fnorm == 0.) 546c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::RelativeErrorTooSmall; 547c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 548c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* tests for termination and stringent tolerances. */ 549c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (nfev >= parameters.maxfev) 550c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::TooManyFunctionEvaluation; 551c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (Scalar(.1) * (std::max)(Scalar(.1) * delta, pnorm) <= NumTraits<Scalar>::epsilon() * xnorm) 552c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::TolTooSmall; 553c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (nslow2 == 5) 554c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::NotMakingProgressJacobian; 555c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (nslow1 == 10) 556c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::NotMakingProgressIterations; 557c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 558c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* criterion for recalculating jacobian. */ 559c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ncfail == 2) 560c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; // leave inner loop and go for the next outer loop iteration 561c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 562c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* calculate the rank one modification to the jacobian */ 563c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* and update qtf if necessary. */ 564c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa1 = diag.cwiseProduct( diag.cwiseProduct(wa1)/pnorm ); 565c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa2 = fjac.transpose() * wa4; 566c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ratio >= Scalar(1e-4)) 567c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath qtf = wa2; 568c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath wa2 = (wa2-wa3)/pnorm; 569c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 570c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* compute the qr factorization of the updated jacobian. */ 571c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::r1updt<Scalar>(R, wa1, v_givens, w_givens, wa2, wa3, &sing); 572c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::r1mpyq<Scalar>(n, n, fjac.data(), v_givens, w_givens); 573c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::r1mpyq<Scalar>(1, n, qtf.data(), v_givens, w_givens); 574c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 575c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jeval = false; 576c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 577c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HybridNonLinearSolverSpace::Running; 578c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 579c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 580c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename FunctorType, typename Scalar> 581c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathHybridNonLinearSolverSpace::Status 582c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathHybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiff(FVectorType &x) 583c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 584c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolverSpace::Status status = solveNumericalDiffInit(x); 585c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (status==HybridNonLinearSolverSpace::ImproperInputParameters) 586c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return status; 587c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while (status==HybridNonLinearSolverSpace::Running) 588c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath status = solveNumericalDiffOneStep(x); 589c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return status; 590c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 591c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 592c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 593c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 594c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_HYBRIDNONLINEARSOLVER_H 595c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 596c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath//vim: ai ts=4 sts=4 et sw=4 597