xref: /petsc/src/tao/constrained/tutorials/ex1.c (revision 5f80ce2ab25dff0f4601e710601cbbcecf323266)
1aad13602SShrirang Abhyankar /* Program usage: mpiexec -n 2 ./ex1 [-help] [all TAO options] */
2aad13602SShrirang Abhyankar 
3aad13602SShrirang Abhyankar /* ----------------------------------------------------------------------
4628da978Sresundermann 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
8628da978Sresundermann -->
9628da978Sresundermann       g(x)  = 0
10628da978Sresundermann       h(x) >= 0
11628da978Sresundermann       -1 <= x0, x1 <= 2
12628da978Sresundermann where
13628da978Sresundermann       g(x) = x0^2 + x1 - 2
14628da978Sresundermann       h(x) = [x0^2 - x1
15628da978Sresundermann               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\
24f560b561SHong Zhang   -init_view         : view initial object setup\n\
25aad13602SShrirang Abhyankar   -snes_fd           : snes with finite difference Jacobian (needed for pdipm)\n\
26628da978Sresundermann   -snes_compare_explicit : compare user Jacobian with finite difference Jacobian \n\
27aad13602SShrirang Abhyankar   -tao_cmonitor      : convergence monitor with constraint norm \n\
28a5b23f4aSJose E. Roman   -tao_view_solution : view exact solution at each iteration\n\
29f560b561SHong Zhang   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";
30aad13602SShrirang Abhyankar 
31aad13602SShrirang Abhyankar /*
32628da978Sresundermann    User-defined application context - contains data needed by the application
33aad13602SShrirang Abhyankar */
34aad13602SShrirang Abhyankar typedef struct {
35aad13602SShrirang Abhyankar   PetscInt   n;  /* Global length of x */
36aad13602SShrirang Abhyankar   PetscInt   ne; /* Global number of equality constraints */
37aad13602SShrirang Abhyankar   PetscInt   ni; /* Global number of inequality constraints */
38f560b561SHong Zhang   PetscBool  noeqflag, initview;
39aad13602SShrirang Abhyankar   Vec        x,xl,xu;
40aad13602SShrirang Abhyankar   Vec        ce,ci,bl,bu,Xseq;
41aad13602SShrirang Abhyankar   Mat        Ae,Ai,H;
42aad13602SShrirang Abhyankar   VecScatter scat;
43aad13602SShrirang Abhyankar } AppCtx;
44aad13602SShrirang Abhyankar 
45aad13602SShrirang Abhyankar /* -------- User-defined Routines --------- */
46aad13602SShrirang Abhyankar PetscErrorCode InitializeProblem(AppCtx *);
47aad13602SShrirang Abhyankar PetscErrorCode DestroyProblem(AppCtx *);
48aad13602SShrirang Abhyankar PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
49aad13602SShrirang Abhyankar PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
50aad13602SShrirang Abhyankar PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*);
51aad13602SShrirang Abhyankar PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*);
52aad13602SShrirang Abhyankar PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*);
53aad13602SShrirang Abhyankar PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*);
54aad13602SShrirang Abhyankar 
55aad13602SShrirang Abhyankar PetscErrorCode main(int argc,char **argv)
56aad13602SShrirang Abhyankar {
57aad13602SShrirang Abhyankar   PetscErrorCode ierr;  /* used to check for functions returning nonzeros */
58aad13602SShrirang Abhyankar   Tao            tao;
59aad13602SShrirang Abhyankar   KSP            ksp;
60aad13602SShrirang Abhyankar   PC             pc;
61aad13602SShrirang Abhyankar   AppCtx         user;  /* application context */
62f560b561SHong Zhang   Vec            x,G,CI,CE;
63f560b561SHong Zhang   PetscMPIInt    size;
64661095bbSAlp Dener   TaoType        type;
65f560b561SHong Zhang   PetscReal      f;
66661095bbSAlp Dener   PetscBool      pdipm;
67aad13602SShrirang Abhyankar 
68aad13602SShrirang Abhyankar   ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
69*5f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
703c859ba3SBarry Smith   PetscCheck(size <= 2,PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,"More than 2 processors detected. Example written to use max of 2 processors.");
71aad13602SShrirang Abhyankar 
72*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"---- Constrained Problem -----\n"));
73*5f80ce2aSJacob Faibussowitsch   CHKERRQ(InitializeProblem(&user)); /* sets up problem, function below */
74*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoCreate(PETSC_COMM_WORLD,&tao));
75*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetType(tao,TAOPDIPM));
76*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetSolution(tao,user.x)); /* gets solution vector from problem */
77*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetVariableBounds(tao,user.xl,user.xu)); /* sets lower upper bounds from given solution */
78*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void*)&user));
79aad13602SShrirang Abhyankar 
808e58fa1dSresundermann   if (!user.noeqflag) {
81*5f80ce2aSJacob Faibussowitsch     CHKERRQ(TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user));
828e58fa1dSresundermann   }
83*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user));
848e58fa1dSresundermann   if (!user.noeqflag) {
85*5f80ce2aSJacob Faibussowitsch     CHKERRQ(TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user)); /* equality jacobian */
868e58fa1dSresundermann   }
87*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user)); /* inequality jacobian */
88*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetTolerances(tao,1.e-6,1.e-6,1.e-6));
89*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetConstraintTolerances(tao,1.e-6,1.e-6));
90aad13602SShrirang Abhyankar 
91*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoGetKSP(tao,&ksp));
92*5f80ce2aSJacob Faibussowitsch   CHKERRQ(KSPGetPC(ksp,&pc));
93*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PCSetType(pc,PCCHOLESKY));
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   */
98f560b561SHong Zhang #if defined(PETSC_HAVE_MUMPS)
99*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PCFactorSetMatSolverType(pc,MATSOLVERMUMPS));  /* requires mumps to solve pdipm */
100f560b561SHong Zhang #else
1013c859ba3SBarry Smith   PetscCheck(size == 1,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Requires an external package that supports parallel PCCHOLESKY, e.g., MUMPS.");
102f560b561SHong Zhang #endif
103*5f80ce2aSJacob Faibussowitsch   CHKERRQ(KSPSetType(ksp,KSPPREONLY));
104*5f80ce2aSJacob Faibussowitsch   CHKERRQ(KSPSetFromOptions(ksp));
105aad13602SShrirang Abhyankar 
106*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetFromOptions(tao));
107*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoGetType(tao,&type));
108*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectTypeCompare((PetscObject)tao,TAOPDIPM,&pdipm));
109661095bbSAlp Dener   if (pdipm) {
110*5f80ce2aSJacob Faibussowitsch     CHKERRQ(TaoSetHessian(tao,user.H,user.H,FormHessian,(void*)&user));
111f560b561SHong Zhang     if (user.initview) {
112*5f80ce2aSJacob Faibussowitsch       CHKERRQ(TaoSetUp(tao));
113*5f80ce2aSJacob Faibussowitsch       CHKERRQ(VecDuplicate(user.x, &G));
114*5f80ce2aSJacob Faibussowitsch       CHKERRQ(FormFunctionGradient(tao, user.x, &f, G, (void*)&user));
115*5f80ce2aSJacob Faibussowitsch       CHKERRQ(FormHessian(tao, user.x, user.H, user.H, (void*)&user));
116*5f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD));
117*5f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nInitial point X:\n",f));
118*5f80ce2aSJacob Faibussowitsch       CHKERRQ(VecView(user.x, PETSC_VIEWER_STDOUT_WORLD));
119*5f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nInitial objective f(x) = %g\n",f));
120*5f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nInitial gradient and Hessian:\n",f));
121*5f80ce2aSJacob Faibussowitsch       CHKERRQ(VecView(G, PETSC_VIEWER_STDOUT_WORLD));
122*5f80ce2aSJacob Faibussowitsch       CHKERRQ(MatView(user.H, PETSC_VIEWER_STDOUT_WORLD));
123*5f80ce2aSJacob Faibussowitsch       CHKERRQ(VecDestroy(&G));
124*5f80ce2aSJacob Faibussowitsch       CHKERRQ(FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void*)&user));
125*5f80ce2aSJacob Faibussowitsch       CHKERRQ(MatCreateVecs(user.Ai, NULL, &CI));
126*5f80ce2aSJacob Faibussowitsch       CHKERRQ(FormInequalityConstraints(tao, user.x, CI, (void*)&user));
127*5f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nInitial inequality constraints and Jacobian:\n",f));
128*5f80ce2aSJacob Faibussowitsch       CHKERRQ(VecView(CI, PETSC_VIEWER_STDOUT_WORLD));
129*5f80ce2aSJacob Faibussowitsch       CHKERRQ(MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD));
130*5f80ce2aSJacob Faibussowitsch       CHKERRQ(VecDestroy(&CI));
131f560b561SHong Zhang       if (!user.noeqflag) {
132*5f80ce2aSJacob Faibussowitsch         CHKERRQ(FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void*)&user));
133*5f80ce2aSJacob Faibussowitsch         CHKERRQ(MatCreateVecs(user.Ae, NULL, &CE));
134*5f80ce2aSJacob Faibussowitsch         CHKERRQ(FormEqualityConstraints(tao, user.x, CE, (void*)&user));
135*5f80ce2aSJacob Faibussowitsch         CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nInitial equality constraints and Jacobian:\n",f));
136*5f80ce2aSJacob Faibussowitsch         CHKERRQ(VecView(CE, PETSC_VIEWER_STDOUT_WORLD));
137*5f80ce2aSJacob Faibussowitsch         CHKERRQ(MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD));
138*5f80ce2aSJacob Faibussowitsch         CHKERRQ(VecDestroy(&CE));
139f560b561SHong Zhang       }
140*5f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "\n"));
141*5f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD));
142f560b561SHong Zhang     }
143661095bbSAlp Dener   }
144aad13602SShrirang Abhyankar 
145*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSolve(tao));
146*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoGetSolution(tao,&x));
147*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecView(x,PETSC_VIEWER_STDOUT_WORLD));
148aad13602SShrirang Abhyankar 
149aad13602SShrirang Abhyankar   /* Free objects */
150*5f80ce2aSJacob Faibussowitsch   CHKERRQ(DestroyProblem(&user));
151*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoDestroy(&tao));
152aad13602SShrirang Abhyankar   ierr = PetscFinalize();
153aad13602SShrirang Abhyankar   return ierr;
154aad13602SShrirang Abhyankar }
155aad13602SShrirang Abhyankar 
156aad13602SShrirang Abhyankar PetscErrorCode InitializeProblem(AppCtx *user)
157aad13602SShrirang Abhyankar {
158aad13602SShrirang Abhyankar   PetscMPIInt    size;
159aad13602SShrirang Abhyankar   PetscMPIInt    rank;
160aad13602SShrirang Abhyankar   PetscInt       nloc,neloc,niloc;
161aad13602SShrirang Abhyankar 
162aad13602SShrirang Abhyankar   PetscFunctionBegin;
163*5f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
164*5f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
1658e58fa1dSresundermann   user->noeqflag = PETSC_FALSE;
166*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-no_eq",&user->noeqflag,NULL));
167f560b561SHong Zhang   user->initview = PETSC_FALSE;
168*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-init_view",&user->initview,NULL));
169628da978Sresundermann 
1708e58fa1dSresundermann   if (!user->noeqflag) {
171f560b561SHong Zhang     /* Tell user the correct solution, not an error checking */
172*5f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n"));
1738e58fa1dSresundermann   }
174aad13602SShrirang Abhyankar 
1752d4ee042Sprj-   /* create vector x and set initial values */
176aad13602SShrirang Abhyankar   user->n = 2; /* global length */
177f560b561SHong Zhang   nloc = (size==1)?user->n:1;
178*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user->x));
179*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSetSizes(user->x,nloc,user->n));
180*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSetFromOptions(user->x));
181*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSet(user->x,0));
182aad13602SShrirang Abhyankar 
183aad13602SShrirang Abhyankar   /* create and set lower and upper bound vectors */
184*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDuplicate(user->x,&user->xl));
185*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDuplicate(user->x,&user->xu));
186*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSet(user->xl,-1.0));
187*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSet(user->xu,2.0));
188aad13602SShrirang Abhyankar 
189aad13602SShrirang Abhyankar   /* create scater to zero */
190*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterCreateToZero(user->x,&user->scat,&user->Xseq));
191aad13602SShrirang Abhyankar 
192aad13602SShrirang Abhyankar   user->ne = 1;
193628da978Sresundermann   user->ni = 2;
194aad13602SShrirang Abhyankar   neloc = (rank==0)?user->ne:0;
195f560b561SHong Zhang   niloc = (size==1)?user->ni:1;
196aad13602SShrirang Abhyankar 
1978e58fa1dSresundermann   if (!user->noeqflag) {
198*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user->ce)); /* a 1x1 vec for equality constraints */
199*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetSizes(user->ce,neloc,user->ne));
200*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetFromOptions(user->ce));
201*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetUp(user->ce));
2028e58fa1dSresundermann   }
203628da978Sresundermann 
204*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user->ci)); /* a 2x1 vec for inequality constraints */
205*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSetSizes(user->ci,niloc,user->ni));
206*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSetFromOptions(user->ci));
207*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSetUp(user->ci));
208aad13602SShrirang Abhyankar 
209a5b23f4aSJose E. Roman   /* nexn & nixn matricies for equally and inequalty constraints */
2108e58fa1dSresundermann   if (!user->noeqflag) {
211*5f80ce2aSJacob Faibussowitsch     CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user->Ae));
212*5f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n));
213*5f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetFromOptions(user->Ae));
214*5f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetUp(user->Ae));
2158e58fa1dSresundermann   }
2168e58fa1dSresundermann 
217*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user->Ai));
218*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n));
219*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetFromOptions(user->Ai));
220*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetUp(user->Ai));
221628da978Sresundermann 
222*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user->H));
223*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetSizes(user->H,nloc,nloc,user->n,user->n));
224*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetFromOptions(user->H));
225*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetUp(user->H));
226aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
227aad13602SShrirang Abhyankar }
228aad13602SShrirang Abhyankar 
229aad13602SShrirang Abhyankar PetscErrorCode DestroyProblem(AppCtx *user)
230aad13602SShrirang Abhyankar {
231aad13602SShrirang Abhyankar   PetscFunctionBegin;
2328e58fa1dSresundermann   if (!user->noeqflag) {
233*5f80ce2aSJacob Faibussowitsch     CHKERRQ(MatDestroy(&user->Ae));
2348e58fa1dSresundermann   }
235*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatDestroy(&user->Ai));
236*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatDestroy(&user->H));
237aad13602SShrirang Abhyankar 
238*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&user->x));
2398e58fa1dSresundermann   if (!user->noeqflag) {
240*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecDestroy(&user->ce));
2418e58fa1dSresundermann   }
242*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&user->ci));
243*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&user->xl));
244*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&user->xu));
245*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&user->Xseq));
246*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterDestroy(&user->scat));
247aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
248aad13602SShrirang Abhyankar }
249aad13602SShrirang Abhyankar 
250628da978Sresundermann /* Evaluate
251628da978Sresundermann    f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
252628da978Sresundermann    G = grad f = [2*(x0 - 2) - 2;
253628da978Sresundermann                  2*(x1 - 2) - 2]
254aad13602SShrirang Abhyankar */
255aad13602SShrirang Abhyankar PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
256aad13602SShrirang Abhyankar {
257aad13602SShrirang Abhyankar   PetscScalar       g;
258aad13602SShrirang Abhyankar   const PetscScalar *x;
259aad13602SShrirang Abhyankar   MPI_Comm          comm;
260aad13602SShrirang Abhyankar   PetscMPIInt       rank;
261aad13602SShrirang Abhyankar   PetscReal         fin;
262aad13602SShrirang Abhyankar   AppCtx            *user=(AppCtx*)ctx;
263aad13602SShrirang Abhyankar   Vec               Xseq=user->Xseq;
264aad13602SShrirang Abhyankar   VecScatter        scat=user->scat;
265aad13602SShrirang Abhyankar 
266aad13602SShrirang Abhyankar   PetscFunctionBegin;
267*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm));
268*5f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm,&rank));
269aad13602SShrirang Abhyankar 
270*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
271*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
272aad13602SShrirang Abhyankar 
273aad13602SShrirang Abhyankar   fin = 0.0;
274dd400576SPatrick Sanan   if (rank == 0) {
275*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecGetArrayRead(Xseq,&x));
276aad13602SShrirang Abhyankar     fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]);
277aad13602SShrirang Abhyankar     g = 2.0*(x[0]-2.0) - 2.0;
278*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetValue(G,0,g,INSERT_VALUES));
279aad13602SShrirang Abhyankar     g = 2.0*(x[1]-2.0) - 2.0;
280*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetValue(G,1,g,INSERT_VALUES));
281*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecRestoreArrayRead(Xseq,&x));
282aad13602SShrirang Abhyankar   }
283*5f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm));
284*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyBegin(G));
285*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyEnd(G));
286aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
287aad13602SShrirang Abhyankar }
288aad13602SShrirang Abhyankar 
289628da978Sresundermann /* Evaluate
290628da978Sresundermann    H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)]
291628da978Sresundermann      = [ 2*(1+de[0]-di[0]+di[1]), 0;
292628da978Sresundermann                    0,             2]
293628da978Sresundermann */
294aad13602SShrirang Abhyankar PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx)
295aad13602SShrirang Abhyankar {
296aad13602SShrirang Abhyankar   AppCtx            *user=(AppCtx*)ctx;
297aad13602SShrirang Abhyankar   Vec               DE,DI;
298aad13602SShrirang Abhyankar   const PetscScalar *de,*di;
299aad13602SShrirang Abhyankar   PetscInt          zero=0,one=1;
300aad13602SShrirang Abhyankar   PetscScalar       two=2.0;
301aad13602SShrirang Abhyankar   PetscScalar       val=0.0;
302f560b561SHong Zhang   Vec               Deseq,Diseq;
303f560b561SHong Zhang   VecScatter        Descat,Discat;
304aad13602SShrirang Abhyankar   PetscMPIInt       rank;
305aad13602SShrirang Abhyankar   MPI_Comm          comm;
306aad13602SShrirang Abhyankar 
307aad13602SShrirang Abhyankar   PetscFunctionBegin;
308*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoGetDualVariables(tao,&DE,&DI));
309aad13602SShrirang Abhyankar 
310*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm));
311*5f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm,&rank));
312aad13602SShrirang Abhyankar 
3138e58fa1dSresundermann   if (!user->noeqflag) {
314*5f80ce2aSJacob Faibussowitsch    CHKERRQ(VecScatterCreateToZero(DE,&Descat,&Deseq));
315*5f80ce2aSJacob Faibussowitsch    CHKERRQ(VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD));
316*5f80ce2aSJacob Faibussowitsch    CHKERRQ(VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD));
3178e58fa1dSresundermann   }
318*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterCreateToZero(DI,&Discat,&Diseq));
319*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD));
320*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD));
321aad13602SShrirang Abhyankar 
322dd400576SPatrick Sanan   if (rank == 0) {
3238e58fa1dSresundermann     if (!user->noeqflag) {
324*5f80ce2aSJacob Faibussowitsch       CHKERRQ(VecGetArrayRead(Deseq,&de));  /* places equality constraint dual into array */
3258e58fa1dSresundermann     }
326*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecGetArrayRead(Diseq,&di));  /* places inequality constraint dual into array */
327628da978Sresundermann 
3288e58fa1dSresundermann     if (!user->noeqflag) {
329628da978Sresundermann       val = 2.0 * (1 + de[0] - di[0] + di[1]);
330*5f80ce2aSJacob Faibussowitsch       CHKERRQ(VecRestoreArrayRead(Deseq,&de));
331*5f80ce2aSJacob Faibussowitsch       CHKERRQ(VecRestoreArrayRead(Diseq,&di));
3328e58fa1dSresundermann     } else {
333628da978Sresundermann       val = 2.0 * (1 - di[0] + di[1]);
3348e58fa1dSresundermann     }
335*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecRestoreArrayRead(Diseq,&di));
336*5f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES));
337*5f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES));
338aad13602SShrirang Abhyankar   }
339*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY));
340*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY));
3418e58fa1dSresundermann   if (!user->noeqflag) {
342*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecScatterDestroy(&Descat));
343*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecDestroy(&Deseq));
3448e58fa1dSresundermann   }
345*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterDestroy(&Discat));
346*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&Diseq));
347aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
348aad13602SShrirang Abhyankar }
349aad13602SShrirang Abhyankar 
350628da978Sresundermann /* Evaluate
351628da978Sresundermann    h = [ x0^2 - x1;
352628da978Sresundermann          1 -(x0^2 - x1)]
353628da978Sresundermann */
354aad13602SShrirang Abhyankar PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx)
355aad13602SShrirang Abhyankar {
356aad13602SShrirang Abhyankar   const PetscScalar *x;
357aad13602SShrirang Abhyankar   PetscScalar       ci;
358aad13602SShrirang Abhyankar   MPI_Comm          comm;
359aad13602SShrirang Abhyankar   PetscMPIInt       rank;
360aad13602SShrirang Abhyankar   AppCtx            *user=(AppCtx*)ctx;
361aad13602SShrirang Abhyankar   Vec               Xseq=user->Xseq;
362aad13602SShrirang Abhyankar   VecScatter        scat=user->scat;
363aad13602SShrirang Abhyankar 
364aad13602SShrirang Abhyankar   PetscFunctionBegin;
365*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm));
366*5f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm,&rank));
367aad13602SShrirang Abhyankar 
368*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
369*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
370aad13602SShrirang Abhyankar 
371dd400576SPatrick Sanan   if (rank == 0) {
372*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecGetArrayRead(Xseq,&x));
373aad13602SShrirang Abhyankar     ci = x[0]*x[0] - x[1];
374*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetValue(CI,0,ci,INSERT_VALUES));
375aad13602SShrirang Abhyankar     ci = -x[0]*x[0] + x[1] + 1.0;
376*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetValue(CI,1,ci,INSERT_VALUES));
377*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecRestoreArrayRead(Xseq,&x));
378aad13602SShrirang Abhyankar   }
379*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyBegin(CI));
380*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyEnd(CI));
381aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
382aad13602SShrirang Abhyankar }
383aad13602SShrirang Abhyankar 
384628da978Sresundermann /* Evaluate
385628da978Sresundermann    g = [ x0^2 + x1 - 2]
386628da978Sresundermann */
387aad13602SShrirang Abhyankar PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx)
388aad13602SShrirang Abhyankar {
389aad13602SShrirang Abhyankar   const PetscScalar *x;
390aad13602SShrirang Abhyankar   PetscScalar       ce;
391aad13602SShrirang Abhyankar   MPI_Comm          comm;
392aad13602SShrirang Abhyankar   PetscMPIInt       rank;
393aad13602SShrirang Abhyankar   AppCtx            *user=(AppCtx*)ctx;
394aad13602SShrirang Abhyankar   Vec               Xseq=user->Xseq;
395aad13602SShrirang Abhyankar   VecScatter        scat=user->scat;
396aad13602SShrirang Abhyankar 
397aad13602SShrirang Abhyankar   PetscFunctionBegin;
398*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm));
399*5f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm,&rank));
400aad13602SShrirang Abhyankar 
401*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
402*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
403aad13602SShrirang Abhyankar 
404dd400576SPatrick Sanan   if (rank == 0) {
405*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecGetArrayRead(Xseq,&x));
406aad13602SShrirang Abhyankar     ce = x[0]*x[0] + x[1] - 2.0;
407*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetValue(CE,0,ce,INSERT_VALUES));
408*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecRestoreArrayRead(Xseq,&x));
409aad13602SShrirang Abhyankar   }
410*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyBegin(CE));
411*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyEnd(CE));
412aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
413aad13602SShrirang Abhyankar }
414aad13602SShrirang Abhyankar 
415628da978Sresundermann /*
416628da978Sresundermann   grad h = [  2*x0, -1;
417628da978Sresundermann              -2*x0,  1]
418628da978Sresundermann */
419aad13602SShrirang Abhyankar PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre,  void *ctx)
420aad13602SShrirang Abhyankar {
421aad13602SShrirang Abhyankar   AppCtx            *user=(AppCtx*)ctx;
422f560b561SHong Zhang   PetscInt          zero=0,one=1,cols[2];
423aad13602SShrirang Abhyankar   PetscScalar       vals[2];
424aad13602SShrirang Abhyankar   const PetscScalar *x;
425aad13602SShrirang Abhyankar   Vec               Xseq=user->Xseq;
426aad13602SShrirang Abhyankar   VecScatter        scat=user->scat;
427aad13602SShrirang Abhyankar   MPI_Comm          comm;
428aad13602SShrirang Abhyankar   PetscMPIInt       rank;
429aad13602SShrirang Abhyankar 
430aad13602SShrirang Abhyankar   PetscFunctionBegin;
431*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm));
432*5f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm,&rank));
433*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
434*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
435aad13602SShrirang Abhyankar 
436*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArrayRead(Xseq,&x));
437f560b561SHong Zhang   if (!rank) {
438aad13602SShrirang Abhyankar     cols[0] = 0; cols[1] = 1;
439628da978Sresundermann     vals[0] = 2*x[0]; vals[1] = -1.0;
440*5f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetValues(JI,1,&zero,2,cols,vals,INSERT_VALUES));
441628da978Sresundermann     vals[0] = -2*x[0]; vals[1] = 1.0;
442*5f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetValues(JI,1,&one,2,cols,vals,INSERT_VALUES));
443aad13602SShrirang Abhyankar   }
444*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArrayRead(Xseq,&x));
445*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY));
446*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY));
447aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
448aad13602SShrirang Abhyankar }
449aad13602SShrirang Abhyankar 
450628da978Sresundermann /*
451628da978Sresundermann   grad g = [2*x0
452628da978Sresundermann              1.0 ]
453628da978Sresundermann */
454aad13602SShrirang Abhyankar PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx)
455aad13602SShrirang Abhyankar {
456f560b561SHong Zhang   PetscInt          zero=0,cols[2];
457aad13602SShrirang Abhyankar   PetscScalar       vals[2];
458aad13602SShrirang Abhyankar   const PetscScalar *x;
459aad13602SShrirang Abhyankar   PetscMPIInt       rank;
460aad13602SShrirang Abhyankar   MPI_Comm          comm;
461aad13602SShrirang Abhyankar 
462aad13602SShrirang Abhyankar   PetscFunctionBegin;
463*5f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm));
464*5f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm,&rank));
465aad13602SShrirang Abhyankar 
466dd400576SPatrick Sanan   if (rank == 0) {
467*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecGetArrayRead(X,&x));
468f560b561SHong Zhang     cols[0] = 0;       cols[1] = 1;
469aad13602SShrirang Abhyankar     vals[0] = 2*x[0];  vals[1] = 1.0;
470*5f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetValues(JE,1,&zero,2,cols,vals,INSERT_VALUES));
471*5f80ce2aSJacob Faibussowitsch     CHKERRQ(VecRestoreArrayRead(X,&x));
472aad13602SShrirang Abhyankar   }
473*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY));
474*5f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY));
475aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
476aad13602SShrirang Abhyankar }
477aad13602SShrirang Abhyankar 
478aad13602SShrirang Abhyankar /*TEST
479aad13602SShrirang Abhyankar 
480aad13602SShrirang Abhyankar    build:
481f560b561SHong Zhang       requires: !complex !defined(PETSC_USE_CXX)
482aad13602SShrirang Abhyankar 
483aad13602SShrirang Abhyankar    test:
48409ee8bb0SRylee Sundermann       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
485f560b561SHong Zhang       requires: mumps
486aad13602SShrirang Abhyankar 
487aad13602SShrirang Abhyankar    test:
488aad13602SShrirang Abhyankar       suffix: 2
489aad13602SShrirang Abhyankar       nsize: 2
49009ee8bb0SRylee Sundermann       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
491f560b561SHong Zhang       requires: mumps
492aad13602SShrirang Abhyankar 
4938e58fa1dSresundermann    test:
4948e58fa1dSresundermann       suffix: 3
4958e58fa1dSresundermann       args: -tao_converged_reason -no_eq
496f560b561SHong Zhang       requires: mumps
4978e58fa1dSresundermann 
4988e58fa1dSresundermann    test:
4998e58fa1dSresundermann       suffix: 4
5008e58fa1dSresundermann       nsize: 2
5018e58fa1dSresundermann       args: -tao_converged_reason -no_eq
502f560b561SHong Zhang       requires: mumps
5038e58fa1dSresundermann 
504661095bbSAlp Dener    test:
505661095bbSAlp Dener       suffix: 5
506661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm
507f560b561SHong Zhang       requires: mumps
508661095bbSAlp Dener 
509661095bbSAlp Dener    test:
510661095bbSAlp Dener       suffix: 6
511661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm -tao_almm_type phr
512f560b561SHong Zhang       requires: mumps
513661095bbSAlp Dener 
514661095bbSAlp Dener    test:
515661095bbSAlp Dener       suffix: 7
516661095bbSAlp Dener       nsize: 2
517f560b561SHong Zhang       requires: mumps
518661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm
519661095bbSAlp Dener 
520661095bbSAlp Dener    test:
521661095bbSAlp Dener       suffix: 8
522661095bbSAlp Dener       nsize: 2
523f560b561SHong Zhang       requires: cuda mumps
524661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm -vec_type cuda -mat_type aijcusparse
525661095bbSAlp Dener 
526661095bbSAlp Dener    test:
527661095bbSAlp Dener       suffix: 9
528661095bbSAlp Dener       nsize: 2
529661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm -no_eq
530f560b561SHong Zhang       requires: mumps
531661095bbSAlp Dener 
532661095bbSAlp Dener    test:
533661095bbSAlp Dener       suffix: 10
534661095bbSAlp Dener       nsize: 2
535661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm -tao_almm_type phr -no_eq
536f560b561SHong Zhang       requires: mumps
537661095bbSAlp Dener 
538aad13602SShrirang Abhyankar TEST*/
539