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