xref: /petsc/src/tao/constrained/tutorials/ex1.c (revision 628da978d1ab16495f6ac818252fe8b2f4352128)
1aad13602SShrirang Abhyankar /* Program usage: mpiexec -n 2 ./ex1 [-help] [all TAO options] */
2aad13602SShrirang Abhyankar 
3aad13602SShrirang Abhyankar /* ----------------------------------------------------------------------
4*628da978Sresundermann min f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
5aad13602SShrirang Abhyankar s.t.  x0^2 + x1 - 2 = 0
6aad13602SShrirang Abhyankar       0  <= x0^2 - x1 <= 1
7aad13602SShrirang Abhyankar       -1 <= x0, x1 <= 2
8*628da978Sresundermann -->
9*628da978Sresundermann       g(x)  = 0
10*628da978Sresundermann       h(x) >= 0
11*628da978Sresundermann       -1 <= x0, x1 <= 2
12*628da978Sresundermann where
13*628da978Sresundermann       g(x) = x0^2 + x1 - 2
14*628da978Sresundermann       h(x) = [x0^2 - x1
15*628da978Sresundermann               1 -(x0^2 - x1)]
16aad13602SShrirang Abhyankar ---------------------------------------------------------------------- */
17aad13602SShrirang Abhyankar 
18aad13602SShrirang Abhyankar #include <petsctao.h>
19aad13602SShrirang Abhyankar 
20aad13602SShrirang Abhyankar static  char help[]= "Solves constrained optimiztion problem using pdipm.\n\
21aad13602SShrirang Abhyankar Input parameters include:\n\
22aad13602SShrirang Abhyankar   -tao_type pdipm    : sets Tao solver\n\
238e58fa1dSresundermann   -no_eq             : removes the equaility constraints from the problem\n\
24aad13602SShrirang Abhyankar   -snes_fd           : snes with finite difference Jacobian (needed for pdipm)\n\
25*628da978Sresundermann   -snes_compare_explicit : compare user Jacobian with finite difference Jacobian \n\
26aad13602SShrirang Abhyankar   -tao_cmonitor      : convergence monitor with constraint norm \n\
27aad13602SShrirang Abhyankar   -tao_view_solution : view exact solution at each itteration\n\
28*628da978Sresundermann   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";
29aad13602SShrirang Abhyankar 
30aad13602SShrirang Abhyankar /*
31*628da978Sresundermann    User-defined application context - contains data needed by the application
32aad13602SShrirang Abhyankar */
33aad13602SShrirang Abhyankar typedef struct {
34aad13602SShrirang Abhyankar   PetscInt   n;  /* Global length of x */
35aad13602SShrirang Abhyankar   PetscInt   ne; /* Global number of equality constraints */
36aad13602SShrirang Abhyankar   PetscInt   ni; /* Global number of inequality constraints */
378e58fa1dSresundermann   PetscBool  noeqflag;
38aad13602SShrirang Abhyankar   Vec        x,xl,xu;
39aad13602SShrirang Abhyankar   Vec        ce,ci,bl,bu,Xseq;
40aad13602SShrirang Abhyankar   Mat        Ae,Ai,H;
41aad13602SShrirang Abhyankar   VecScatter scat;
42aad13602SShrirang Abhyankar } AppCtx;
43aad13602SShrirang Abhyankar 
44aad13602SShrirang Abhyankar /* -------- User-defined Routines --------- */
45aad13602SShrirang Abhyankar PetscErrorCode InitializeProblem(AppCtx *);
46aad13602SShrirang Abhyankar PetscErrorCode DestroyProblem(AppCtx *);
47aad13602SShrirang Abhyankar PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
48aad13602SShrirang Abhyankar PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
49aad13602SShrirang Abhyankar PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*);
50aad13602SShrirang Abhyankar PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*);
51aad13602SShrirang Abhyankar PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*);
52aad13602SShrirang Abhyankar PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*);
53aad13602SShrirang Abhyankar 
54aad13602SShrirang Abhyankar PetscErrorCode main(int argc,char **argv)
55aad13602SShrirang Abhyankar {
56aad13602SShrirang Abhyankar   PetscErrorCode ierr;  /* used to check for functions returning nonzeros */
57aad13602SShrirang Abhyankar   Tao            tao;
58aad13602SShrirang Abhyankar   KSP            ksp;
59aad13602SShrirang Abhyankar   PC             pc;
60aad13602SShrirang Abhyankar   AppCtx         user;  /* application context */
61aad13602SShrirang Abhyankar   Vec            x;
62aad13602SShrirang Abhyankar   PetscMPIInt    rank;
63661095bbSAlp Dener   TaoType        type;
64661095bbSAlp Dener   PetscBool      pdipm;
65aad13602SShrirang Abhyankar 
66aad13602SShrirang Abhyankar   ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
67ffc4695bSBarry Smith   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr);
68aad13602SShrirang Abhyankar   if (rank>1){
69aad13602SShrirang Abhyankar     SETERRQ(PETSC_COMM_SELF,1,"More than 2 processors detected. Example written to use max of 2 processors.\n");
70aad13602SShrirang Abhyankar   }
71aad13602SShrirang Abhyankar 
72aad13602SShrirang Abhyankar   ierr = PetscPrintf(PETSC_COMM_WORLD,"---- Constrained Problem -----\n");CHKERRQ(ierr);
73aad13602SShrirang Abhyankar   ierr = InitializeProblem(&user);CHKERRQ(ierr); /* sets up problem, function below */
74aad13602SShrirang Abhyankar   ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr);
75aad13602SShrirang Abhyankar   ierr = TaoSetType(tao,TAOPDIPM);CHKERRQ(ierr);
76aad13602SShrirang Abhyankar   ierr = TaoSetInitialVector(tao,user.x);CHKERRQ(ierr); /* gets solution vector from problem */
77aad13602SShrirang Abhyankar   ierr = TaoSetVariableBounds(tao,user.xl,user.xu);CHKERRQ(ierr); /* sets lower upper bounds from given solution */
78aad13602SShrirang Abhyankar   ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);CHKERRQ(ierr);
79aad13602SShrirang Abhyankar 
808e58fa1dSresundermann   if (!user.noeqflag){
81aad13602SShrirang Abhyankar     ierr = TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user);CHKERRQ(ierr);
828e58fa1dSresundermann   }
83aad13602SShrirang Abhyankar   ierr = TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user);CHKERRQ(ierr);
848e58fa1dSresundermann   if (!user.noeqflag){
85aad13602SShrirang Abhyankar     ierr = TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user);CHKERRQ(ierr); /* equality jacobian */
868e58fa1dSresundermann   }
87aad13602SShrirang Abhyankar   ierr = TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user);CHKERRQ(ierr); /* inequality jacobian */
88aad13602SShrirang Abhyankar   ierr = TaoSetTolerances(tao,1.e-6,1.e-6,1.e-6);CHKERRQ(ierr);
89aad13602SShrirang Abhyankar   ierr = TaoSetConstraintTolerances(tao,1.e-6,1.e-6);CHKERRQ(ierr);
90aad13602SShrirang Abhyankar 
91aad13602SShrirang Abhyankar   ierr = TaoGetKSP(tao,&ksp);CHKERRQ(ierr);
92aad13602SShrirang Abhyankar   ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
9312d688e0SRylee Sundermann   ierr = PCSetType(pc,PCCHOLESKY);CHKERRQ(ierr);
94aad13602SShrirang Abhyankar   /*
95aad13602SShrirang Abhyankar       This algorithm produces matrices with zeros along the diagonal therefore we use
9612d688e0SRylee Sundermann     MUMPS which provides solver for indefinite matrices
97aad13602SShrirang Abhyankar   */
9812d688e0SRylee Sundermann   ierr = PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);CHKERRQ(ierr);  /* requires mumps to solve pdipm */
99aad13602SShrirang Abhyankar   ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr);
100aad13602SShrirang Abhyankar   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
101aad13602SShrirang Abhyankar 
102aad13602SShrirang Abhyankar   ierr = TaoSetFromOptions(tao);CHKERRQ(ierr);
103661095bbSAlp Dener   ierr = TaoGetType(tao,&type);CHKERRQ(ierr);
104661095bbSAlp Dener   ierr = PetscObjectTypeCompare((PetscObject)tao,TAOPDIPM,&pdipm);CHKERRQ(ierr);
105661095bbSAlp Dener   if (pdipm) {
106661095bbSAlp Dener     ierr = TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);CHKERRQ(ierr);
107661095bbSAlp Dener   }
108aad13602SShrirang Abhyankar 
109aad13602SShrirang Abhyankar   ierr = TaoSolve(tao);CHKERRQ(ierr);
110aad13602SShrirang Abhyankar   ierr = TaoGetSolutionVector(tao,&x);CHKERRQ(ierr);
111aad13602SShrirang Abhyankar   ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
112aad13602SShrirang Abhyankar 
113aad13602SShrirang Abhyankar   /* Free objects */
114aad13602SShrirang Abhyankar   ierr = DestroyProblem(&user);CHKERRQ(ierr);
115aad13602SShrirang Abhyankar   ierr = TaoDestroy(&tao);CHKERRQ(ierr);
116aad13602SShrirang Abhyankar   ierr = PetscFinalize();
117aad13602SShrirang Abhyankar   return ierr;
118aad13602SShrirang Abhyankar }
119aad13602SShrirang Abhyankar 
120aad13602SShrirang Abhyankar PetscErrorCode InitializeProblem(AppCtx *user)
121aad13602SShrirang Abhyankar {
122aad13602SShrirang Abhyankar   PetscErrorCode ierr;
123aad13602SShrirang Abhyankar   PetscMPIInt    size;
124aad13602SShrirang Abhyankar   PetscMPIInt    rank;
125aad13602SShrirang Abhyankar   PetscInt       nloc,neloc,niloc;
126aad13602SShrirang Abhyankar 
127aad13602SShrirang Abhyankar   PetscFunctionBegin;
128ffc4695bSBarry Smith   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr);
129ffc4695bSBarry Smith   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr);
1308e58fa1dSresundermann   user->noeqflag = PETSC_FALSE;
1318e58fa1dSresundermann   ierr = PetscOptionsGetBool(NULL,NULL,"-no_eq",&user->noeqflag,NULL);CHKERRQ(ierr);
132*628da978Sresundermann 
1338e58fa1dSresundermann   if (!user->noeqflag) {
1348e58fa1dSresundermann     ierr = PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n");CHKERRQ(ierr);
1358e58fa1dSresundermann   }
136aad13602SShrirang Abhyankar 
1372d4ee042Sprj-   /* create vector x and set initial values */
138aad13602SShrirang Abhyankar   user->n = 2; /* global length */
139aad13602SShrirang Abhyankar   nloc = (rank==0)?user->n:0;
140661095bbSAlp Dener   ierr = VecCreate(PETSC_COMM_WORLD,&user->x);CHKERRQ(ierr);
141661095bbSAlp Dener   ierr = VecSetSizes(user->x,nloc,user->n);CHKERRQ(ierr);
142661095bbSAlp Dener   ierr = VecSetFromOptions(user->x);CHKERRQ(ierr);
143aad13602SShrirang Abhyankar   ierr = VecSet(user->x,0);CHKERRQ(ierr);
144aad13602SShrirang Abhyankar 
145aad13602SShrirang Abhyankar   /* create and set lower and upper bound vectors */
146aad13602SShrirang Abhyankar   ierr = VecDuplicate(user->x,&user->xl);CHKERRQ(ierr);
147aad13602SShrirang Abhyankar   ierr = VecDuplicate(user->x,&user->xu);CHKERRQ(ierr);
148aad13602SShrirang Abhyankar   ierr = VecSet(user->xl,-1.0);CHKERRQ(ierr);
149aad13602SShrirang Abhyankar   ierr = VecSet(user->xu,2.0);CHKERRQ(ierr);
150aad13602SShrirang Abhyankar 
151aad13602SShrirang Abhyankar   /* create scater to zero */
152aad13602SShrirang Abhyankar   ierr = VecScatterCreateToZero(user->x,&user->scat,&user->Xseq);CHKERRQ(ierr);
153aad13602SShrirang Abhyankar 
154aad13602SShrirang Abhyankar     user->ne = 1;
155*628da978Sresundermann     user->ni = 2;
156aad13602SShrirang Abhyankar     neloc = (rank==0)?user->ne:0;
157*628da978Sresundermann     niloc = (rank==0)?user->ni:0;
158aad13602SShrirang Abhyankar 
1598e58fa1dSresundermann   if (!user->noeqflag){
160661095bbSAlp Dener     ierr = VecCreate(PETSC_COMM_WORLD,&user->ce);CHKERRQ(ierr); /* a 1x1 vec for equality constraints */
161661095bbSAlp Dener     ierr = VecSetSizes(user->ce,neloc,user->ne);CHKERRQ(ierr);
162661095bbSAlp Dener     ierr = VecSetFromOptions(user->ce);CHKERRQ(ierr);
163661095bbSAlp Dener     ierr = VecSetUp(user->ce);CHKERRQ(ierr);
1648e58fa1dSresundermann   }
165*628da978Sresundermann 
166661095bbSAlp Dener   ierr = VecCreate(PETSC_COMM_WORLD,&user->ci);CHKERRQ(ierr); /* a 2x1 vec for inequality constraints */
167661095bbSAlp Dener   ierr = VecSetSizes(user->ci,niloc,user->ni);CHKERRQ(ierr);
168661095bbSAlp Dener   ierr = VecSetFromOptions(user->ci);CHKERRQ(ierr);
169661095bbSAlp Dener   ierr = VecSetUp(user->ci);CHKERRQ(ierr);
170aad13602SShrirang Abhyankar 
171aad13602SShrirang Abhyankar   /* nexn & nixn matricies for equaly and inequalty constriants */
1728e58fa1dSresundermann   if (!user->noeqflag){
173aad13602SShrirang Abhyankar     ierr = MatCreate(PETSC_COMM_WORLD,&user->Ae);CHKERRQ(ierr);
1748e58fa1dSresundermann     ierr = MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n);CHKERRQ(ierr);
1758e58fa1dSresundermann     ierr = MatSetFromOptions(user->Ae);CHKERRQ(ierr);
176661095bbSAlp Dener     ierr = MatSetUp(user->Ae);CHKERRQ(ierr);
1778e58fa1dSresundermann   }
1788e58fa1dSresundermann 
179aad13602SShrirang Abhyankar   ierr = MatCreate(PETSC_COMM_WORLD,&user->Ai);CHKERRQ(ierr);
180aad13602SShrirang Abhyankar   ierr = MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n);CHKERRQ(ierr);
181aad13602SShrirang Abhyankar   ierr = MatSetFromOptions(user->Ai);CHKERRQ(ierr);
182661095bbSAlp Dener   ierr = MatSetUp(user->Ai);CHKERRQ(ierr);
183*628da978Sresundermann 
184*628da978Sresundermann   ierr = MatCreate(PETSC_COMM_WORLD,&user->H);CHKERRQ(ierr);
185*628da978Sresundermann   ierr = MatSetSizes(user->H,nloc,nloc,user->n,user->n);CHKERRQ(ierr);
186*628da978Sresundermann   ierr = MatSetFromOptions(user->H);CHKERRQ(ierr);
187661095bbSAlp Dener   ierr = MatSetUp(user->H);CHKERRQ(ierr);
188aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
189aad13602SShrirang Abhyankar }
190aad13602SShrirang Abhyankar 
191aad13602SShrirang Abhyankar PetscErrorCode DestroyProblem(AppCtx *user)
192aad13602SShrirang Abhyankar {
193aad13602SShrirang Abhyankar   PetscErrorCode ierr;
194aad13602SShrirang Abhyankar 
195aad13602SShrirang Abhyankar   PetscFunctionBegin;
1968e58fa1dSresundermann   if (!user->noeqflag){
197aad13602SShrirang Abhyankar    ierr = MatDestroy(&user->Ae);CHKERRQ(ierr);
1988e58fa1dSresundermann   }
199aad13602SShrirang Abhyankar   ierr = MatDestroy(&user->Ai);CHKERRQ(ierr);
200aad13602SShrirang Abhyankar   ierr = MatDestroy(&user->H);CHKERRQ(ierr);
201aad13602SShrirang Abhyankar 
202aad13602SShrirang Abhyankar   ierr = VecDestroy(&user->x);CHKERRQ(ierr);
2038e58fa1dSresundermann   if (!user->noeqflag){
204aad13602SShrirang Abhyankar     ierr = VecDestroy(&user->ce);CHKERRQ(ierr);
2058e58fa1dSresundermann   }
206aad13602SShrirang Abhyankar   ierr = VecDestroy(&user->ci);CHKERRQ(ierr);
207aad13602SShrirang Abhyankar   ierr = VecDestroy(&user->xl);CHKERRQ(ierr);
208aad13602SShrirang Abhyankar   ierr = VecDestroy(&user->xu);CHKERRQ(ierr);
209aad13602SShrirang Abhyankar   ierr = VecDestroy(&user->Xseq);CHKERRQ(ierr);
210aad13602SShrirang Abhyankar   ierr = VecScatterDestroy(&user->scat);CHKERRQ(ierr);
211aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
212aad13602SShrirang Abhyankar }
213aad13602SShrirang Abhyankar 
214*628da978Sresundermann /* Evaluate
215*628da978Sresundermann    f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
216*628da978Sresundermann    G = grad f = [2*(x0 - 2) - 2;
217*628da978Sresundermann                  2*(x1 - 2) - 2]
218aad13602SShrirang Abhyankar */
219aad13602SShrirang Abhyankar PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
220aad13602SShrirang Abhyankar {
221aad13602SShrirang Abhyankar   PetscScalar       g;
222aad13602SShrirang Abhyankar   const PetscScalar *x;
223aad13602SShrirang Abhyankar   MPI_Comm          comm;
224aad13602SShrirang Abhyankar   PetscMPIInt       rank;
225aad13602SShrirang Abhyankar   PetscErrorCode    ierr;
226aad13602SShrirang Abhyankar   PetscReal         fin;
227aad13602SShrirang Abhyankar   AppCtx            *user=(AppCtx*)ctx;
228aad13602SShrirang Abhyankar   Vec               Xseq=user->Xseq;
229aad13602SShrirang Abhyankar   VecScatter        scat=user->scat;
230aad13602SShrirang Abhyankar 
231aad13602SShrirang Abhyankar   PetscFunctionBegin;
232aad13602SShrirang Abhyankar   ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr);
233ffc4695bSBarry Smith   ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr);
234aad13602SShrirang Abhyankar 
235aad13602SShrirang Abhyankar   ierr = VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
236aad13602SShrirang Abhyankar   ierr = VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
237aad13602SShrirang Abhyankar 
238aad13602SShrirang Abhyankar   fin = 0.0;
239aad13602SShrirang Abhyankar   if (!rank) {
240aad13602SShrirang Abhyankar     ierr = VecGetArrayRead(Xseq,&x);CHKERRQ(ierr);
241aad13602SShrirang Abhyankar     fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]);
242aad13602SShrirang Abhyankar     g = 2.0*(x[0]-2.0) - 2.0;
243aad13602SShrirang Abhyankar     ierr = VecSetValue(G,0,g,INSERT_VALUES);CHKERRQ(ierr);
244aad13602SShrirang Abhyankar     g = 2.0*(x[1]-2.0) - 2.0;
245aad13602SShrirang Abhyankar     ierr = VecSetValue(G,1,g,INSERT_VALUES);CHKERRQ(ierr);
246aad13602SShrirang Abhyankar     ierr = VecRestoreArrayRead(Xseq,&x);CHKERRQ(ierr);
247aad13602SShrirang Abhyankar   }
248ffc4695bSBarry Smith   ierr = MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm);CHKERRMPI(ierr);
249aad13602SShrirang Abhyankar   ierr = VecAssemblyBegin(G);CHKERRQ(ierr);
250aad13602SShrirang Abhyankar   ierr = VecAssemblyEnd(G);CHKERRQ(ierr);
251aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
252aad13602SShrirang Abhyankar }
253aad13602SShrirang Abhyankar 
254*628da978Sresundermann /* Evaluate
255*628da978Sresundermann    H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)]
256*628da978Sresundermann      = [ 2*(1+de[0]-di[0]+di[1]), 0;
257*628da978Sresundermann                    0,             2]
258*628da978Sresundermann */
259aad13602SShrirang Abhyankar PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx)
260aad13602SShrirang Abhyankar {
261aad13602SShrirang Abhyankar   AppCtx            *user=(AppCtx*)ctx;
262aad13602SShrirang Abhyankar   Vec               DE,DI;
263aad13602SShrirang Abhyankar   const PetscScalar *de,*di;
264aad13602SShrirang Abhyankar   PetscInt          zero=0,one=1;
265aad13602SShrirang Abhyankar   PetscScalar       two=2.0;
266aad13602SShrirang Abhyankar   PetscScalar       val=0.0;
267aad13602SShrirang Abhyankar   Vec               Deseq,Diseq=user->Xseq;
268aad13602SShrirang Abhyankar   VecScatter        Descat,Discat=user->scat;
269aad13602SShrirang Abhyankar   PetscMPIInt       rank;
270aad13602SShrirang Abhyankar   MPI_Comm          comm;
271aad13602SShrirang Abhyankar   PetscErrorCode    ierr;
272aad13602SShrirang Abhyankar 
273aad13602SShrirang Abhyankar   PetscFunctionBegin;
274aad13602SShrirang Abhyankar   ierr = TaoGetDualVariables(tao,&DE,&DI);CHKERRQ(ierr);
275aad13602SShrirang Abhyankar 
276aad13602SShrirang Abhyankar   ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr);
277ffc4695bSBarry Smith   ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr);
278aad13602SShrirang Abhyankar 
2798e58fa1dSresundermann   if (!user->noeqflag){
280aad13602SShrirang Abhyankar    ierr = VecScatterCreateToZero(DE,&Descat,&Deseq);CHKERRQ(ierr);
281aad13602SShrirang Abhyankar    ierr = VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
282aad13602SShrirang Abhyankar    ierr = VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
2838e58fa1dSresundermann   }
2848e58fa1dSresundermann   ierr = VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
285aad13602SShrirang Abhyankar   ierr = VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
286aad13602SShrirang Abhyankar 
287aad13602SShrirang Abhyankar   if (!rank){
2888e58fa1dSresundermann     if (!user->noeqflag){
289aad13602SShrirang Abhyankar       ierr = VecGetArrayRead(Deseq,&de);CHKERRQ(ierr);  /* places equality constraint dual into array */
2908e58fa1dSresundermann     }
291aad13602SShrirang Abhyankar     ierr = VecGetArrayRead(Diseq,&di);CHKERRQ(ierr);  /* places inequality constraint dual into array */
292*628da978Sresundermann 
2938e58fa1dSresundermann     if (!user->noeqflag){
294*628da978Sresundermann       val = 2.0 * (1 + de[0] - di[0] + di[1]);
295aad13602SShrirang Abhyankar       ierr = VecRestoreArrayRead(Deseq,&de);CHKERRQ(ierr);
296*628da978Sresundermann       ierr = VecRestoreArrayRead(Diseq,&di);CHKERRQ(ierr);
2978e58fa1dSresundermann     } else {
298*628da978Sresundermann       val = 2.0 * (1 - di[0] + di[1]);
2998e58fa1dSresundermann     }
300aad13602SShrirang Abhyankar     ierr = VecRestoreArrayRead(Diseq,&di);CHKERRQ(ierr);
301aad13602SShrirang Abhyankar     ierr = MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES);CHKERRQ(ierr);
302aad13602SShrirang Abhyankar     ierr = MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES);CHKERRQ(ierr);
303aad13602SShrirang Abhyankar   }
304aad13602SShrirang Abhyankar   ierr = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
305aad13602SShrirang Abhyankar   ierr = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
3068e58fa1dSresundermann   if (!user->noeqflag){
307aad13602SShrirang Abhyankar     ierr = VecScatterDestroy(&Descat);CHKERRQ(ierr);
308aad13602SShrirang Abhyankar     ierr = VecDestroy(&Deseq);CHKERRQ(ierr);
3098e58fa1dSresundermann   }
310aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
311aad13602SShrirang Abhyankar }
312aad13602SShrirang Abhyankar 
313*628da978Sresundermann /* Evaluate
314*628da978Sresundermann    h = [ x0^2 - x1;
315*628da978Sresundermann          1 -(x0^2 - x1)]
316*628da978Sresundermann */
317aad13602SShrirang Abhyankar PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx)
318aad13602SShrirang Abhyankar {
319aad13602SShrirang Abhyankar   const PetscScalar *x;
320aad13602SShrirang Abhyankar   PetscScalar       ci;
321aad13602SShrirang Abhyankar   PetscErrorCode    ierr;
322aad13602SShrirang Abhyankar   MPI_Comm          comm;
323aad13602SShrirang Abhyankar   PetscMPIInt       rank;
324aad13602SShrirang Abhyankar   AppCtx            *user=(AppCtx*)ctx;
325aad13602SShrirang Abhyankar   Vec               Xseq=user->Xseq;
326aad13602SShrirang Abhyankar   VecScatter        scat=user->scat;
327aad13602SShrirang Abhyankar 
328aad13602SShrirang Abhyankar   PetscFunctionBegin;
329aad13602SShrirang Abhyankar   ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr);
330ffc4695bSBarry Smith   ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr);
331aad13602SShrirang Abhyankar 
332aad13602SShrirang Abhyankar   ierr = VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
333aad13602SShrirang Abhyankar   ierr = VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
334aad13602SShrirang Abhyankar 
335aad13602SShrirang Abhyankar   if (!rank) {
336aad13602SShrirang Abhyankar     ierr = VecGetArrayRead(Xseq,&x);CHKERRQ(ierr);
337aad13602SShrirang Abhyankar     ci = x[0]*x[0] - x[1];
338aad13602SShrirang Abhyankar     ierr = VecSetValue(CI,0,ci,INSERT_VALUES);CHKERRQ(ierr);
339aad13602SShrirang Abhyankar     ci = -x[0]*x[0] + x[1] + 1.0;
340aad13602SShrirang Abhyankar     ierr = VecSetValue(CI,1,ci,INSERT_VALUES);CHKERRQ(ierr);
341aad13602SShrirang Abhyankar     ierr = VecRestoreArrayRead(Xseq,&x);CHKERRQ(ierr);
342aad13602SShrirang Abhyankar   }
343aad13602SShrirang Abhyankar   ierr = VecAssemblyBegin(CI);CHKERRQ(ierr);
344aad13602SShrirang Abhyankar   ierr = VecAssemblyEnd(CI);CHKERRQ(ierr);
345aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
346aad13602SShrirang Abhyankar }
347aad13602SShrirang Abhyankar 
348*628da978Sresundermann /* Evaluate
349*628da978Sresundermann    g = [ x0^2 + x1 - 2]
350*628da978Sresundermann */
351aad13602SShrirang Abhyankar PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx)
352aad13602SShrirang Abhyankar {
353aad13602SShrirang Abhyankar   const PetscScalar *x;
354aad13602SShrirang Abhyankar   PetscScalar       ce;
355aad13602SShrirang Abhyankar   PetscErrorCode    ierr;
356aad13602SShrirang Abhyankar   MPI_Comm          comm;
357aad13602SShrirang Abhyankar   PetscMPIInt       rank;
358aad13602SShrirang Abhyankar   AppCtx            *user=(AppCtx*)ctx;
359aad13602SShrirang Abhyankar   Vec               Xseq=user->Xseq;
360aad13602SShrirang Abhyankar   VecScatter        scat=user->scat;
361aad13602SShrirang Abhyankar 
362aad13602SShrirang Abhyankar   PetscFunctionBegin;
363aad13602SShrirang Abhyankar   ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr);
364ffc4695bSBarry Smith   ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr);
365aad13602SShrirang Abhyankar 
366aad13602SShrirang Abhyankar   ierr = VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
367aad13602SShrirang Abhyankar   ierr = VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
368aad13602SShrirang Abhyankar 
369aad13602SShrirang Abhyankar   if (!rank) {
370aad13602SShrirang Abhyankar     ierr = VecGetArrayRead(Xseq,&x);CHKERRQ(ierr);
371aad13602SShrirang Abhyankar     ce = x[0]*x[0] + x[1] - 2.0;
372aad13602SShrirang Abhyankar     ierr = VecSetValue(CE,0,ce,INSERT_VALUES);CHKERRQ(ierr);
373aad13602SShrirang Abhyankar     ierr = VecRestoreArrayRead(Xseq,&x);CHKERRQ(ierr);
374aad13602SShrirang Abhyankar   }
375aad13602SShrirang Abhyankar   ierr = VecAssemblyBegin(CE);CHKERRQ(ierr);
376aad13602SShrirang Abhyankar   ierr = VecAssemblyEnd(CE);CHKERRQ(ierr);
377aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
378aad13602SShrirang Abhyankar }
379aad13602SShrirang Abhyankar 
380*628da978Sresundermann /*
381*628da978Sresundermann   grad h = [  2*x0, -1;
382*628da978Sresundermann              -2*x0,  1]
383*628da978Sresundermann */
384aad13602SShrirang Abhyankar PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre,  void *ctx)
385aad13602SShrirang Abhyankar {
386aad13602SShrirang Abhyankar   AppCtx            *user=(AppCtx*)ctx;
387aad13602SShrirang Abhyankar   PetscInt          cols[2],min,max,i;
388aad13602SShrirang Abhyankar   PetscScalar       vals[2];
389aad13602SShrirang Abhyankar   const PetscScalar *x;
390aad13602SShrirang Abhyankar   PetscErrorCode    ierr;
391aad13602SShrirang Abhyankar   Vec               Xseq=user->Xseq;
392aad13602SShrirang Abhyankar   VecScatter        scat=user->scat;
393aad13602SShrirang Abhyankar   MPI_Comm          comm;
394aad13602SShrirang Abhyankar   PetscMPIInt       rank;
395aad13602SShrirang Abhyankar 
396aad13602SShrirang Abhyankar   PetscFunctionBegin;
397aad13602SShrirang Abhyankar   ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr);
398ffc4695bSBarry Smith   ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr);
399aad13602SShrirang Abhyankar   ierr = VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
400aad13602SShrirang Abhyankar   ierr = VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
401aad13602SShrirang Abhyankar 
402aad13602SShrirang Abhyankar   ierr = VecGetArrayRead(Xseq,&x);CHKERRQ(ierr);
403aad13602SShrirang Abhyankar   ierr = MatGetOwnershipRange(JI,&min,&max);CHKERRQ(ierr);
404aad13602SShrirang Abhyankar 
405aad13602SShrirang Abhyankar   cols[0] = 0; cols[1] = 1;
406aad13602SShrirang Abhyankar   for (i=min;i<max;i++) {
407aad13602SShrirang Abhyankar     if (i==0){
408*628da978Sresundermann       vals[0] = 2*x[0]; vals[1] = -1.0;
409aad13602SShrirang Abhyankar       ierr = MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);CHKERRQ(ierr);
410aad13602SShrirang Abhyankar     }
411aad13602SShrirang Abhyankar     if (i==1) {
412*628da978Sresundermann       vals[0] = -2*x[0]; vals[1] = 1.0;
413aad13602SShrirang Abhyankar       ierr = MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);CHKERRQ(ierr);
414aad13602SShrirang Abhyankar     }
415aad13602SShrirang Abhyankar   }
416aad13602SShrirang Abhyankar   ierr = VecRestoreArrayRead(Xseq,&x);CHKERRQ(ierr);
417aad13602SShrirang Abhyankar 
418aad13602SShrirang Abhyankar   ierr = MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
419aad13602SShrirang Abhyankar   ierr = MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
420aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
421aad13602SShrirang Abhyankar }
422aad13602SShrirang Abhyankar 
423*628da978Sresundermann /*
424*628da978Sresundermann   grad g = [2*x0
425*628da978Sresundermann              1.0 ]
426*628da978Sresundermann */
427aad13602SShrirang Abhyankar PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx)
428aad13602SShrirang Abhyankar {
429aad13602SShrirang Abhyankar   PetscInt          rows[2];
430aad13602SShrirang Abhyankar   PetscScalar       vals[2];
431aad13602SShrirang Abhyankar   const PetscScalar *x;
432aad13602SShrirang Abhyankar   PetscMPIInt       rank;
433aad13602SShrirang Abhyankar   MPI_Comm          comm;
434aad13602SShrirang Abhyankar   PetscErrorCode    ierr;
435aad13602SShrirang Abhyankar 
436aad13602SShrirang Abhyankar   PetscFunctionBegin;
437aad13602SShrirang Abhyankar   ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr);
438ffc4695bSBarry Smith   ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr);
439aad13602SShrirang Abhyankar 
440aad13602SShrirang Abhyankar   if (!rank) {
441aad13602SShrirang Abhyankar     ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
442aad13602SShrirang Abhyankar     rows[0] = 0;       rows[1] = 1;
443aad13602SShrirang Abhyankar     vals[0] = 2*x[0];  vals[1] = 1.0;
444aad13602SShrirang Abhyankar     ierr = MatSetValues(JE,1,rows,2,rows,vals,INSERT_VALUES);CHKERRQ(ierr);
445aad13602SShrirang Abhyankar     ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
446aad13602SShrirang Abhyankar   }
447aad13602SShrirang Abhyankar   ierr = MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
448aad13602SShrirang Abhyankar   ierr = MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
449aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
450aad13602SShrirang Abhyankar }
451aad13602SShrirang Abhyankar 
452aad13602SShrirang Abhyankar 
453aad13602SShrirang Abhyankar /*TEST
454aad13602SShrirang Abhyankar 
455aad13602SShrirang Abhyankar    build:
45612d688e0SRylee Sundermann       requires: !complex !define(PETSC_USE_CXX) mumps
457aad13602SShrirang Abhyankar 
458aad13602SShrirang Abhyankar    test:
45909ee8bb0SRylee Sundermann       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
460aad13602SShrirang Abhyankar 
461aad13602SShrirang Abhyankar    test:
462aad13602SShrirang Abhyankar       suffix: 2
463aad13602SShrirang Abhyankar       nsize: 2
46409ee8bb0SRylee Sundermann       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
465aad13602SShrirang Abhyankar 
4668e58fa1dSresundermann    test:
4678e58fa1dSresundermann       suffix: 3
4688e58fa1dSresundermann       args: -tao_converged_reason -no_eq
4698e58fa1dSresundermann 
4708e58fa1dSresundermann    test:
4718e58fa1dSresundermann       suffix: 4
4728e58fa1dSresundermann       nsize: 2
4738e58fa1dSresundermann       args: -tao_converged_reason -no_eq
4748e58fa1dSresundermann 
475661095bbSAlp Dener    test:
476661095bbSAlp Dener       suffix: 5
477661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm
478661095bbSAlp Dener 
479661095bbSAlp Dener    test:
480661095bbSAlp Dener       suffix: 6
481661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm -tao_almm_type phr
482661095bbSAlp Dener 
483661095bbSAlp Dener    test:
484661095bbSAlp Dener       suffix: 7
485661095bbSAlp Dener       nsize: 2
486661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm
487661095bbSAlp Dener 
488661095bbSAlp Dener    test:
489661095bbSAlp Dener       suffix: 8
490661095bbSAlp Dener       nsize: 2
491661095bbSAlp Dener       requires: cuda
492661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm -vec_type cuda -mat_type aijcusparse
493661095bbSAlp Dener 
494661095bbSAlp Dener    test:
495661095bbSAlp Dener       suffix: 9
496661095bbSAlp Dener       nsize: 2
497661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm -no_eq
498661095bbSAlp Dener 
499661095bbSAlp Dener    test:
500661095bbSAlp Dener       suffix: 10
501661095bbSAlp Dener       nsize: 2
502661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm -tao_almm_type phr -no_eq
503661095bbSAlp Dener 
504aad13602SShrirang Abhyankar TEST*/
505