1 /*
2 Include "petsctao.h" so that we can use TAO solvers. Note that this
3 file automatically includes libraries such as:
4 petsc.h - base PETSc routines petscvec.h - vectors
5 petscsys.h - system routines petscmat.h - matrices
6 petscis.h - index sets petscksp.h - Krylov subspace methods
7 petscviewer.h - viewers petscpc.h - preconditioners
8
9 */
10
11 #include <petsctao.h>
12
13 /*
14 Description: These data are the result of a NIST study involving
15 ultrasonic calibration. The response variable is
16 ultrasonic response, and the predictor variable is
17 metal distance.
18
19 Reference: Chwirut, D., NIST (197?).
20 Ultrasonic Reference Block Study.
21 */
22
23 static char help[] = "Finds the nonlinear least-squares solution to the model \n\
24 y = exp[-b1*x]/(b2+b3*x) + e \n";
25
26 #define NOBSERVATIONS 214
27 #define NPARAMETERS 3
28
29 /* User-defined application context */
30 typedef struct {
31 /* Working space */
32 PetscReal t[NOBSERVATIONS]; /* array of independent variables of observation */
33 PetscReal y[NOBSERVATIONS]; /* array of dependent variables */
34 PetscReal j[NOBSERVATIONS][NPARAMETERS]; /* dense jacobian matrix array*/
35 PetscInt idm[NOBSERVATIONS]; /* Matrix indices for jacobian */
36 PetscInt idn[NPARAMETERS];
37 } AppCtx;
38
39 /* User provided Routines */
40 PetscErrorCode InitializeData(AppCtx *user);
41 PetscErrorCode FormStartingPoint(Vec);
42 PetscErrorCode EvaluateFunction(Tao, Vec, Vec, void *);
43 PetscErrorCode EvaluateJacobian(Tao, Vec, Mat, Mat, void *);
44
45 /*--------------------------------------------------------------------*/
main(int argc,char ** argv)46 int main(int argc, char **argv)
47 {
48 Vec x, f; /* solution, function */
49 Mat J; /* Jacobian matrix */
50 Tao tao; /* Tao solver context */
51 PetscInt i; /* iteration information */
52 PetscReal hist[100], resid[100];
53 PetscInt lits[100];
54 AppCtx user; /* user-defined work context */
55
56 PetscFunctionBeginUser;
57 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
58 /* Allocate vectors */
59 PetscCall(VecCreateSeq(MPI_COMM_SELF, NPARAMETERS, &x));
60 PetscCall(VecCreateSeq(MPI_COMM_SELF, NOBSERVATIONS, &f));
61
62 /* Create the Jacobian matrix. */
63 PetscCall(MatCreateSeqDense(MPI_COMM_SELF, NOBSERVATIONS, NPARAMETERS, NULL, &J));
64
65 for (i = 0; i < NOBSERVATIONS; i++) user.idm[i] = i;
66
67 for (i = 0; i < NPARAMETERS; i++) user.idn[i] = i;
68
69 /* Create TAO solver and set desired solution method */
70 PetscCall(TaoCreate(PETSC_COMM_SELF, &tao));
71 PetscCall(TaoSetType(tao, TAOPOUNDERS));
72
73 /* Set the function and Jacobian routines. */
74 PetscCall(InitializeData(&user));
75 PetscCall(FormStartingPoint(x));
76 PetscCall(TaoSetSolution(tao, x));
77 PetscCall(TaoSetResidualRoutine(tao, f, EvaluateFunction, (void *)&user));
78 PetscCall(TaoSetJacobianResidualRoutine(tao, J, J, EvaluateJacobian, (void *)&user));
79
80 /* Check for any TAO command line arguments */
81 PetscCall(TaoSetFromOptions(tao));
82
83 PetscCall(TaoSetConvergenceHistory(tao, hist, resid, 0, lits, 100, PETSC_TRUE));
84 /* Perform the Solve */
85 PetscCall(TaoSolve(tao));
86
87 /* View the vector; then destroy it. */
88 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_SELF));
89
90 /* Free TAO data structures */
91 PetscCall(TaoDestroy(&tao));
92
93 /* Free PETSc data structures */
94 PetscCall(VecDestroy(&x));
95 PetscCall(VecDestroy(&f));
96 PetscCall(MatDestroy(&J));
97
98 PetscCall(PetscFinalize());
99 return 0;
100 }
101
102 /*--------------------------------------------------------------------*/
EvaluateFunction(Tao tao,Vec X,Vec F,void * ptr)103 PetscErrorCode EvaluateFunction(Tao tao, Vec X, Vec F, void *ptr)
104 {
105 AppCtx *user = (AppCtx *)ptr;
106 PetscInt i;
107 const PetscReal *x;
108 PetscReal *y = user->y, *f, *t = user->t;
109
110 PetscFunctionBegin;
111 PetscCall(VecGetArrayRead(X, &x));
112 PetscCall(VecGetArray(F, &f));
113
114 for (i = 0; i < NOBSERVATIONS; i++) f[i] = y[i] - PetscExpScalar(-x[0] * t[i]) / (x[1] + x[2] * t[i]);
115 PetscCall(VecRestoreArrayRead(X, &x));
116 PetscCall(VecRestoreArray(F, &f));
117 PetscCall(PetscLogFlops(6 * NOBSERVATIONS));
118 PetscFunctionReturn(PETSC_SUCCESS);
119 }
120
121 /*------------------------------------------------------------*/
122 /* J[i][j] = df[i]/dt[j] */
EvaluateJacobian(Tao tao,Vec X,Mat J,Mat Jpre,void * ptr)123 PetscErrorCode EvaluateJacobian(Tao tao, Vec X, Mat J, Mat Jpre, void *ptr)
124 {
125 AppCtx *user = (AppCtx *)ptr;
126 PetscInt i;
127 const PetscReal *x;
128 PetscReal *t = user->t;
129 PetscReal base;
130
131 PetscFunctionBegin;
132 PetscCall(VecGetArrayRead(X, &x));
133 for (i = 0; i < NOBSERVATIONS; i++) {
134 base = PetscExpScalar(-x[0] * t[i]) / (x[1] + x[2] * t[i]);
135
136 user->j[i][0] = t[i] * base;
137 user->j[i][1] = base / (x[1] + x[2] * t[i]);
138 user->j[i][2] = base * t[i] / (x[1] + x[2] * t[i]);
139 }
140
141 /* Assemble the matrix */
142 PetscCall(MatSetValues(J, NOBSERVATIONS, user->idm, NPARAMETERS, user->idn, (PetscReal *)user->j, INSERT_VALUES));
143 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
144 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
145
146 PetscCall(VecRestoreArrayRead(X, &x));
147 PetscCall(PetscLogFlops(NOBSERVATIONS * 13));
148 PetscFunctionReturn(PETSC_SUCCESS);
149 }
150
151 /* ------------------------------------------------------------ */
FormStartingPoint(Vec X)152 PetscErrorCode FormStartingPoint(Vec X)
153 {
154 PetscReal *x;
155
156 PetscFunctionBegin;
157 PetscCall(VecGetArray(X, &x));
158 x[0] = 0.15;
159 x[1] = 0.008;
160 x[2] = 0.010;
161 PetscCall(VecRestoreArray(X, &x));
162 PetscFunctionReturn(PETSC_SUCCESS);
163 }
164
165 /* ---------------------------------------------------------------------- */
InitializeData(AppCtx * user)166 PetscErrorCode InitializeData(AppCtx *user)
167 {
168 PetscReal *t = user->t, *y = user->y;
169 PetscInt i = 0;
170
171 PetscFunctionBegin;
172 y[i] = 92.9000;
173 t[i++] = 0.5000;
174 y[i] = 78.7000;
175 t[i++] = 0.6250;
176 y[i] = 64.2000;
177 t[i++] = 0.7500;
178 y[i] = 64.9000;
179 t[i++] = 0.8750;
180 y[i] = 57.1000;
181 t[i++] = 1.0000;
182 y[i] = 43.3000;
183 t[i++] = 1.2500;
184 y[i] = 31.1000;
185 t[i++] = 1.7500;
186 y[i] = 23.6000;
187 t[i++] = 2.2500;
188 y[i] = 31.0500;
189 t[i++] = 1.7500;
190 y[i] = 23.7750;
191 t[i++] = 2.2500;
192 y[i] = 17.7375;
193 t[i++] = 2.7500;
194 y[i] = 13.8000;
195 t[i++] = 3.2500;
196 y[i] = 11.5875;
197 t[i++] = 3.7500;
198 y[i] = 9.4125;
199 t[i++] = 4.2500;
200 y[i] = 7.7250;
201 t[i++] = 4.7500;
202 y[i] = 7.3500;
203 t[i++] = 5.2500;
204 y[i] = 8.0250;
205 t[i++] = 5.7500;
206 y[i] = 90.6000;
207 t[i++] = 0.5000;
208 y[i] = 76.9000;
209 t[i++] = 0.6250;
210 y[i] = 71.6000;
211 t[i++] = 0.7500;
212 y[i] = 63.6000;
213 t[i++] = 0.8750;
214 y[i] = 54.0000;
215 t[i++] = 1.0000;
216 y[i] = 39.2000;
217 t[i++] = 1.2500;
218 y[i] = 29.3000;
219 t[i++] = 1.7500;
220 y[i] = 21.4000;
221 t[i++] = 2.2500;
222 y[i] = 29.1750;
223 t[i++] = 1.7500;
224 y[i] = 22.1250;
225 t[i++] = 2.2500;
226 y[i] = 17.5125;
227 t[i++] = 2.7500;
228 y[i] = 14.2500;
229 t[i++] = 3.2500;
230 y[i] = 9.4500;
231 t[i++] = 3.7500;
232 y[i] = 9.1500;
233 t[i++] = 4.2500;
234 y[i] = 7.9125;
235 t[i++] = 4.7500;
236 y[i] = 8.4750;
237 t[i++] = 5.2500;
238 y[i] = 6.1125;
239 t[i++] = 5.7500;
240 y[i] = 80.0000;
241 t[i++] = 0.5000;
242 y[i] = 79.0000;
243 t[i++] = 0.6250;
244 y[i] = 63.8000;
245 t[i++] = 0.7500;
246 y[i] = 57.2000;
247 t[i++] = 0.8750;
248 y[i] = 53.2000;
249 t[i++] = 1.0000;
250 y[i] = 42.5000;
251 t[i++] = 1.2500;
252 y[i] = 26.8000;
253 t[i++] = 1.7500;
254 y[i] = 20.4000;
255 t[i++] = 2.2500;
256 y[i] = 26.8500;
257 t[i++] = 1.7500;
258 y[i] = 21.0000;
259 t[i++] = 2.2500;
260 y[i] = 16.4625;
261 t[i++] = 2.7500;
262 y[i] = 12.5250;
263 t[i++] = 3.2500;
264 y[i] = 10.5375;
265 t[i++] = 3.7500;
266 y[i] = 8.5875;
267 t[i++] = 4.2500;
268 y[i] = 7.1250;
269 t[i++] = 4.7500;
270 y[i] = 6.1125;
271 t[i++] = 5.2500;
272 y[i] = 5.9625;
273 t[i++] = 5.7500;
274 y[i] = 74.1000;
275 t[i++] = 0.5000;
276 y[i] = 67.3000;
277 t[i++] = 0.6250;
278 y[i] = 60.8000;
279 t[i++] = 0.7500;
280 y[i] = 55.5000;
281 t[i++] = 0.8750;
282 y[i] = 50.3000;
283 t[i++] = 1.0000;
284 y[i] = 41.0000;
285 t[i++] = 1.2500;
286 y[i] = 29.4000;
287 t[i++] = 1.7500;
288 y[i] = 20.4000;
289 t[i++] = 2.2500;
290 y[i] = 29.3625;
291 t[i++] = 1.7500;
292 y[i] = 21.1500;
293 t[i++] = 2.2500;
294 y[i] = 16.7625;
295 t[i++] = 2.7500;
296 y[i] = 13.2000;
297 t[i++] = 3.2500;
298 y[i] = 10.8750;
299 t[i++] = 3.7500;
300 y[i] = 8.1750;
301 t[i++] = 4.2500;
302 y[i] = 7.3500;
303 t[i++] = 4.7500;
304 y[i] = 5.9625;
305 t[i++] = 5.2500;
306 y[i] = 5.6250;
307 t[i++] = 5.7500;
308 y[i] = 81.5000;
309 t[i++] = .5000;
310 y[i] = 62.4000;
311 t[i++] = .7500;
312 y[i] = 32.5000;
313 t[i++] = 1.5000;
314 y[i] = 12.4100;
315 t[i++] = 3.0000;
316 y[i] = 13.1200;
317 t[i++] = 3.0000;
318 y[i] = 15.5600;
319 t[i++] = 3.0000;
320 y[i] = 5.6300;
321 t[i++] = 6.0000;
322 y[i] = 78.0000;
323 t[i++] = .5000;
324 y[i] = 59.9000;
325 t[i++] = .7500;
326 y[i] = 33.2000;
327 t[i++] = 1.5000;
328 y[i] = 13.8400;
329 t[i++] = 3.0000;
330 y[i] = 12.7500;
331 t[i++] = 3.0000;
332 y[i] = 14.6200;
333 t[i++] = 3.0000;
334 y[i] = 3.9400;
335 t[i++] = 6.0000;
336 y[i] = 76.8000;
337 t[i++] = .5000;
338 y[i] = 61.0000;
339 t[i++] = .7500;
340 y[i] = 32.9000;
341 t[i++] = 1.5000;
342 y[i] = 13.8700;
343 t[i++] = 3.0000;
344 y[i] = 11.8100;
345 t[i++] = 3.0000;
346 y[i] = 13.3100;
347 t[i++] = 3.0000;
348 y[i] = 5.4400;
349 t[i++] = 6.0000;
350 y[i] = 78.0000;
351 t[i++] = .5000;
352 y[i] = 63.5000;
353 t[i++] = .7500;
354 y[i] = 33.8000;
355 t[i++] = 1.5000;
356 y[i] = 12.5600;
357 t[i++] = 3.0000;
358 y[i] = 5.6300;
359 t[i++] = 6.0000;
360 y[i] = 12.7500;
361 t[i++] = 3.0000;
362 y[i] = 13.1200;
363 t[i++] = 3.0000;
364 y[i] = 5.4400;
365 t[i++] = 6.0000;
366 y[i] = 76.8000;
367 t[i++] = .5000;
368 y[i] = 60.0000;
369 t[i++] = .7500;
370 y[i] = 47.8000;
371 t[i++] = 1.0000;
372 y[i] = 32.0000;
373 t[i++] = 1.5000;
374 y[i] = 22.2000;
375 t[i++] = 2.0000;
376 y[i] = 22.5700;
377 t[i++] = 2.0000;
378 y[i] = 18.8200;
379 t[i++] = 2.5000;
380 y[i] = 13.9500;
381 t[i++] = 3.0000;
382 y[i] = 11.2500;
383 t[i++] = 4.0000;
384 y[i] = 9.0000;
385 t[i++] = 5.0000;
386 y[i] = 6.6700;
387 t[i++] = 6.0000;
388 y[i] = 75.8000;
389 t[i++] = .5000;
390 y[i] = 62.0000;
391 t[i++] = .7500;
392 y[i] = 48.8000;
393 t[i++] = 1.0000;
394 y[i] = 35.2000;
395 t[i++] = 1.5000;
396 y[i] = 20.0000;
397 t[i++] = 2.0000;
398 y[i] = 20.3200;
399 t[i++] = 2.0000;
400 y[i] = 19.3100;
401 t[i++] = 2.5000;
402 y[i] = 12.7500;
403 t[i++] = 3.0000;
404 y[i] = 10.4200;
405 t[i++] = 4.0000;
406 y[i] = 7.3100;
407 t[i++] = 5.0000;
408 y[i] = 7.4200;
409 t[i++] = 6.0000;
410 y[i] = 70.5000;
411 t[i++] = .5000;
412 y[i] = 59.5000;
413 t[i++] = .7500;
414 y[i] = 48.5000;
415 t[i++] = 1.0000;
416 y[i] = 35.8000;
417 t[i++] = 1.5000;
418 y[i] = 21.0000;
419 t[i++] = 2.0000;
420 y[i] = 21.6700;
421 t[i++] = 2.0000;
422 y[i] = 21.0000;
423 t[i++] = 2.5000;
424 y[i] = 15.6400;
425 t[i++] = 3.0000;
426 y[i] = 8.1700;
427 t[i++] = 4.0000;
428 y[i] = 8.5500;
429 t[i++] = 5.0000;
430 y[i] = 10.1200;
431 t[i++] = 6.0000;
432 y[i] = 78.0000;
433 t[i++] = .5000;
434 y[i] = 66.0000;
435 t[i++] = .6250;
436 y[i] = 62.0000;
437 t[i++] = .7500;
438 y[i] = 58.0000;
439 t[i++] = .8750;
440 y[i] = 47.7000;
441 t[i++] = 1.0000;
442 y[i] = 37.8000;
443 t[i++] = 1.2500;
444 y[i] = 20.2000;
445 t[i++] = 2.2500;
446 y[i] = 21.0700;
447 t[i++] = 2.2500;
448 y[i] = 13.8700;
449 t[i++] = 2.7500;
450 y[i] = 9.6700;
451 t[i++] = 3.2500;
452 y[i] = 7.7600;
453 t[i++] = 3.7500;
454 y[i] = 5.4400;
455 t[i++] = 4.2500;
456 y[i] = 4.8700;
457 t[i++] = 4.7500;
458 y[i] = 4.0100;
459 t[i++] = 5.2500;
460 y[i] = 3.7500;
461 t[i++] = 5.7500;
462 y[i] = 24.1900;
463 t[i++] = 3.0000;
464 y[i] = 25.7600;
465 t[i++] = 3.0000;
466 y[i] = 18.0700;
467 t[i++] = 3.0000;
468 y[i] = 11.8100;
469 t[i++] = 3.0000;
470 y[i] = 12.0700;
471 t[i++] = 3.0000;
472 y[i] = 16.1200;
473 t[i++] = 3.0000;
474 y[i] = 70.8000;
475 t[i++] = .5000;
476 y[i] = 54.7000;
477 t[i++] = .7500;
478 y[i] = 48.0000;
479 t[i++] = 1.0000;
480 y[i] = 39.8000;
481 t[i++] = 1.5000;
482 y[i] = 29.8000;
483 t[i++] = 2.0000;
484 y[i] = 23.7000;
485 t[i++] = 2.5000;
486 y[i] = 29.6200;
487 t[i++] = 2.0000;
488 y[i] = 23.8100;
489 t[i++] = 2.5000;
490 y[i] = 17.7000;
491 t[i++] = 3.0000;
492 y[i] = 11.5500;
493 t[i++] = 4.0000;
494 y[i] = 12.0700;
495 t[i++] = 5.0000;
496 y[i] = 8.7400;
497 t[i++] = 6.0000;
498 y[i] = 80.7000;
499 t[i++] = .5000;
500 y[i] = 61.3000;
501 t[i++] = .7500;
502 y[i] = 47.5000;
503 t[i++] = 1.0000;
504 y[i] = 29.0000;
505 t[i++] = 1.5000;
506 y[i] = 24.0000;
507 t[i++] = 2.0000;
508 y[i] = 17.7000;
509 t[i++] = 2.5000;
510 y[i] = 24.5600;
511 t[i++] = 2.0000;
512 y[i] = 18.6700;
513 t[i++] = 2.5000;
514 y[i] = 16.2400;
515 t[i++] = 3.0000;
516 y[i] = 8.7400;
517 t[i++] = 4.0000;
518 y[i] = 7.8700;
519 t[i++] = 5.0000;
520 y[i] = 8.5100;
521 t[i++] = 6.0000;
522 y[i] = 66.7000;
523 t[i++] = .5000;
524 y[i] = 59.2000;
525 t[i++] = .7500;
526 y[i] = 40.8000;
527 t[i++] = 1.0000;
528 y[i] = 30.7000;
529 t[i++] = 1.5000;
530 y[i] = 25.7000;
531 t[i++] = 2.0000;
532 y[i] = 16.3000;
533 t[i++] = 2.5000;
534 y[i] = 25.9900;
535 t[i++] = 2.0000;
536 y[i] = 16.9500;
537 t[i++] = 2.5000;
538 y[i] = 13.3500;
539 t[i++] = 3.0000;
540 y[i] = 8.6200;
541 t[i++] = 4.0000;
542 y[i] = 7.2000;
543 t[i++] = 5.0000;
544 y[i] = 6.6400;
545 t[i++] = 6.0000;
546 y[i] = 13.6900;
547 t[i++] = 3.0000;
548 y[i] = 81.0000;
549 t[i++] = .5000;
550 y[i] = 64.5000;
551 t[i++] = .7500;
552 y[i] = 35.5000;
553 t[i++] = 1.5000;
554 y[i] = 13.3100;
555 t[i++] = 3.0000;
556 y[i] = 4.8700;
557 t[i++] = 6.0000;
558 y[i] = 12.9400;
559 t[i++] = 3.0000;
560 y[i] = 5.0600;
561 t[i++] = 6.0000;
562 y[i] = 15.1900;
563 t[i++] = 3.0000;
564 y[i] = 14.6200;
565 t[i++] = 3.0000;
566 y[i] = 15.6400;
567 t[i++] = 3.0000;
568 y[i] = 25.5000;
569 t[i++] = 1.7500;
570 y[i] = 25.9500;
571 t[i++] = 1.7500;
572 y[i] = 81.7000;
573 t[i++] = .5000;
574 y[i] = 61.6000;
575 t[i++] = .7500;
576 y[i] = 29.8000;
577 t[i++] = 1.7500;
578 y[i] = 29.8100;
579 t[i++] = 1.7500;
580 y[i] = 17.1700;
581 t[i++] = 2.7500;
582 y[i] = 10.3900;
583 t[i++] = 3.7500;
584 y[i] = 28.4000;
585 t[i++] = 1.7500;
586 y[i] = 28.6900;
587 t[i++] = 1.7500;
588 y[i] = 81.3000;
589 t[i++] = .5000;
590 y[i] = 60.9000;
591 t[i++] = .7500;
592 y[i] = 16.6500;
593 t[i++] = 2.7500;
594 y[i] = 10.0500;
595 t[i++] = 3.7500;
596 y[i] = 28.9000;
597 t[i++] = 1.7500;
598 y[i] = 28.9500;
599 t[i++] = 1.7500;
600 PetscFunctionReturn(PETSC_SUCCESS);
601 }
602
603 /*TEST
604
605 build:
606 requires: !complex !single
607
608 test:
609 args: -tao_monitor_short -tao_max_it 100 -tao_type pounders -tao_gatol 1.e-5
610
611 test:
612 suffix: 2
613 args: -tao_monitor_short -tao_max_it 100 -tao_type brgn -tao_brgn_regularization_type l2prox -tao_brgn_regularizer_weight 1e-4 -tao_gatol 1.e-5 -tao_brgn_subsolver_tao_monitor_short
614
615 test:
616 suffix: 3
617 args: -tao_monitor_short -tao_max_it 100 -tao_type brgn -tao_brgn_regularization_type l1dict -tao_brgn_regularizer_weight 1e-4 -tao_brgn_l1_smooth_epsilon 1e-6 -tao_gatol 1.e-5 -tao_brgn_subsolver_tao_monitor_short
618
619 test:
620 suffix: 4
621 args: -tao_monitor_short -tao_max_it 100 -tao_type brgn -tao_brgn_regularization_type lm -tao_gatol 1.e-5 -tao_brgn_subsolver_tao_type bnls -tao_brgn_subsolver_tao_monitor_short
622
623 TEST*/
624