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 #include <assert.h> // for assert
19 #include <float.h> // for DBL_EPSILON
20 #include <stdio.h> // for fprintf
21 #include <stdlib.h> // for exit, EXIT_...
22 #include "NonSmoothDrivers.h" // for soclcp_driver
23 #include "NumericsFwd.h" // for SolverOptions
24 #include "SOCLCP_Solvers.h" // for soclcp_VI_E...
25 #include "SOCLCP_cst.h" // for SICONOS_SOC...
26 #include "SecondOrderConeLinearComplementarityProblem.h" // for SecondOrder...
27 #include "SolverOptions.h" // for SolverOptions
28 #include "numerics_verbose.h" // for numerics_pr...
29
30 const char* const SICONOS_SOCLCP_NSGS_STR = "SOCLCP_NSGS";
31 const char* const SICONOS_SOCLCP_NSGSV_STR = "SOCLCP_NSGSV";
32 const char* const SICONOS_SOCLCP_TFP_STR = "SOCLCP_TFP";
33 const char* const SICONOS_SOCLCP_NSN_AC_STR = "SOCLCP_NSN_AC";
34 const char* const SICONOS_SOCLCP_NSN_FB_STR = "SOCLCP_NSN_FB";
35 const char* const SICONOS_SOCLCP_DSFP_STR = "SOCLCP_DeSaxceFixedPoint";
36 const char* const SICONOS_SOCLCP_NCPGlockerFBFixedPoint_STR = "SOCLCP_NCPGlockerFBFixedPoint";
37 const char* const SICONOS_SOCLCP_AlartCurnierNewton_STR = "SOCLCP_AlartCurnierNewton";
38 const char* const SICONOS_SOCLCP_DampedAlartCurnierNewton_STR = "SOCLCP_DampedAlartCurnierNewton";
39 const char* const SICONOS_SOCLCP_NCPGlockerFBNewton_STR = "SOCLCP_NCPGlockerFBNewton";
40 const char* const SICONOS_SOCLCP_ProjectionOnConeWithDiagonalization_STR = "SOCLCP_ProjectionOnConeWithDiagonalization";
41 const char* const SICONOS_SOCLCP_ProjectionOnCone_STR = "SOCLCP_ProjectionOnCone";
42 const char* const SICONOS_SOCLCP_ProjectionOnConeWithLocalIteration_STR = "SOCLCP_ProjectionOnConeWithLocalIteration";
43 const char* const SICONOS_SOCLCP_ProjectionOnConeWithRegularization_STR = "SOCLCP_ProjectionOnConeWithRegularization";
44 const char* const SICONOS_SOCLCP_NCPGlockerFBPATH_STR = "SOCLCP_NCPGlockerFBPATH";
45 const char* const SICONOS_SOCLCP_projectionOnCylinder_STR = "SOCLCP_projectionOnCylinder";
46 const char* const SICONOS_SOCLCP_ProjectionOnCone_velocity_STR = "SOCLCP_ProjectionOnCone_velocity";
47 const char* const SICONOS_SOCLCP_PGoC_STR = "SOCLCP_PGoC";
48 const char* const SICONOS_SOCLCP_DeSaxceFixedPoint_STR = "SOCLCP_DeSaxceFixedPoint";
49 const char* const SICONOS_SOCLCP_EG_STR = "SOCLCP_ExtraGradient";
50 const char* const SICONOS_SOCLCP_FPP_STR = "SOCLCP_FixedPointProjection";
51 const char* const SICONOS_SOCLCP_VI_EG_STR = "SOCLCP_VI_ExtraGradient";
52 const char* const SICONOS_SOCLCP_VI_FPP_STR = "SOCLCP_VI_FixedPointProjection";
53 const char* const SICONOS_SOCLCP_HP_STR = "SOCLCP_HyperplaneProjection";
54 const char* const SICONOS_SOCLCP_PROX_STR = "SOCLCP_PROX";
55 const char* const SICONOS_SOCLCP_QUARTIC_STR = "SOCLCP_QUARTIC";
56 const char* const SICONOS_SOCLCP_QUARTIC_NU_STR = "SOCLCP_QUARTIC_NU";
57
58
soclcp_driver(SecondOrderConeLinearComplementarityProblem * problem,double * r,double * v,SolverOptions * options)59 int soclcp_driver(SecondOrderConeLinearComplementarityProblem* problem,
60 double *r, double *v,
61 SolverOptions* options)
62 {
63 if(options == NULL)
64 numerics_error("soclcp_driver", "null input for solver and/or global options");
65
66 assert(options->isSet);
67
68 if(verbose > 0)
69 solver_options_print(options);
70
71 /* Solver name */
72 /*const char* const name = options->solverName;*/
73
74 int info = -1 ;
75
76 /* Check for trivial case */
77 info = soclcp_checkTrivialCase(problem, v, r, options);
78
79
80 if(info == 0)
81 return info;
82
83
84 switch(options->solverId)
85 {
86 /* Non Smooth Gauss Seidel (NSGS) */
87 case SICONOS_SOCLCP_NSGS:
88 {
89 numerics_printf(" ========================== Call NSGS solver for Second Order Cone LCP problem ==========================\n");
90 soclcp_nsgs(problem, r, v, &info, options);
91 break;
92 }
93 /* case SICONOS_SOCLCP_NSGSV: */
94 /* { */
95 /* numerics_printf(numerics_printf(" ========================== Call NSGSV solver for Second Order Cone LCP problem ==========================\n"); */
96 /* soclcp_nsgs_v(problem, r , v , &info , options); */
97 /* break; */
98 /* } */
99 /* /\* Proximal point algorithm *\/ */
100 /* case SICONOS_SOCLCP_PROX: */
101 /* { */
102 /* numerics_printf(numerics_printf(" ========================== Call PROX (Proximal Point) solver for Second Order Cone LCP problem ==========================\n"); */
103 /* soclcp_proximal(problem, r , v , &info , options); */
104 /* break; */
105 /* } */
106 /* /\* Tresca Fixed point algorithm *\/ */
107 /* case SICONOS_SOCLCP_TFP: */
108 /* { */
109 /* numerics_printf(" ========================== Call TFP (Tresca Fixed Point) solver for Second Order Cone LCP problem ==========================\n"); */
110 /* soclcp_TrescaFixedPoint(problem, r , v , &info , options); */
111 /* break; */
112 /* } */
113 case SICONOS_SOCLCP_VI_FPP:
114 {
115 numerics_printf(" ========================== Call VI_FixedPointProjection (VI_FPP) solver for Second Order Cone LCP problem ==========================\n");
116 soclcp_VI_FixedPointProjection(problem, r, v, &info, options);
117 break;
118 }
119 /* VI Extra Gradient algorithm */
120 case SICONOS_SOCLCP_VI_EG:
121 {
122 numerics_printf(" ========================== Call VI_ExtraGradient (VI_EG) solver for Second Order Cone LCP problem ==========================\n");
123 soclcp_VI_ExtraGradient(problem, r, v, &info, options);
124 break;
125 }
126 /* /\* Hyperplane Projection algorithm *\/ */
127 /* case SICONOS_SOCLCP_HP: */
128 /* { */
129 /* numerics_printf(" ========================== Call Hyperplane Projection (HP) solver for Second Order Cone LCP problem ==========================\n"); */
130 /* soclcp_HyperplaneProjection(problem, r , v , &info , options); */
131 /* break; */
132 /* } */
133 /* /\* Alart Curnier in local coordinates *\/ */
134 /* case SICONOS_SOCLCP_NSN_AC: */
135 /* { */
136 /* numerics_printf(" ========================== Call Alart Curnier solver for Second Order Cone LCP problem ==========================\n"); */
137 /* if (problem->M->matrix0) */
138 /* { */
139 /* soclcp_nonsmooth_Newton_AlartCurnier(problem, r , v , &info , options); */
140 /* } */
141 /* else */
142 /* { */
143 /* soclcp_nonsmooth_Newton_AlartCurnier(problem, r , v , &info , options); */
144 /* } */
145 /* break; */
146 /* } */
147 /* /\* Fischer Burmeister in local coordinates *\/ */
148 /* case SICONOS_SOCLCP_NSN_FB: */
149 /* { */
150 /* numerics_printf(" ========================== Call Fischer Burmeister solver for Second Order Cone LCP problem ==========================\n"); */
151 /* soclcp_nonsmooth_Newton_FischerBurmeister(problem, r , v , &info , options); */
152 /* break; */
153 /* } */
154 /* case SICONOS_SOCLCP_QUARTIC_NU: */
155 /* case SICONOS_SOCLCP_QUARTIC: */
156 /* { */
157 /* numerics_printf(" ========================== Call Quartic solver for Second Order Cone LCP problem ==========================\n"); */
158 /* soclcp_unitary_enumerative(problem, r , v , &info , options); */
159 /* break; */
160 /* } */
161 /* case SICONOS_SOCLCP_AlartCurnierNewton: */
162 /* case SICONOS_SOCLCP_DampedAlartCurnierNewton: */
163 /* { */
164 /* numerics_printf(" ========================== Call Quartic solver for Second Order Cone LCP problem ==========================\n"); */
165 /* info =soclcp_Newton_solve(problem, r , options); */
166 /* break; */
167 /* } */
168 default:
169 {
170 fprintf(stderr, "Numerics, SecondOrderConeLinearComplementarity_driver failed. Unknown solver.\n");
171 exit(EXIT_FAILURE);
172
173 }
174 }
175
176 return info;
177
178 }
179
soclcp_checkTrivialCase(SecondOrderConeLinearComplementarityProblem * problem,double * v,double * r,SolverOptions * options)180 int soclcp_checkTrivialCase(SecondOrderConeLinearComplementarityProblem* problem, double* v,
181 double* r, SolverOptions* options)
182 {
183 /* Number of contacts */
184 int nc = problem->nc;
185 double* q = problem->q;
186
187 int i = 0;
188 /* Dimension of the problem */
189 int n = problem->n;
190
191 /* r=0 ?*/
192 int normal=0;
193 for(i = 0; i < nc; i++)
194 {
195 normal = problem->coneIndex[i];
196 if(q[normal] < -DBL_EPSILON)
197 return -1;
198
199 }
200 for(i = 0 ; i < n ; ++i)
201 {
202 v[i] = q[i];
203 r[i] = 0.;
204 }
205 options->iparam[2] = 0;
206 options->iparam[4] = 0;
207 options->dparam[SICONOS_IPARAM_ITER_DONE] = 0.0;
208 options->dparam[3] = 0.0;
209 if(verbose == 1)
210 printf("SecondOrderConeLinearComplementarity driver, trivial solution r = 0, v = q.\n");
211 return 0;
212 }
213