1 #include <petsc/private/taoimpl.h> /*I "petsctao.h" I*/
2
3
4 /*@C
5 TaoSetHessianRoutine - Sets the function to compute the Hessian as well as the location to store the matrix.
6
7 Logically collective on Tao
8
9 Input Parameters:
10 + tao - the Tao context
11 . H - Matrix used for the hessian
12 . Hpre - Matrix that will be used operated on by preconditioner, can be same as H
13 . func - Hessian evaluation routine
14 - ctx - [optional] user-defined context for private data for the
15 Hessian evaluation routine (may be NULL)
16
17 Calling sequence of func:
18 $ func(Tao tao,Vec x,Mat H,Mat Hpre,void *ctx);
19
20 + tao - the Tao context
21 . x - input vector
22 . H - Hessian matrix
23 . Hpre - preconditioner matrix, usually the same as H
24 - ctx - [optional] user-defined Hessian context
25
26 Level: beginner
27 @*/
TaoSetHessianRoutine(Tao tao,Mat H,Mat Hpre,PetscErrorCode (* func)(Tao,Vec,Mat,Mat,void *),void * ctx)28 PetscErrorCode TaoSetHessianRoutine(Tao tao, Mat H, Mat Hpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void*), void *ctx)
29 {
30 PetscErrorCode ierr;
31
32 PetscFunctionBegin;
33 PetscValidHeaderSpecific(tao,TAO_CLASSID,1);
34 if (H) {
35 PetscValidHeaderSpecific(H,MAT_CLASSID,2);
36 PetscCheckSameComm(tao,1,H,2);
37 }
38 if (Hpre) {
39 PetscValidHeaderSpecific(Hpre,MAT_CLASSID,3);
40 PetscCheckSameComm(tao,1,Hpre,3);
41 }
42 if (ctx) {
43 tao->user_hessP = ctx;
44 }
45 if (func) {
46 tao->ops->computehessian = func;
47 }
48 if (H) {
49 ierr = PetscObjectReference((PetscObject)H);CHKERRQ(ierr);
50 ierr = MatDestroy(&tao->hessian);CHKERRQ(ierr);
51 tao->hessian = H;
52 }
53 if (Hpre) {
54 ierr = PetscObjectReference((PetscObject)Hpre);CHKERRQ(ierr);
55 ierr = MatDestroy(&tao->hessian_pre);CHKERRQ(ierr);
56 tao->hessian_pre = Hpre;
57 }
58 PetscFunctionReturn(0);
59 }
60
TaoTestHessian(Tao tao)61 PetscErrorCode TaoTestHessian(Tao tao)
62 {
63 Mat A,B,C,D,hessian;
64 Vec x = tao->solution;
65 PetscErrorCode ierr;
66 PetscReal nrm,gnorm;
67 PetscReal threshold = 1.e-5;
68 PetscInt m,n,M,N;
69 PetscBool complete_print = PETSC_FALSE,test = PETSC_FALSE,flg;
70 PetscViewer viewer,mviewer;
71 MPI_Comm comm;
72 PetscInt tabs;
73 static PetscBool directionsprinted = PETSC_FALSE;
74 PetscViewerFormat format;
75
76 PetscFunctionBegin;
77 ierr = PetscObjectOptionsBegin((PetscObject)tao);CHKERRQ(ierr);
78 ierr = PetscOptionsName("-tao_test_hessian","Compare hand-coded and finite difference Hessians","None",&test);CHKERRQ(ierr);
79 ierr = PetscOptionsReal("-tao_test_hessian", "Threshold for element difference between hand-coded and finite difference being meaningful","None",threshold,&threshold,NULL);CHKERRQ(ierr);
80 ierr = PetscOptionsViewer("-tao_test_hessian_view","View difference between hand-coded and finite difference Hessians element entries","None",&mviewer,&format,&complete_print);CHKERRQ(ierr);
81 ierr = PetscOptionsEnd();CHKERRQ(ierr);
82 if (!test) PetscFunctionReturn(0);
83
84 ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr);
85 ierr = PetscViewerASCIIGetStdout(comm,&viewer);CHKERRQ(ierr);
86 ierr = PetscViewerASCIIGetTab(viewer, &tabs);CHKERRQ(ierr);
87 ierr = PetscViewerASCIISetTab(viewer, ((PetscObject)tao)->tablevel);CHKERRQ(ierr);
88 ierr = PetscViewerASCIIPrintf(viewer," ---------- Testing Hessian -------------\n");CHKERRQ(ierr);
89 if (!complete_print && !directionsprinted) {
90 ierr = PetscViewerASCIIPrintf(viewer," Run with -tao_test_hessian_view and optionally -tao_test_hessian <threshold> to show difference\n");CHKERRQ(ierr);
91 ierr = PetscViewerASCIIPrintf(viewer," of hand-coded and finite difference Hessian entries greater than <threshold>.\n");CHKERRQ(ierr);
92 }
93 if (!directionsprinted) {
94 ierr = PetscViewerASCIIPrintf(viewer," Testing hand-coded Hessian, if (for double precision runs) ||J - Jfd||_F/||J||_F is\n");CHKERRQ(ierr);
95 ierr = PetscViewerASCIIPrintf(viewer," O(1.e-8), the hand-coded Hessian is probably correct.\n");CHKERRQ(ierr);
96 directionsprinted = PETSC_TRUE;
97 }
98 if (complete_print) {
99 ierr = PetscViewerPushFormat(mviewer,format);CHKERRQ(ierr);
100 }
101
102 ierr = PetscObjectTypeCompare((PetscObject)tao->hessian,MATMFFD,&flg);CHKERRQ(ierr);
103 if (!flg) hessian = tao->hessian;
104 else hessian = tao->hessian_pre;
105
106 while (hessian) {
107 ierr = PetscObjectBaseTypeCompareAny((PetscObject)hessian,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");CHKERRQ(ierr);
108 if (flg) {
109 A = hessian;
110 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
111 } else {
112 ierr = MatComputeOperator(hessian,MATAIJ,&A);CHKERRQ(ierr);
113 }
114
115 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr);
116 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
117 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
118 ierr = MatSetSizes(B,m,n,M,N);CHKERRQ(ierr);
119 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
120 ierr = MatSetUp(B);CHKERRQ(ierr);
121 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
122
123 ierr = TaoDefaultComputeHessian(tao,x,B,B,NULL);CHKERRQ(ierr);
124
125 ierr = MatDuplicate(B,MAT_COPY_VALUES,&D);CHKERRQ(ierr);
126 ierr = MatAYPX(D,-1.0,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
127 ierr = MatNorm(D,NORM_FROBENIUS,&nrm);CHKERRQ(ierr);
128 ierr = MatNorm(A,NORM_FROBENIUS,&gnorm);CHKERRQ(ierr);
129 ierr = MatDestroy(&D);CHKERRQ(ierr);
130 if (!gnorm) gnorm = 1; /* just in case */
131 ierr = PetscViewerASCIIPrintf(viewer," ||H - Hfd||_F/||H||_F = %g, ||H - Hfd||_F = %g\n",(double)(nrm/gnorm),(double)nrm);CHKERRQ(ierr);
132
133 if (complete_print) {
134 ierr = PetscViewerASCIIPrintf(viewer," Hand-coded Hessian ----------\n");CHKERRQ(ierr);
135 ierr = MatView(A,mviewer);CHKERRQ(ierr);
136 ierr = PetscViewerASCIIPrintf(viewer," Finite difference Hessian ----------\n");CHKERRQ(ierr);
137 ierr = MatView(B,mviewer);CHKERRQ(ierr);
138 }
139
140 if (complete_print) {
141 PetscInt Istart, Iend, *ccols, bncols, cncols, j, row;
142 PetscScalar *cvals;
143 const PetscInt *bcols;
144 const PetscScalar *bvals;
145
146 ierr = MatAYPX(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
147 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr);
148 ierr = MatSetSizes(C,m,n,M,N);CHKERRQ(ierr);
149 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr);
150 ierr = MatSetUp(C);CHKERRQ(ierr);
151 ierr = MatSetOption(C,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
152 ierr = MatGetOwnershipRange(B,&Istart,&Iend);CHKERRQ(ierr);
153
154 for (row = Istart; row < Iend; row++) {
155 ierr = MatGetRow(B,row,&bncols,&bcols,&bvals);CHKERRQ(ierr);
156 ierr = PetscMalloc2(bncols,&ccols,bncols,&cvals);CHKERRQ(ierr);
157 for (j = 0, cncols = 0; j < bncols; j++) {
158 if (PetscAbsScalar(bvals[j]) > threshold) {
159 ccols[cncols] = bcols[j];
160 cvals[cncols] = bvals[j];
161 cncols += 1;
162 }
163 }
164 if (cncols) {
165 ierr = MatSetValues(C,1,&row,cncols,ccols,cvals,INSERT_VALUES);CHKERRQ(ierr);
166 }
167 ierr = MatRestoreRow(B,row,&bncols,&bcols,&bvals);CHKERRQ(ierr);
168 ierr = PetscFree2(ccols,cvals);CHKERRQ(ierr);
169 }
170 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
171 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
172 ierr = PetscViewerASCIIPrintf(viewer," Finite-difference minus hand-coded Hessian with tolerance %g ----------\n",(double)threshold);CHKERRQ(ierr);
173 ierr = MatView(C,mviewer);CHKERRQ(ierr);
174 ierr = MatDestroy(&C);CHKERRQ(ierr);
175 }
176 ierr = MatDestroy(&A);CHKERRQ(ierr);
177 ierr = MatDestroy(&B);CHKERRQ(ierr);
178
179 if (hessian != tao->hessian_pre) {
180 hessian = tao->hessian_pre;
181 ierr = PetscViewerASCIIPrintf(viewer," ---------- Testing Hessian for preconditioner -------------\n");CHKERRQ(ierr);
182 } else hessian = NULL;
183 }
184 if (complete_print) {
185 ierr = PetscViewerPopFormat(mviewer);CHKERRQ(ierr);
186 ierr = PetscViewerDestroy(&mviewer);CHKERRQ(ierr);
187 }
188 ierr = PetscViewerASCIISetTab(viewer,tabs);CHKERRQ(ierr);
189 PetscFunctionReturn(0);
190 }
191
192 /*@C
193 TaoComputeHessian - Computes the Hessian matrix that has been
194 set with TaoSetHessianRoutine().
195
196 Collective on Tao
197
198 Input Parameters:
199 + tao - the Tao solver context
200 - X - input vector
201
202 Output Parameters:
203 + H - Hessian matrix
204 - Hpre - Preconditioning matrix
205
206 Options Database Keys:
207 + -tao_test_hessian - compare the user provided Hessian with one compute via finite differences to check for errors
208 . -tao_test_hessian <numerical value> - display entries in the difference between the user provided Hessian and finite difference Hessian that are greater than a certain value to help users detect errors
209 - -tao_test_hessian_view - display the user provided Hessian, the finite difference Hessian and the difference between them to help users detect the location of errors in the user provided Hessian
210
211 Notes:
212 Most users should not need to explicitly call this routine, as it
213 is used internally within the minimization solvers.
214
215 TaoComputeHessian() is typically used within minimization
216 implementations, so most users would not generally call this routine
217 themselves.
218
219 Developer Notes:
220 The Hessian test mechanism follows SNESTestJacobian().
221
222 Level: developer
223
224 .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetHessianRoutine()
225 @*/
TaoComputeHessian(Tao tao,Vec X,Mat H,Mat Hpre)226 PetscErrorCode TaoComputeHessian(Tao tao, Vec X, Mat H, Mat Hpre)
227 {
228 PetscErrorCode ierr;
229
230 PetscFunctionBegin;
231 PetscValidHeaderSpecific(tao,TAO_CLASSID,1);
232 PetscValidHeaderSpecific(X, VEC_CLASSID,2);
233 PetscCheckSameComm(tao,1,X,2);
234 if (!tao->ops->computehessian) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetHessianRoutine() first");
235 ++tao->nhess;
236 ierr = VecLockReadPush(X);CHKERRQ(ierr);
237 ierr = PetscLogEventBegin(TAO_HessianEval,tao,X,H,Hpre);CHKERRQ(ierr);
238 PetscStackPush("Tao user Hessian function");
239 ierr = (*tao->ops->computehessian)(tao,X,H,Hpre,tao->user_hessP);CHKERRQ(ierr);
240 PetscStackPop;
241 ierr = PetscLogEventEnd(TAO_HessianEval,tao,X,H,Hpre);CHKERRQ(ierr);
242 ierr = VecLockReadPop(X);CHKERRQ(ierr);
243
244 ierr = TaoTestHessian(tao);CHKERRQ(ierr);
245 PetscFunctionReturn(0);
246 }
247
248 /*@C
249 TaoComputeJacobian - Computes the Jacobian matrix that has been
250 set with TaoSetJacobianRoutine().
251
252 Collective on Tao
253
254 Input Parameters:
255 + tao - the Tao solver context
256 - X - input vector
257
258 Output Parameters:
259 + J - Jacobian matrix
260 - Jpre - Preconditioning matrix
261
262 Notes:
263 Most users should not need to explicitly call this routine, as it
264 is used internally within the minimization solvers.
265
266 TaoComputeJacobian() is typically used within minimization
267 implementations, so most users would not generally call this routine
268 themselves.
269
270 Level: developer
271
272 .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianRoutine()
273 @*/
TaoComputeJacobian(Tao tao,Vec X,Mat J,Mat Jpre)274 PetscErrorCode TaoComputeJacobian(Tao tao, Vec X, Mat J, Mat Jpre)
275 {
276 PetscErrorCode ierr;
277
278 PetscFunctionBegin;
279 PetscValidHeaderSpecific(tao,TAO_CLASSID,1);
280 PetscValidHeaderSpecific(X, VEC_CLASSID,2);
281 PetscCheckSameComm(tao,1,X,2);
282 if (!tao->ops->computejacobian) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetJacobian() first");
283 ++tao->njac;
284 ierr = VecLockReadPush(X);CHKERRQ(ierr);
285 ierr = PetscLogEventBegin(TAO_JacobianEval,tao,X,J,Jpre);CHKERRQ(ierr);
286 PetscStackPush("Tao user Jacobian function");
287 ierr = (*tao->ops->computejacobian)(tao,X,J,Jpre,tao->user_jacP);CHKERRQ(ierr);
288 PetscStackPop;
289 ierr = PetscLogEventEnd(TAO_JacobianEval,tao,X,J,Jpre);CHKERRQ(ierr);
290 ierr = VecLockReadPop(X);CHKERRQ(ierr);
291 PetscFunctionReturn(0);
292 }
293
294 /*@C
295 TaoComputeResidualJacobian - Computes the least-squares residual Jacobian matrix that has been
296 set with TaoSetJacobianResidual().
297
298 Collective on Tao
299
300 Input Parameters:
301 + tao - the Tao solver context
302 - X - input vector
303
304 Output Parameters:
305 + J - Jacobian matrix
306 - Jpre - Preconditioning matrix
307
308 Notes:
309 Most users should not need to explicitly call this routine, as it
310 is used internally within the minimization solvers.
311
312 TaoComputeResidualJacobian() is typically used within least-squares
313 implementations, so most users would not generally call this routine
314 themselves.
315
316 Level: developer
317
318 .seealso: TaoComputeResidual(), TaoSetJacobianResidual()
319 @*/
TaoComputeResidualJacobian(Tao tao,Vec X,Mat J,Mat Jpre)320 PetscErrorCode TaoComputeResidualJacobian(Tao tao, Vec X, Mat J, Mat Jpre)
321 {
322 PetscErrorCode ierr;
323
324 PetscFunctionBegin;
325 PetscValidHeaderSpecific(tao,TAO_CLASSID,1);
326 PetscValidHeaderSpecific(X, VEC_CLASSID,2);
327 PetscCheckSameComm(tao,1,X,2);
328 if (!tao->ops->computeresidualjacobian) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetResidualJacobian() first");
329 ++tao->njac;
330 ierr = VecLockReadPush(X);CHKERRQ(ierr);
331 ierr = PetscLogEventBegin(TAO_JacobianEval,tao,X,J,Jpre);CHKERRQ(ierr);
332 PetscStackPush("Tao user least-squares residual Jacobian function");
333 ierr = (*tao->ops->computeresidualjacobian)(tao,X,J,Jpre,tao->user_lsjacP);CHKERRQ(ierr);
334 PetscStackPop;
335 ierr = PetscLogEventEnd(TAO_JacobianEval,tao,X,J,Jpre);CHKERRQ(ierr);
336 ierr = VecLockReadPop(X);CHKERRQ(ierr);
337 PetscFunctionReturn(0);
338 }
339
340 /*@C
341 TaoComputeJacobianState - Computes the Jacobian matrix that has been
342 set with TaoSetJacobianStateRoutine().
343
344 Collective on Tao
345
346 Input Parameters:
347 + tao - the Tao solver context
348 - X - input vector
349
350 Output Parameters:
351 + Jpre - Jacobian matrix
352 - Jinv - Preconditioning matrix
353
354 Notes:
355 Most users should not need to explicitly call this routine, as it
356 is used internally within the minimization solvers.
357
358 TaoComputeJacobianState() is typically used within minimization
359 implementations, so most users would not generally call this routine
360 themselves.
361
362 Level: developer
363
364 .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianStateRoutine(), TaoComputeJacobianDesign(), TaoSetStateDesignIS()
365 @*/
TaoComputeJacobianState(Tao tao,Vec X,Mat J,Mat Jpre,Mat Jinv)366 PetscErrorCode TaoComputeJacobianState(Tao tao, Vec X, Mat J, Mat Jpre, Mat Jinv)
367 {
368 PetscErrorCode ierr;
369
370 PetscFunctionBegin;
371 PetscValidHeaderSpecific(tao,TAO_CLASSID,1);
372 PetscValidHeaderSpecific(X, VEC_CLASSID,2);
373 PetscCheckSameComm(tao,1,X,2);
374 if (!tao->ops->computejacobianstate) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetJacobianState() first");
375 ++tao->njac_state;
376 ierr = VecLockReadPush(X);CHKERRQ(ierr);
377 ierr = PetscLogEventBegin(TAO_JacobianEval,tao,X,J,Jpre);CHKERRQ(ierr);
378 PetscStackPush("Tao user Jacobian(state) function");
379 ierr = (*tao->ops->computejacobianstate)(tao,X,J,Jpre,Jinv,tao->user_jac_stateP);CHKERRQ(ierr);
380 PetscStackPop;
381 ierr = PetscLogEventEnd(TAO_JacobianEval,tao,X,J,Jpre);CHKERRQ(ierr);
382 ierr = VecLockReadPop(X);CHKERRQ(ierr);
383 PetscFunctionReturn(0);
384 }
385
386 /*@C
387 TaoComputeJacobianDesign - Computes the Jacobian matrix that has been
388 set with TaoSetJacobianDesignRoutine().
389
390 Collective on Tao
391
392 Input Parameters:
393 + tao - the Tao solver context
394 - X - input vector
395
396 Output Parameters:
397 . J - Jacobian matrix
398
399 Notes:
400 Most users should not need to explicitly call this routine, as it
401 is used internally within the minimization solvers.
402
403 TaoComputeJacobianDesign() is typically used within minimization
404 implementations, so most users would not generally call this routine
405 themselves.
406
407 Level: developer
408
409 .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianDesignRoutine(), TaoComputeJacobianDesign(), TaoSetStateDesignIS()
410 @*/
TaoComputeJacobianDesign(Tao tao,Vec X,Mat J)411 PetscErrorCode TaoComputeJacobianDesign(Tao tao, Vec X, Mat J)
412 {
413 PetscErrorCode ierr;
414
415 PetscFunctionBegin;
416 PetscValidHeaderSpecific(tao,TAO_CLASSID,1);
417 PetscValidHeaderSpecific(X, VEC_CLASSID,2);
418 PetscCheckSameComm(tao,1,X,2);
419 if (!tao->ops->computejacobiandesign) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetJacobianDesign() first");
420 ++tao->njac_design;
421 ierr = VecLockReadPush(X);CHKERRQ(ierr);
422 ierr = PetscLogEventBegin(TAO_JacobianEval,tao,X,J,NULL);CHKERRQ(ierr);
423 PetscStackPush("Tao user Jacobian(design) function");
424 ierr = (*tao->ops->computejacobiandesign)(tao,X,J,tao->user_jac_designP);CHKERRQ(ierr);
425 PetscStackPop;
426 ierr = PetscLogEventEnd(TAO_JacobianEval,tao,X,J,NULL);CHKERRQ(ierr);
427 ierr = VecLockReadPop(X);CHKERRQ(ierr);
428 PetscFunctionReturn(0);
429 }
430
431 /*@C
432 TaoSetJacobianRoutine - Sets the function to compute the Jacobian as well as the location to store the matrix.
433
434 Logically collective on Tao
435
436 Input Parameters:
437 + tao - the Tao context
438 . J - Matrix used for the jacobian
439 . Jpre - Matrix that will be used operated on by preconditioner, can be same as J
440 . func - Jacobian evaluation routine
441 - ctx - [optional] user-defined context for private data for the
442 Jacobian evaluation routine (may be NULL)
443
444 Calling sequence of func:
445 $ func(Tao tao,Vec x,Mat J,Mat Jpre,void *ctx);
446
447 + tao - the Tao context
448 . x - input vector
449 . J - Jacobian matrix
450 . Jpre - preconditioning matrix, usually the same as J
451 - ctx - [optional] user-defined Jacobian context
452
453 Level: intermediate
454 @*/
TaoSetJacobianRoutine(Tao tao,Mat J,Mat Jpre,PetscErrorCode (* func)(Tao,Vec,Mat,Mat,void *),void * ctx)455 PetscErrorCode TaoSetJacobianRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void*), void *ctx)
456 {
457 PetscErrorCode ierr;
458
459 PetscFunctionBegin;
460 PetscValidHeaderSpecific(tao,TAO_CLASSID,1);
461 if (J) {
462 PetscValidHeaderSpecific(J,MAT_CLASSID,2);
463 PetscCheckSameComm(tao,1,J,2);
464 }
465 if (Jpre) {
466 PetscValidHeaderSpecific(Jpre,MAT_CLASSID,3);
467 PetscCheckSameComm(tao,1,Jpre,3);
468 }
469 if (ctx) {
470 tao->user_jacP = ctx;
471 }
472 if (func) {
473 tao->ops->computejacobian = func;
474 }
475 if (J) {
476 ierr = PetscObjectReference((PetscObject)J);CHKERRQ(ierr);
477 ierr = MatDestroy(&tao->jacobian);CHKERRQ(ierr);
478 tao->jacobian = J;
479 }
480 if (Jpre) {
481 ierr = PetscObjectReference((PetscObject)Jpre);CHKERRQ(ierr);
482 ierr = MatDestroy(&tao->jacobian_pre);CHKERRQ(ierr);
483 tao->jacobian_pre=Jpre;
484 }
485 PetscFunctionReturn(0);
486 }
487
488 /*@C
489 TaoSetJacobianResidualRoutine - Sets the function to compute the least-squares residual Jacobian as well as the
490 location to store the matrix.
491
492 Logically collective on Tao
493
494 Input Parameters:
495 + tao - the Tao context
496 . J - Matrix used for the jacobian
497 . Jpre - Matrix that will be used operated on by preconditioner, can be same as J
498 . func - Jacobian evaluation routine
499 - ctx - [optional] user-defined context for private data for the
500 Jacobian evaluation routine (may be NULL)
501
502 Calling sequence of func:
503 $ func(Tao tao,Vec x,Mat J,Mat Jpre,void *ctx);
504
505 + tao - the Tao context
506 . x - input vector
507 . J - Jacobian matrix
508 . Jpre - preconditioning matrix, usually the same as J
509 - ctx - [optional] user-defined Jacobian context
510
511 Level: intermediate
512 @*/
TaoSetJacobianResidualRoutine(Tao tao,Mat J,Mat Jpre,PetscErrorCode (* func)(Tao,Vec,Mat,Mat,void *),void * ctx)513 PetscErrorCode TaoSetJacobianResidualRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void*), void *ctx)
514 {
515 PetscErrorCode ierr;
516
517 PetscFunctionBegin;
518 PetscValidHeaderSpecific(tao,TAO_CLASSID,1);
519 if (J) {
520 PetscValidHeaderSpecific(J,MAT_CLASSID,2);
521 PetscCheckSameComm(tao,1,J,2);
522 }
523 if (Jpre) {
524 PetscValidHeaderSpecific(Jpre,MAT_CLASSID,3);
525 PetscCheckSameComm(tao,1,Jpre,3);
526 }
527 if (ctx) {
528 tao->user_lsjacP = ctx;
529 }
530 if (func) {
531 tao->ops->computeresidualjacobian = func;
532 }
533 if (J) {
534 ierr = PetscObjectReference((PetscObject)J);CHKERRQ(ierr);
535 ierr = MatDestroy(&tao->ls_jac);CHKERRQ(ierr);
536 tao->ls_jac = J;
537 }
538 if (Jpre) {
539 ierr = PetscObjectReference((PetscObject)Jpre);CHKERRQ(ierr);
540 ierr = MatDestroy(&tao->ls_jac_pre);CHKERRQ(ierr);
541 tao->ls_jac_pre=Jpre;
542 }
543 PetscFunctionReturn(0);
544 }
545
546 /*@C
547 TaoSetJacobianStateRoutine - Sets the function to compute the Jacobian
548 (and its inverse) of the constraint function with respect to the state variables.
549 Used only for pde-constrained optimization.
550
551 Logically collective on Tao
552
553 Input Parameters:
554 + tao - the Tao context
555 . J - Matrix used for the jacobian
556 . Jpre - Matrix that will be used operated on by PETSc preconditioner, can be same as J. Only used if Jinv is NULL
557 . Jinv - [optional] Matrix used to apply the inverse of the state jacobian. Use NULL to default to PETSc KSP solvers to apply the inverse.
558 . func - Jacobian evaluation routine
559 - ctx - [optional] user-defined context for private data for the
560 Jacobian evaluation routine (may be NULL)
561
562 Calling sequence of func:
563 $ func(Tao tao,Vec x,Mat J,Mat Jpre,Mat Jinv,void *ctx);
564
565 + tao - the Tao context
566 . x - input vector
567 . J - Jacobian matrix
568 . Jpre - preconditioner matrix, usually the same as J
569 . Jinv - inverse of J
570 - ctx - [optional] user-defined Jacobian context
571
572 Level: intermediate
573 .seealso: TaoComputeJacobianState(), TaoSetJacobianDesignRoutine(), TaoSetStateDesignIS()
574 @*/
TaoSetJacobianStateRoutine(Tao tao,Mat J,Mat Jpre,Mat Jinv,PetscErrorCode (* func)(Tao,Vec,Mat,Mat,Mat,void *),void * ctx)575 PetscErrorCode TaoSetJacobianStateRoutine(Tao tao, Mat J, Mat Jpre, Mat Jinv, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, Mat, void*), void *ctx)
576 {
577 PetscErrorCode ierr;
578
579 PetscFunctionBegin;
580 PetscValidHeaderSpecific(tao,TAO_CLASSID,1);
581 if (J) {
582 PetscValidHeaderSpecific(J,MAT_CLASSID,2);
583 PetscCheckSameComm(tao,1,J,2);
584 }
585 if (Jpre) {
586 PetscValidHeaderSpecific(Jpre,MAT_CLASSID,3);
587 PetscCheckSameComm(tao,1,Jpre,3);
588 }
589 if (Jinv) {
590 PetscValidHeaderSpecific(Jinv,MAT_CLASSID,4);
591 PetscCheckSameComm(tao,1,Jinv,4);
592 }
593 if (ctx) {
594 tao->user_jac_stateP = ctx;
595 }
596 if (func) {
597 tao->ops->computejacobianstate = func;
598 }
599 if (J) {
600 ierr = PetscObjectReference((PetscObject)J);CHKERRQ(ierr);
601 ierr = MatDestroy(&tao->jacobian_state);CHKERRQ(ierr);
602 tao->jacobian_state = J;
603 }
604 if (Jpre) {
605 ierr = PetscObjectReference((PetscObject)Jpre);CHKERRQ(ierr);
606 ierr = MatDestroy(&tao->jacobian_state_pre);CHKERRQ(ierr);
607 tao->jacobian_state_pre=Jpre;
608 }
609 if (Jinv) {
610 ierr = PetscObjectReference((PetscObject)Jinv);CHKERRQ(ierr);
611 ierr = MatDestroy(&tao->jacobian_state_inv);CHKERRQ(ierr);
612 tao->jacobian_state_inv=Jinv;
613 }
614 PetscFunctionReturn(0);
615 }
616
617 /*@C
618 TaoSetJacobianDesignRoutine - Sets the function to compute the Jacobian of
619 the constraint function with respect to the design variables. Used only for
620 pde-constrained optimization.
621
622 Logically collective on Tao
623
624 Input Parameters:
625 + tao - the Tao context
626 . J - Matrix used for the jacobian
627 . func - Jacobian evaluation routine
628 - ctx - [optional] user-defined context for private data for the
629 Jacobian evaluation routine (may be NULL)
630
631 Calling sequence of func:
632 $ func(Tao tao,Vec x,Mat J,void *ctx);
633
634 + tao - the Tao context
635 . x - input vector
636 . J - Jacobian matrix
637 - ctx - [optional] user-defined Jacobian context
638
639 Level: intermediate
640
641 .seealso: TaoComputeJacobianDesign(), TaoSetJacobianStateRoutine(), TaoSetStateDesignIS()
642 @*/
TaoSetJacobianDesignRoutine(Tao tao,Mat J,PetscErrorCode (* func)(Tao,Vec,Mat,void *),void * ctx)643 PetscErrorCode TaoSetJacobianDesignRoutine(Tao tao, Mat J, PetscErrorCode (*func)(Tao, Vec, Mat, void*), void *ctx)
644 {
645 PetscErrorCode ierr;
646
647 PetscFunctionBegin;
648 PetscValidHeaderSpecific(tao,TAO_CLASSID,1);
649 if (J) {
650 PetscValidHeaderSpecific(J,MAT_CLASSID,2);
651 PetscCheckSameComm(tao,1,J,2);
652 }
653 if (ctx) {
654 tao->user_jac_designP = ctx;
655 }
656 if (func) {
657 tao->ops->computejacobiandesign = func;
658 }
659 if (J) {
660 ierr = PetscObjectReference((PetscObject)J);CHKERRQ(ierr);
661 ierr = MatDestroy(&tao->jacobian_design);CHKERRQ(ierr);
662 tao->jacobian_design = J;
663 }
664 PetscFunctionReturn(0);
665 }
666
667 /*@
668 TaoSetStateDesignIS - Indicate to the Tao which variables in the
669 solution vector are state variables and which are design. Only applies to
670 pde-constrained optimization.
671
672 Logically Collective on Tao
673
674 Input Parameters:
675 + tao - The Tao context
676 . s_is - the index set corresponding to the state variables
677 - d_is - the index set corresponding to the design variables
678
679 Level: intermediate
680
681 .seealso: TaoSetJacobianStateRoutine(), TaoSetJacobianDesignRoutine()
682 @*/
TaoSetStateDesignIS(Tao tao,IS s_is,IS d_is)683 PetscErrorCode TaoSetStateDesignIS(Tao tao, IS s_is, IS d_is)
684 {
685 PetscErrorCode ierr;
686
687 PetscFunctionBegin;
688 ierr = PetscObjectReference((PetscObject)s_is);CHKERRQ(ierr);
689 ierr = ISDestroy(&tao->state_is);CHKERRQ(ierr);
690 tao->state_is = s_is;
691 ierr = PetscObjectReference((PetscObject)(d_is));CHKERRQ(ierr);
692 ierr = ISDestroy(&tao->design_is);CHKERRQ(ierr);
693 tao->design_is = d_is;
694 PetscFunctionReturn(0);
695 }
696
697 /*@C
698 TaoComputeJacobianEquality - Computes the Jacobian matrix that has been
699 set with TaoSetJacobianEqualityRoutine().
700
701 Collective on Tao
702
703 Input Parameters:
704 + tao - the Tao solver context
705 - X - input vector
706
707 Output Parameters:
708 + J - Jacobian matrix
709 - Jpre - Preconditioning matrix
710
711 Notes:
712 Most users should not need to explicitly call this routine, as it
713 is used internally within the minimization solvers.
714
715 Level: developer
716
717 .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianStateRoutine(), TaoComputeJacobianDesign(), TaoSetStateDesignIS()
718 @*/
TaoComputeJacobianEquality(Tao tao,Vec X,Mat J,Mat Jpre)719 PetscErrorCode TaoComputeJacobianEquality(Tao tao, Vec X, Mat J, Mat Jpre)
720 {
721 PetscErrorCode ierr;
722
723 PetscFunctionBegin;
724 PetscValidHeaderSpecific(tao,TAO_CLASSID,1);
725 PetscValidHeaderSpecific(X, VEC_CLASSID,2);
726 PetscCheckSameComm(tao,1,X,2);
727 if (!tao->ops->computejacobianequality) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetJacobianEquality() first");
728 ++tao->njac_equality;
729 ierr = VecLockReadPush(X);CHKERRQ(ierr);
730 ierr = PetscLogEventBegin(TAO_JacobianEval,tao,X,J,Jpre);CHKERRQ(ierr);
731 PetscStackPush("Tao user Jacobian(equality) function");
732 ierr = (*tao->ops->computejacobianequality)(tao,X,J,Jpre,tao->user_jac_equalityP);CHKERRQ(ierr);
733 PetscStackPop;
734 ierr = PetscLogEventEnd(TAO_JacobianEval,tao,X,J,Jpre);CHKERRQ(ierr);
735 ierr = VecLockReadPop(X);CHKERRQ(ierr);
736 PetscFunctionReturn(0);
737 }
738
739 /*@C
740 TaoComputeJacobianInequality - Computes the Jacobian matrix that has been
741 set with TaoSetJacobianInequalityRoutine().
742
743 Collective on Tao
744
745 Input Parameters:
746 + tao - the Tao solver context
747 - X - input vector
748
749 Output Parameters:
750 + J - Jacobian matrix
751 - Jpre - Preconditioning matrix
752
753 Notes:
754 Most users should not need to explicitly call this routine, as it
755 is used internally within the minimization solvers.
756
757 Level: developer
758
759 .seealso: TaoComputeObjective(), TaoComputeObjectiveAndGradient(), TaoSetJacobianStateRoutine(), TaoComputeJacobianDesign(), TaoSetStateDesignIS()
760 @*/
TaoComputeJacobianInequality(Tao tao,Vec X,Mat J,Mat Jpre)761 PetscErrorCode TaoComputeJacobianInequality(Tao tao, Vec X, Mat J, Mat Jpre)
762 {
763 PetscErrorCode ierr;
764
765 PetscFunctionBegin;
766 PetscValidHeaderSpecific(tao,TAO_CLASSID,1);
767 PetscValidHeaderSpecific(X, VEC_CLASSID,2);
768 PetscCheckSameComm(tao,1,X,2);
769 if (!tao->ops->computejacobianinequality) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TaoSetJacobianInequality() first");
770 ++tao->njac_inequality;
771 ierr = VecLockReadPush(X);CHKERRQ(ierr);
772 ierr = PetscLogEventBegin(TAO_JacobianEval,tao,X,J,Jpre);CHKERRQ(ierr);
773 PetscStackPush("Tao user Jacobian(inequality) function");
774 ierr = (*tao->ops->computejacobianinequality)(tao,X,J,Jpre,tao->user_jac_inequalityP);CHKERRQ(ierr);
775 PetscStackPop;
776 ierr = PetscLogEventEnd(TAO_JacobianEval,tao,X,J,Jpre);CHKERRQ(ierr);
777 ierr = VecLockReadPop(X);CHKERRQ(ierr);
778 PetscFunctionReturn(0);
779 }
780
781 /*@C
782 TaoSetJacobianEqualityRoutine - Sets the function to compute the Jacobian
783 (and its inverse) of the constraint function with respect to the equality variables.
784 Used only for pde-constrained optimization.
785
786 Logically collective on Tao
787
788 Input Parameters:
789 + tao - the Tao context
790 . J - Matrix used for the jacobian
791 . Jpre - Matrix that will be used operated on by PETSc preconditioner, can be same as J.
792 . func - Jacobian evaluation routine
793 - ctx - [optional] user-defined context for private data for the
794 Jacobian evaluation routine (may be NULL)
795
796 Calling sequence of func:
797 $ func(Tao tao,Vec x,Mat J,Mat Jpre,void *ctx);
798
799 + tao - the Tao context
800 . x - input vector
801 . J - Jacobian matrix
802 . Jpre - preconditioner matrix, usually the same as J
803 - ctx - [optional] user-defined Jacobian context
804
805 Level: intermediate
806
807 .seealso: TaoComputeJacobianEquality(), TaoSetJacobianDesignRoutine(), TaoSetEqualityDesignIS()
808 @*/
TaoSetJacobianEqualityRoutine(Tao tao,Mat J,Mat Jpre,PetscErrorCode (* func)(Tao,Vec,Mat,Mat,void *),void * ctx)809 PetscErrorCode TaoSetJacobianEqualityRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat,void*), void *ctx)
810 {
811 PetscErrorCode ierr;
812
813 PetscFunctionBegin;
814 PetscValidHeaderSpecific(tao,TAO_CLASSID,1);
815 if (J) {
816 PetscValidHeaderSpecific(J,MAT_CLASSID,2);
817 PetscCheckSameComm(tao,1,J,2);
818 }
819 if (Jpre) {
820 PetscValidHeaderSpecific(Jpre,MAT_CLASSID,3);
821 PetscCheckSameComm(tao,1,Jpre,3);
822 }
823 if (ctx) {
824 tao->user_jac_equalityP = ctx;
825 }
826 if (func) {
827 tao->ops->computejacobianequality = func;
828 }
829 if (J) {
830 ierr = PetscObjectReference((PetscObject)J);CHKERRQ(ierr);
831 ierr = MatDestroy(&tao->jacobian_equality);CHKERRQ(ierr);
832 tao->jacobian_equality = J;
833 }
834 if (Jpre) {
835 ierr = PetscObjectReference((PetscObject)Jpre);CHKERRQ(ierr);
836 ierr = MatDestroy(&tao->jacobian_equality_pre);CHKERRQ(ierr);
837 tao->jacobian_equality_pre=Jpre;
838 }
839 PetscFunctionReturn(0);
840 }
841
842 /*@C
843 TaoSetJacobianInequalityRoutine - Sets the function to compute the Jacobian
844 (and its inverse) of the constraint function with respect to the inequality variables.
845 Used only for pde-constrained optimization.
846
847 Logically collective on Tao
848
849 Input Parameters:
850 + tao - the Tao context
851 . J - Matrix used for the jacobian
852 . Jpre - Matrix that will be used operated on by PETSc preconditioner, can be same as J.
853 . func - Jacobian evaluation routine
854 - ctx - [optional] user-defined context for private data for the
855 Jacobian evaluation routine (may be NULL)
856
857 Calling sequence of func:
858 $ func(Tao tao,Vec x,Mat J,Mat Jpre,void *ctx);
859
860 + tao - the Tao context
861 . x - input vector
862 . J - Jacobian matrix
863 . Jpre - preconditioner matrix, usually the same as J
864 - ctx - [optional] user-defined Jacobian context
865
866 Level: intermediate
867
868 .seealso: TaoComputeJacobianInequality(), TaoSetJacobianDesignRoutine(), TaoSetInequalityDesignIS()
869 @*/
TaoSetJacobianInequalityRoutine(Tao tao,Mat J,Mat Jpre,PetscErrorCode (* func)(Tao,Vec,Mat,Mat,void *),void * ctx)870 PetscErrorCode TaoSetJacobianInequalityRoutine(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat,void*), void *ctx)
871 {
872 PetscErrorCode ierr;
873
874 PetscFunctionBegin;
875 PetscValidHeaderSpecific(tao,TAO_CLASSID,1);
876 if (J) {
877 PetscValidHeaderSpecific(J,MAT_CLASSID,2);
878 PetscCheckSameComm(tao,1,J,2);
879 }
880 if (Jpre) {
881 PetscValidHeaderSpecific(Jpre,MAT_CLASSID,3);
882 PetscCheckSameComm(tao,1,Jpre,3);
883 }
884 if (ctx) {
885 tao->user_jac_inequalityP = ctx;
886 }
887 if (func) {
888 tao->ops->computejacobianinequality = func;
889 }
890 if (J) {
891 ierr = PetscObjectReference((PetscObject)J);CHKERRQ(ierr);
892 ierr = MatDestroy(&tao->jacobian_inequality);CHKERRQ(ierr);
893 tao->jacobian_inequality = J;
894 }
895 if (Jpre) {
896 ierr = PetscObjectReference((PetscObject)Jpre);CHKERRQ(ierr);
897 ierr = MatDestroy(&tao->jacobian_inequality_pre);CHKERRQ(ierr);
898 tao->jacobian_inequality_pre=Jpre;
899 }
900 PetscFunctionReturn(0);
901 }
902