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\ 2710978b7dSBarry Smith -tao_monitor_constraint_norm : convergence monitor with constraint norm \n\ 2810978b7dSBarry Smith -tao_monitor_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 *); 4983f04556SHong Zhang PetscErrorCode FormPDIPMHessian(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 553ba16761SJacob Faibussowitsch int 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)); 10583f04556SHong Zhang PetscCall(PCSetType(pc, PCCHOLESKY)); 1069a09327dSAlp Dener PetscCall(KSPSetFromOptions(ksp)); 10783f04556SHong Zhang PetscCall(TaoSetHessian(tao, user.H, user.H, FormPDIPMHessian, (void *)&user)); 1089a09327dSAlp Dener } 1099a09327dSAlp Dener 1109a09327dSAlp Dener /* Print out an initial view of the problem */ 111f560b561SHong Zhang if (user.initview) { 1129566063dSJacob Faibussowitsch PetscCall(TaoSetUp(tao)); 1139566063dSJacob Faibussowitsch PetscCall(VecDuplicate(user.x, &G)); 1149566063dSJacob Faibussowitsch PetscCall(FormFunctionGradient(tao, user.x, &f, G, (void *)&user)); 11583f04556SHong Zhang PetscCall(FormPDIPMHessian(tao, user.x, user.H, user.H, (void *)&user)); 1169566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD)); 11763a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial point X:\n")); 1189566063dSJacob Faibussowitsch PetscCall(VecView(user.x, PETSC_VIEWER_STDOUT_WORLD)); 11963a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial objective f(x) = %g\n", (double)f)); 12063a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial gradient and Hessian:\n")); 1219566063dSJacob Faibussowitsch PetscCall(VecView(G, PETSC_VIEWER_STDOUT_WORLD)); 1229566063dSJacob Faibussowitsch PetscCall(MatView(user.H, PETSC_VIEWER_STDOUT_WORLD)); 1239566063dSJacob Faibussowitsch PetscCall(VecDestroy(&G)); 1249566063dSJacob Faibussowitsch PetscCall(FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void *)&user)); 1259566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(user.Ai, NULL, &CI)); 1269566063dSJacob Faibussowitsch PetscCall(FormInequalityConstraints(tao, user.x, CI, (void *)&user)); 12763a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial inequality constraints and Jacobian:\n")); 1289566063dSJacob Faibussowitsch PetscCall(VecView(CI, PETSC_VIEWER_STDOUT_WORLD)); 1299566063dSJacob Faibussowitsch PetscCall(MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD)); 1309566063dSJacob Faibussowitsch PetscCall(VecDestroy(&CI)); 131f560b561SHong Zhang if (!user.noeqflag) { 1329566063dSJacob Faibussowitsch PetscCall(FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void *)&user)); 1339566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(user.Ae, NULL, &CE)); 1349566063dSJacob Faibussowitsch PetscCall(FormEqualityConstraints(tao, user.x, CE, (void *)&user)); 13563a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial equality constraints and Jacobian:\n")); 1369566063dSJacob Faibussowitsch PetscCall(VecView(CE, PETSC_VIEWER_STDOUT_WORLD)); 1379566063dSJacob Faibussowitsch PetscCall(MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD)); 1389566063dSJacob Faibussowitsch PetscCall(VecDestroy(&CE)); 139f560b561SHong Zhang } 1409566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n")); 1419566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD)); 142f560b561SHong Zhang } 143aad13602SShrirang Abhyankar 1449566063dSJacob Faibussowitsch PetscCall(TaoSolve(tao)); 1459566063dSJacob Faibussowitsch PetscCall(TaoGetSolution(tao, &x)); 1469a09327dSAlp Dener PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Found solution:\n")); 1479566063dSJacob Faibussowitsch PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD)); 148aad13602SShrirang Abhyankar 149aad13602SShrirang Abhyankar /* Free objects */ 1509566063dSJacob Faibussowitsch PetscCall(DestroyProblem(&user)); 1519566063dSJacob Faibussowitsch PetscCall(TaoDestroy(&tao)); 1529566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 1533ba16761SJacob Faibussowitsch return PETSC_SUCCESS; 154aad13602SShrirang Abhyankar } 155aad13602SShrirang Abhyankar 156d71ae5a4SJacob Faibussowitsch PetscErrorCode InitializeProblem(AppCtx *user) 157d71ae5a4SJacob Faibussowitsch { 158aad13602SShrirang Abhyankar PetscMPIInt size; 159aad13602SShrirang Abhyankar PetscMPIInt rank; 160aad13602SShrirang Abhyankar PetscInt nloc, neloc, niloc; 161aad13602SShrirang Abhyankar 162aad13602SShrirang Abhyankar PetscFunctionBegin; 1639566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 1649566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 1658e58fa1dSresundermann user->noeqflag = PETSC_FALSE; 1669566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_eq", &user->noeqflag, NULL)); 167f560b561SHong Zhang user->initview = PETSC_FALSE; 1689566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL, NULL, "-init_view", &user->initview, NULL)); 169628da978Sresundermann 170f560b561SHong Zhang /* Tell user the correct solution, not an error checking */ 1719a09327dSAlp Dener if (!user->noeqflag) { 1729a09327dSAlp Dener PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution: f(1, 1) = -2\n")); 1739a09327dSAlp Dener } else { 1749a09327dSAlp Dener PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution (-no_eq): f(1.73205, 2) = -7.3923\n")); 1758e58fa1dSresundermann } 176aad13602SShrirang Abhyankar 1772d4ee042Sprj- /* create vector x and set initial values */ 178aad13602SShrirang Abhyankar user->n = 2; /* global length */ 179f560b561SHong Zhang nloc = (size == 1) ? user->n : 1; 1809566063dSJacob Faibussowitsch PetscCall(VecCreate(PETSC_COMM_WORLD, &user->x)); 1819566063dSJacob Faibussowitsch PetscCall(VecSetSizes(user->x, nloc, user->n)); 1829566063dSJacob Faibussowitsch PetscCall(VecSetFromOptions(user->x)); 1839a09327dSAlp Dener PetscCall(VecSet(user->x, 0.0)); 184aad13602SShrirang Abhyankar 185aad13602SShrirang Abhyankar /* create and set lower and upper bound vectors */ 1869566063dSJacob Faibussowitsch PetscCall(VecDuplicate(user->x, &user->xl)); 1879566063dSJacob Faibussowitsch PetscCall(VecDuplicate(user->x, &user->xu)); 1889566063dSJacob Faibussowitsch PetscCall(VecSet(user->xl, -1.0)); 1899566063dSJacob Faibussowitsch PetscCall(VecSet(user->xu, 2.0)); 190aad13602SShrirang Abhyankar 191aad13602SShrirang Abhyankar /* create scater to zero */ 1929566063dSJacob Faibussowitsch PetscCall(VecScatterCreateToZero(user->x, &user->scat, &user->Xseq)); 193aad13602SShrirang Abhyankar 194aad13602SShrirang Abhyankar user->ne = 1; 195628da978Sresundermann user->ni = 2; 196aad13602SShrirang Abhyankar neloc = (rank == 0) ? user->ne : 0; 197f560b561SHong Zhang niloc = (size == 1) ? user->ni : 1; 198aad13602SShrirang Abhyankar 1998e58fa1dSresundermann if (!user->noeqflag) { 2009566063dSJacob Faibussowitsch PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ce)); /* a 1x1 vec for equality constraints */ 2019566063dSJacob Faibussowitsch PetscCall(VecSetSizes(user->ce, neloc, user->ne)); 2029566063dSJacob Faibussowitsch PetscCall(VecSetFromOptions(user->ce)); 2039566063dSJacob Faibussowitsch PetscCall(VecSetUp(user->ce)); 2048e58fa1dSresundermann } 205628da978Sresundermann 2069566063dSJacob Faibussowitsch PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ci)); /* a 2x1 vec for inequality constraints */ 2079566063dSJacob Faibussowitsch PetscCall(VecSetSizes(user->ci, niloc, user->ni)); 2089566063dSJacob Faibussowitsch PetscCall(VecSetFromOptions(user->ci)); 2099566063dSJacob Faibussowitsch PetscCall(VecSetUp(user->ci)); 210aad13602SShrirang Abhyankar 211baca6076SPierre Jolivet /* nexn & nixn matrices for equally and inequalty constraints */ 2128e58fa1dSresundermann if (!user->noeqflag) { 2139566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ae)); 2149566063dSJacob Faibussowitsch PetscCall(MatSetSizes(user->Ae, neloc, nloc, user->ne, user->n)); 2159566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(user->Ae)); 2169566063dSJacob Faibussowitsch PetscCall(MatSetUp(user->Ae)); 2178e58fa1dSresundermann } 2188e58fa1dSresundermann 2199566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ai)); 2209566063dSJacob Faibussowitsch PetscCall(MatSetSizes(user->Ai, niloc, nloc, user->ni, user->n)); 2219566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(user->Ai)); 2229566063dSJacob Faibussowitsch PetscCall(MatSetUp(user->Ai)); 223628da978Sresundermann 2249566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD, &user->H)); 2259566063dSJacob Faibussowitsch PetscCall(MatSetSizes(user->H, nloc, nloc, user->n, user->n)); 2269566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(user->H)); 2279566063dSJacob Faibussowitsch PetscCall(MatSetUp(user->H)); 2283ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 229aad13602SShrirang Abhyankar } 230aad13602SShrirang Abhyankar 231d71ae5a4SJacob Faibussowitsch PetscErrorCode DestroyProblem(AppCtx *user) 232d71ae5a4SJacob Faibussowitsch { 233aad13602SShrirang Abhyankar PetscFunctionBegin; 23448a46eb9SPierre Jolivet if (!user->noeqflag) PetscCall(MatDestroy(&user->Ae)); 2359566063dSJacob Faibussowitsch PetscCall(MatDestroy(&user->Ai)); 2369566063dSJacob Faibussowitsch PetscCall(MatDestroy(&user->H)); 237aad13602SShrirang Abhyankar 2389566063dSJacob Faibussowitsch PetscCall(VecDestroy(&user->x)); 23948a46eb9SPierre Jolivet if (!user->noeqflag) PetscCall(VecDestroy(&user->ce)); 2409566063dSJacob Faibussowitsch PetscCall(VecDestroy(&user->ci)); 2419566063dSJacob Faibussowitsch PetscCall(VecDestroy(&user->xl)); 2429566063dSJacob Faibussowitsch PetscCall(VecDestroy(&user->xu)); 2439566063dSJacob Faibussowitsch PetscCall(VecDestroy(&user->Xseq)); 2449566063dSJacob Faibussowitsch PetscCall(VecScatterDestroy(&user->scat)); 2453ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 246aad13602SShrirang Abhyankar } 247aad13602SShrirang Abhyankar 248628da978Sresundermann /* Evaluate 249628da978Sresundermann f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1) 25083f04556SHong Zhang G = grad f = [2*(x0 - 2) - 2, 2*(x1 - 2) - 2] = 2*X - 6 251aad13602SShrirang Abhyankar */ 252d71ae5a4SJacob Faibussowitsch PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx) 253d71ae5a4SJacob Faibussowitsch { 254aad13602SShrirang Abhyankar const PetscScalar *x; 255aad13602SShrirang Abhyankar MPI_Comm comm; 256aad13602SShrirang Abhyankar PetscMPIInt rank; 257aad13602SShrirang Abhyankar PetscReal fin; 258aad13602SShrirang Abhyankar AppCtx *user = (AppCtx *)ctx; 259aad13602SShrirang Abhyankar Vec Xseq = user->Xseq; 260aad13602SShrirang Abhyankar VecScatter scat = user->scat; 261aad13602SShrirang Abhyankar 262aad13602SShrirang Abhyankar PetscFunctionBegin; 26383f04556SHong Zhang /* f = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1) */ 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]); 2749566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Xseq, &x)); 275aad13602SShrirang Abhyankar } 276712fec58SPierre Jolivet PetscCall(MPIU_Allreduce(&fin, f, 1, MPIU_REAL, MPIU_SUM, comm)); 27783f04556SHong Zhang 27883f04556SHong Zhang /* G = 2*X - 6 */ 27983f04556SHong Zhang PetscCall(VecSet(G, -6.0)); 28083f04556SHong Zhang PetscCall(VecAXPY(G, 2.0, X)); 2813ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 282aad13602SShrirang Abhyankar } 283aad13602SShrirang Abhyankar 28483f04556SHong Zhang /* Evaluate PDIPM Hessian, see Eqn(22) in http://doi.org/10.1049/gtd2.12708 28583f04556SHong Zhang H = Wxx = fxx + Jacobian(grad g^T*DE) - Jacobian(grad h^T*DI)] 28683f04556SHong Zhang = fxx + Jacobin([2*x0; 1]de[0]) - Jacobian([2*x0, -2*x0; -1, 1][di[0] di[1]]^T) 28783f04556SHong Zhang = [2 0; 0 2] + [2*de[0] 0; 0 0] - [2*di[0]-2*di[1], 0; 0, 0] 288628da978Sresundermann = [ 2*(1+de[0]-di[0]+di[1]), 0; 289628da978Sresundermann 0, 2] 290628da978Sresundermann */ 29183f04556SHong Zhang PetscErrorCode FormPDIPMHessian(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)); 3423ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 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)); 3763ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 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)); 4073ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 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)); 4453ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 446aad13602SShrirang Abhyankar } 447aad13602SShrirang Abhyankar 448628da978Sresundermann /* 44983f04556SHong Zhang grad g = [2*x0, 1.0] 450628da978Sresundermann */ 451d71ae5a4SJacob Faibussowitsch PetscErrorCode FormEqualityJacobian(Tao tao, Vec X, Mat JE, Mat JEpre, void *ctx) 452d71ae5a4SJacob Faibussowitsch { 453f560b561SHong Zhang PetscInt zero = 0, cols[2]; 454aad13602SShrirang Abhyankar PetscScalar vals[2]; 455aad13602SShrirang Abhyankar const PetscScalar *x; 456aad13602SShrirang Abhyankar PetscMPIInt rank; 457aad13602SShrirang Abhyankar MPI_Comm comm; 458aad13602SShrirang Abhyankar 459aad13602SShrirang Abhyankar PetscFunctionBegin; 4609566063dSJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 4619566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 462aad13602SShrirang Abhyankar 463dd400576SPatrick Sanan if (rank == 0) { 4649566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(X, &x)); 4659371c9d4SSatish Balay cols[0] = 0; 4669371c9d4SSatish Balay cols[1] = 1; 4679371c9d4SSatish Balay vals[0] = 2 * x[0]; 4689371c9d4SSatish Balay vals[1] = 1.0; 4699566063dSJacob Faibussowitsch PetscCall(MatSetValues(JE, 1, &zero, 2, cols, vals, INSERT_VALUES)); 4709566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(X, &x)); 471aad13602SShrirang Abhyankar } 4729566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(JE, MAT_FINAL_ASSEMBLY)); 4739566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(JE, MAT_FINAL_ASSEMBLY)); 4743ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 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: 4839a09327dSAlp Dener args: -tao_converged_reason -tao_gatol 1.e-6 -tao_type pdipm -tao_pdipm_kkt_shift_pd 484f560b561SHong Zhang requires: mumps 4859a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 486aad13602SShrirang Abhyankar 487aad13602SShrirang Abhyankar test: 48883f04556SHong Zhang suffix: pdipm_2 48983f04556SHong Zhang requires: mumps 49083f04556SHong Zhang nsize: 2 49183f04556SHong Zhang args: -tao_converged_reason -tao_gatol 1.e-6 -tao_type pdipm 49283f04556SHong Zhang filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 49383f04556SHong Zhang 49483f04556SHong Zhang test: 495aad13602SShrirang Abhyankar suffix: 2 4969a09327dSAlp Dener args: -tao_converged_reason 4979a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 498aad13602SShrirang Abhyankar 4998e58fa1dSresundermann test: 5008e58fa1dSresundermann suffix: 3 5018e58fa1dSresundermann args: -tao_converged_reason -no_eq 5029a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 5038e58fa1dSresundermann 5048e58fa1dSresundermann test: 5058e58fa1dSresundermann suffix: 4 5069a09327dSAlp Dener args: -tao_converged_reason -tao_almm_type classic 5079a09327dSAlp Dener requires: !single 5089a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 5098e58fa1dSresundermann 510661095bbSAlp Dener test: 511661095bbSAlp Dener suffix: 5 512*910a5171SToby Isaac args: -tao_converged_reason -tao_almm_type classic -no_eq -tao_almm_subsolver_tao_max_it 100 513b9e5a4c7SSatish Balay requires: !single !defined(PETSCTEST_VALGRIND) 5149a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 515661095bbSAlp Dener 516661095bbSAlp Dener test: 517661095bbSAlp Dener suffix: 6 5189a09327dSAlp Dener args: -tao_converged_reason -tao_almm_subsolver_tao_type bqnktr 5199a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 520661095bbSAlp Dener 521661095bbSAlp Dener test: 522661095bbSAlp Dener suffix: 7 5239a09327dSAlp Dener args: -tao_converged_reason -tao_almm_subsolver_tao_type bncg 5249a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 525661095bbSAlp Dener 526661095bbSAlp Dener test: 527661095bbSAlp Dener suffix: 8 528661095bbSAlp Dener nsize: 2 5299a09327dSAlp Dener args: -tao_converged_reason 5309a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 531661095bbSAlp Dener 532661095bbSAlp Dener test: 533661095bbSAlp Dener suffix: 9 534661095bbSAlp Dener nsize: 2 5359a09327dSAlp Dener args: -tao_converged_reason -vec_type cuda -mat_type aijcusparse 5369a09327dSAlp Dener requires: cuda 5379a09327dSAlp Dener filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 538661095bbSAlp Dener 539aad13602SShrirang Abhyankar TEST*/ 540