xref: /petsc/src/tao/constrained/tutorials/ex1.c (revision b122ec5aa1bd4469eb4e0673542fb7de3f411254)
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   Tao            tao;
58aad13602SShrirang Abhyankar   KSP            ksp;
59aad13602SShrirang Abhyankar   PC             pc;
60aad13602SShrirang Abhyankar   AppCtx         user;  /* application context */
61f560b561SHong Zhang   Vec            x,G,CI,CE;
62f560b561SHong Zhang   PetscMPIInt    size;
63661095bbSAlp Dener   TaoType        type;
64f560b561SHong Zhang   PetscReal      f;
65661095bbSAlp Dener   PetscBool      pdipm;
66aad13602SShrirang Abhyankar 
67*b122ec5aSJacob Faibussowitsch   CHKERRQ(PetscInitialize(&argc,&argv,(char*)0,help));
685f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
693c859ba3SBarry 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.");
70aad13602SShrirang Abhyankar 
715f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"---- Constrained Problem -----\n"));
725f80ce2aSJacob Faibussowitsch   CHKERRQ(InitializeProblem(&user)); /* sets up problem, function below */
735f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoCreate(PETSC_COMM_WORLD,&tao));
745f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetType(tao,TAOPDIPM));
755f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetSolution(tao,user.x)); /* gets solution vector from problem */
765f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetVariableBounds(tao,user.xl,user.xu)); /* sets lower upper bounds from given solution */
775f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void*)&user));
78aad13602SShrirang Abhyankar 
798e58fa1dSresundermann   if (!user.noeqflag) {
805f80ce2aSJacob Faibussowitsch     CHKERRQ(TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user));
818e58fa1dSresundermann   }
825f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user));
838e58fa1dSresundermann   if (!user.noeqflag) {
845f80ce2aSJacob Faibussowitsch     CHKERRQ(TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user)); /* equality jacobian */
858e58fa1dSresundermann   }
865f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user)); /* inequality jacobian */
875f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetTolerances(tao,1.e-6,1.e-6,1.e-6));
885f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetConstraintTolerances(tao,1.e-6,1.e-6));
89aad13602SShrirang Abhyankar 
905f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoGetKSP(tao,&ksp));
915f80ce2aSJacob Faibussowitsch   CHKERRQ(KSPGetPC(ksp,&pc));
925f80ce2aSJacob Faibussowitsch   CHKERRQ(PCSetType(pc,PCCHOLESKY));
93aad13602SShrirang Abhyankar   /*
94aad13602SShrirang Abhyankar       This algorithm produces matrices with zeros along the diagonal therefore we use
9512d688e0SRylee Sundermann     MUMPS which provides solver for indefinite matrices
96aad13602SShrirang Abhyankar   */
97f560b561SHong Zhang #if defined(PETSC_HAVE_MUMPS)
985f80ce2aSJacob Faibussowitsch   CHKERRQ(PCFactorSetMatSolverType(pc,MATSOLVERMUMPS));  /* requires mumps to solve pdipm */
99f560b561SHong Zhang #else
1003c859ba3SBarry Smith   PetscCheck(size == 1,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Requires an external package that supports parallel PCCHOLESKY, e.g., MUMPS.");
101f560b561SHong Zhang #endif
1025f80ce2aSJacob Faibussowitsch   CHKERRQ(KSPSetType(ksp,KSPPREONLY));
1035f80ce2aSJacob Faibussowitsch   CHKERRQ(KSPSetFromOptions(ksp));
104aad13602SShrirang Abhyankar 
1055f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetFromOptions(tao));
1065f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoGetType(tao,&type));
1075f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectTypeCompare((PetscObject)tao,TAOPDIPM,&pdipm));
108661095bbSAlp Dener   if (pdipm) {
1095f80ce2aSJacob Faibussowitsch     CHKERRQ(TaoSetHessian(tao,user.H,user.H,FormHessian,(void*)&user));
110f560b561SHong Zhang     if (user.initview) {
1115f80ce2aSJacob Faibussowitsch       CHKERRQ(TaoSetUp(tao));
1125f80ce2aSJacob Faibussowitsch       CHKERRQ(VecDuplicate(user.x, &G));
1135f80ce2aSJacob Faibussowitsch       CHKERRQ(FormFunctionGradient(tao, user.x, &f, G, (void*)&user));
1145f80ce2aSJacob Faibussowitsch       CHKERRQ(FormHessian(tao, user.x, user.H, user.H, (void*)&user));
1155f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD));
1165f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nInitial point X:\n",f));
1175f80ce2aSJacob Faibussowitsch       CHKERRQ(VecView(user.x, PETSC_VIEWER_STDOUT_WORLD));
1185f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nInitial objective f(x) = %g\n",f));
1195f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nInitial gradient and Hessian:\n",f));
1205f80ce2aSJacob Faibussowitsch       CHKERRQ(VecView(G, PETSC_VIEWER_STDOUT_WORLD));
1215f80ce2aSJacob Faibussowitsch       CHKERRQ(MatView(user.H, PETSC_VIEWER_STDOUT_WORLD));
1225f80ce2aSJacob Faibussowitsch       CHKERRQ(VecDestroy(&G));
1235f80ce2aSJacob Faibussowitsch       CHKERRQ(FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void*)&user));
1245f80ce2aSJacob Faibussowitsch       CHKERRQ(MatCreateVecs(user.Ai, NULL, &CI));
1255f80ce2aSJacob Faibussowitsch       CHKERRQ(FormInequalityConstraints(tao, user.x, CI, (void*)&user));
1265f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nInitial inequality constraints and Jacobian:\n",f));
1275f80ce2aSJacob Faibussowitsch       CHKERRQ(VecView(CI, PETSC_VIEWER_STDOUT_WORLD));
1285f80ce2aSJacob Faibussowitsch       CHKERRQ(MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD));
1295f80ce2aSJacob Faibussowitsch       CHKERRQ(VecDestroy(&CI));
130f560b561SHong Zhang       if (!user.noeqflag) {
1315f80ce2aSJacob Faibussowitsch         CHKERRQ(FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void*)&user));
1325f80ce2aSJacob Faibussowitsch         CHKERRQ(MatCreateVecs(user.Ae, NULL, &CE));
1335f80ce2aSJacob Faibussowitsch         CHKERRQ(FormEqualityConstraints(tao, user.x, CE, (void*)&user));
1345f80ce2aSJacob Faibussowitsch         CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nInitial equality constraints and Jacobian:\n",f));
1355f80ce2aSJacob Faibussowitsch         CHKERRQ(VecView(CE, PETSC_VIEWER_STDOUT_WORLD));
1365f80ce2aSJacob Faibussowitsch         CHKERRQ(MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD));
1375f80ce2aSJacob Faibussowitsch         CHKERRQ(VecDestroy(&CE));
138f560b561SHong Zhang       }
1395f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "\n"));
1405f80ce2aSJacob Faibussowitsch       CHKERRQ(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD));
141f560b561SHong Zhang     }
142661095bbSAlp Dener   }
143aad13602SShrirang Abhyankar 
1445f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSolve(tao));
1455f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoGetSolution(tao,&x));
1465f80ce2aSJacob Faibussowitsch   CHKERRQ(VecView(x,PETSC_VIEWER_STDOUT_WORLD));
147aad13602SShrirang Abhyankar 
148aad13602SShrirang Abhyankar   /* Free objects */
1495f80ce2aSJacob Faibussowitsch   CHKERRQ(DestroyProblem(&user));
1505f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoDestroy(&tao));
151*b122ec5aSJacob Faibussowitsch   CHKERRQ(PetscFinalize());
152*b122ec5aSJacob Faibussowitsch   return 0;
153aad13602SShrirang Abhyankar }
154aad13602SShrirang Abhyankar 
155aad13602SShrirang Abhyankar PetscErrorCode InitializeProblem(AppCtx *user)
156aad13602SShrirang Abhyankar {
157aad13602SShrirang Abhyankar   PetscMPIInt    size;
158aad13602SShrirang Abhyankar   PetscMPIInt    rank;
159aad13602SShrirang Abhyankar   PetscInt       nloc,neloc,niloc;
160aad13602SShrirang Abhyankar 
161aad13602SShrirang Abhyankar   PetscFunctionBegin;
1625f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
1635f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
1648e58fa1dSresundermann   user->noeqflag = PETSC_FALSE;
1655f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-no_eq",&user->noeqflag,NULL));
166f560b561SHong Zhang   user->initview = PETSC_FALSE;
1675f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-init_view",&user->initview,NULL));
168628da978Sresundermann 
1698e58fa1dSresundermann   if (!user->noeqflag) {
170f560b561SHong Zhang     /* Tell user the correct solution, not an error checking */
1715f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n"));
1728e58fa1dSresundermann   }
173aad13602SShrirang Abhyankar 
1742d4ee042Sprj-   /* create vector x and set initial values */
175aad13602SShrirang Abhyankar   user->n = 2; /* global length */
176f560b561SHong Zhang   nloc = (size==1)?user->n:1;
1775f80ce2aSJacob Faibussowitsch   CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user->x));
1785f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSetSizes(user->x,nloc,user->n));
1795f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSetFromOptions(user->x));
1805f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSet(user->x,0));
181aad13602SShrirang Abhyankar 
182aad13602SShrirang Abhyankar   /* create and set lower and upper bound vectors */
1835f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDuplicate(user->x,&user->xl));
1845f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDuplicate(user->x,&user->xu));
1855f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSet(user->xl,-1.0));
1865f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSet(user->xu,2.0));
187aad13602SShrirang Abhyankar 
188aad13602SShrirang Abhyankar   /* create scater to zero */
1895f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterCreateToZero(user->x,&user->scat,&user->Xseq));
190aad13602SShrirang Abhyankar 
191aad13602SShrirang Abhyankar   user->ne = 1;
192628da978Sresundermann   user->ni = 2;
193aad13602SShrirang Abhyankar   neloc = (rank==0)?user->ne:0;
194f560b561SHong Zhang   niloc = (size==1)?user->ni:1;
195aad13602SShrirang Abhyankar 
1968e58fa1dSresundermann   if (!user->noeqflag) {
1975f80ce2aSJacob Faibussowitsch     CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user->ce)); /* a 1x1 vec for equality constraints */
1985f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetSizes(user->ce,neloc,user->ne));
1995f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetFromOptions(user->ce));
2005f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetUp(user->ce));
2018e58fa1dSresundermann   }
202628da978Sresundermann 
2035f80ce2aSJacob Faibussowitsch   CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user->ci)); /* a 2x1 vec for inequality constraints */
2045f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSetSizes(user->ci,niloc,user->ni));
2055f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSetFromOptions(user->ci));
2065f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSetUp(user->ci));
207aad13602SShrirang Abhyankar 
208a5b23f4aSJose E. Roman   /* nexn & nixn matricies for equally and inequalty constraints */
2098e58fa1dSresundermann   if (!user->noeqflag) {
2105f80ce2aSJacob Faibussowitsch     CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user->Ae));
2115f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n));
2125f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetFromOptions(user->Ae));
2135f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetUp(user->Ae));
2148e58fa1dSresundermann   }
2158e58fa1dSresundermann 
2165f80ce2aSJacob Faibussowitsch   CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user->Ai));
2175f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n));
2185f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetFromOptions(user->Ai));
2195f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetUp(user->Ai));
220628da978Sresundermann 
2215f80ce2aSJacob Faibussowitsch   CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user->H));
2225f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetSizes(user->H,nloc,nloc,user->n,user->n));
2235f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetFromOptions(user->H));
2245f80ce2aSJacob Faibussowitsch   CHKERRQ(MatSetUp(user->H));
225aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
226aad13602SShrirang Abhyankar }
227aad13602SShrirang Abhyankar 
228aad13602SShrirang Abhyankar PetscErrorCode DestroyProblem(AppCtx *user)
229aad13602SShrirang Abhyankar {
230aad13602SShrirang Abhyankar   PetscFunctionBegin;
2318e58fa1dSresundermann   if (!user->noeqflag) {
2325f80ce2aSJacob Faibussowitsch     CHKERRQ(MatDestroy(&user->Ae));
2338e58fa1dSresundermann   }
2345f80ce2aSJacob Faibussowitsch   CHKERRQ(MatDestroy(&user->Ai));
2355f80ce2aSJacob Faibussowitsch   CHKERRQ(MatDestroy(&user->H));
236aad13602SShrirang Abhyankar 
2375f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&user->x));
2388e58fa1dSresundermann   if (!user->noeqflag) {
2395f80ce2aSJacob Faibussowitsch     CHKERRQ(VecDestroy(&user->ce));
2408e58fa1dSresundermann   }
2415f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&user->ci));
2425f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&user->xl));
2435f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&user->xu));
2445f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&user->Xseq));
2455f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterDestroy(&user->scat));
246aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
247aad13602SShrirang Abhyankar }
248aad13602SShrirang Abhyankar 
249628da978Sresundermann /* Evaluate
250628da978Sresundermann    f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
251628da978Sresundermann    G = grad f = [2*(x0 - 2) - 2;
252628da978Sresundermann                  2*(x1 - 2) - 2]
253aad13602SShrirang Abhyankar */
254aad13602SShrirang Abhyankar PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
255aad13602SShrirang Abhyankar {
256aad13602SShrirang Abhyankar   PetscScalar       g;
257aad13602SShrirang Abhyankar   const PetscScalar *x;
258aad13602SShrirang Abhyankar   MPI_Comm          comm;
259aad13602SShrirang Abhyankar   PetscMPIInt       rank;
260aad13602SShrirang Abhyankar   PetscReal         fin;
261aad13602SShrirang Abhyankar   AppCtx            *user=(AppCtx*)ctx;
262aad13602SShrirang Abhyankar   Vec               Xseq=user->Xseq;
263aad13602SShrirang Abhyankar   VecScatter        scat=user->scat;
264aad13602SShrirang Abhyankar 
265aad13602SShrirang Abhyankar   PetscFunctionBegin;
2665f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm));
2675f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm,&rank));
268aad13602SShrirang Abhyankar 
2695f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
2705f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
271aad13602SShrirang Abhyankar 
272aad13602SShrirang Abhyankar   fin = 0.0;
273dd400576SPatrick Sanan   if (rank == 0) {
2745f80ce2aSJacob Faibussowitsch     CHKERRQ(VecGetArrayRead(Xseq,&x));
275aad13602SShrirang Abhyankar     fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]);
276aad13602SShrirang Abhyankar     g = 2.0*(x[0]-2.0) - 2.0;
2775f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetValue(G,0,g,INSERT_VALUES));
278aad13602SShrirang Abhyankar     g = 2.0*(x[1]-2.0) - 2.0;
2795f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetValue(G,1,g,INSERT_VALUES));
2805f80ce2aSJacob Faibussowitsch     CHKERRQ(VecRestoreArrayRead(Xseq,&x));
281aad13602SShrirang Abhyankar   }
2825f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm));
2835f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyBegin(G));
2845f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyEnd(G));
285aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
286aad13602SShrirang Abhyankar }
287aad13602SShrirang Abhyankar 
288628da978Sresundermann /* Evaluate
289628da978Sresundermann    H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)]
290628da978Sresundermann      = [ 2*(1+de[0]-di[0]+di[1]), 0;
291628da978Sresundermann                    0,             2]
292628da978Sresundermann */
293aad13602SShrirang Abhyankar PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx)
294aad13602SShrirang Abhyankar {
295aad13602SShrirang Abhyankar   AppCtx            *user=(AppCtx*)ctx;
296aad13602SShrirang Abhyankar   Vec               DE,DI;
297aad13602SShrirang Abhyankar   const PetscScalar *de,*di;
298aad13602SShrirang Abhyankar   PetscInt          zero=0,one=1;
299aad13602SShrirang Abhyankar   PetscScalar       two=2.0;
300aad13602SShrirang Abhyankar   PetscScalar       val=0.0;
301f560b561SHong Zhang   Vec               Deseq,Diseq;
302f560b561SHong Zhang   VecScatter        Descat,Discat;
303aad13602SShrirang Abhyankar   PetscMPIInt       rank;
304aad13602SShrirang Abhyankar   MPI_Comm          comm;
305aad13602SShrirang Abhyankar 
306aad13602SShrirang Abhyankar   PetscFunctionBegin;
3075f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoGetDualVariables(tao,&DE,&DI));
308aad13602SShrirang Abhyankar 
3095f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm));
3105f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm,&rank));
311aad13602SShrirang Abhyankar 
3128e58fa1dSresundermann   if (!user->noeqflag) {
3135f80ce2aSJacob Faibussowitsch    CHKERRQ(VecScatterCreateToZero(DE,&Descat,&Deseq));
3145f80ce2aSJacob Faibussowitsch    CHKERRQ(VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD));
3155f80ce2aSJacob Faibussowitsch    CHKERRQ(VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD));
3168e58fa1dSresundermann   }
3175f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterCreateToZero(DI,&Discat,&Diseq));
3185f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD));
3195f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD));
320aad13602SShrirang Abhyankar 
321dd400576SPatrick Sanan   if (rank == 0) {
3228e58fa1dSresundermann     if (!user->noeqflag) {
3235f80ce2aSJacob Faibussowitsch       CHKERRQ(VecGetArrayRead(Deseq,&de));  /* places equality constraint dual into array */
3248e58fa1dSresundermann     }
3255f80ce2aSJacob Faibussowitsch     CHKERRQ(VecGetArrayRead(Diseq,&di));  /* places inequality constraint dual into array */
326628da978Sresundermann 
3278e58fa1dSresundermann     if (!user->noeqflag) {
328628da978Sresundermann       val = 2.0 * (1 + de[0] - di[0] + di[1]);
3295f80ce2aSJacob Faibussowitsch       CHKERRQ(VecRestoreArrayRead(Deseq,&de));
3305f80ce2aSJacob Faibussowitsch       CHKERRQ(VecRestoreArrayRead(Diseq,&di));
3318e58fa1dSresundermann     } else {
332628da978Sresundermann       val = 2.0 * (1 - di[0] + di[1]);
3338e58fa1dSresundermann     }
3345f80ce2aSJacob Faibussowitsch     CHKERRQ(VecRestoreArrayRead(Diseq,&di));
3355f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES));
3365f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES));
337aad13602SShrirang Abhyankar   }
3385f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY));
3395f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY));
3408e58fa1dSresundermann   if (!user->noeqflag) {
3415f80ce2aSJacob Faibussowitsch     CHKERRQ(VecScatterDestroy(&Descat));
3425f80ce2aSJacob Faibussowitsch     CHKERRQ(VecDestroy(&Deseq));
3438e58fa1dSresundermann   }
3445f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterDestroy(&Discat));
3455f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&Diseq));
346aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
347aad13602SShrirang Abhyankar }
348aad13602SShrirang Abhyankar 
349628da978Sresundermann /* Evaluate
350628da978Sresundermann    h = [ x0^2 - x1;
351628da978Sresundermann          1 -(x0^2 - x1)]
352628da978Sresundermann */
353aad13602SShrirang Abhyankar PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx)
354aad13602SShrirang Abhyankar {
355aad13602SShrirang Abhyankar   const PetscScalar *x;
356aad13602SShrirang Abhyankar   PetscScalar       ci;
357aad13602SShrirang Abhyankar   MPI_Comm          comm;
358aad13602SShrirang Abhyankar   PetscMPIInt       rank;
359aad13602SShrirang Abhyankar   AppCtx            *user=(AppCtx*)ctx;
360aad13602SShrirang Abhyankar   Vec               Xseq=user->Xseq;
361aad13602SShrirang Abhyankar   VecScatter        scat=user->scat;
362aad13602SShrirang Abhyankar 
363aad13602SShrirang Abhyankar   PetscFunctionBegin;
3645f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm));
3655f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm,&rank));
366aad13602SShrirang Abhyankar 
3675f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
3685f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
369aad13602SShrirang Abhyankar 
370dd400576SPatrick Sanan   if (rank == 0) {
3715f80ce2aSJacob Faibussowitsch     CHKERRQ(VecGetArrayRead(Xseq,&x));
372aad13602SShrirang Abhyankar     ci = x[0]*x[0] - x[1];
3735f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetValue(CI,0,ci,INSERT_VALUES));
374aad13602SShrirang Abhyankar     ci = -x[0]*x[0] + x[1] + 1.0;
3755f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetValue(CI,1,ci,INSERT_VALUES));
3765f80ce2aSJacob Faibussowitsch     CHKERRQ(VecRestoreArrayRead(Xseq,&x));
377aad13602SShrirang Abhyankar   }
3785f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyBegin(CI));
3795f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyEnd(CI));
380aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
381aad13602SShrirang Abhyankar }
382aad13602SShrirang Abhyankar 
383628da978Sresundermann /* Evaluate
384628da978Sresundermann    g = [ x0^2 + x1 - 2]
385628da978Sresundermann */
386aad13602SShrirang Abhyankar PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx)
387aad13602SShrirang Abhyankar {
388aad13602SShrirang Abhyankar   const PetscScalar *x;
389aad13602SShrirang Abhyankar   PetscScalar       ce;
390aad13602SShrirang Abhyankar   MPI_Comm          comm;
391aad13602SShrirang Abhyankar   PetscMPIInt       rank;
392aad13602SShrirang Abhyankar   AppCtx            *user=(AppCtx*)ctx;
393aad13602SShrirang Abhyankar   Vec               Xseq=user->Xseq;
394aad13602SShrirang Abhyankar   VecScatter        scat=user->scat;
395aad13602SShrirang Abhyankar 
396aad13602SShrirang Abhyankar   PetscFunctionBegin;
3975f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm));
3985f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm,&rank));
399aad13602SShrirang Abhyankar 
4005f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
4015f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
402aad13602SShrirang Abhyankar 
403dd400576SPatrick Sanan   if (rank == 0) {
4045f80ce2aSJacob Faibussowitsch     CHKERRQ(VecGetArrayRead(Xseq,&x));
405aad13602SShrirang Abhyankar     ce = x[0]*x[0] + x[1] - 2.0;
4065f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetValue(CE,0,ce,INSERT_VALUES));
4075f80ce2aSJacob Faibussowitsch     CHKERRQ(VecRestoreArrayRead(Xseq,&x));
408aad13602SShrirang Abhyankar   }
4095f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyBegin(CE));
4105f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyEnd(CE));
411aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
412aad13602SShrirang Abhyankar }
413aad13602SShrirang Abhyankar 
414628da978Sresundermann /*
415628da978Sresundermann   grad h = [  2*x0, -1;
416628da978Sresundermann              -2*x0,  1]
417628da978Sresundermann */
418aad13602SShrirang Abhyankar PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre,  void *ctx)
419aad13602SShrirang Abhyankar {
420aad13602SShrirang Abhyankar   AppCtx            *user=(AppCtx*)ctx;
421f560b561SHong Zhang   PetscInt          zero=0,one=1,cols[2];
422aad13602SShrirang Abhyankar   PetscScalar       vals[2];
423aad13602SShrirang Abhyankar   const PetscScalar *x;
424aad13602SShrirang Abhyankar   Vec               Xseq=user->Xseq;
425aad13602SShrirang Abhyankar   VecScatter        scat=user->scat;
426aad13602SShrirang Abhyankar   MPI_Comm          comm;
427aad13602SShrirang Abhyankar   PetscMPIInt       rank;
428aad13602SShrirang Abhyankar 
429aad13602SShrirang Abhyankar   PetscFunctionBegin;
4305f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm));
4315f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm,&rank));
4325f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
4335f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD));
434aad13602SShrirang Abhyankar 
4355f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArrayRead(Xseq,&x));
436f560b561SHong Zhang   if (!rank) {
437aad13602SShrirang Abhyankar     cols[0] = 0; cols[1] = 1;
438628da978Sresundermann     vals[0] = 2*x[0]; vals[1] = -1.0;
4395f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetValues(JI,1,&zero,2,cols,vals,INSERT_VALUES));
440628da978Sresundermann     vals[0] = -2*x[0]; vals[1] = 1.0;
4415f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetValues(JI,1,&one,2,cols,vals,INSERT_VALUES));
442aad13602SShrirang Abhyankar   }
4435f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArrayRead(Xseq,&x));
4445f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY));
4455f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY));
446aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
447aad13602SShrirang Abhyankar }
448aad13602SShrirang Abhyankar 
449628da978Sresundermann /*
450628da978Sresundermann   grad g = [2*x0
451628da978Sresundermann              1.0 ]
452628da978Sresundermann */
453aad13602SShrirang Abhyankar PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx)
454aad13602SShrirang Abhyankar {
455f560b561SHong Zhang   PetscInt          zero=0,cols[2];
456aad13602SShrirang Abhyankar   PetscScalar       vals[2];
457aad13602SShrirang Abhyankar   const PetscScalar *x;
458aad13602SShrirang Abhyankar   PetscMPIInt       rank;
459aad13602SShrirang Abhyankar   MPI_Comm          comm;
460aad13602SShrirang Abhyankar 
461aad13602SShrirang Abhyankar   PetscFunctionBegin;
4625f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm));
4635f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_rank(comm,&rank));
464aad13602SShrirang Abhyankar 
465dd400576SPatrick Sanan   if (rank == 0) {
4665f80ce2aSJacob Faibussowitsch     CHKERRQ(VecGetArrayRead(X,&x));
467f560b561SHong Zhang     cols[0] = 0;       cols[1] = 1;
468aad13602SShrirang Abhyankar     vals[0] = 2*x[0];  vals[1] = 1.0;
4695f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetValues(JE,1,&zero,2,cols,vals,INSERT_VALUES));
4705f80ce2aSJacob Faibussowitsch     CHKERRQ(VecRestoreArrayRead(X,&x));
471aad13602SShrirang Abhyankar   }
4725f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY));
4735f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY));
474aad13602SShrirang Abhyankar   PetscFunctionReturn(0);
475aad13602SShrirang Abhyankar }
476aad13602SShrirang Abhyankar 
477aad13602SShrirang Abhyankar /*TEST
478aad13602SShrirang Abhyankar 
479aad13602SShrirang Abhyankar    build:
480f560b561SHong Zhang       requires: !complex !defined(PETSC_USE_CXX)
481aad13602SShrirang Abhyankar 
482aad13602SShrirang Abhyankar    test:
48309ee8bb0SRylee Sundermann       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
484f560b561SHong Zhang       requires: mumps
485aad13602SShrirang Abhyankar 
486aad13602SShrirang Abhyankar    test:
487aad13602SShrirang Abhyankar       suffix: 2
488aad13602SShrirang Abhyankar       nsize: 2
48909ee8bb0SRylee Sundermann       args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
490f560b561SHong Zhang       requires: mumps
491aad13602SShrirang Abhyankar 
4928e58fa1dSresundermann    test:
4938e58fa1dSresundermann       suffix: 3
4948e58fa1dSresundermann       args: -tao_converged_reason -no_eq
495f560b561SHong Zhang       requires: mumps
4968e58fa1dSresundermann 
4978e58fa1dSresundermann    test:
4988e58fa1dSresundermann       suffix: 4
4998e58fa1dSresundermann       nsize: 2
5008e58fa1dSresundermann       args: -tao_converged_reason -no_eq
501f560b561SHong Zhang       requires: mumps
5028e58fa1dSresundermann 
503661095bbSAlp Dener    test:
504661095bbSAlp Dener       suffix: 5
505661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm
506f560b561SHong Zhang       requires: mumps
507661095bbSAlp Dener 
508661095bbSAlp Dener    test:
509661095bbSAlp Dener       suffix: 6
510661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm -tao_almm_type phr
511f560b561SHong Zhang       requires: mumps
512661095bbSAlp Dener 
513661095bbSAlp Dener    test:
514661095bbSAlp Dener       suffix: 7
515661095bbSAlp Dener       nsize: 2
516f560b561SHong Zhang       requires: mumps
517661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm
518661095bbSAlp Dener 
519661095bbSAlp Dener    test:
520661095bbSAlp Dener       suffix: 8
521661095bbSAlp Dener       nsize: 2
522f560b561SHong Zhang       requires: cuda mumps
523661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm -vec_type cuda -mat_type aijcusparse
524661095bbSAlp Dener 
525661095bbSAlp Dener    test:
526661095bbSAlp Dener       suffix: 9
527661095bbSAlp Dener       nsize: 2
528661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm -no_eq
529f560b561SHong Zhang       requires: mumps
530661095bbSAlp Dener 
531661095bbSAlp Dener    test:
532661095bbSAlp Dener       suffix: 10
533661095bbSAlp Dener       nsize: 2
534661095bbSAlp Dener       args: -tao_cmonitor -tao_type almm -tao_almm_type phr -no_eq
535f560b561SHong Zhang       requires: mumps
536661095bbSAlp Dener 
537aad13602SShrirang Abhyankar TEST*/
538