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