1 /* -----------------------------------------------------------------
2 * Programmer(s): Radu Serban @ LLNL
3 * -----------------------------------------------------------------
4 * Example (serial):
5 *
6 * This example solves a nonlinear system from.
7 *
8 * Source: "Handbook of Test Problems in Local and Global Optimization",
9 * C.A. Floudas, P.M. Pardalos et al.
10 * Kluwer Academic Publishers, 1999.
11 * Test problem 4 from Section 14.1, Chapter 14: Ferraris and Tronconi
12 *
13 * This problem involves a blend of trigonometric and exponential terms.
14 * 0.5 sin(x1 x2) - 0.25 x2/pi - 0.5 x1 = 0
15 * (1-0.25/pi) ( exp(2 x1)-e ) + e x2 / pi - 2 e x1 = 0
16 * such that
17 * 0.25 <= x1 <=1.0
18 * 1.5 <= x2 <= 2 pi
19 *
20 * The treatment of the bound constraints on x1 and x2 is done using
21 * the additional variables
22 * l1 = x1 - x1_min >= 0
23 * L1 = x1 - x1_max <= 0
24 * l2 = x2 - x2_min >= 0
25 * L2 = x2 - x2_max >= 0
26 *
27 * and using the constraint feature in KINSOL to impose
28 * l1 >= 0 l2 >= 0
29 * L1 <= 0 L2 <= 0
30 *
31 * The Ferraris-Tronconi test problem has two known solutions.
32 * The nonlinear system is solved by KINSOL using different
33 * combinations of globalization and Jacobian update strategies
34 * and with different initial guesses (leading to one or the other
35 * of the known solutions).
36 *
37 * Constraints are imposed to make all components of the solution
38 * positive.
39 * -----------------------------------------------------------------
40 */
41
42 #include <stdio.h>
43 #include <stdlib.h>
44 #include <math.h>
45
46 #include <kinsol/kinsol.h> /* access to KINSOL func., consts. */
47 #include <nvector/nvector_serial.h> /* access to serial N_Vector */
48 #include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix */
49 #include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver */
50 #include <sundials/sundials_types.h> /* defs. of realtype, sunindextype */
51
52 /* Problem Constants */
53
54 #define NVAR 2
55 #define NEQ 3*NVAR
56
57 #define FTOL RCONST(1.e-5) /* function tolerance */
58 #define STOL RCONST(1.e-5) /* step tolerance */
59
60 #define ZERO RCONST(0.0)
61 #define PT25 RCONST(0.25)
62 #define PT5 RCONST(0.5)
63 #define ONE RCONST(1.0)
64 #define ONEPT5 RCONST(1.5)
65 #define TWO RCONST(2.0)
66
67 #define PI RCONST(3.1415926)
68 #define E RCONST(2.7182818)
69
70 typedef struct {
71 realtype lb[NVAR];
72 realtype ub[NVAR];
73 } *UserData;
74
75 /* Accessor macro */
76 #define Ith(v,i) NV_Ith_S(v,i-1)
77
78 /* Functions Called by the KINSOL Solver */
79 static int func(N_Vector u, N_Vector f, void *user_data);
80
81 /* Private Helper Functions */
82 static void SetInitialGuess1(N_Vector u, UserData data);
83 static void SetInitialGuess2(N_Vector u, UserData data);
84 static int SolveIt(void *kmem, N_Vector u, N_Vector s, int glstr, int mset);
85 static void PrintHeader(realtype fnormtol, realtype scsteptol);
86 static void PrintOutput(N_Vector u);
87 static void PrintFinalStats(void *kmem);
88 static int check_flag(void *flagvalue, const char *funcname, int opt);
89
90 /*
91 *--------------------------------------------------------------------
92 * MAIN PROGRAM
93 *--------------------------------------------------------------------
94 */
95
main()96 int main()
97 {
98 UserData data;
99 realtype fnormtol, scsteptol;
100 N_Vector u1, u2, u, s, c;
101 int glstr, mset, flag;
102 void *kmem;
103 SUNMatrix J;
104 SUNLinearSolver LS;
105
106 u1 = u2 = u = NULL;
107 s = c = NULL;
108 kmem = NULL;
109 J = NULL;
110 LS = NULL;
111 data = NULL;
112
113 /* User data */
114
115 data = (UserData)malloc(sizeof *data);
116 data->lb[0] = PT25; data->ub[0] = ONE;
117 data->lb[1] = ONEPT5; data->ub[1] = TWO*PI;
118
119 /* Create serial vectors of length NEQ */
120 u1 = N_VNew_Serial(NEQ);
121 if (check_flag((void *)u1, "N_VNew_Serial", 0)) return(1);
122
123 u2 = N_VNew_Serial(NEQ);
124 if (check_flag((void *)u2, "N_VNew_Serial", 0)) return(1);
125
126 u = N_VNew_Serial(NEQ);
127 if (check_flag((void *)u, "N_VNew_Serial", 0)) return(1);
128
129 s = N_VNew_Serial(NEQ);
130 if (check_flag((void *)s, "N_VNew_Serial", 0)) return(1);
131
132 c = N_VNew_Serial(NEQ);
133 if (check_flag((void *)c, "N_VNew_Serial", 0)) return(1);
134
135 SetInitialGuess1(u1,data);
136 SetInitialGuess2(u2,data);
137
138 N_VConst(ONE,s); /* no scaling */
139
140 Ith(c,1) = ZERO; /* no constraint on x1 */
141 Ith(c,2) = ZERO; /* no constraint on x2 */
142 Ith(c,3) = ONE; /* l1 = x1 - x1_min >= 0 */
143 Ith(c,4) = -ONE; /* L1 = x1 - x1_max <= 0 */
144 Ith(c,5) = ONE; /* l2 = x2 - x2_min >= 0 */
145 Ith(c,6) = -ONE; /* L2 = x2 - x22_min <= 0 */
146
147 fnormtol=FTOL; scsteptol=STOL;
148
149
150 kmem = KINCreate();
151 if (check_flag((void *)kmem, "KINCreate", 0)) return(1);
152
153 flag = KINSetUserData(kmem, data);
154 if (check_flag(&flag, "KINSetUserData", 1)) return(1);
155 flag = KINSetConstraints(kmem, c);
156 if (check_flag(&flag, "KINSetConstraints", 1)) return(1);
157 flag = KINSetFuncNormTol(kmem, fnormtol);
158 if (check_flag(&flag, "KINSetFuncNormTol", 1)) return(1);
159 flag = KINSetScaledStepTol(kmem, scsteptol);
160 if (check_flag(&flag, "KINSetScaledStepTol", 1)) return(1);
161
162 flag = KINInit(kmem, func, u);
163 if (check_flag(&flag, "KINInit", 1)) return(1);
164
165 /* Create dense SUNMatrix */
166 J = SUNDenseMatrix(NEQ, NEQ);
167 if(check_flag((void *)J, "SUNDenseMatrix", 0)) return(1);
168
169 /* Create dense SUNLinearSolver object */
170 LS = SUNLinSol_Dense(u, J);
171 if(check_flag((void *)LS, "SUNLinSol_Dense", 0)) return(1);
172
173 /* Attach the matrix and linear solver to KINSOL */
174 flag = KINSetLinearSolver(kmem, LS, J);
175 if(check_flag(&flag, "KINSetLinearSolver", 1)) return(1);
176
177 /* Print out the problem size, solution parameters, initial guess. */
178 PrintHeader(fnormtol, scsteptol);
179
180 /* --------------------------- */
181
182 printf("\n------------------------------------------\n");
183 printf("\nInitial guess on lower bounds\n");
184 printf(" [x1,x2] = ");
185 PrintOutput(u1);
186
187 N_VScale(ONE,u1,u);
188 glstr = KIN_NONE;
189 mset = 1;
190 SolveIt(kmem, u, s, glstr, mset);
191
192 /* --------------------------- */
193
194 N_VScale(ONE,u1,u);
195 glstr = KIN_LINESEARCH;
196 mset = 1;
197 SolveIt(kmem, u, s, glstr, mset);
198
199 /* --------------------------- */
200
201 N_VScale(ONE,u1,u);
202 glstr = KIN_NONE;
203 mset = 0;
204 SolveIt(kmem, u, s, glstr, mset);
205
206 /* --------------------------- */
207
208 N_VScale(ONE,u1,u);
209 glstr = KIN_LINESEARCH;
210 mset = 0;
211 SolveIt(kmem, u, s, glstr, mset);
212
213
214
215 /* --------------------------- */
216
217 printf("\n------------------------------------------\n");
218 printf("\nInitial guess in middle of feasible region\n");
219 printf(" [x1,x2] = ");
220 PrintOutput(u2);
221
222 N_VScale(ONE,u2,u);
223 glstr = KIN_NONE;
224 mset = 1;
225 SolveIt(kmem, u, s, glstr, mset);
226
227 /* --------------------------- */
228
229 N_VScale(ONE,u2,u);
230 glstr = KIN_LINESEARCH;
231 mset = 1;
232 SolveIt(kmem, u, s, glstr, mset);
233
234 /* --------------------------- */
235
236 N_VScale(ONE,u2,u);
237 glstr = KIN_NONE;
238 mset = 0;
239 SolveIt(kmem, u, s, glstr, mset);
240
241 /* --------------------------- */
242
243 N_VScale(ONE,u2,u);
244 glstr = KIN_LINESEARCH;
245 mset = 0;
246 SolveIt(kmem, u, s, glstr, mset);
247
248
249
250
251 /* Free memory */
252
253 N_VDestroy(u1);
254 N_VDestroy(u2);
255 N_VDestroy(u);
256 N_VDestroy(s);
257 N_VDestroy(c);
258 KINFree(&kmem);
259 SUNLinSolFree(LS);
260 SUNMatDestroy(J);
261 free(data);
262
263 return(0);
264 }
265
266
SolveIt(void * kmem,N_Vector u,N_Vector s,int glstr,int mset)267 static int SolveIt(void *kmem, N_Vector u, N_Vector s, int glstr, int mset)
268 {
269 int flag;
270
271 printf("\n");
272
273 if (mset==1)
274 printf("Exact Newton");
275 else
276 printf("Modified Newton");
277
278 if (glstr == KIN_NONE)
279 printf("\n");
280 else
281 printf(" with line search\n");
282
283 flag = KINSetMaxSetupCalls(kmem, mset);
284 if (check_flag(&flag, "KINSetMaxSetupCalls", 1)) return(1);
285
286 flag = KINSol(kmem, u, glstr, s, s);
287 if (check_flag(&flag, "KINSol", 1)) return(1);
288
289 printf("Solution:\n [x1,x2] = ");
290 PrintOutput(u);
291
292 PrintFinalStats(kmem);
293
294 return(0);
295
296 }
297
298
299 /*
300 *--------------------------------------------------------------------
301 * FUNCTIONS CALLED BY KINSOL
302 *--------------------------------------------------------------------
303 */
304
305 /*
306 * System function for predator-prey system
307 */
308
func(N_Vector u,N_Vector f,void * user_data)309 static int func(N_Vector u, N_Vector f, void *user_data)
310 {
311 realtype *udata, *fdata;
312 realtype x1, l1, L1, x2, l2, L2;
313 realtype *lb, *ub;
314 UserData data;
315
316 data = (UserData)user_data;
317 lb = data->lb;
318 ub = data->ub;
319
320 udata = N_VGetArrayPointer(u);
321 fdata = N_VGetArrayPointer(f);
322
323 x1 = udata[0];
324 x2 = udata[1];
325 l1 = udata[2];
326 L1 = udata[3];
327 l2 = udata[4];
328 L2 = udata[5];
329
330 fdata[0] = PT5 * sin(x1*x2) - PT25 * x2 / PI - PT5 * x1;
331 fdata[1] = (ONE - PT25/PI)*(exp(TWO*x1)-E) + E*x2/PI - TWO*E*x1;
332 fdata[2] = l1 - x1 + lb[0];
333 fdata[3] = L1 - x1 + ub[0];
334 fdata[4] = l2 - x2 + lb[1];
335 fdata[5] = L2 - x2 + ub[1];
336
337 return(0);
338 }
339
340 /*
341 *--------------------------------------------------------------------
342 * PRIVATE FUNCTIONS
343 *--------------------------------------------------------------------
344 */
345
346 /*
347 * Initial guesses
348 */
349
SetInitialGuess1(N_Vector u,UserData data)350 static void SetInitialGuess1(N_Vector u, UserData data)
351 {
352 realtype x1, x2;
353 realtype *udata;
354 realtype *lb, *ub;
355
356 udata = N_VGetArrayPointer(u);
357
358 lb = data->lb;
359 ub = data->ub;
360
361 /* There are two known solutions for this problem */
362
363 /* this init. guess should take us to (0.29945; 2.83693) */
364 x1 = lb[0];
365 x2 = lb[1];
366
367 udata[0] = x1;
368 udata[1] = x2;
369 udata[2] = x1 - lb[0];
370 udata[3] = x1 - ub[0];
371 udata[4] = x2 - lb[1];
372 udata[5] = x2 - ub[1];
373 }
374
SetInitialGuess2(N_Vector u,UserData data)375 static void SetInitialGuess2(N_Vector u, UserData data)
376 {
377 realtype x1, x2;
378 realtype *udata;
379 realtype *lb, *ub;
380
381 udata = N_VGetArrayPointer(u);
382
383 lb = data->lb;
384 ub = data->ub;
385
386 /* There are two known solutions for this problem */
387
388 /* this init. guess should take us to (0.5; 3.1415926) */
389 x1 = PT5 * (lb[0] + ub[0]);
390 x2 = PT5 * (lb[1] + ub[1]);
391
392 udata[0] = x1;
393 udata[1] = x2;
394 udata[2] = x1 - lb[0];
395 udata[3] = x1 - ub[0];
396 udata[4] = x2 - lb[1];
397 udata[5] = x2 - ub[1];
398 }
399
400 /*
401 * Print first lines of output (problem description)
402 */
403
PrintHeader(realtype fnormtol,realtype scsteptol)404 static void PrintHeader(realtype fnormtol, realtype scsteptol)
405 {
406 printf("\nFerraris and Tronconi test problem\n");
407 printf("Tolerance parameters:\n");
408 #if defined(SUNDIALS_EXTENDED_PRECISION)
409 printf(" fnormtol = %10.6Lg\n scsteptol = %10.6Lg\n",
410 fnormtol, scsteptol);
411 #elif defined(SUNDIALS_DOUBLE_PRECISION)
412 printf(" fnormtol = %10.6g\n scsteptol = %10.6g\n",
413 fnormtol, scsteptol);
414 #else
415 printf(" fnormtol = %10.6g\n scsteptol = %10.6g\n",
416 fnormtol, scsteptol);
417 #endif
418 }
419
420 /*
421 * Print solution
422 */
423
PrintOutput(N_Vector u)424 static void PrintOutput(N_Vector u)
425 {
426 #if defined(SUNDIALS_EXTENDED_PRECISION)
427 printf(" %8.6Lg %8.6Lg\n", Ith(u,1), Ith(u,2));
428 #elif defined(SUNDIALS_DOUBLE_PRECISION)
429 printf(" %8.6g %8.6g\n", Ith(u,1), Ith(u,2));
430 #else
431 printf(" %8.6g %8.6g\n", Ith(u,1), Ith(u,2));
432 #endif
433 }
434
435 /*
436 * Print final statistics contained in iopt
437 */
438
PrintFinalStats(void * kmem)439 static void PrintFinalStats(void *kmem)
440 {
441 long int nni, nfe, nje, nfeD;
442 int flag;
443
444 flag = KINGetNumNonlinSolvIters(kmem, &nni);
445 check_flag(&flag, "KINGetNumNonlinSolvIters", 1);
446 flag = KINGetNumFuncEvals(kmem, &nfe);
447 check_flag(&flag, "KINGetNumFuncEvals", 1);
448
449 flag = KINGetNumJacEvals(kmem, &nje);
450 check_flag(&flag, "KINGetNumJacEvals", 1);
451 flag = KINGetNumLinFuncEvals(kmem, &nfeD);
452 check_flag(&flag, "KINGetNumLinFuncEvals", 1);
453
454 printf("Final Statistics:\n");
455 printf(" nni = %5ld nfe = %5ld \n", nni, nfe);
456 printf(" nje = %5ld nfeD = %5ld \n", nje, nfeD);
457 }
458
459 /*
460 * Check function return value...
461 * opt == 0 means SUNDIALS function allocates memory so check if
462 * returned NULL pointer
463 * opt == 1 means SUNDIALS function returns a flag so check if
464 * flag >= 0
465 * opt == 2 means function allocates memory so check if returned
466 * NULL pointer
467 */
468
check_flag(void * flagvalue,const char * funcname,int opt)469 static int check_flag(void *flagvalue, const char *funcname, int opt)
470 {
471 int *errflag;
472
473 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
474 if (opt == 0 && flagvalue == NULL) {
475 fprintf(stderr,
476 "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
477 funcname);
478 return(1);
479 }
480
481 /* Check if flag < 0 */
482 else if (opt == 1) {
483 errflag = (int *) flagvalue;
484 if (*errflag < 0) {
485 fprintf(stderr,
486 "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n",
487 funcname, *errflag);
488 return(1);
489 }
490 }
491
492 /* Check if function returned NULL pointer - no memory allocated */
493 else if (opt == 2 && flagvalue == NULL) {
494 fprintf(stderr,
495 "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
496 funcname);
497 return(1);
498 }
499
500 return(0);
501 }
502