1 /* Siconos is a program dedicated to modeling, simulation and control 2 * of non smooth dynamical systems. 3 * 4 * Copyright 2021 INRIA. 5 * 6 * Licensed under the Apache License, Version 2.0 (the "License"); 7 * you may not use this file except in compliance with the License. 8 * You may obtain a copy of the License at 9 * 10 * http://www.apache.org/licenses/LICENSE-2.0 11 * 12 * Unless required by applicable law or agreed to in writing, software 13 * distributed under the License is distributed on an "AS IS" BASIS, 14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 15 * See the License for the specific language governing permissions and 16 * limitations under the License. 17 */ 18 19 #ifndef NEWTON_METHODS_H 20 #define NEWTON_METHODS_H 21 22 /*!\file Newton_methods.h 23 * \brief Data structure and function for using Newton based solvers 24 * 25 * The reference used for the implementation is "Finite-Dimensional Variational 26 * Inequalities and Complementarity Problems" by Facchinei and Pang. 27 * 28 * More precisely, the function newton_LSA() is algorithm VFBLSA. 29 * 30 */ 31 32 #include <stddef.h> // for NULL, size_t 33 #include "NumericsFwd.h" // for SolverOptions, NumericsMatrix 34 #include "SiconosConfig.h" // for BUILD_AS_CPP // IWYU pragma: keep 35 #ifndef __cplusplus 36 #include <stdbool.h> // for bool 37 #endif 38 39 typedef void (*compute_F_ptr) (void* data_opaque, double* z, double* F); 40 typedef void (*compute_F_merit_ptr) (void* data_opaque, double* z, double* F, double* F_merit); 41 42 /** \struct functions_LSA Newton_methods.h 43 * Struct holding the necessary pointers to functions needed by the 44 * newton_LSA() procedure. 45 */ 46 typedef struct { 47 compute_F_ptr compute_F; /**< function to evaluate w = F(z) */ 48 compute_F_merit_ptr compute_F_merit; /**< function to evaluate F_merit(z) (e.g. F_FB, F_{min}, ...) */ 49 void (*compute_H)(void* data_opaque, double* z, double* w, double* workV1, double* workV2, NumericsMatrix* H); /**< function to get an element H of T */ 50 void (*compute_error)(void* data_opaque, double* z, double* w, double* nabla_theta, double tol, double* err); /**< function to compute the error */ 51 void (*compute_RHS_desc)(void* data_opaque, double* z, double* w, double* F_desc); /**< function to evaluate F_desc(z) (e.g. F_FB, F_{min}, ...), optional */ 52 void (*compute_H_desc)(void* data_opaque, double* z, double* w, double* workV1, double* workV2, NumericsMatrix* H_desc); /**< function to get an element H_desc of T_desc, optional */ 53 int (*compute_descent_direction)(void* data_opaque, double* z, double* w, double* descent_dir, SolverOptions* options); /**< function to get the descent direction, used for instance in the Newton-Josephy method */ 54 void (*compute_JacTheta_merit)(void* data_opaque, double* z, double* w, double* F_merit, double* workV, double* JacThetaF_merit, SolverOptions* options); /**< function to get the descent direction, used for instance in the Newton-Josephy method */ 55 void* (*get_set_from_problem_data)(void* problem); /**< Function returning the set description from the */ 56 int (*ls_failure_fn)(void* problem, double* z, double* w, double* descent_dir, double err, size_t status); /**< Function to call when the line search fails */ 57 } functions_LSA; 58 59 // id of the stat structure 60 #define NEWTON_STATS_ITERATION 1 61 62 /** \struct newton_stats Newton_methods.h */ 63 typedef struct { 64 int id; /**< id of this structure */ 65 double merit_value; /**< value of the merit function at the end of the iteration */ 66 double alpha; /**< value of the LS parameter */ 67 unsigned int status; /**< status of this newton iteration */ 68 } newton_stats; 69 70 /** \struct newton_LSA_param Newton_methods.h*/ 71 typedef struct { 72 double p; /**< p value for the acceptance test of the direction solution of the linear system */ 73 double sigma; /**< ratio for the decrease in norm of the C-function (\f$gamma'\f$ in VFBLSA)*/ 74 double rho; /**< coefficient for the direction check*/ 75 bool keep_H; /**< keep the matrix H untouched. Only used in the dense case, where a copy of the matrix is factorized */ 76 bool check_dir_quality; /**< Check the quality of the descent direction (Eqn 9.1.6 p. 805 in Facchinei & Pang)*/ 77 } newton_LSA_param; 78 79 /** \struct newton_LSA_data Newton_methods.h*/ 80 typedef struct { 81 NumericsMatrix* H; /**< matrix */ 82 } newton_LSA_data; 83 84 85 enum NEWTON_SOLVER 86 { 87 SICONOS_NEWTON_LSA = 10000 88 }; 89 90 extern const char* const SICONOS_NEWTON_LSA_STR ; 91 92 enum SICONOS_NEWTON_IPARAM 93 { 94 /** line search based algo use this */ 95 SICONOS_IPARAM_LSA_NONMONOTONE_LS = 3, 96 SICONOS_IPARAM_LSA_NONMONOTONE_LS_M = 4, 97 SICONOS_IPARAM_LSA_FORCE_ARCSEARCH = 5, 98 SICONOS_IPARAM_LSA_SEARCH_CRITERION=6, 99 SICONOS_IPARAM_STOPPING_CRITERION=10 100 }; 101 102 enum SICONOS_STOPPING_CRITERION 103 { 104 SICONOS_STOPPING_CRITERION_RESIDU=0, 105 SICONOS_STOPPING_CRITERION_STATIONARITY=1, 106 SICONOS_STOPPING_CRITERION_RESIDU_AND_STATIONARITY=2, 107 SICONOS_STOPPING_CRITERION_USER_ROUTINE=3 108 }; 109 110 111 enum SICONOS_GOLDSTEIN_IPARAM 112 { 113 SICONOS_IPARAM_GOLDSTEIN_ITERMAX=7 114 }; 115 116 enum SICONOS_NMS_IPARAM 117 { 118 /** non-monotone specific part */ 119 SICONOS_IPARAM_NMS_WATCHDOG_TYPE=7, 120 SICONOS_IPARAM_NMS_PROJECTED_GRADIENT_TYPE=8, 121 SICONOS_IPARAM_NMS_N_MAX=9 122 }; 123 124 enum SICONOS_NEWTON_DPARAM{ 125 /** line-search */ 126 SICONOS_DPARAM_LSA_ALPHA_MIN=2, 127 SICONOS_DPARAM_GOLDSTEIN_C=3, 128 SICONOS_DPARAM_GOLDSTEIN_ALPHAMAX=4 129 }; 130 131 enum SICONOS_NMS_DPARAM 132 { 133 /** non-monotone specific part */ 134 SICONOS_DPARAM_NMS_DELTA = 2, 135 SICONOS_DPARAM_NMS_DELTA_VAR = 3, 136 SICONOS_DPARAM_NMS_SIGMA = 4, 137 SICONOS_DPARAM_NMS_ALPHA_MIN_WATCHDOG = 5, 138 SICONOS_DPARAM_NMS_ALPHA_MIN_PGRAD = 6, 139 SICONOS_DPARAM_NMS_MERIT_INCR=7 140 }; 141 142 143 144 145 // status of the newton step 146 #define NEWTON_STATS_NEWTON_STEP 1 147 #define NEWTON_STATS_DESC_DIR 2 148 149 #if defined(__cplusplus) && !defined(BUILD_AS_CPP) 150 extern "C" 151 { 152 #endif 153 154 /** Newton algorithm for finding the zero of a function with a line search. 155 * Mainly used for equation-based reformulation of CP or VI. 156 * \param n size of the problem 157 * \param z variable 158 * \param w value of F(z) 159 * \param info solver-specific values 160 * \param data opaque problem definition 161 * \param options options for this solver 162 * \param functions struct of function pointers to compute F, H and the error 163 */ 164 void newton_LSA(unsigned n, double *z, double *w, int *info, void* data, SolverOptions* options, functions_LSA* functions); 165 166 /** Set the functions to compute F and F_merit and all the other pointers to NULL 167 * \param functions structure to fill 168 * \param compute_F function to compute F 169 * \param merit_function function to compute F_merit 170 */ init_lsa_functions(functions_LSA * functions,compute_F_ptr compute_F,compute_F_merit_ptr merit_function)171 static inline void init_lsa_functions(functions_LSA* functions, compute_F_ptr compute_F, compute_F_merit_ptr merit_function) 172 { 173 functions->compute_F = compute_F; 174 functions->compute_F_merit = merit_function; 175 functions->compute_H = NULL; 176 functions->compute_error = NULL; 177 functions->compute_RHS_desc = NULL; 178 functions->compute_H_desc = NULL; 179 functions->compute_descent_direction = NULL; 180 functions->compute_JacTheta_merit = NULL; 181 functions->get_set_from_problem_data = NULL; 182 functions->ls_failure_fn = NULL; 183 } 184 185 /** Set the parameters and data for newton_LSA 186 * \param options the solver option 187 * \param mat the */ 188 void set_lsa_params_data(SolverOptions* options, NumericsMatrix* mat); 189 190 /** clear the solver-specific data 191 * \param options the SolverOption structure 192 */ 193 void newton_LSA_free_solverOptions(SolverOptions* options); 194 195 /** @addtogroup SetSolverOptions 196 @{ 197 */ 198 void newton_lsa_set_default(SolverOptions* options); 199 /** @} */ 200 201 #if defined(__cplusplus) && !defined(BUILD_AS_CPP) 202 } 203 #endif 204 205 #endif 206