Actual source code: chwirut1.c

  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

  9: */

 11: #include <petsctao.h>

 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.

 19: Reference:     Chwirut, D., NIST (197?).
 20:                Ultrasonic Reference Block Study.
 21: */

 23: static char help[] = "Finds the nonlinear least-squares solution to the model \n\
 24:             y = exp[-b1*x]/(b2+b3*x)  +  e \n";

 26: #define NOBSERVATIONS 214
 27: #define NPARAMETERS   3

 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;

 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 *);

 45: /*--------------------------------------------------------------------*/
 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 */

 57:   PetscInitialize(&argc, &argv, (char *)0, help);
 58:   /* Allocate vectors */
 59:   VecCreateSeq(MPI_COMM_SELF, NPARAMETERS, &x);
 60:   VecCreateSeq(MPI_COMM_SELF, NOBSERVATIONS, &f);

 62:   /* Create the Jacobian matrix. */
 63:   MatCreateSeqDense(MPI_COMM_SELF, NOBSERVATIONS, NPARAMETERS, NULL, &J);

 65:   for (i = 0; i < NOBSERVATIONS; i++) user.idm[i] = i;

 67:   for (i = 0; i < NPARAMETERS; i++) user.idn[i] = i;

 69:   /* Create TAO solver and set desired solution method */
 70:   TaoCreate(PETSC_COMM_SELF, &tao);
 71:   TaoSetType(tao, TAOPOUNDERS);

 73:   /* Set the function and Jacobian routines. */
 74:   InitializeData(&user);
 75:   FormStartingPoint(x);
 76:   TaoSetSolution(tao, x);
 77:   TaoSetResidualRoutine(tao, f, EvaluateFunction, (void *)&user);
 78:   TaoSetJacobianResidualRoutine(tao, J, J, EvaluateJacobian, (void *)&user);

 80:   /* Check for any TAO command line arguments */
 81:   TaoSetFromOptions(tao);

 83:   TaoSetConvergenceHistory(tao, hist, resid, 0, lits, 100, PETSC_TRUE);
 84:   /* Perform the Solve */
 85:   TaoSolve(tao);

 87:   /* View the vector; then destroy it.  */
 88:   VecView(x, PETSC_VIEWER_STDOUT_SELF);

 90:   /* Free TAO data structures */
 91:   TaoDestroy(&tao);

 93:   /* Free PETSc data structures */
 94:   VecDestroy(&x);
 95:   VecDestroy(&f);
 96:   MatDestroy(&J);

 98:   PetscFinalize();
 99:   return 0;
100: }

102: /*--------------------------------------------------------------------*/
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;

110:   VecGetArrayRead(X, &x);
111:   VecGetArray(F, &f);

113:   for (i = 0; i < NOBSERVATIONS; i++) f[i] = y[i] - PetscExpScalar(-x[0] * t[i]) / (x[1] + x[2] * t[i]);
114:   VecRestoreArrayRead(X, &x);
115:   VecRestoreArray(F, &f);
116:   PetscLogFlops(6 * NOBSERVATIONS);
117:   return 0;
118: }

120: /*------------------------------------------------------------*/
121: /* J[i][j] = df[i]/dt[j] */
122: PetscErrorCode EvaluateJacobian(Tao tao, Vec X, Mat J, Mat Jpre, void *ptr)
123: {
124:   AppCtx          *user = (AppCtx *)ptr;
125:   PetscInt         i;
126:   const PetscReal *x;
127:   PetscReal       *t = user->t;
128:   PetscReal        base;

130:   VecGetArrayRead(X, &x);
131:   for (i = 0; i < NOBSERVATIONS; i++) {
132:     base = PetscExpScalar(-x[0] * t[i]) / (x[1] + x[2] * t[i]);

134:     user->j[i][0] = t[i] * base;
135:     user->j[i][1] = base / (x[1] + x[2] * t[i]);
136:     user->j[i][2] = base * t[i] / (x[1] + x[2] * t[i]);
137:   }

139:   /* Assemble the matrix */
140:   MatSetValues(J, NOBSERVATIONS, user->idm, NPARAMETERS, user->idn, (PetscReal *)user->j, INSERT_VALUES);
141:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
142:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);

144:   VecRestoreArrayRead(X, &x);
145:   PetscLogFlops(NOBSERVATIONS * 13);
146:   return 0;
147: }

149: /* ------------------------------------------------------------ */
150: PetscErrorCode FormStartingPoint(Vec X)
151: {
152:   PetscReal *x;

154:   VecGetArray(X, &x);
155:   x[0] = 0.15;
156:   x[1] = 0.008;
157:   x[2] = 0.010;
158:   VecRestoreArray(X, &x);
159:   return 0;
160: }

