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