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 20da81f932SPierre Jolivet static char help[] = "Solves constrained optimization 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 55d71ae5a4SJacob Faibussowitsch PetscErrorCode main(int argc, char **argv) 56d71ae5a4SJacob Faibussowitsch { 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; 659a09327dSAlp Dener PetscBool pdipm, mumps; 66aad13602SShrirang Abhyankar 67327415f7SBarry Smith PetscFunctionBeginUser; 689566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 699566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 709a09327dSAlp Dener PetscCheck(size <= 2, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "More than 2 processors detected. Example written to use max of 2 processors."); 71aad13602SShrirang Abhyankar 729566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "---- Constrained Problem -----\n")); 739566063dSJacob Faibussowitsch PetscCall(InitializeProblem(&user)); /* sets up problem, function below */ 749a09327dSAlp Dener 759566063dSJacob Faibussowitsch PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao)); 769a09327dSAlp Dener PetscCall(TaoSetType(tao, TAOALMM)); 779a09327dSAlp Dener PetscCall(TaoSetSolution(tao, user.x)); 789a09327dSAlp Dener PetscCall(TaoSetVariableBounds(tao, user.xl, user.xu)); 799566063dSJacob Faibussowitsch PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user)); 809a09327dSAlp Dener PetscCall(TaoSetTolerances(tao, 1.e-4, 0.0, 0.0)); 819a09327dSAlp Dener PetscCall(TaoSetConstraintTolerances(tao, 1.e-4, 0.0)); 829a09327dSAlp Dener PetscCall(TaoSetMaximumFunctionEvaluations(tao, 1e6)); 839566063dSJacob Faibussowitsch PetscCall(TaoSetFromOptions(tao)); 849a09327dSAlp Dener 859a09327dSAlp Dener if (!user.noeqflag) { 869a09327dSAlp Dener PetscCall(TaoSetEqualityConstraintsRoutine(tao, user.ce, FormEqualityConstraints, (void *)&user)); 879a09327dSAlp Dener PetscCall(TaoSetJacobianEqualityRoutine(tao, user.Ae, user.Ae, FormEqualityJacobian, (void *)&user)); 889a09327dSAlp Dener } 899a09327dSAlp Dener PetscCall(TaoSetInequalityConstraintsRoutine(tao, user.ci, FormInequalityConstraints, (void *)&user)); 909a09327dSAlp Dener PetscCall(TaoSetJacobianInequalityRoutine(tao, user.Ai, user.Ai, FormInequalityJacobian, (void *)&user)); 919a09327dSAlp Dener 929566063dSJacob Faibussowitsch PetscCall(TaoGetType(tao, &type)); 939566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOPDIPM, &pdipm)); 94661095bbSAlp Dener if (pdipm) { 959a09327dSAlp Dener /* 969a09327dSAlp Dener PDIPM produces an indefinite KKT matrix with some zeros along the diagonal 979a09327dSAlp Dener Inverting this indefinite matrix requires PETSc to be configured with MUMPS 989a09327dSAlp Dener */ 999a09327dSAlp Dener PetscCall(PetscHasExternalPackage("mumps", &mumps)); 1009a09327dSAlp Dener PetscCheck(mumps, PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "TAOPDIPM requires PETSc to be configured with MUMPS (--download-mumps)"); 1019a09327dSAlp Dener PetscCall(TaoGetKSP(tao, &ksp)); 1029a09327dSAlp Dener PetscCall(KSPGetPC(ksp, &pc)); 1039a09327dSAlp Dener PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS)); 1049a09327dSAlp Dener PetscCall(KSPSetType(ksp, KSPPREONLY)); 1059a09327dSAlp Dener PetscCall(KSPSetFromOptions(ksp)); 1069566063dSJacob Faibussowitsch PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user)); 1079a09327dSAlp Dener } 1089a09327dSAlp Dener 1099a09327dSAlp Dener /* Print out an initial view of the problem */ 110f560b561SHong Zhang if (user.initview) { 1119566063dSJacob Faibussowitsch PetscCall(TaoSetUp(tao)); 1129566063dSJacob Faibussowitsch PetscCall(VecDuplicate(user.x, &G)); 1139566063dSJacob Faibussowitsch PetscCall(FormFunctionGradient(tao, user.x, &f, G, (void *)&user)); 1149566063dSJacob Faibussowitsch PetscCall(FormHessian(tao, user.x, user.H, user.H, (void *)&user)); 1159566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD)); 11663a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial point X:\n")); 1179566063dSJacob Faibussowitsch PetscCall(VecView(user.x, PETSC_VIEWER_STDOUT_WORLD)); 11863a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial objective f(x) = %g\n", (double)f)); 11963a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial gradient and Hessian:\n")); 1209566063dSJacob Faibussowitsch PetscCall(VecView(G, PETSC_VIEWER_STDOUT_WORLD)); 1219566063dSJacob Faibussowitsch PetscCall(MatView(user.H, PETSC_VIEWER_STDOUT_WORLD)); 1229566063dSJacob Faibussowitsch PetscCall(VecDestroy(&G)); 1239566063dSJacob Faibussowitsch PetscCall(FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void *)&user)); 1249566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(user.Ai, NULL, &CI)); 1259566063dSJacob Faibussowitsch PetscCall(FormInequalityConstraints(tao, user.x, CI, (void *)&user)); 12663a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial inequality constraints and Jacobian:\n")); 1279566063dSJacob Faibussowitsch PetscCall(VecView(CI, PETSC_VIEWER_STDOUT_WORLD)); 1289566063dSJacob Faibussowitsch PetscCall(MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD)); 1299566063dSJacob Faibussowitsch PetscCall(VecDestroy(&CI)); 130f560b561SHong Zhang if (!user.noeqflag) { 1319566063dSJacob Faibussowitsch PetscCall(FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void *)&user)); 1329566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(user.Ae, NULL, &CE)); 1339566063dSJacob Faibussowitsch PetscCall(FormEqualityConstraints(tao, user.x, CE, (void *)&user)); 13463a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial equality constraints and Jacobian:\n")); 1359566063dSJacob Faibussowitsch PetscCall(VecView(CE, PETSC_VIEWER_STDOUT_WORLD)); 1369566063dSJacob Faibussowitsch PetscCall(MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD)); 1379566063dSJacob Faibussowitsch PetscCall(VecDestroy(&CE)); 138f560b561SHong Zhang } 1399566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n")); 1409566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD)); 141f560b561SHong Zhang } 142aad13602SShrirang Abhyankar 1439566063dSJacob Faibussowitsch PetscCall(TaoSolve(tao)); 1449566063dSJacob Faibussowitsch PetscCall(TaoGetSolution(tao, &x)); 1459a09327dSAlp Dener PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Found solution:\n")); 1469566063dSJacob Faibussowitsch PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD)); 147aad13602SShrirang Abhyankar 148aad13602SShrirang Abhyankar /* Free objects */ 1499566063dSJacob Faibussowitsch PetscCall(DestroyProblem(&user)); 1509566063dSJacob Faibussowitsch PetscCall(TaoDestroy(&tao)); 1519566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 152b122ec5aSJacob Faibussowitsch return 0; 153aad13602SShrirang Abhyankar } 154aad13602SShrirang Abhyankar 155d71ae5a4SJacob Faibussowitsch PetscErrorCode InitializeProblem(AppCtx *user) 156d71ae5a4SJacob Faibussowitsch { 157aad13602SShrirang Abhyankar PetscMPIInt size; 158aad13602SShrirang Abhyankar PetscMPIInt rank; 159aad13602SShrirang Abhyankar PetscInt nloc, neloc, niloc; 160aad13602SShrirang Abhyankar 161aad13602SShrirang Abhyankar PetscFunctionBegin; 1629566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 1639566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 1648e58fa1dSresundermann user->noeqflag = PETSC_FALSE; 1659566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_eq", &user->noeqflag, NULL)); 166f560b561SHong Zhang user->initview = PETSC_FALSE; 1679566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL, NULL, "-init_view", &user->initview, NULL)); 168628da978Sresundermann 169f560b561SHong Zhang /* Tell user the correct solution, not an error checking */ 1709a09327dSAlp Dener if (!user->noeqflag) { 1719a09327dSAlp Dener PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution: f(1, 1) = -2\n")); 1729a09327dSAlp Dener } else { 1739a09327dSAlp Dener PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution (-no_eq): f(1.73205, 2) = -7.3923\n")); 1748e58fa1dSresundermann } 175aad13602SShrirang Abhyankar 1762d4ee042Sprj- /* create vector x and set initial values */ 177aad13602SShrirang Abhyankar user->n = 2; /* global length */ 178f560b561SHong Zhang nloc = (size == 1) ? user->n : 1; 1799566063dSJacob Faibussowitsch PetscCall(VecCreate(PETSC_COMM_WORLD, &user->x)); 1809566063dSJacob Faibussowitsch PetscCall(VecSetSizes(user->x, nloc, user->n)); 1819566063dSJacob Faibussowitsch PetscCall(VecSetFromOptions(user->x)); 1829a09327dSAlp Dener PetscCall(VecSet(user->x, 0.0)); 183aad13602SShrirang Abhyankar 184aad13602SShrirang Abhyankar /* create and set lower and upper bound vectors */ 1859566063dSJacob Faibussowitsch PetscCall(VecDuplicate(user->x, &user->xl)); 1869566063dSJacob Faibussowitsch PetscCall(VecDuplicate(user->x, &user->xu)); 1879566063dSJacob Faibussowitsch PetscCall(VecSet(user->xl, -1.0)); 1889566063dSJacob Faibussowitsch PetscCall(VecSet(user->xu, 2.0)); 189aad13602SShrirang Abhyankar 190aad13602SShrirang Abhyankar /* create scater to zero */ 1919566063dSJacob Faibussowitsch PetscCall(VecScatterCreateToZero(user->x, &user->scat, &user->Xseq)); 192aad13602SShrirang Abhyankar 193aad13602SShrirang Abhyankar user->ne = 1; 194628da978Sresundermann user->ni = 2; 195aad13602SShrirang Abhyankar neloc = (rank == 0) ? user->ne : 0; 196f560b561SHong Zhang niloc = (size == 1) ? user->ni : 1; 197aad13602SShrirang Abhyankar 1988e58fa1dSresundermann if (!user->noeqflag) { 1999566063dSJacob Faibussowitsch PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ce)); /* a 1x1 vec for equality constraints */ 2009566063dSJacob Faibussowitsch PetscCall(VecSetSizes(user->ce, neloc, user->ne)); 2019566063dSJacob Faibussowitsch PetscCall(VecSetFromOptions(user->ce)); 2029566063dSJacob Faibussowitsch PetscCall(VecSetUp(user->ce)); 2038e58fa1dSresundermann } 204628da978Sresundermann 2059566063dSJacob Faibussowitsch PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ci)); /* a 2x1 vec for inequality constraints */ 2069566063dSJacob Faibussowitsch PetscCall(VecSetSizes(user->ci, niloc, user->ni)); 2079566063dSJacob Faibussowitsch PetscCall(VecSetFromOptions(user->ci)); 2089566063dSJacob Faibussowitsch PetscCall(VecSetUp(user->ci)); 209aad13602SShrirang Abhyankar 210a5b23f4aSJose E. Roman /* nexn & nixn matricies for equally and inequalty constraints */ 2118e58fa1dSresundermann if (!user->noeqflag) { 2129566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ae)); 2139566063dSJacob Faibussowitsch PetscCall(MatSetSizes(user->Ae, neloc, nloc, user->ne, user->n)); 2149566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(user->Ae)); 2159566063dSJacob Faibussowitsch PetscCall(MatSetUp(user->Ae)); 2168e58fa1dSresundermann } 2178e58fa1dSresundermann 2189566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ai)); 2199566063dSJacob Faibussowitsch PetscCall(MatSetSizes(user->Ai, niloc, nloc, user->ni, user->n)); 2209566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(user->Ai)); 2219566063dSJacob Faibussowitsch PetscCall(MatSetUp(user->Ai)); 222628da978Sresundermann 2239566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD, &user->H)); 2249566063dSJacob Faibussowitsch PetscCall(MatSetSizes(user->H, nloc, nloc, user->n, user->n)); 2259566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(user->H)); 2269566063dSJacob Faibussowitsch PetscCall(MatSetUp(user->H)); 227aad13602SShrirang Abhyankar PetscFunctionReturn(0); 228aad13602SShrirang Abhyankar } 229aad13602SShrirang Abhyankar 230d71ae5a4SJacob Faibussowitsch PetscErrorCode DestroyProblem(AppCtx *user) 231d71ae5a4SJacob Faibussowitsch { 232aad13602SShrirang Abhyankar PetscFunctionBegin; 23348a46eb9SPierre Jolivet if (!user->noeqflag) PetscCall(MatDestroy(&user->Ae)); 2349566063dSJacob Faibussowitsch PetscCall(MatDestroy(&user->Ai)); 2359566063dSJacob Faibussowitsch PetscCall(MatDestroy(&user->H)); 236aad13602SShrirang Abhyankar 2379566063dSJacob Faibussowitsch PetscCall(VecDestroy(&user->x)); 23848a46eb9SPierre Jolivet if (!user->noeqflag) PetscCall(VecDestroy(&user->ce)); 2399566063dSJacob Faibussowitsch PetscCall(VecDestroy(&user->ci)); 2409566063dSJacob Faibussowitsch PetscCall(VecDestroy(&user->xl)); 2419566063dSJacob Faibussowitsch PetscCall(VecDestroy(&user->xu)); 2429566063dSJacob Faibussowitsch PetscCall(VecDestroy(&user->Xseq)); 2439566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&user->scat)); 244aad13602SShrirang Abhyankar PetscFunctionReturn(0); 245aad13602SShrirang Abhyankar } 246aad13602SShrirang Abhyankar 247628da978Sresundermann /* Evaluate 248628da978Sresundermann f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1) 249628da978Sresundermann G = grad f = [2*(x0 - 2) - 2; 250628da978Sresundermann 2*(x1 - 2) - 2] 251aad13602SShrirang Abhyankar */ 252d71ae5a4SJacob Faibussowitsch PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx) 253d71ae5a4SJacob Faibussowitsch { 254aad13602SShrirang Abhyankar PetscScalar g; 255aad13602SShrirang Abhyankar const PetscScalar *x; 256aad13602SShrirang Abhyankar MPI_Comm comm; 257aad13602SShrirang Abhyankar PetscMPIInt rank; 258aad13602SShrirang Abhyankar PetscReal fin; 259aad13602SShrirang Abhyankar AppCtx *user = (AppCtx *)ctx; 260aad13602SShrirang Abhyankar Vec Xseq = user->Xseq; 261aad13602SShrirang Abhyankar VecScatter scat = user->scat; 262aad13602SShrirang Abhyankar 263aad13602SShrirang Abhyankar PetscFunctionBegin; 2649566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 2659566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 266aad13602SShrirang Abhyankar 2679566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 2689566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 269aad13602SShrirang Abhyankar 270aad13602SShrirang Abhyankar fin = 0.0; 271dd400576SPatrick Sanan if (rank == 0) { 2729566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(Xseq, &x)); 273aad13602SShrirang Abhyankar fin = (x[0] - 2.0) * (x[0] - 2.0) + (x[1] - 2.0) * (x[1] - 2.0) - 2.0 * (x[0] + x[1]); 274aad13602SShrirang Abhyankar g = 2.0 * (x[0] - 2.0) - 2.0; 2759566063dSJacob Faibussowitsch PetscCall(VecSetValue(G, 0, g, INSERT_VALUES)); 276aad13602SShrirang Abhyankar g = 2.0 * (x[1] - 2.0) - 2.0; 2779566063dSJacob Faibussowitsch PetscCall(VecSetValue(G, 1, g, INSERT_VALUES)); 2789566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Xseq, &x)); 279aad13602SShrirang Abhyankar } 2809566063dSJacob Faibussowitsch PetscCallMPI(MPI_Allreduce(&fin, f, 1, MPIU_REAL, MPIU_SUM, comm)); 2819566063dSJacob Faibussowitsch PetscCall(VecAssemblyBegin(G)); 2829566063dSJacob Faibussowitsch PetscCall(VecAssemblyEnd(G)); 283aad13602SShrirang Abhyankar PetscFunctionReturn(0); 284aad13602SShrirang Abhyankar } 285aad13602SShrirang Abhyankar 286628da978Sresundermann /* Evaluate 287628da978Sresundermann H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)] 288628da978Sresundermann = [ 2*(1+de[0]-di[0]+di[1]), 0; 289628da978Sresundermann 0, 2] 290628da978Sresundermann */ 291d71ae5a4SJacob Faibussowitsch PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx) 292d71ae5a4SJacob Faibussowitsch { 293aad13602SShrirang Abhyankar AppCtx *user = (AppCtx *)ctx; 294aad13602SShrirang Abhyankar Vec DE, DI; 295aad13602SShrirang Abhyankar const PetscScalar *de, *di; 296aad13602SShrirang Abhyankar PetscInt zero = 0, one = 1; 297aad13602SShrirang Abhyankar PetscScalar two = 2.0; 298aad13602SShrirang Abhyankar PetscScalar val = 0.0; 299f560b561SHong Zhang Vec Deseq, Diseq; 300f560b561SHong Zhang VecScatter Descat, Discat; 301aad13602SShrirang Abhyankar PetscMPIInt rank; 302aad13602SShrirang Abhyankar MPI_Comm comm; 303aad13602SShrirang Abhyankar 304aad13602SShrirang Abhyankar PetscFunctionBegin; 3059566063dSJacob Faibussowitsch PetscCall(TaoGetDualVariables(tao, &DE, &DI)); 306aad13602SShrirang Abhyankar 3079566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 3089566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 309aad13602SShrirang Abhyankar 3108e58fa1dSresundermann if (!user->noeqflag) { 3119566063dSJacob Faibussowitsch PetscCall(VecScatterCreateToZero(DE, &Descat, &Deseq)); 3129566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD)); 3139566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD)); 3148e58fa1dSresundermann } 3159566063dSJacob Faibussowitsch PetscCall(VecScatterCreateToZero(DI, &Discat, &Diseq)); 3169566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD)); 3179566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD)); 318aad13602SShrirang Abhyankar 319dd400576SPatrick Sanan if (rank == 0) { 3209371c9d4SSatish Balay if (!user->noeqflag) { PetscCall(VecGetArrayRead(Deseq, &de)); /* places equality constraint dual into array */ } 3219566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(Diseq, &di)); /* places inequality constraint dual into array */ 322628da978Sresundermann 3238e58fa1dSresundermann if (!user->noeqflag) { 324628da978Sresundermann val = 2.0 * (1 + de[0] - di[0] + di[1]); 3259566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Deseq, &de)); 3269566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Diseq, &di)); 3278e58fa1dSresundermann } else { 328628da978Sresundermann val = 2.0 * (1 - di[0] + di[1]); 3298e58fa1dSresundermann } 3309566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Diseq, &di)); 3319566063dSJacob Faibussowitsch PetscCall(MatSetValues(H, 1, &zero, 1, &zero, &val, INSERT_VALUES)); 3329566063dSJacob Faibussowitsch PetscCall(MatSetValues(H, 1, &one, 1, &one, &two, INSERT_VALUES)); 333aad13602SShrirang Abhyankar } 3349566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY)); 3359566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY)); 3368e58fa1dSresundermann if (!user->noeqflag) { 3379566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&Descat)); 3389566063dSJacob Faibussowitsch PetscCall(VecDestroy(&Deseq)); 3398e58fa1dSresundermann } 3409566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&Discat)); 3419566063dSJacob Faibussowitsch PetscCall(VecDestroy(&Diseq)); 342aad13602SShrirang Abhyankar PetscFunctionReturn(0); 343aad13602SShrirang Abhyankar } 344aad13602SShrirang Abhyankar 345628da978Sresundermann /* Evaluate 346628da978Sresundermann h = [ x0^2 - x1; 347628da978Sresundermann 1 -(x0^2 - x1)] 348628da978Sresundermann */ 349d71ae5a4SJacob Faibussowitsch PetscErrorCode FormInequalityConstraints(Tao tao, Vec X, Vec CI, void *ctx) 350d71ae5a4SJacob Faibussowitsch { 351aad13602SShrirang Abhyankar const PetscScalar *x; 352aad13602SShrirang Abhyankar PetscScalar ci; 353aad13602SShrirang Abhyankar MPI_Comm comm; 354aad13602SShrirang Abhyankar PetscMPIInt rank; 355aad13602SShrirang Abhyankar AppCtx *user = (AppCtx *)ctx; 356aad13602SShrirang Abhyankar Vec Xseq = user->Xseq; 357aad13602SShrirang Abhyankar VecScatter scat = user->scat; 358aad13602SShrirang Abhyankar 359aad13602SShrirang Abhyankar PetscFunctionBegin; 3609566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 3619566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 362aad13602SShrirang Abhyankar 3639566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 3649566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 365aad13602SShrirang Abhyankar 366dd400576SPatrick Sanan if (rank == 0) { 3679566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(Xseq, &x)); 368aad13602SShrirang Abhyankar ci = x[0] * x[0] - x[1]; 3699566063dSJacob Faibussowitsch PetscCall(VecSetValue(CI, 0, ci, INSERT_VALUES)); 370aad13602SShrirang Abhyankar ci = -x[0] * x[0] + x[1] + 1.0; 3719566063dSJacob Faibussowitsch PetscCall(VecSetValue(CI, 1, ci, INSERT_VALUES)); 3729566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Xseq, &x)); 373aad13602SShrirang Abhyankar } 3749566063dSJacob Faibussowitsch PetscCall(VecAssemblyBegin(CI)); 3759566063dSJacob Faibussowitsch PetscCall(VecAssemblyEnd(CI)); 376aad13602SShrirang Abhyankar PetscFunctionReturn(0); 377aad13602SShrirang Abhyankar } 378aad13602SShrirang Abhyankar 379628da978Sresundermann /* Evaluate 380628da978Sresundermann g = [ x0^2 + x1 - 2] 381628da978Sresundermann */ 382d71ae5a4SJacob Faibussowitsch PetscErrorCode FormEqualityConstraints(Tao tao, Vec X, Vec CE, void *ctx) 383d71ae5a4SJacob Faibussowitsch { 384aad13602SShrirang Abhyankar const PetscScalar *x; 385aad13602SShrirang Abhyankar PetscScalar ce; 386aad13602SShrirang Abhyankar MPI_Comm comm; 387aad13602SShrirang Abhyankar PetscMPIInt rank; 388aad13602SShrirang Abhyankar AppCtx *user = (AppCtx *)ctx; 389aad13602SShrirang Abhyankar Vec Xseq = user->Xseq; 390aad13602SShrirang Abhyankar VecScatter scat = user->scat; 391aad13602SShrirang Abhyankar 392aad13602SShrirang Abhyankar PetscFunctionBegin; 3939566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 3949566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 395aad13602SShrirang Abhyankar 3969566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 3979566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 398aad13602SShrirang Abhyankar 399dd400576SPatrick Sanan if (rank == 0) { 4009566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(Xseq, &x)); 401aad13602SShrirang Abhyankar ce = x[0] * x[0] + x[1] - 2.0; 4029566063dSJacob Faibussowitsch PetscCall(VecSetValue(CE, 0, ce, INSERT_VALUES)); 4039566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Xseq, &x)); 404aad13602SShrirang Abhyankar } 4059566063dSJacob Faibussowitsch PetscCall(VecAssemblyBegin(CE)); 4069566063dSJacob Faibussowitsch PetscCall(VecAssemblyEnd(CE)); 407aad13602SShrirang Abhyankar PetscFunctionReturn(0); 408aad13602SShrirang Abhyankar } 409aad13602SShrirang Abhyankar 410628da978Sresundermann /* 411628da978Sresundermann grad h = [ 2*x0, -1; 412628da978Sresundermann -2*x0, 1] 413628da978Sresundermann */ 414d71ae5a4SJacob Faibussowitsch PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx) 415d71ae5a4SJacob Faibussowitsch { 416aad13602SShrirang Abhyankar AppCtx *user = (AppCtx *)ctx; 417f560b561SHong Zhang PetscInt zero = 0, one = 1, cols[2]; 418aad13602SShrirang Abhyankar PetscScalar vals[2]; 419aad13602SShrirang Abhyankar const PetscScalar *x; 420aad13602SShrirang Abhyankar Vec Xseq = user->Xseq; 421aad13602SShrirang Abhyankar VecScatter scat = user->scat; 422aad13602SShrirang Abhyankar MPI_Comm comm; 423aad13602SShrirang Abhyankar PetscMPIInt rank; 424aad13602SShrirang Abhyankar 425aad13602SShrirang Abhyankar PetscFunctionBegin; 4269566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 4279566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 4289566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 4299566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 430aad13602SShrirang Abhyankar 4319566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(Xseq, &x)); 432c5853193SPierre Jolivet if (rank == 0) { 4339371c9d4SSatish Balay cols[0] = 0; 4349371c9d4SSatish Balay cols[1] = 1; 4359371c9d4SSatish Balay vals[0] = 2 * x[0]; 4369371c9d4SSatish Balay vals[1] = -1.0; 4379566063dSJacob Faibussowitsch PetscCall(MatSetValues(JI, 1, &zero, 2, cols, vals, INSERT_VALUES)); 4389371c9d4SSatish Balay vals[0] = -2 * x[0]; 4399371c9d4SSatish Balay vals[1] = 1.0; 4409566063dSJacob Faibussowitsch PetscCall(MatSetValues(JI, 1, &one, 2, cols, vals, INSERT_VALUES)); 441aad13602SShrirang Abhyankar } 4429566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Xseq, &x)); 4439566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(JI, MAT_FINAL_ASSEMBLY)); 4449566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(JI, MAT_FINAL_ASSEMBLY)); 445aad13602SShrirang Abhyankar PetscFunctionReturn(0); 446aad13602SShrirang Abhyankar } 447aad13602SShrirang Abhyankar 448628da978Sresundermann /* 449628da978Sresundermann grad g = [2*x0 450628da978Sresundermann 1.0 ] 451628da978Sresundermann */ 452d71ae5a4SJacob Faibussowitsch PetscErrorCode FormEqualityJacobian(Tao tao, Vec X, Mat JE, Mat JEpre, void *ctx) 453d71ae5a4SJacob Faibussowitsch { 454f560b561SHong Zhang PetscInt zero = 0, cols[2]; 455aad13602SShrirang Abhyankar PetscScalar vals[2]; 456aad13602SShrirang Abhyankar const PetscScalar *x; 457aad13602SShrirang Abhyankar PetscMPIInt rank; 458aad13602SShrirang Abhyankar MPI_Comm comm; 459aad13602SShrirang Abhyankar 460aad13602SShrirang Abhyankar PetscFunctionBegin; 4619566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 4629566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 463aad13602SShrirang Abhyankar 464dd400576SPatrick Sanan if (rank == 0) { 4659566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(X, &x)); 4669371c9d4SSatish Balay cols[0] = 0; 4679371c9d4SSatish Balay cols[1] = 1; 4689371c9d4SSatish Balay vals[0] = 2 * x[0]; 4699371c9d4SSatish Balay vals[1] = 1.0; 4709566063dSJacob Faibussowitsch PetscCall(MatSetValues(JE, 1, &zero, 2, cols, vals, INSERT_VALUES)); 4719566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(X, &x)); 472aad13602SShrirang Abhyankar } 4739566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(JE, MAT_FINAL_ASSEMBLY)); 4749566063dSJacob Faibussowitsch PetscCall(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: 4849a09327dSAlp Dener args: -tao_converged_reason -tao_gatol 1.e-6 -tao_type pdipm -tao_pdipm_kkt_shift_pd 485f560b561SHong Zhang requires: mumps 4869a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 487aad13602SShrirang Abhyankar 488aad13602SShrirang Abhyankar test: 489aad13602SShrirang Abhyankar suffix: 2 4909a09327dSAlp Dener args: -tao_converged_reason 4919a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 492aad13602SShrirang Abhyankar 4938e58fa1dSresundermann test: 4948e58fa1dSresundermann suffix: 3 4958e58fa1dSresundermann args: -tao_converged_reason -no_eq 4969a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 4978e58fa1dSresundermann 4988e58fa1dSresundermann test: 4998e58fa1dSresundermann suffix: 4 5009a09327dSAlp Dener args: -tao_converged_reason -tao_almm_type classic 5019a09327dSAlp Dener requires: !single 5029a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 5038e58fa1dSresundermann 504661095bbSAlp Dener test: 505661095bbSAlp Dener suffix: 5 5069a09327dSAlp Dener args: -tao_converged_reason -tao_almm_type classic -no_eq 507*b9e5a4c7SSatish Balay requires: !single !defined(PETSCTEST_VALGRIND) 5089a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 509661095bbSAlp Dener 510661095bbSAlp Dener test: 511661095bbSAlp Dener suffix: 6 5129a09327dSAlp Dener args: -tao_converged_reason -tao_almm_subsolver_tao_type bqnktr 5139a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 514661095bbSAlp Dener 515661095bbSAlp Dener test: 516661095bbSAlp Dener suffix: 7 5179a09327dSAlp Dener args: -tao_converged_reason -tao_almm_subsolver_tao_type bncg 5189a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 519661095bbSAlp Dener 520661095bbSAlp Dener test: 521661095bbSAlp Dener suffix: 8 522661095bbSAlp Dener nsize: 2 5239a09327dSAlp Dener args: -tao_converged_reason 5249a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 525661095bbSAlp Dener 526661095bbSAlp Dener test: 527661095bbSAlp Dener suffix: 9 528661095bbSAlp Dener nsize: 2 5299a09327dSAlp Dener args: -tao_converged_reason -vec_type cuda -mat_type aijcusparse 5309a09327dSAlp Dener requires: cuda 5319a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 532661095bbSAlp Dener 533aad13602SShrirang Abhyankar TEST*/ 534