162: /* ---------------------------------------------------------------------- */
163: PetscErrorCode InitializeData(AppCtx *user)
164: {
165:   PetscReal *t = user->t, *y = user->y;
166:   PetscInt   i = 0;

168:   y[i]   = 92.9000;
169:   t[i++] = 0.5000;
170:   y[i]   = 78.7000;
171:   t[i++] = 0.6250;
172:   y[i]   = 64.2000;
173:   t[i++] = 0.7500;
174:   y[i]   = 64.9000;
175:   t[i++] = 0.8750;
176:   y[i]   = 57.1000;
177:   t[i++] = 1.0000;
178:   y[i]   = 43.3000;
179:   t[i++] = 1.2500;
180:   y[i]   = 31.1000;
181:   t[i++] = 1.7500;
182:   y[i]   = 23.6000;
183:   t[i++] = 2.2500;
184:   y[i]   = 31.0500;
185:   t[i++] = 1.7500;
186:   y[i]   = 23.7750;
187:   t[i++] = 2.2500;
188:   y[i]   = 17.7375;
189:   t[i++] = 2.7500;
190:   y[i]   = 13.8000;
191:   t[i++] = 3.2500;
192:   y[i]   = 11.5875;
193:   t[i++] = 3.7500;
194:   y[i]   = 9.4125;
195:   t[i++] = 4.2500;
196:   y[i]   = 7.7250;
197:   t[i++] = 4.7500;
198:   y[i]   = 7.3500;
199:   t[i++] = 5.2500;
200:   y[i]   = 8.0250;
201:   t[i++] = 5.7500;
202:   y[i]   = 90.6000;
203:   t[i++] = 0.5000;
204:   y[i]   = 76.9000;
205:   t[i++] = 0.6250;
206:   y[i]   = 71.6000;
207:   t[i++] = 0.7500;
208:   y[i]   = 63.6000;
209:   t[i++] = 0.8750;
210:   y[i]   = 54.0000;
211:   t[i++] = 1.0000;
212:   y[i]   = 39.2000;
213:   t[i++] = 1.2500;
214:   y[i]   = 29.3000;
215:   t[i++] = 1.7500;
216:   y[i]   = 21.4000;
217:   t[i++] = 2.2500;
218:   y[i]   = 29.1750;
219:   t[i++] = 1.7500;
220:   y[i]   = 22.1250;
221:   t[i++] = 2.2500;
222:   y[i]   = 17.5125;
223:   t[i++] = 2.7500;
224:   y[i]   = 14.2500;
225:   t[i++] = 3.2500;
226:   y[i]   = 9.4500;
227:   t[i++] = 3.7500;
228:   y[i]   = 9.1500;
229:   t[i++] = 4.2500;
230:   y[i]   = 7.9125;
231:   t[i++] = 4.7500;
232:   y[i]   = 8.4750;
233:   t[i++] = 5.2500;
234:   y[i]   = 6.1125;
235:   t[i++] = 5.7500;
236:   y[i]   = 80.0000;
237:   t[i++] = 0.5000;
238:   y[i]   = 79.0000;
239:   t[i++] = 0.6250;
240:   y[i]   = 63.8000;
241:   t[i++] = 0.7500;
242:   y[i]   = 57.2000;
243:   t[i++] = 0.8750;
244:   y[i]   = 53.2000;
245:   t[i++] = 1.0000;
246:   y[i]   = 42.5000;
247:   t[i++] = 1.2500;
248:   y[i]   = 26.8000;
249:   t[i++] = 1.7500;
250:   y[i]   = 20.4000;
251:   t[i++] = 2.2500;
252:   y[i]   = 26.8500;
253:   t[i++] = 1.7500;
254:   y[i]   = 21.0000;
255:   t[i++] = 2.2500;
256:   y[i]   = 16.4625;
257:   t[i++] = 2.7500;
258:   y[i]   = 12.5250;
259:   t[i++] = 3.2500;
260:   y[i]   = 10.5375;
261:   t[i++] = 3.7500;
262:   y[i]   = 8.5875;
263:   t[i++] = 4.2500;
264:   y[i]   = 7.1250;
265:   t[i++] = 4.7500;
266:   y[i]   = 6.1125;
267:   t[i++] = 5.2500;
268:   y[i]   = 5.9625;
269:   t[i++] = 5.7500;
270:   y[i]   = 74.1000;
271:   t[i++] = 0.5000;
272:   y[i]   = 67.3000;
273:   t[i++] = 0.6250;
274:   y[i]   = 60.8000;
275:   t[i++] = 0.7500;
276:   y[i]   = 55.5000;
277:   t[i++] = 0.8750;
278:   y[i]   = 50.3000;
279:   t[i++] = 1.0000;
280:   y[i]   = 41.0000;
281:   t[i++] = 1.2500;
282:   y[i]   = 29.4000;
283:   t[i++] = 1.7500;
284:   y[i]   = 20.4000;
285:   t[i++] = 2.2500;
286:   y[i]   = 29.3625;
287:   t[i++] = 1.7500;
288:   y[i]   = 21.1500;
289:   t[i++] = 2.2500;
290:   y[i]   = 16.7625;
291:   t[i++] = 2.7500;
292:   y[i]   = 13.2000;
293:   t[i++] = 3.2500;
294:   y[i]   = 10.8750;
295:   t[i++] = 3.7500;
296:   y[i]   = 8.1750;
297:   t[i++] = 4.2500;
298:   y[i]   = 7.3500;
299:   t[i++] = 4.7500;
300:   y[i]   = 5.9625;
301:   t[i++] = 5.2500;
302:   y[i]   = 5.6250;
303:   t[i++] = 5.7500;
304:   y[i]   = 81.5000;
305:   t[i++] = .5000;
306:   y[i]   = 62.4000;
307:   t[i++] = .7500;
308:   y[i]   = 32.5000;
309:   t[i++] = 1.5000;
310:   y[i]   = 12.4100;
311:   t[i++] = 3.0000;
312:   y[i]   = 13.1200;
313:   t[i++] = 3.0000;
314:   y[i]   = 15.5600;
315:   t[i++] = 3.0000;
316:   y[i]   = 5.6300;
317:   t[i++] = 6.0000;
318:   y[i]   = 78.0000;
319:   t[i++] = .5000;
320:   y[i]   = 59.9000;
321:   t[i++] = .7500;
322:   y[i]   = 33.2000;
323:   t[i++] = 1.5000;
324:   y[i]   = 13.8400;
325:   t[i++] = 3.0000;
326:   y[i]   = 12.7500;
327:   t[i++] = 3.0000;
328:   y[i]   = 14.6200;
329:   t[i++] = 3.0000;
330:   y[i]   = 3.9400;
331:   t[i++] = 6.0000;
332:   y[i]   = 76.8000;
333:   t[i++] = .5000;
334:   y[i]   = 61.0000;
335:   t[i++] = .7500;
336:   y[i]   = 32.9000;
337:   t[i++] = 1.5000;
338:   y[i]   = 13.8700;
339:   t[i++] = 3.0000;
340:   y[i]   = 11.8100;
341:   t[i++] = 3.0000;
342:   y[i]   = 13.3100;
343:   t[i++] = 3.0000;
344:   y[i]   = 5.4400;
345:   t[i++] = 6.0000;
346:   y[i]   = 78.0000;
347:   t[i++] = .5000;
348:   y[i]   = 63.5000;
349:   t[i++] = .7500;
350:   y[i]   = 33.8000;
351:   t[i++] = 1.5000;
352:   y[i]   = 12.5600;
353:   t[i++] = 3.0000;
354:   y[i]   = 5.6300;
355:   t[i++] = 6.0000;
356:   y[i]   = 12.7500;
357:   t[i++] = 3.0000;
358:   y[i]   = 13.1200;
359:   t[i++] = 3.0000;
360:   y[i]   = 5.4400;
361:   t[i++] = 6.0000;
362:   y[i]   = 76.8000;
363:   t[i++] = .5000;
364:   y[i]   = 60.0000;
365:   t[i++] = .7500;
366:   y[i]   = 47.8000;
367:   t[i++] = 1.0000;
368:   y[i]   = 32.0000;
369:   t[i++] = 1.5000;
370:   y[i]   = 22.2000;
371:   t[i++] = 2.0000;
372:   y[i]   = 22.5700;
373:   t[i++] = 2.0000;
374:   y[i]   = 18.8200;
375:   t[i++] = 2.5000;
376:   y[i]   = 13.9500;
377:   t[i++] = 3.0000;
378:   y[i]   = 11.2500;
379:   t[i++] = 4.0000;
380:   y[i]   = 9.0000;
381:   t[i++] = 5.0000;
382:   y[i]   = 6.6700;
383:   t[i++] = 6.0000;
384:   y[i]   = 75.8000;
385:   t[i++] = .5000;
386:   y[i]   = 62.0000;
387:   t[i++] = .7500;
388:   y[i]   = 48.8000;
389:   t[i++] = 1.0000;
390:   y[i]   = 35.2000;
391:   t[i++] = 1.5000;
392:   y[i]   = 20.0000;
393:   t[i++] = 2.0000;
394:   y[i]   = 20.3200;
395:   t[i++] = 2.0000;
396:   y[i]   = 19.3100;
397:   t[i++] = 2.5000;
398:   y[i]   = 12.7500;
399:   t[i++] = 3.0000;
400:   y[i]   = 10.4200;
401:   t[i++] = 4.0000;
402:   y[i]   = 7.3100;
403:   t[i++] = 5.0000;
404:   y[i]   = 7.4200;
405:   t[i++] = 6.0000;
406:   y[i]   = 70.5000;
407:   t[i++] = .5000;
408:   y[i]   = 59.5000;
409:   t[i++] = .7500;
410:   y[i]   = 48.5000;
411:   t[i++] = 1.0000;
412:   y[i]   = 35.8000;
413:   t[i++] = 1.5000;
414:   y[i]   = 21.0000;
415:   t[i++] = 2.0000;
416:   y[i]   = 21.6700;
417:   t[i++] = 2.0000;
418:   y[i]   = 21.0000;
419:   t[i++] = 2.5000;
420:   y[i]   = 15.6400;
421:   t[i++] = 3.0000;
422:   y[i]   = 8.1700;
423:   t[i++] = 4.0000;
424:   y[i]   = 8.5500;
425:   t[i++] = 5.0000;
426:   y[i]   = 10.1200;
427:   t[i++] = 6.0000;
428:   y[i]   = 78.0000;
429:   t[i++] = .5000;
430:   y[i]   = 66.0000;
431:   t[i++] = .6250;
432:   y[i]   = 62.0000;
433:   t[i++] = .7500;
434:   y[i]   = 58.0000;
435:   t[i++] = .8750;
436:   y[i]   = 47.7000;
437:   t[i++] = 1.0000;
438:   y[i]   = 37.8000;
439:   t[i++] = 1.2500;
440:   y[i]   = 20.2000;
441:   t[i++] = 2.2500;
442:   y[i]   = 21.0700;
443:   t[i++] = 2.2500;
444:   y[i]   = 13.8700;
445:   t[i++] = 2.7500;
446:   y[i]   = 9.6700;
447:   t[i++] = 3.2500;
448:   y[i]   = 7.7600;
449:   t[i++] = 3.7500;
450:   y[i]   = 5.4400;
451:   t[i++] = 4.2500;
452:   y[i]   = 4.8700;
453:   t[i++] = 4.7500;
454:   y[i]   = 4.0100;
455:   t[i++] = 5.2500;
456:   y[i]   = 3.7500;
457:   t[i++] = 5.7500;
458:   y[i]   = 24.1900;
459:   t[i++] = 3.0000;
460:   y[i]   = 25.7600;
461:   t[i++] = 3.0000;
462:   y[i]   = 18.0700;
463:   t[i++] = 3.0000;
464:   y[i]   = 11.8100;
465:   t[i++] = 3.0000;
466:   y[i]   = 12.0700;
467:   t[i++] = 3.0000;
468:   y[i]   = 16.1200;
469:   t[i++] = 3.0000;
470:   y[i]   = 70.8000;
471:   t[i++] = .5000;
472:   y[i]   = 54.7000;
473:   t[i++] = .7500;
474:   y[i]   = 48.0000;
475:   t[i++] = 1.0000;
476:   y[i]   = 39.8000;
477:   t[i++] = 1.5000;
478:   y[i]   = 29.8000;
479:   t[i++] = 2.0000;
480:   y[i]   = 23.7000;
481:   t[i++] = 2.5000;
482:   y[i]   = 29.6200;
483:   t[i++] = 2.0000;
484:   y[i]   = 23.8100;
485:   t[i++] = 2.5000;
486:   y[i]   = 17.7000;
487:   t[i++] = 3.0000;
488:   y[i]   = 11.5500;
489:   t[i++] = 4.0000;
490:   y[i]   = 12.0700;
491:   t[i++] = 5.0000;
492:   y[i]   = 8.7400;
493:   t[i++] = 6.0000;
494:   y[i]   = 80.7000;
495:   t[i++] = .5000;
496:   y[i]   = 61.3000;
497:   t[i++] = .7500;
498:   y[i]   = 47.5000;
499:   t[i++] = 1.0000;
500:   y[i]   = 29.0000;
501:   t[i++] = 1.5000;
502:   y[i]   = 24.0000;
503:   t[i++] = 2.0000;
504:   y[i]   = 17.7000;
505:   t[i++] = 2.5000;
506:   y[i]   = 24.5600;
507:   t[i++] = 2.0000;
508:   y[i]   = 18.6700;
509:   t[i++] = 2.5000;
510:   y[i]   = 16.2400;
511:   t[i++] = 3.0000;
512:   y[i]   = 8.7400;
513:   t[i++] = 4.0000;
514:   y[i]   = 7.8700;
515:   t[i++] = 5.0000;
516:   y[i]   = 8.5100;
517:   t[i++] = 6.0000;
518:   y[i]   = 66.7000;
519:   t[i++] = .5000;
520:   y[i]   = 59.2000;
521:   t[i++] = .7500;
522:   y[i]   = 40.8000;
523:   t[i++] = 1.0000;
524:   y[i]   = 30.7000;
525:   t[i++] = 1.5000;
526:   y[i]   = 25.7000;
527:   t[i++] = 2.0000;
528:   y[i]   = 16.3000;
529:   t[i++] = 2.5000;
530:   y[i]   = 25.9900;
531:   t[i++] = 2.0000;
532:   y[i]   = 16.9500;
533:   t[i++] = 2.5000;
534:   y[i]   = 13.3500;
535:   t[i++] = 3.0000;
536:   y[i]   = 8.6200;
537:   t[i++] = 4.0000;
538:   y[i]   = 7.2000;
539:   t[i++] = 5.0000;
540:   y[i]   = 6.6400;
541:   t[i++] = 6.0000;
542:   y[i]   = 13.6900;
543:   t[i++] = 3.0000;
544:   y[i]   = 81.0000;
545:   t[i++] = .5000;
546:   y[i]   = 64.5000;
547:   t[i++] = .7500;
548:   y[i]   = 35.5000;
549:   t[i++] = 1.5000;
550:   y[i]   = 13.3100;
551:   t[i++] = 3.0000;
552:   y[i]   = 4.8700;
553:   t[i++] = 6.0000;
554:   y[i]   = 12.9400;
555:   t[i++] = 3.0000;
556:   y[i]   = 5.0600;
557:   t[i++] = 6.0000;
558:   y[i]   = 15.1900;
559:   t[i++] = 3.0000;
560:   y[i]   = 14.6200;
561:   t[i++] = 3.0000;
562:   y[i]   = 15.6400;
563:   t[i++] = 3.0000;
564:   y[i]   = 25.5000;
565:   t[i++] = 1.7500;
566:   y[i]   = 25.9500;
567:   t[i++] = 1.7500;
568:   y[i]   = 81.7000;
569:   t[i++] = .5000;
570:   y[i]   = 61.6000;
571:   t[i++] = .7500;
572:   y[i]   = 29.8000;
573:   t[i++] = 1.7500;
574:   y[i]   = 29.8100;
575:   t[i++] = 1.7500;
576:   y[i]   = 17.1700;
577:   t[i++] = 2.7500;
578:   y[i]   = 10.3900;
579:   t[i++] = 3.7500;
580:   y[i]   = 28.4000;
581:   t[i++] = 1.7500;
582:   y[i]   = 28.6900;
583:   t[i++] = 1.7500;
584:   y[i]   = 81.3000;
585:   t[i++] = .5000;
586:   y[i]   = 60.9000;
587:   t[i++] = .7500;
588:   y[i]   = 16.6500;
589:   t[i++] = 2.7500;
590:   y[i]   = 10.0500;
591:   t[i++] = 3.7500;
592:   y[i]   = 28.9000;
593:   t[i++] = 1.7500;
594:   y[i]   = 28.9500;
595:   t[i++] = 1.7500;
596:   return 0;
597: }

599: /*TEST

601:    build:
602:       requires: !complex !single

604:    test:
605:       args: -tao_smonitor -tao_max_it 100 -tao_type pounders -tao_gatol 1.e-5

607:    test:
608:       suffix: 2
609:       args: -tao_smonitor -tao_max_it 100 -tao_type brgn -tao_brgn_regularization_type l2prox -tao_brgn_regularizer_weight 1e-4 -tao_gatol 1.e-5

611:    test:
612:       suffix: 3
613:       args: -tao_smonitor -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

615:    test:
616:       suffix: 4
617:       args: -tao_smonitor -tao_max_it 100 -tao_type brgn -tao_brgn_regularization_type lm -tao_gatol 1.e-5 -tao_brgn_subsolver_tao_type bnls

619: TEST*/