xref: /petsc/src/tao/constrained/tutorials/ex1.c (revision 800f99ff9e85495c69e9e5819c0be0dbd8cbc57c)
1 /* Program usage: mpiexec -n 2 ./ex1 [-help] [all TAO options] */
2 
3 /* ----------------------------------------------------------------------
4 min f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
5 s.t.  x0^2 + x1 - 2 = 0
6       0  <= x0^2 - x1 <= 1
7       -1 <= x0, x1 <= 2
8 -->
9       g(x)  = 0
10       h(x) >= 0
11       -1 <= x0, x1 <= 2
12 where
13       g(x) = x0^2 + x1 - 2
14       h(x) = [x0^2 - x1
15               1 -(x0^2 - x1)]
16 ---------------------------------------------------------------------- */
17 
18 #include <petsctao.h>
19 
20 static  char help[]= "Solves constrained optimiztion problem using pdipm.\n\
21 Input parameters include:\n\
22   -tao_type pdipm    : sets Tao solver\n\
23   -no_eq             : removes the equaility constraints from the problem\n\
24   -init_view         : view initial object setup\n\
25   -snes_fd           : snes with finite difference Jacobian (needed for pdipm)\n\
26   -snes_compare_explicit : compare user Jacobian with finite difference Jacobian \n\
27   -tao_cmonitor      : convergence monitor with constraint norm \n\
28   -tao_view_solution : view exact solution at each iteration\n\
29   Note: external package MUMPS is required to run pdipm in parallel. This is designed for a maximum of 2 processors, the code will error if size > 2.\n";
30 
31 /*
32    User-defined application context - contains data needed by the application
33 */
34 typedef struct {
35   PetscInt   n;  /* Global length of x */
36   PetscInt   ne; /* Global number of equality constraints */
37   PetscInt   ni; /* Global number of inequality constraints */
38   PetscBool  noeqflag, initview;
39   Vec        x,xl,xu;
40   Vec        ce,ci,bl,bu,Xseq;
41   Mat        Ae,Ai,H;
42   VecScatter scat;
43 } AppCtx;
44 
45 /* -------- User-defined Routines --------- */
46 PetscErrorCode InitializeProblem(AppCtx *);
47 PetscErrorCode DestroyProblem(AppCtx *);
48 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
49 PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
50 PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*);
51 PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*);
52 PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*);
53 PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*);
54 
55 PetscErrorCode main(int argc,char **argv)
56 {
57   Tao            tao;
58   KSP            ksp;
59   PC             pc;
60   AppCtx         user;  /* application context */
61   Vec            x,G,CI,CE;
62   PetscMPIInt    size;
63   TaoType        type;
64   PetscReal      f;
65   PetscBool      pdipm;
66 
67   PetscFunctionBeginUser;
68   PetscCall(PetscInitialize(&argc,&argv,(char*)0,help));
69   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
70   PetscCheck(size <= 2,PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,"More than 2 processors detected. Example written to use max of 2 processors.");
71 
72   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"---- Constrained Problem -----\n"));
73   PetscCall(InitializeProblem(&user)); /* sets up problem, function below */
74   PetscCall(TaoCreate(PETSC_COMM_WORLD,&tao));
75   PetscCall(TaoSetType(tao,TAOPDIPM));
76   PetscCall(TaoSetSolution(tao,user.x)); /* gets solution vector from problem */
77   PetscCall(TaoSetVariableBounds(tao,user.xl,user.xu)); /* sets lower upper bounds from given solution */
78   PetscCall(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void*)&user));
79 
80   if (!user.noeqflag) {
81     PetscCall(TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user));
82   }
83   PetscCall(TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user));
84   if (!user.noeqflag) {
85     PetscCall(TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user)); /* equality jacobian */
86   }
87   PetscCall(TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user)); /* inequality jacobian */
88   PetscCall(TaoSetTolerances(tao,1.e-6,1.e-6,1.e-6));
89   PetscCall(TaoSetConstraintTolerances(tao,1.e-6,1.e-6));
90 
91   PetscCall(TaoGetKSP(tao,&ksp));
92   PetscCall(KSPGetPC(ksp,&pc));
93   PetscCall(PCSetType(pc,PCCHOLESKY));
94   /*
95       This algorithm produces matrices with zeros along the diagonal therefore we use
96     MUMPS which provides solver for indefinite matrices
97   */
98 #if defined(PETSC_HAVE_MUMPS)
99   PetscCall(PCFactorSetMatSolverType(pc,MATSOLVERMUMPS));  /* requires mumps to solve pdipm */
100 #else
101   PetscCheck(size == 1,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Requires an external package that supports parallel PCCHOLESKY, e.g., MUMPS.");
102 #endif
103   PetscCall(KSPSetType(ksp,KSPPREONLY));
104   PetscCall(KSPSetFromOptions(ksp));
105 
106   PetscCall(TaoSetFromOptions(tao));
107   PetscCall(TaoGetType(tao,&type));
108   PetscCall(PetscObjectTypeCompare((PetscObject)tao,TAOPDIPM,&pdipm));
109   if (pdipm) {
110     PetscCall(TaoSetHessian(tao,user.H,user.H,FormHessian,(void*)&user));
111     if (user.initview) {
112       PetscCall(TaoSetUp(tao));
113       PetscCall(VecDuplicate(user.x, &G));
114       PetscCall(FormFunctionGradient(tao, user.x, &f, G, (void*)&user));
115       PetscCall(FormHessian(tao, user.x, user.H, user.H, (void*)&user));
116       PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD));
117       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nInitial point X:\n"));
118       PetscCall(VecView(user.x, PETSC_VIEWER_STDOUT_WORLD));
119       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nInitial objective f(x) = %g\n",(double)f));
120       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nInitial gradient and Hessian:\n"));
121       PetscCall(VecView(G, PETSC_VIEWER_STDOUT_WORLD));
122       PetscCall(MatView(user.H, PETSC_VIEWER_STDOUT_WORLD));
123       PetscCall(VecDestroy(&G));
124       PetscCall(FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void*)&user));
125       PetscCall(MatCreateVecs(user.Ai, NULL, &CI));
126       PetscCall(FormInequalityConstraints(tao, user.x, CI, (void*)&user));
127       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nInitial inequality constraints and Jacobian:\n"));
128       PetscCall(VecView(CI, PETSC_VIEWER_STDOUT_WORLD));
129       PetscCall(MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD));
130       PetscCall(VecDestroy(&CI));
131       if (!user.noeqflag) {
132         PetscCall(FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void*)&user));
133         PetscCall(MatCreateVecs(user.Ae, NULL, &CE));
134         PetscCall(FormEqualityConstraints(tao, user.x, CE, (void*)&user));
135         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nInitial equality constraints and Jacobian:\n"));
136         PetscCall(VecView(CE, PETSC_VIEWER_STDOUT_WORLD));
137         PetscCall(MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD));
138         PetscCall(VecDestroy(&CE));
139       }
140       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n"));
141       PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD));
142     }
143   }
144 
145   PetscCall(TaoSolve(tao));
146   PetscCall(TaoGetSolution(tao,&x));
147   PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD));
148 
149   /* Free objects */
150   PetscCall(DestroyProblem(&user));
151   PetscCall(TaoDestroy(&tao));
152   PetscCall(PetscFinalize());
153   return 0;
154 }
155 
156 PetscErrorCode InitializeProblem(AppCtx *user)
157 {
158   PetscMPIInt    size;
159   PetscMPIInt    rank;
160   PetscInt       nloc,neloc,niloc;
161 
162   PetscFunctionBegin;
163   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
164   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
165   user->noeqflag = PETSC_FALSE;
166   PetscCall(PetscOptionsGetBool(NULL,NULL,"-no_eq",&user->noeqflag,NULL));
167   user->initview = PETSC_FALSE;
168   PetscCall(PetscOptionsGetBool(NULL,NULL,"-init_view",&user->initview,NULL));
169 
170   if (!user->noeqflag) {
171     /* Tell user the correct solution, not an error checking */
172     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n"));
173   }
174 
175   /* create vector x and set initial values */
176   user->n = 2; /* global length */
177   nloc = (size==1)?user->n:1;
178   PetscCall(VecCreate(PETSC_COMM_WORLD,&user->x));
179   PetscCall(VecSetSizes(user->x,nloc,user->n));
180   PetscCall(VecSetFromOptions(user->x));
181   PetscCall(VecSet(user->x,0));
182 
183   /* create and set lower and upper bound vectors */
184   PetscCall(VecDuplicate(user->x,&user->xl));
185   PetscCall(VecDuplicate(user->x,&user->xu));
186   PetscCall(VecSet(user->xl,-1.0));
187   PetscCall(VecSet(user->xu,2.0));
188 
189   /* create scater to zero */
190   PetscCall(VecScatterCreateToZero(user->x,&user->scat,&user->Xseq));
191 
192   user->ne = 1;
193   user->ni = 2;
194   neloc = (rank==0)?user->ne:0;
195   niloc = (size==1)?user->ni:1;
196 
197   if (!user->noeqflag) {
198     PetscCall(VecCreate(PETSC_COMM_WORLD,&user->ce)); /* a 1x1 vec for equality constraints */
199     PetscCall(VecSetSizes(user->ce,neloc,user->ne));
200     PetscCall(VecSetFromOptions(user->ce));
201     PetscCall(VecSetUp(user->ce));
202   }
203 
204   PetscCall(VecCreate(PETSC_COMM_WORLD,&user->ci)); /* a 2x1 vec for inequality constraints */
205   PetscCall(VecSetSizes(user->ci,niloc,user->ni));
206   PetscCall(VecSetFromOptions(user->ci));
207   PetscCall(VecSetUp(user->ci));
208 
209   /* nexn & nixn matricies for equally and inequalty constraints */
210   if (!user->noeqflag) {
211     PetscCall(MatCreate(PETSC_COMM_WORLD,&user->Ae));
212     PetscCall(MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n));
213     PetscCall(MatSetFromOptions(user->Ae));
214     PetscCall(MatSetUp(user->Ae));
215   }
216 
217   PetscCall(MatCreate(PETSC_COMM_WORLD,&user->Ai));
218   PetscCall(MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n));
219   PetscCall(MatSetFromOptions(user->Ai));
220   PetscCall(MatSetUp(user->Ai));
221 
222   PetscCall(MatCreate(PETSC_COMM_WORLD,&user->H));
223   PetscCall(MatSetSizes(user->H,nloc,nloc,user->n,user->n));
224   PetscCall(MatSetFromOptions(user->H));
225   PetscCall(MatSetUp(user->H));
226   PetscFunctionReturn(0);
227 }
228 
229 PetscErrorCode DestroyProblem(AppCtx *user)
230 {
231   PetscFunctionBegin;
232   if (!user->noeqflag) {
233     PetscCall(MatDestroy(&user->Ae));
234   }
235   PetscCall(MatDestroy(&user->Ai));
236   PetscCall(MatDestroy(&user->H));
237 
238   PetscCall(VecDestroy(&user->x));
239   if (!user->noeqflag) {
240     PetscCall(VecDestroy(&user->ce));
241   }
242   PetscCall(VecDestroy(&user->ci));
243   PetscCall(VecDestroy(&user->xl));
244   PetscCall(VecDestroy(&user->xu));
245   PetscCall(VecDestroy(&user->Xseq));
246   PetscCall(VecScatterDestroy(&user->scat));
247   PetscFunctionReturn(0);
248 }
249 
250 /* Evaluate
251    f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
252    G = grad f = [2*(x0 - 2) - 2;
253                  2*(x1 - 2) - 2]
254 */
255 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
256 {
257   PetscScalar       g;
258   const PetscScalar *x;
259   MPI_Comm          comm;
260   PetscMPIInt       rank;
261   PetscReal         fin;
262   AppCtx            *user=(AppCtx*)ctx;
263   Vec               Xseq=user->Xseq;
264   VecScatter        scat=user->scat;
265 
266   PetscFunctionBegin;
267   PetscCall(PetscObjectGetComm((PetscObject)tao,&comm));
268   PetscCallMPI(MPI_Comm_rank(comm,&rank));
269 
270   PetscCall(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
271   PetscCall(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
272 
273   fin = 0.0;
274   if (rank == 0) {
275     PetscCall(VecGetArrayRead(Xseq,&x));
276     fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]);
277     g = 2.0*(x[0]-2.0) - 2.0;
278     PetscCall(VecSetValue(G,0,g,INSERT_VALUES));
279     g = 2.0*(x[1]-2.0) - 2.0;
280     PetscCall(VecSetValue(G,1,g,INSERT_VALUES));
281     PetscCall(VecRestoreArrayRead(Xseq,&x));
282   }
283   PetscCallMPI(MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm));
284   PetscCall(VecAssemblyBegin(G));
285   PetscCall(VecAssemblyEnd(G));
286   PetscFunctionReturn(0);
287 }
288 
289 /* Evaluate
290    H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)]
291      = [ 2*(1+de[0]-di[0]+di[1]), 0;
292                    0,             2]
293 */
294 PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx)
295 {
296   AppCtx            *user=(AppCtx*)ctx;
297   Vec               DE,DI;
298   const PetscScalar *de,*di;
299   PetscInt          zero=0,one=1;
300   PetscScalar       two=2.0;
301   PetscScalar       val=0.0;
302   Vec               Deseq,Diseq;
303   VecScatter        Descat,Discat;
304   PetscMPIInt       rank;
305   MPI_Comm          comm;
306 
307   PetscFunctionBegin;
308   PetscCall(TaoGetDualVariables(tao,&DE,&DI));
309 
310   PetscCall(PetscObjectGetComm((PetscObject)tao,&comm));
311   PetscCallMPI(MPI_Comm_rank(comm,&rank));
312 
313   if (!user->noeqflag) {
314    PetscCall(VecScatterCreateToZero(DE,&Descat,&Deseq));
315    PetscCall(VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD));
316    PetscCall(VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD));
317   }
318   PetscCall(VecScatterCreateToZero(DI,&Discat,&Diseq));
319   PetscCall(VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD));
320   PetscCall(VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD));
321 
322   if (rank == 0) {
323     if (!user->noeqflag) {
324       PetscCall(VecGetArrayRead(Deseq,&de));  /* places equality constraint dual into array */
325     }
326     PetscCall(VecGetArrayRead(Diseq,&di));  /* places inequality constraint dual into array */
327 
328     if (!user->noeqflag) {
329       val = 2.0 * (1 + de[0] - di[0] + di[1]);
330       PetscCall(VecRestoreArrayRead(Deseq,&de));
331       PetscCall(VecRestoreArrayRead(Diseq,&di));
332     } else {
333       val = 2.0 * (1 - di[0] + di[1]);
334     }
335     PetscCall(VecRestoreArrayRead(Diseq,&di));
336     PetscCall(MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES));
337     PetscCall(MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES));
338   }
339   PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY));
340   PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY));
341   if (!user->noeqflag) {
342     PetscCall(VecScatterDestroy(&Descat));
343     PetscCall(VecDestroy(&Deseq));
344   }
345   PetscCall(VecScatterDestroy(&Discat));
346   PetscCall(VecDestroy(&Diseq));
347   PetscFunctionReturn(0);
348 }
349 
350 /* Evaluate
351    h = [ x0^2 - x1;
352          1 -(x0^2 - x1)]
353 */
354 PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx)
355 {
356   const PetscScalar *x;
357   PetscScalar       ci;
358   MPI_Comm          comm;
359   PetscMPIInt       rank;
360   AppCtx            *user=(AppCtx*)ctx;
361   Vec               Xseq=user->Xseq;
362   VecScatter        scat=user->scat;
363 
364   PetscFunctionBegin;
365   PetscCall(PetscObjectGetComm((PetscObject)tao,&comm));
366   PetscCallMPI(MPI_Comm_rank(comm,&rank));
367 
368   PetscCall(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
369   PetscCall(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
370 
371   if (rank == 0) {
372     PetscCall(VecGetArrayRead(Xseq,&x));
373     ci = x[0]*x[0] - x[1];
374     PetscCall(VecSetValue(CI,0,ci,INSERT_VALUES));
375     ci = -x[0]*x[0] + x[1] + 1.0;
376     PetscCall(VecSetValue(CI,1,ci,INSERT_VALUES));
377     PetscCall(VecRestoreArrayRead(Xseq,&x));
378   }
379   PetscCall(VecAssemblyBegin(CI));
380   PetscCall(VecAssemblyEnd(CI));
381   PetscFunctionReturn(0);
382 }
383 
384 /* Evaluate
385    g = [ x0^2 + x1 - 2]
386 */
387 PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx)
388 {
389   const PetscScalar *x;
390   PetscScalar       ce;
391   MPI_Comm          comm;
392   PetscMPIInt       rank;
393   AppCtx            *user=(AppCtx*)ctx;
394   Vec               Xseq=user->Xseq;
395   VecScatter        scat=user->scat;
396 
397   PetscFunctionBegin;
398   PetscCall(PetscObjectGetComm((PetscObject)tao,&comm));
399   PetscCallMPI(MPI_Comm_rank(comm,&rank));
400 
401   PetscCall(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
402   PetscCall(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
403 
404   if (rank == 0) {
405     PetscCall(VecGetArrayRead(Xseq,&x));
406     ce = x[0]*x[0] + x[1] - 2.0;
407     PetscCall(VecSetValue(CE,0,ce,INSERT_VALUES));
408     PetscCall(VecRestoreArrayRead(Xseq,&x));
409   }
410   PetscCall(VecAssemblyBegin(CE));
411   PetscCall(VecAssemblyEnd(CE));
412   PetscFunctionReturn(0);
413 }
414 
415 /*
416   grad h = [  2*x0, -1;
417              -2*x0,  1]
418 */
419 PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre,  void *ctx)
420 {
421   AppCtx            *user=(AppCtx*)ctx;
422   PetscInt          zero=0,one=1,cols[2];
423   PetscScalar       vals[2];
424   const PetscScalar *x;
425   Vec               Xseq=user->Xseq;
426   VecScatter        scat=user->scat;
427   MPI_Comm          comm;
428   PetscMPIInt       rank;
429 
430   PetscFunctionBegin;
431   PetscCall(PetscObjectGetComm((PetscObject)tao,&comm));
432   PetscCallMPI(MPI_Comm_rank(comm,&rank));
433   PetscCall(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
434   PetscCall(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
435 
436   PetscCall(VecGetArrayRead(Xseq,&x));
437   if (rank == 0) {
438     cols[0] = 0; cols[1] = 1;
439     vals[0] = 2*x[0]; vals[1] = -1.0;
440     PetscCall(MatSetValues(JI,1,&zero,2,cols,vals,INSERT_VALUES));
441     vals[0] = -2*x[0]; vals[1] = 1.0;
442     PetscCall(MatSetValues(JI,1,&one,2,cols,vals,INSERT_VALUES));
443   }
444   PetscCall(VecRestoreArrayRead(Xseq,&x));
445   PetscCall(MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY));
446   PetscCall(MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY));
447   PetscFunctionReturn(0);
448 }
449 
450 /*
451   grad g = [2*x0
452              1.0 ]
453 */
454 PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx)
455 {
456   PetscInt          zero=0,cols[2];
457   PetscScalar       vals[2];
458   const PetscScalar *x;
459   PetscMPIInt       rank;
460   MPI_Comm          comm;
461 
462   PetscFunctionBegin;
463   PetscCall(PetscObjectGetComm((PetscObject)tao,&comm));
464   PetscCallMPI(MPI_Comm_rank(comm,&rank));
465 
466   if (rank == 0) {
467     PetscCall(VecGetArrayRead(X,&x));
468     cols[0] = 0;       cols[1] = 1;
469     vals[0] = 2*x[0];  vals[1] = 1.0;
470     PetscCall(MatSetValues(JE,1,&zero,2,cols,vals,INSERT_VALUES));
471     PetscCall(VecRestoreArrayRead(X,&x));
472   }
473   PetscCall(MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY));
474   PetscCall(MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY));
475   PetscFunctionReturn(0);
476 }
477 
478 /*TEST
479 
480    build:
481       requires: !complex !defined(PETSC_USE_CXX)
482 
483    test:
484       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
485       requires: mumps
486 
487    test:
488       suffix: 2
489       nsize: 2
490       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
491       requires: mumps
492 
493    test:
494       suffix: 3
495       args: -tao_converged_reason -no_eq
496       requires: mumps
497 
498    test:
499       suffix: 4
500       nsize: 2
501       args: -tao_converged_reason -no_eq
502       requires: mumps
503 
504    test:
505       suffix: 5
506       args: -tao_cmonitor -tao_type almm
507       requires: mumps
508 
509    test:
510       suffix: 6
511       args: -tao_cmonitor -tao_type almm -tao_almm_type phr
512       requires: mumps
513 
514    test:
515       suffix: 7
516       nsize: 2
517       requires: mumps
518       args: -tao_cmonitor -tao_type almm
519 
520    test:
521       suffix: 8
522       nsize: 2
523       requires: cuda mumps
524       args: -tao_cmonitor -tao_type almm -vec_type cuda -mat_type aijcusparse
525 
526    test:
527       suffix: 9
528       nsize: 2
529       args: -tao_cmonitor -tao_type almm -no_eq
530       requires: mumps
531 
532    test:
533       suffix: 10
534       nsize: 2
535       args: -tao_cmonitor -tao_type almm -tao_almm_type phr -no_eq
536       requires: mumps
537 
538 TEST*/
539