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 559371c9d4SSatish Balay PetscErrorCode main(int argc, char **argv) { 56aad13602SShrirang Abhyankar Tao tao; 57aad13602SShrirang Abhyankar KSP ksp; 58aad13602SShrirang Abhyankar PC pc; 59aad13602SShrirang Abhyankar AppCtx user; /* application context */ 60f560b561SHong Zhang Vec x, G, CI, CE; 61f560b561SHong Zhang PetscMPIInt size; 62661095bbSAlp Dener TaoType type; 63f560b561SHong Zhang PetscReal f; 64661095bbSAlp Dener PetscBool pdipm; 65aad13602SShrirang Abhyankar 66327415f7SBarry Smith PetscFunctionBeginUser; 679566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 689566063dSJacob Faibussowitsch PetscCallMPI(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 719566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "---- Constrained Problem -----\n")); 729566063dSJacob Faibussowitsch PetscCall(InitializeProblem(&user)); /* sets up problem, function below */ 739566063dSJacob Faibussowitsch PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao)); 749566063dSJacob Faibussowitsch PetscCall(TaoSetType(tao, TAOPDIPM)); 759566063dSJacob Faibussowitsch PetscCall(TaoSetSolution(tao, user.x)); /* gets solution vector from problem */ 769566063dSJacob Faibussowitsch PetscCall(TaoSetVariableBounds(tao, user.xl, user.xu)); /* sets lower upper bounds from given solution */ 779566063dSJacob Faibussowitsch PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user)); 78aad13602SShrirang Abhyankar 79*48a46eb9SPierre Jolivet if (!user.noeqflag) PetscCall(TaoSetEqualityConstraintsRoutine(tao, user.ce, FormEqualityConstraints, (void *)&user)); 809566063dSJacob Faibussowitsch PetscCall(TaoSetInequalityConstraintsRoutine(tao, user.ci, FormInequalityConstraints, (void *)&user)); 819371c9d4SSatish Balay if (!user.noeqflag) { PetscCall(TaoSetJacobianEqualityRoutine(tao, user.Ae, user.Ae, FormEqualityJacobian, (void *)&user)); /* equality jacobian */ } 829566063dSJacob Faibussowitsch PetscCall(TaoSetJacobianInequalityRoutine(tao, user.Ai, user.Ai, FormInequalityJacobian, (void *)&user)); /* inequality jacobian */ 839566063dSJacob Faibussowitsch PetscCall(TaoSetTolerances(tao, 1.e-6, 1.e-6, 1.e-6)); 849566063dSJacob Faibussowitsch PetscCall(TaoSetConstraintTolerances(tao, 1.e-6, 1.e-6)); 85aad13602SShrirang Abhyankar 869566063dSJacob Faibussowitsch PetscCall(TaoGetKSP(tao, &ksp)); 879566063dSJacob Faibussowitsch PetscCall(KSPGetPC(ksp, &pc)); 889566063dSJacob Faibussowitsch PetscCall(PCSetType(pc, PCCHOLESKY)); 89aad13602SShrirang Abhyankar /* 90aad13602SShrirang Abhyankar This algorithm produces matrices with zeros along the diagonal therefore we use 9112d688e0SRylee Sundermann MUMPS which provides solver for indefinite matrices 92aad13602SShrirang Abhyankar */ 93f560b561SHong Zhang #if defined(PETSC_HAVE_MUMPS) 949566063dSJacob Faibussowitsch PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS)); /* requires mumps to solve pdipm */ 95f560b561SHong Zhang #else 963c859ba3SBarry Smith PetscCheck(size == 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Requires an external package that supports parallel PCCHOLESKY, e.g., MUMPS."); 97f560b561SHong Zhang #endif 989566063dSJacob Faibussowitsch PetscCall(KSPSetType(ksp, KSPPREONLY)); 999566063dSJacob Faibussowitsch PetscCall(KSPSetFromOptions(ksp)); 100aad13602SShrirang Abhyankar 1019566063dSJacob Faibussowitsch PetscCall(TaoSetFromOptions(tao)); 1029566063dSJacob Faibussowitsch PetscCall(TaoGetType(tao, &type)); 1039566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOPDIPM, &pdipm)); 104661095bbSAlp Dener if (pdipm) { 1059566063dSJacob Faibussowitsch PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user)); 106f560b561SHong Zhang if (user.initview) { 1079566063dSJacob Faibussowitsch PetscCall(TaoSetUp(tao)); 1089566063dSJacob Faibussowitsch PetscCall(VecDuplicate(user.x, &G)); 1099566063dSJacob Faibussowitsch PetscCall(FormFunctionGradient(tao, user.x, &f, G, (void *)&user)); 1109566063dSJacob Faibussowitsch PetscCall(FormHessian(tao, user.x, user.H, user.H, (void *)&user)); 1119566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD)); 11263a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial point X:\n")); 1139566063dSJacob Faibussowitsch PetscCall(VecView(user.x, PETSC_VIEWER_STDOUT_WORLD)); 11463a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial objective f(x) = %g\n", (double)f)); 11563a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial gradient and Hessian:\n")); 1169566063dSJacob Faibussowitsch PetscCall(VecView(G, PETSC_VIEWER_STDOUT_WORLD)); 1179566063dSJacob Faibussowitsch PetscCall(MatView(user.H, PETSC_VIEWER_STDOUT_WORLD)); 1189566063dSJacob Faibussowitsch PetscCall(VecDestroy(&G)); 1199566063dSJacob Faibussowitsch PetscCall(FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void *)&user)); 1209566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(user.Ai, NULL, &CI)); 1219566063dSJacob Faibussowitsch PetscCall(FormInequalityConstraints(tao, user.x, CI, (void *)&user)); 12263a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial inequality constraints and Jacobian:\n")); 1239566063dSJacob Faibussowitsch PetscCall(VecView(CI, PETSC_VIEWER_STDOUT_WORLD)); 1249566063dSJacob Faibussowitsch PetscCall(MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD)); 1259566063dSJacob Faibussowitsch PetscCall(VecDestroy(&CI)); 126f560b561SHong Zhang if (!user.noeqflag) { 1279566063dSJacob Faibussowitsch PetscCall(FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void *)&user)); 1289566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(user.Ae, NULL, &CE)); 1299566063dSJacob Faibussowitsch PetscCall(FormEqualityConstraints(tao, user.x, CE, (void *)&user)); 13063a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial equality constraints and Jacobian:\n")); 1319566063dSJacob Faibussowitsch PetscCall(VecView(CE, PETSC_VIEWER_STDOUT_WORLD)); 1329566063dSJacob Faibussowitsch PetscCall(MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD)); 1339566063dSJacob Faibussowitsch PetscCall(VecDestroy(&CE)); 134f560b561SHong Zhang } 1359566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n")); 1369566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD)); 137f560b561SHong Zhang } 138661095bbSAlp Dener } 139aad13602SShrirang Abhyankar 1409566063dSJacob Faibussowitsch PetscCall(TaoSolve(tao)); 1419566063dSJacob Faibussowitsch PetscCall(TaoGetSolution(tao, &x)); 1429566063dSJacob Faibussowitsch PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD)); 143aad13602SShrirang Abhyankar 144aad13602SShrirang Abhyankar /* Free objects */ 1459566063dSJacob Faibussowitsch PetscCall(DestroyProblem(&user)); 1469566063dSJacob Faibussowitsch PetscCall(TaoDestroy(&tao)); 1479566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 148b122ec5aSJacob Faibussowitsch return 0; 149aad13602SShrirang Abhyankar } 150aad13602SShrirang Abhyankar 1519371c9d4SSatish Balay PetscErrorCode InitializeProblem(AppCtx *user) { 152aad13602SShrirang Abhyankar PetscMPIInt size; 153aad13602SShrirang Abhyankar PetscMPIInt rank; 154aad13602SShrirang Abhyankar PetscInt nloc, neloc, niloc; 155aad13602SShrirang Abhyankar 156aad13602SShrirang Abhyankar PetscFunctionBegin; 1579566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 1589566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 1598e58fa1dSresundermann user->noeqflag = PETSC_FALSE; 1609566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_eq", &user->noeqflag, NULL)); 161f560b561SHong Zhang user->initview = PETSC_FALSE; 1629566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL, NULL, "-init_view", &user->initview, NULL)); 163628da978Sresundermann 1648e58fa1dSresundermann if (!user->noeqflag) { 165f560b561SHong Zhang /* Tell user the correct solution, not an error checking */ 1669566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Solution should be f(1,1)=-2\n")); 1678e58fa1dSresundermann } 168aad13602SShrirang Abhyankar 1692d4ee042Sprj- /* create vector x and set initial values */ 170aad13602SShrirang Abhyankar user->n = 2; /* global length */ 171f560b561SHong Zhang nloc = (size == 1) ? user->n : 1; 1729566063dSJacob Faibussowitsch PetscCall(VecCreate(PETSC_COMM_WORLD, &user->x)); 1739566063dSJacob Faibussowitsch PetscCall(VecSetSizes(user->x, nloc, user->n)); 1749566063dSJacob Faibussowitsch PetscCall(VecSetFromOptions(user->x)); 1759566063dSJacob Faibussowitsch PetscCall(VecSet(user->x, 0)); 176aad13602SShrirang Abhyankar 177aad13602SShrirang Abhyankar /* create and set lower and upper bound vectors */ 1789566063dSJacob Faibussowitsch PetscCall(VecDuplicate(user->x, &user->xl)); 1799566063dSJacob Faibussowitsch PetscCall(VecDuplicate(user->x, &user->xu)); 1809566063dSJacob Faibussowitsch PetscCall(VecSet(user->xl, -1.0)); 1819566063dSJacob Faibussowitsch PetscCall(VecSet(user->xu, 2.0)); 182aad13602SShrirang Abhyankar 183aad13602SShrirang Abhyankar /* create scater to zero */ 1849566063dSJacob Faibussowitsch PetscCall(VecScatterCreateToZero(user->x, &user->scat, &user->Xseq)); 185aad13602SShrirang Abhyankar 186aad13602SShrirang Abhyankar user->ne = 1; 187628da978Sresundermann user->ni = 2; 188aad13602SShrirang Abhyankar neloc = (rank == 0) ? user->ne : 0; 189f560b561SHong Zhang niloc = (size == 1) ? user->ni : 1; 190aad13602SShrirang Abhyankar 1918e58fa1dSresundermann if (!user->noeqflag) { 1929566063dSJacob Faibussowitsch PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ce)); /* a 1x1 vec for equality constraints */ 1939566063dSJacob Faibussowitsch PetscCall(VecSetSizes(user->ce, neloc, user->ne)); 1949566063dSJacob Faibussowitsch PetscCall(VecSetFromOptions(user->ce)); 1959566063dSJacob Faibussowitsch PetscCall(VecSetUp(user->ce)); 1968e58fa1dSresundermann } 197628da978Sresundermann 1989566063dSJacob Faibussowitsch PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ci)); /* a 2x1 vec for inequality constraints */ 1999566063dSJacob Faibussowitsch PetscCall(VecSetSizes(user->ci, niloc, user->ni)); 2009566063dSJacob Faibussowitsch PetscCall(VecSetFromOptions(user->ci)); 2019566063dSJacob Faibussowitsch PetscCall(VecSetUp(user->ci)); 202aad13602SShrirang Abhyankar 203a5b23f4aSJose E. Roman /* nexn & nixn matricies for equally and inequalty constraints */ 2048e58fa1dSresundermann if (!user->noeqflag) { 2059566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ae)); 2069566063dSJacob Faibussowitsch PetscCall(MatSetSizes(user->Ae, neloc, nloc, user->ne, user->n)); 2079566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(user->Ae)); 2089566063dSJacob Faibussowitsch PetscCall(MatSetUp(user->Ae)); 2098e58fa1dSresundermann } 2108e58fa1dSresundermann 2119566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ai)); 2129566063dSJacob Faibussowitsch PetscCall(MatSetSizes(user->Ai, niloc, nloc, user->ni, user->n)); 2139566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(user->Ai)); 2149566063dSJacob Faibussowitsch PetscCall(MatSetUp(user->Ai)); 215628da978Sresundermann 2169566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD, &user->H)); 2179566063dSJacob Faibussowitsch PetscCall(MatSetSizes(user->H, nloc, nloc, user->n, user->n)); 2189566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(user->H)); 2199566063dSJacob Faibussowitsch PetscCall(MatSetUp(user->H)); 220aad13602SShrirang Abhyankar PetscFunctionReturn(0); 221aad13602SShrirang Abhyankar } 222aad13602SShrirang Abhyankar 2239371c9d4SSatish Balay PetscErrorCode DestroyProblem(AppCtx *user) { 224aad13602SShrirang Abhyankar PetscFunctionBegin; 225*48a46eb9SPierre Jolivet if (!user->noeqflag) PetscCall(MatDestroy(&user->Ae)); 2269566063dSJacob Faibussowitsch PetscCall(MatDestroy(&user->Ai)); 2279566063dSJacob Faibussowitsch PetscCall(MatDestroy(&user->H)); 228aad13602SShrirang Abhyankar 2299566063dSJacob Faibussowitsch PetscCall(VecDestroy(&user->x)); 230*48a46eb9SPierre Jolivet if (!user->noeqflag) PetscCall(VecDestroy(&user->ce)); 2319566063dSJacob Faibussowitsch PetscCall(VecDestroy(&user->ci)); 2329566063dSJacob Faibussowitsch PetscCall(VecDestroy(&user->xl)); 2339566063dSJacob Faibussowitsch PetscCall(VecDestroy(&user->xu)); 2349566063dSJacob Faibussowitsch PetscCall(VecDestroy(&user->Xseq)); 2359566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&user->scat)); 236aad13602SShrirang Abhyankar PetscFunctionReturn(0); 237aad13602SShrirang Abhyankar } 238aad13602SShrirang Abhyankar 239628da978Sresundermann /* Evaluate 240628da978Sresundermann f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1) 241628da978Sresundermann G = grad f = [2*(x0 - 2) - 2; 242628da978Sresundermann 2*(x1 - 2) - 2] 243aad13602SShrirang Abhyankar */ 2449371c9d4SSatish Balay PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx) { 245aad13602SShrirang Abhyankar PetscScalar g; 246aad13602SShrirang Abhyankar const PetscScalar *x; 247aad13602SShrirang Abhyankar MPI_Comm comm; 248aad13602SShrirang Abhyankar PetscMPIInt rank; 249aad13602SShrirang Abhyankar PetscReal fin; 250aad13602SShrirang Abhyankar AppCtx *user = (AppCtx *)ctx; 251aad13602SShrirang Abhyankar Vec Xseq = user->Xseq; 252aad13602SShrirang Abhyankar VecScatter scat = user->scat; 253aad13602SShrirang Abhyankar 254aad13602SShrirang Abhyankar PetscFunctionBegin; 2559566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 2569566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 257aad13602SShrirang Abhyankar 2589566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 2599566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 260aad13602SShrirang Abhyankar 261aad13602SShrirang Abhyankar fin = 0.0; 262dd400576SPatrick Sanan if (rank == 0) { 2639566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(Xseq, &x)); 264aad13602SShrirang Abhyankar fin = (x[0] - 2.0) * (x[0] - 2.0) + (x[1] - 2.0) * (x[1] - 2.0) - 2.0 * (x[0] + x[1]); 265aad13602SShrirang Abhyankar g = 2.0 * (x[0] - 2.0) - 2.0; 2669566063dSJacob Faibussowitsch PetscCall(VecSetValue(G, 0, g, INSERT_VALUES)); 267aad13602SShrirang Abhyankar g = 2.0 * (x[1] - 2.0) - 2.0; 2689566063dSJacob Faibussowitsch PetscCall(VecSetValue(G, 1, g, INSERT_VALUES)); 2699566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Xseq, &x)); 270aad13602SShrirang Abhyankar } 2719566063dSJacob Faibussowitsch PetscCallMPI(MPI_Allreduce(&fin, f, 1, MPIU_REAL, MPIU_SUM, comm)); 2729566063dSJacob Faibussowitsch PetscCall(VecAssemblyBegin(G)); 2739566063dSJacob Faibussowitsch PetscCall(VecAssemblyEnd(G)); 274aad13602SShrirang Abhyankar PetscFunctionReturn(0); 275aad13602SShrirang Abhyankar } 276aad13602SShrirang Abhyankar 277628da978Sresundermann /* Evaluate 278628da978Sresundermann H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)] 279628da978Sresundermann = [ 2*(1+de[0]-di[0]+di[1]), 0; 280628da978Sresundermann 0, 2] 281628da978Sresundermann */ 2829371c9d4SSatish Balay PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx) { 283aad13602SShrirang Abhyankar AppCtx *user = (AppCtx *)ctx; 284aad13602SShrirang Abhyankar Vec DE, DI; 285aad13602SShrirang Abhyankar const PetscScalar *de, *di; 286aad13602SShrirang Abhyankar PetscInt zero = 0, one = 1; 287aad13602SShrirang Abhyankar PetscScalar two = 2.0; 288aad13602SShrirang Abhyankar PetscScalar val = 0.0; 289f560b561SHong Zhang Vec Deseq, Diseq; 290f560b561SHong Zhang VecScatter Descat, Discat; 291aad13602SShrirang Abhyankar PetscMPIInt rank; 292aad13602SShrirang Abhyankar MPI_Comm comm; 293aad13602SShrirang Abhyankar 294aad13602SShrirang Abhyankar PetscFunctionBegin; 2959566063dSJacob Faibussowitsch PetscCall(TaoGetDualVariables(tao, &DE, &DI)); 296aad13602SShrirang Abhyankar 2979566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 2989566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 299aad13602SShrirang Abhyankar 3008e58fa1dSresundermann if (!user->noeqflag) { 3019566063dSJacob Faibussowitsch PetscCall(VecScatterCreateToZero(DE, &Descat, &Deseq)); 3029566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD)); 3039566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD)); 3048e58fa1dSresundermann } 3059566063dSJacob Faibussowitsch PetscCall(VecScatterCreateToZero(DI, &Discat, &Diseq)); 3069566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD)); 3079566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD)); 308aad13602SShrirang Abhyankar 309dd400576SPatrick Sanan if (rank == 0) { 3109371c9d4SSatish Balay if (!user->noeqflag) { PetscCall(VecGetArrayRead(Deseq, &de)); /* places equality constraint dual into array */ } 3119566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(Diseq, &di)); /* places inequality constraint dual into array */ 312628da978Sresundermann 3138e58fa1dSresundermann if (!user->noeqflag) { 314628da978Sresundermann val = 2.0 * (1 + de[0] - di[0] + di[1]); 3159566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Deseq, &de)); 3169566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Diseq, &di)); 3178e58fa1dSresundermann } else { 318628da978Sresundermann val = 2.0 * (1 - di[0] + di[1]); 3198e58fa1dSresundermann } 3209566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Diseq, &di)); 3219566063dSJacob Faibussowitsch PetscCall(MatSetValues(H, 1, &zero, 1, &zero, &val, INSERT_VALUES)); 3229566063dSJacob Faibussowitsch PetscCall(MatSetValues(H, 1, &one, 1, &one, &two, INSERT_VALUES)); 323aad13602SShrirang Abhyankar } 3249566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY)); 3259566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY)); 3268e58fa1dSresundermann if (!user->noeqflag) { 3279566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&Descat)); 3289566063dSJacob Faibussowitsch PetscCall(VecDestroy(&Deseq)); 3298e58fa1dSresundermann } 3309566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&Discat)); 3319566063dSJacob Faibussowitsch PetscCall(VecDestroy(&Diseq)); 332aad13602SShrirang Abhyankar PetscFunctionReturn(0); 333aad13602SShrirang Abhyankar } 334aad13602SShrirang Abhyankar 335628da978Sresundermann /* Evaluate 336628da978Sresundermann h = [ x0^2 - x1; 337628da978Sresundermann 1 -(x0^2 - x1)] 338628da978Sresundermann */ 3399371c9d4SSatish Balay PetscErrorCode FormInequalityConstraints(Tao tao, Vec X, Vec CI, void *ctx) { 340aad13602SShrirang Abhyankar const PetscScalar *x; 341aad13602SShrirang Abhyankar PetscScalar ci; 342aad13602SShrirang Abhyankar MPI_Comm comm; 343aad13602SShrirang Abhyankar PetscMPIInt rank; 344aad13602SShrirang Abhyankar AppCtx *user = (AppCtx *)ctx; 345aad13602SShrirang Abhyankar Vec Xseq = user->Xseq; 346aad13602SShrirang Abhyankar VecScatter scat = user->scat; 347aad13602SShrirang Abhyankar 348aad13602SShrirang Abhyankar PetscFunctionBegin; 3499566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 3509566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 351aad13602SShrirang Abhyankar 3529566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 3539566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 354aad13602SShrirang Abhyankar 355dd400576SPatrick Sanan if (rank == 0) { 3569566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(Xseq, &x)); 357aad13602SShrirang Abhyankar ci = x[0] * x[0] - x[1]; 3589566063dSJacob Faibussowitsch PetscCall(VecSetValue(CI, 0, ci, INSERT_VALUES)); 359aad13602SShrirang Abhyankar ci = -x[0] * x[0] + x[1] + 1.0; 3609566063dSJacob Faibussowitsch PetscCall(VecSetValue(CI, 1, ci, INSERT_VALUES)); 3619566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Xseq, &x)); 362aad13602SShrirang Abhyankar } 3639566063dSJacob Faibussowitsch PetscCall(VecAssemblyBegin(CI)); 3649566063dSJacob Faibussowitsch PetscCall(VecAssemblyEnd(CI)); 365aad13602SShrirang Abhyankar PetscFunctionReturn(0); 366aad13602SShrirang Abhyankar } 367aad13602SShrirang Abhyankar 368628da978Sresundermann /* Evaluate 369628da978Sresundermann g = [ x0^2 + x1 - 2] 370628da978Sresundermann */ 3719371c9d4SSatish Balay PetscErrorCode FormEqualityConstraints(Tao tao, Vec X, Vec CE, void *ctx) { 372aad13602SShrirang Abhyankar const PetscScalar *x; 373aad13602SShrirang Abhyankar PetscScalar ce; 374aad13602SShrirang Abhyankar MPI_Comm comm; 375aad13602SShrirang Abhyankar PetscMPIInt rank; 376aad13602SShrirang Abhyankar AppCtx *user = (AppCtx *)ctx; 377aad13602SShrirang Abhyankar Vec Xseq = user->Xseq; 378aad13602SShrirang Abhyankar VecScatter scat = user->scat; 379aad13602SShrirang Abhyankar 380aad13602SShrirang Abhyankar PetscFunctionBegin; 3819566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 3829566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 383aad13602SShrirang Abhyankar 3849566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 3859566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 386aad13602SShrirang Abhyankar 387dd400576SPatrick Sanan if (rank == 0) { 3889566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(Xseq, &x)); 389aad13602SShrirang Abhyankar ce = x[0] * x[0] + x[1] - 2.0; 3909566063dSJacob Faibussowitsch PetscCall(VecSetValue(CE, 0, ce, INSERT_VALUES)); 3919566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Xseq, &x)); 392aad13602SShrirang Abhyankar } 3939566063dSJacob Faibussowitsch PetscCall(VecAssemblyBegin(CE)); 3949566063dSJacob Faibussowitsch PetscCall(VecAssemblyEnd(CE)); 395aad13602SShrirang Abhyankar PetscFunctionReturn(0); 396aad13602SShrirang Abhyankar } 397aad13602SShrirang Abhyankar 398628da978Sresundermann /* 399628da978Sresundermann grad h = [ 2*x0, -1; 400628da978Sresundermann -2*x0, 1] 401628da978Sresundermann */ 4029371c9d4SSatish Balay PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx) { 403aad13602SShrirang Abhyankar AppCtx *user = (AppCtx *)ctx; 404f560b561SHong Zhang PetscInt zero = 0, one = 1, cols[2]; 405aad13602SShrirang Abhyankar PetscScalar vals[2]; 406aad13602SShrirang Abhyankar const PetscScalar *x; 407aad13602SShrirang Abhyankar Vec Xseq = user->Xseq; 408aad13602SShrirang Abhyankar VecScatter scat = user->scat; 409aad13602SShrirang Abhyankar MPI_Comm comm; 410aad13602SShrirang Abhyankar PetscMPIInt rank; 411aad13602SShrirang Abhyankar 412aad13602SShrirang Abhyankar PetscFunctionBegin; 4139566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 4149566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 4159566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 4169566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 417aad13602SShrirang Abhyankar 4189566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(Xseq, &x)); 419c5853193SPierre Jolivet if (rank == 0) { 4209371c9d4SSatish Balay cols[0] = 0; 4219371c9d4SSatish Balay cols[1] = 1; 4229371c9d4SSatish Balay vals[0] = 2 * x[0]; 4239371c9d4SSatish Balay vals[1] = -1.0; 4249566063dSJacob Faibussowitsch PetscCall(MatSetValues(JI, 1, &zero, 2, cols, vals, INSERT_VALUES)); 4259371c9d4SSatish Balay vals[0] = -2 * x[0]; 4269371c9d4SSatish Balay vals[1] = 1.0; 4279566063dSJacob Faibussowitsch PetscCall(MatSetValues(JI, 1, &one, 2, cols, vals, INSERT_VALUES)); 428aad13602SShrirang Abhyankar } 4299566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Xseq, &x)); 4309566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(JI, MAT_FINAL_ASSEMBLY)); 4319566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(JI, MAT_FINAL_ASSEMBLY)); 432aad13602SShrirang Abhyankar PetscFunctionReturn(0); 433aad13602SShrirang Abhyankar } 434aad13602SShrirang Abhyankar 435628da978Sresundermann /* 436628da978Sresundermann grad g = [2*x0 437628da978Sresundermann 1.0 ] 438628da978Sresundermann */ 4399371c9d4SSatish Balay PetscErrorCode FormEqualityJacobian(Tao tao, Vec X, Mat JE, Mat JEpre, void *ctx) { 440f560b561SHong Zhang PetscInt zero = 0, cols[2]; 441aad13602SShrirang Abhyankar PetscScalar vals[2]; 442aad13602SShrirang Abhyankar const PetscScalar *x; 443aad13602SShrirang Abhyankar PetscMPIInt rank; 444aad13602SShrirang Abhyankar MPI_Comm comm; 445aad13602SShrirang Abhyankar 446aad13602SShrirang Abhyankar PetscFunctionBegin; 4479566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 4489566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 449aad13602SShrirang Abhyankar 450dd400576SPatrick Sanan if (rank == 0) { 4519566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(X, &x)); 4529371c9d4SSatish Balay cols[0] = 0; 4539371c9d4SSatish Balay cols[1] = 1; 4549371c9d4SSatish Balay vals[0] = 2 * x[0]; 4559371c9d4SSatish Balay vals[1] = 1.0; 4569566063dSJacob Faibussowitsch PetscCall(MatSetValues(JE, 1, &zero, 2, cols, vals, INSERT_VALUES)); 4579566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(X, &x)); 458aad13602SShrirang Abhyankar } 4599566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(JE, MAT_FINAL_ASSEMBLY)); 4609566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(JE, MAT_FINAL_ASSEMBLY)); 461aad13602SShrirang Abhyankar PetscFunctionReturn(0); 462aad13602SShrirang Abhyankar } 463aad13602SShrirang Abhyankar 464aad13602SShrirang Abhyankar /*TEST 465aad13602SShrirang Abhyankar 466aad13602SShrirang Abhyankar build: 467f560b561SHong Zhang requires: !complex !defined(PETSC_USE_CXX) 468aad13602SShrirang Abhyankar 469aad13602SShrirang Abhyankar test: 47009ee8bb0SRylee Sundermann args: -tao_converged_reason -tao_pdipm_kkt_shift_pd 471f560b561SHong Zhang requires: mumps 472aad13602SShrirang Abhyankar 473aad13602SShrirang Abhyankar test: 474aad13602SShrirang Abhyankar suffix: 2 475aad13602SShrirang Abhyankar nsize: 2 47609ee8bb0SRylee Sundermann args: -tao_converged_reason -tao_pdipm_kkt_shift_pd 477f560b561SHong Zhang requires: mumps 478aad13602SShrirang Abhyankar 4798e58fa1dSresundermann test: 4808e58fa1dSresundermann suffix: 3 4818e58fa1dSresundermann args: -tao_converged_reason -no_eq 482f560b561SHong Zhang requires: mumps 4838e58fa1dSresundermann 4848e58fa1dSresundermann test: 4858e58fa1dSresundermann suffix: 4 4868e58fa1dSresundermann nsize: 2 4878e58fa1dSresundermann args: -tao_converged_reason -no_eq 488f560b561SHong Zhang requires: mumps 4898e58fa1dSresundermann 490661095bbSAlp Dener test: 491661095bbSAlp Dener suffix: 5 492661095bbSAlp Dener args: -tao_cmonitor -tao_type almm 493f560b561SHong Zhang requires: mumps 494661095bbSAlp Dener 495661095bbSAlp Dener test: 496661095bbSAlp Dener suffix: 6 497661095bbSAlp Dener args: -tao_cmonitor -tao_type almm -tao_almm_type phr 498f560b561SHong Zhang requires: mumps 499661095bbSAlp Dener 500661095bbSAlp Dener test: 501661095bbSAlp Dener suffix: 7 502661095bbSAlp Dener nsize: 2 503f560b561SHong Zhang requires: mumps 504661095bbSAlp Dener args: -tao_cmonitor -tao_type almm 505661095bbSAlp Dener 506661095bbSAlp Dener test: 507661095bbSAlp Dener suffix: 8 508661095bbSAlp Dener nsize: 2 509f560b561SHong Zhang requires: cuda mumps 510661095bbSAlp Dener args: -tao_cmonitor -tao_type almm -vec_type cuda -mat_type aijcusparse 511661095bbSAlp Dener 512661095bbSAlp Dener test: 513661095bbSAlp Dener suffix: 9 514661095bbSAlp Dener nsize: 2 515661095bbSAlp Dener args: -tao_cmonitor -tao_type almm -no_eq 516f560b561SHong Zhang requires: mumps 517661095bbSAlp Dener 518661095bbSAlp Dener test: 519661095bbSAlp Dener suffix: 10 520661095bbSAlp Dener nsize: 2 521661095bbSAlp Dener args: -tao_cmonitor -tao_type almm -tao_almm_type phr -no_eq 522f560b561SHong Zhang requires: mumps 523661095bbSAlp Dener 524aad13602SShrirang Abhyankar TEST*/ 525