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 PetscErrorCode ierr; /* used to check for functions returning nonzeros */ 58aad13602SShrirang Abhyankar Tao tao; 59aad13602SShrirang Abhyankar KSP ksp; 60aad13602SShrirang Abhyankar PC pc; 61aad13602SShrirang Abhyankar AppCtx user; /* application context */ 62f560b561SHong Zhang Vec x,G,CI,CE; 63f560b561SHong Zhang PetscMPIInt size; 64661095bbSAlp Dener TaoType type; 65f560b561SHong Zhang PetscReal f; 66661095bbSAlp Dener PetscBool pdipm; 67aad13602SShrirang Abhyankar 68aad13602SShrirang Abhyankar ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 69*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 703c859ba3SBarry 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."); 71aad13602SShrirang Abhyankar 72*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"---- Constrained Problem -----\n")); 73*5f80ce2aSJacob Faibussowitsch CHKERRQ(InitializeProblem(&user)); /* sets up problem, function below */ 74*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoCreate(PETSC_COMM_WORLD,&tao)); 75*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoSetType(tao,TAOPDIPM)); 76*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoSetSolution(tao,user.x)); /* gets solution vector from problem */ 77*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoSetVariableBounds(tao,user.xl,user.xu)); /* sets lower upper bounds from given solution */ 78*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void*)&user)); 79aad13602SShrirang Abhyankar 808e58fa1dSresundermann if (!user.noeqflag) { 81*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user)); 828e58fa1dSresundermann } 83*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user)); 848e58fa1dSresundermann if (!user.noeqflag) { 85*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user)); /* equality jacobian */ 868e58fa1dSresundermann } 87*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user)); /* inequality jacobian */ 88*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoSetTolerances(tao,1.e-6,1.e-6,1.e-6)); 89*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoSetConstraintTolerances(tao,1.e-6,1.e-6)); 90aad13602SShrirang Abhyankar 91*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoGetKSP(tao,&ksp)); 92*5f80ce2aSJacob Faibussowitsch CHKERRQ(KSPGetPC(ksp,&pc)); 93*5f80ce2aSJacob Faibussowitsch CHKERRQ(PCSetType(pc,PCCHOLESKY)); 94aad13602SShrirang Abhyankar /* 95aad13602SShrirang Abhyankar This algorithm produces matrices with zeros along the diagonal therefore we use 9612d688e0SRylee Sundermann MUMPS which provides solver for indefinite matrices 97aad13602SShrirang Abhyankar */ 98f560b561SHong Zhang #if defined(PETSC_HAVE_MUMPS) 99*5f80ce2aSJacob Faibussowitsch CHKERRQ(PCFactorSetMatSolverType(pc,MATSOLVERMUMPS)); /* requires mumps to solve pdipm */ 100f560b561SHong Zhang #else 1013c859ba3SBarry Smith PetscCheck(size == 1,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Requires an external package that supports parallel PCCHOLESKY, e.g., MUMPS."); 102f560b561SHong Zhang #endif 103*5f80ce2aSJacob Faibussowitsch CHKERRQ(KSPSetType(ksp,KSPPREONLY)); 104*5f80ce2aSJacob Faibussowitsch CHKERRQ(KSPSetFromOptions(ksp)); 105aad13602SShrirang Abhyankar 106*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoSetFromOptions(tao)); 107*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoGetType(tao,&type)); 108*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectTypeCompare((PetscObject)tao,TAOPDIPM,&pdipm)); 109661095bbSAlp Dener if (pdipm) { 110*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoSetHessian(tao,user.H,user.H,FormHessian,(void*)&user)); 111f560b561SHong Zhang if (user.initview) { 112*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoSetUp(tao)); 113*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDuplicate(user.x, &G)); 114*5f80ce2aSJacob Faibussowitsch CHKERRQ(FormFunctionGradient(tao, user.x, &f, G, (void*)&user)); 115*5f80ce2aSJacob Faibussowitsch CHKERRQ(FormHessian(tao, user.x, user.H, user.H, (void*)&user)); 116*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD)); 117*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nInitial point X:\n",f)); 118*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecView(user.x, PETSC_VIEWER_STDOUT_WORLD)); 119*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nInitial objective f(x) = %g\n",f)); 120*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nInitial gradient and Hessian:\n",f)); 121*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecView(G, PETSC_VIEWER_STDOUT_WORLD)); 122*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatView(user.H, PETSC_VIEWER_STDOUT_WORLD)); 123*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&G)); 124*5f80ce2aSJacob Faibussowitsch CHKERRQ(FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void*)&user)); 125*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreateVecs(user.Ai, NULL, &CI)); 126*5f80ce2aSJacob Faibussowitsch CHKERRQ(FormInequalityConstraints(tao, user.x, CI, (void*)&user)); 127*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nInitial inequality constraints and Jacobian:\n",f)); 128*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecView(CI, PETSC_VIEWER_STDOUT_WORLD)); 129*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD)); 130*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&CI)); 131f560b561SHong Zhang if (!user.noeqflag) { 132*5f80ce2aSJacob Faibussowitsch CHKERRQ(FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void*)&user)); 133*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreateVecs(user.Ae, NULL, &CE)); 134*5f80ce2aSJacob Faibussowitsch CHKERRQ(FormEqualityConstraints(tao, user.x, CE, (void*)&user)); 135*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nInitial equality constraints and Jacobian:\n",f)); 136*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecView(CE, PETSC_VIEWER_STDOUT_WORLD)); 137*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD)); 138*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&CE)); 139f560b561SHong Zhang } 140*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "\n")); 141*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD)); 142f560b561SHong Zhang } 143661095bbSAlp Dener } 144aad13602SShrirang Abhyankar 145*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoSolve(tao)); 146*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoGetSolution(tao,&x)); 147*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); 148aad13602SShrirang Abhyankar 149aad13602SShrirang Abhyankar /* Free objects */ 150*5f80ce2aSJacob Faibussowitsch CHKERRQ(DestroyProblem(&user)); 151*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoDestroy(&tao)); 152aad13602SShrirang Abhyankar ierr = PetscFinalize(); 153aad13602SShrirang Abhyankar return ierr; 154aad13602SShrirang Abhyankar } 155aad13602SShrirang Abhyankar 156aad13602SShrirang Abhyankar PetscErrorCode InitializeProblem(AppCtx *user) 157aad13602SShrirang Abhyankar { 158aad13602SShrirang Abhyankar PetscMPIInt size; 159aad13602SShrirang Abhyankar PetscMPIInt rank; 160aad13602SShrirang Abhyankar PetscInt nloc,neloc,niloc; 161aad13602SShrirang Abhyankar 162aad13602SShrirang Abhyankar PetscFunctionBegin; 163*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 164*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 1658e58fa1dSresundermann user->noeqflag = PETSC_FALSE; 166*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-no_eq",&user->noeqflag,NULL)); 167f560b561SHong Zhang user->initview = PETSC_FALSE; 168*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-init_view",&user->initview,NULL)); 169628da978Sresundermann 1708e58fa1dSresundermann if (!user->noeqflag) { 171f560b561SHong Zhang /* Tell user the correct solution, not an error checking */ 172*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n")); 1738e58fa1dSresundermann } 174aad13602SShrirang Abhyankar 1752d4ee042Sprj- /* create vector x and set initial values */ 176aad13602SShrirang Abhyankar user->n = 2; /* global length */ 177f560b561SHong Zhang nloc = (size==1)?user->n:1; 178*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user->x)); 179*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSetSizes(user->x,nloc,user->n)); 180*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSetFromOptions(user->x)); 181*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSet(user->x,0)); 182aad13602SShrirang Abhyankar 183aad13602SShrirang Abhyankar /* create and set lower and upper bound vectors */ 184*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDuplicate(user->x,&user->xl)); 185*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDuplicate(user->x,&user->xu)); 186*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSet(user->xl,-1.0)); 187*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSet(user->xu,2.0)); 188aad13602SShrirang Abhyankar 189aad13602SShrirang Abhyankar /* create scater to zero */ 190*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterCreateToZero(user->x,&user->scat,&user->Xseq)); 191aad13602SShrirang Abhyankar 192aad13602SShrirang Abhyankar user->ne = 1; 193628da978Sresundermann user->ni = 2; 194aad13602SShrirang Abhyankar neloc = (rank==0)?user->ne:0; 195f560b561SHong Zhang niloc = (size==1)?user->ni:1; 196aad13602SShrirang Abhyankar 1978e58fa1dSresundermann if (!user->noeqflag) { 198*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user->ce)); /* a 1x1 vec for equality constraints */ 199*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSetSizes(user->ce,neloc,user->ne)); 200*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSetFromOptions(user->ce)); 201*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSetUp(user->ce)); 2028e58fa1dSresundermann } 203628da978Sresundermann 204*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user->ci)); /* a 2x1 vec for inequality constraints */ 205*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSetSizes(user->ci,niloc,user->ni)); 206*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSetFromOptions(user->ci)); 207*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSetUp(user->ci)); 208aad13602SShrirang Abhyankar 209a5b23f4aSJose E. Roman /* nexn & nixn matricies for equally and inequalty constraints */ 2108e58fa1dSresundermann if (!user->noeqflag) { 211*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user->Ae)); 212*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n)); 213*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetFromOptions(user->Ae)); 214*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetUp(user->Ae)); 2158e58fa1dSresundermann } 2168e58fa1dSresundermann 217*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user->Ai)); 218*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n)); 219*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetFromOptions(user->Ai)); 220*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetUp(user->Ai)); 221628da978Sresundermann 222*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user->H)); 223*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetSizes(user->H,nloc,nloc,user->n,user->n)); 224*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetFromOptions(user->H)); 225*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetUp(user->H)); 226aad13602SShrirang Abhyankar PetscFunctionReturn(0); 227aad13602SShrirang Abhyankar } 228aad13602SShrirang Abhyankar 229aad13602SShrirang Abhyankar PetscErrorCode DestroyProblem(AppCtx *user) 230aad13602SShrirang Abhyankar { 231aad13602SShrirang Abhyankar PetscFunctionBegin; 2328e58fa1dSresundermann if (!user->noeqflag) { 233*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&user->Ae)); 2348e58fa1dSresundermann } 235*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&user->Ai)); 236*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&user->H)); 237aad13602SShrirang Abhyankar 238*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&user->x)); 2398e58fa1dSresundermann if (!user->noeqflag) { 240*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&user->ce)); 2418e58fa1dSresundermann } 242*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&user->ci)); 243*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&user->xl)); 244*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&user->xu)); 245*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&user->Xseq)); 246*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterDestroy(&user->scat)); 247aad13602SShrirang Abhyankar PetscFunctionReturn(0); 248aad13602SShrirang Abhyankar } 249aad13602SShrirang Abhyankar 250628da978Sresundermann /* Evaluate 251628da978Sresundermann f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1) 252628da978Sresundermann G = grad f = [2*(x0 - 2) - 2; 253628da978Sresundermann 2*(x1 - 2) - 2] 254aad13602SShrirang Abhyankar */ 255aad13602SShrirang Abhyankar PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx) 256aad13602SShrirang Abhyankar { 257aad13602SShrirang Abhyankar PetscScalar g; 258aad13602SShrirang Abhyankar const PetscScalar *x; 259aad13602SShrirang Abhyankar MPI_Comm comm; 260aad13602SShrirang Abhyankar PetscMPIInt rank; 261aad13602SShrirang Abhyankar PetscReal fin; 262aad13602SShrirang Abhyankar AppCtx *user=(AppCtx*)ctx; 263aad13602SShrirang Abhyankar Vec Xseq=user->Xseq; 264aad13602SShrirang Abhyankar VecScatter scat=user->scat; 265aad13602SShrirang Abhyankar 266aad13602SShrirang Abhyankar PetscFunctionBegin; 267*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm)); 268*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_rank(comm,&rank)); 269aad13602SShrirang Abhyankar 270*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD)); 271*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD)); 272aad13602SShrirang Abhyankar 273aad13602SShrirang Abhyankar fin = 0.0; 274dd400576SPatrick Sanan if (rank == 0) { 275*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayRead(Xseq,&x)); 276aad13602SShrirang Abhyankar fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]); 277aad13602SShrirang Abhyankar g = 2.0*(x[0]-2.0) - 2.0; 278*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSetValue(G,0,g,INSERT_VALUES)); 279aad13602SShrirang Abhyankar g = 2.0*(x[1]-2.0) - 2.0; 280*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSetValue(G,1,g,INSERT_VALUES)); 281*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(Xseq,&x)); 282aad13602SShrirang Abhyankar } 283*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm)); 284*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecAssemblyBegin(G)); 285*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecAssemblyEnd(G)); 286aad13602SShrirang Abhyankar PetscFunctionReturn(0); 287aad13602SShrirang Abhyankar } 288aad13602SShrirang Abhyankar 289628da978Sresundermann /* Evaluate 290628da978Sresundermann H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)] 291628da978Sresundermann = [ 2*(1+de[0]-di[0]+di[1]), 0; 292628da978Sresundermann 0, 2] 293628da978Sresundermann */ 294aad13602SShrirang Abhyankar PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx) 295aad13602SShrirang Abhyankar { 296aad13602SShrirang Abhyankar AppCtx *user=(AppCtx*)ctx; 297aad13602SShrirang Abhyankar Vec DE,DI; 298aad13602SShrirang Abhyankar const PetscScalar *de,*di; 299aad13602SShrirang Abhyankar PetscInt zero=0,one=1; 300aad13602SShrirang Abhyankar PetscScalar two=2.0; 301aad13602SShrirang Abhyankar PetscScalar val=0.0; 302f560b561SHong Zhang Vec Deseq,Diseq; 303f560b561SHong Zhang VecScatter Descat,Discat; 304aad13602SShrirang Abhyankar PetscMPIInt rank; 305aad13602SShrirang Abhyankar MPI_Comm comm; 306aad13602SShrirang Abhyankar 307aad13602SShrirang Abhyankar PetscFunctionBegin; 308*5f80ce2aSJacob Faibussowitsch CHKERRQ(TaoGetDualVariables(tao,&DE,&DI)); 309aad13602SShrirang Abhyankar 310*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm)); 311*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_rank(comm,&rank)); 312aad13602SShrirang Abhyankar 3138e58fa1dSresundermann if (!user->noeqflag) { 314*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterCreateToZero(DE,&Descat,&Deseq)); 315*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD)); 316*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD)); 3178e58fa1dSresundermann } 318*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterCreateToZero(DI,&Discat,&Diseq)); 319*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD)); 320*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD)); 321aad13602SShrirang Abhyankar 322dd400576SPatrick Sanan if (rank == 0) { 3238e58fa1dSresundermann if (!user->noeqflag) { 324*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayRead(Deseq,&de)); /* places equality constraint dual into array */ 3258e58fa1dSresundermann } 326*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayRead(Diseq,&di)); /* places inequality constraint dual into array */ 327628da978Sresundermann 3288e58fa1dSresundermann if (!user->noeqflag) { 329628da978Sresundermann val = 2.0 * (1 + de[0] - di[0] + di[1]); 330*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(Deseq,&de)); 331*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(Diseq,&di)); 3328e58fa1dSresundermann } else { 333628da978Sresundermann val = 2.0 * (1 - di[0] + di[1]); 3348e58fa1dSresundermann } 335*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(Diseq,&di)); 336*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES)); 337*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES)); 338aad13602SShrirang Abhyankar } 339*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY)); 340*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY)); 3418e58fa1dSresundermann if (!user->noeqflag) { 342*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterDestroy(&Descat)); 343*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&Deseq)); 3448e58fa1dSresundermann } 345*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterDestroy(&Discat)); 346*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&Diseq)); 347aad13602SShrirang Abhyankar PetscFunctionReturn(0); 348aad13602SShrirang Abhyankar } 349aad13602SShrirang Abhyankar 350628da978Sresundermann /* Evaluate 351628da978Sresundermann h = [ x0^2 - x1; 352628da978Sresundermann 1 -(x0^2 - x1)] 353628da978Sresundermann */ 354aad13602SShrirang Abhyankar PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx) 355aad13602SShrirang Abhyankar { 356aad13602SShrirang Abhyankar const PetscScalar *x; 357aad13602SShrirang Abhyankar PetscScalar ci; 358aad13602SShrirang Abhyankar MPI_Comm comm; 359aad13602SShrirang Abhyankar PetscMPIInt rank; 360aad13602SShrirang Abhyankar AppCtx *user=(AppCtx*)ctx; 361aad13602SShrirang Abhyankar Vec Xseq=user->Xseq; 362aad13602SShrirang Abhyankar VecScatter scat=user->scat; 363aad13602SShrirang Abhyankar 364aad13602SShrirang Abhyankar PetscFunctionBegin; 365*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm)); 366*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_rank(comm,&rank)); 367aad13602SShrirang Abhyankar 368*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD)); 369*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD)); 370aad13602SShrirang Abhyankar 371dd400576SPatrick Sanan if (rank == 0) { 372*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayRead(Xseq,&x)); 373aad13602SShrirang Abhyankar ci = x[0]*x[0] - x[1]; 374*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSetValue(CI,0,ci,INSERT_VALUES)); 375aad13602SShrirang Abhyankar ci = -x[0]*x[0] + x[1] + 1.0; 376*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSetValue(CI,1,ci,INSERT_VALUES)); 377*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(Xseq,&x)); 378aad13602SShrirang Abhyankar } 379*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecAssemblyBegin(CI)); 380*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecAssemblyEnd(CI)); 381aad13602SShrirang Abhyankar PetscFunctionReturn(0); 382aad13602SShrirang Abhyankar } 383aad13602SShrirang Abhyankar 384628da978Sresundermann /* Evaluate 385628da978Sresundermann g = [ x0^2 + x1 - 2] 386628da978Sresundermann */ 387aad13602SShrirang Abhyankar PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx) 388aad13602SShrirang Abhyankar { 389aad13602SShrirang Abhyankar const PetscScalar *x; 390aad13602SShrirang Abhyankar PetscScalar ce; 391aad13602SShrirang Abhyankar MPI_Comm comm; 392aad13602SShrirang Abhyankar PetscMPIInt rank; 393aad13602SShrirang Abhyankar AppCtx *user=(AppCtx*)ctx; 394aad13602SShrirang Abhyankar Vec Xseq=user->Xseq; 395aad13602SShrirang Abhyankar VecScatter scat=user->scat; 396aad13602SShrirang Abhyankar 397aad13602SShrirang Abhyankar PetscFunctionBegin; 398*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm)); 399*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_rank(comm,&rank)); 400aad13602SShrirang Abhyankar 401*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD)); 402*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD)); 403aad13602SShrirang Abhyankar 404dd400576SPatrick Sanan if (rank == 0) { 405*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayRead(Xseq,&x)); 406aad13602SShrirang Abhyankar ce = x[0]*x[0] + x[1] - 2.0; 407*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSetValue(CE,0,ce,INSERT_VALUES)); 408*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(Xseq,&x)); 409aad13602SShrirang Abhyankar } 410*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecAssemblyBegin(CE)); 411*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecAssemblyEnd(CE)); 412aad13602SShrirang Abhyankar PetscFunctionReturn(0); 413aad13602SShrirang Abhyankar } 414aad13602SShrirang Abhyankar 415628da978Sresundermann /* 416628da978Sresundermann grad h = [ 2*x0, -1; 417628da978Sresundermann -2*x0, 1] 418628da978Sresundermann */ 419aad13602SShrirang Abhyankar PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx) 420aad13602SShrirang Abhyankar { 421aad13602SShrirang Abhyankar AppCtx *user=(AppCtx*)ctx; 422f560b561SHong Zhang PetscInt zero=0,one=1,cols[2]; 423aad13602SShrirang Abhyankar PetscScalar vals[2]; 424aad13602SShrirang Abhyankar const PetscScalar *x; 425aad13602SShrirang Abhyankar Vec Xseq=user->Xseq; 426aad13602SShrirang Abhyankar VecScatter scat=user->scat; 427aad13602SShrirang Abhyankar MPI_Comm comm; 428aad13602SShrirang Abhyankar PetscMPIInt rank; 429aad13602SShrirang Abhyankar 430aad13602SShrirang Abhyankar PetscFunctionBegin; 431*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm)); 432*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_rank(comm,&rank)); 433*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD)); 434*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD)); 435aad13602SShrirang Abhyankar 436*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayRead(Xseq,&x)); 437f560b561SHong Zhang if (!rank) { 438aad13602SShrirang Abhyankar cols[0] = 0; cols[1] = 1; 439628da978Sresundermann vals[0] = 2*x[0]; vals[1] = -1.0; 440*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValues(JI,1,&zero,2,cols,vals,INSERT_VALUES)); 441628da978Sresundermann vals[0] = -2*x[0]; vals[1] = 1.0; 442*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValues(JI,1,&one,2,cols,vals,INSERT_VALUES)); 443aad13602SShrirang Abhyankar } 444*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(Xseq,&x)); 445*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY)); 446*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY)); 447aad13602SShrirang Abhyankar PetscFunctionReturn(0); 448aad13602SShrirang Abhyankar } 449aad13602SShrirang Abhyankar 450628da978Sresundermann /* 451628da978Sresundermann grad g = [2*x0 452628da978Sresundermann 1.0 ] 453628da978Sresundermann */ 454aad13602SShrirang Abhyankar PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx) 455aad13602SShrirang Abhyankar { 456f560b561SHong Zhang PetscInt zero=0,cols[2]; 457aad13602SShrirang Abhyankar PetscScalar vals[2]; 458aad13602SShrirang Abhyankar const PetscScalar *x; 459aad13602SShrirang Abhyankar PetscMPIInt rank; 460aad13602SShrirang Abhyankar MPI_Comm comm; 461aad13602SShrirang Abhyankar 462aad13602SShrirang Abhyankar PetscFunctionBegin; 463*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectGetComm((PetscObject)tao,&comm)); 464*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_rank(comm,&rank)); 465aad13602SShrirang Abhyankar 466dd400576SPatrick Sanan if (rank == 0) { 467*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayRead(X,&x)); 468f560b561SHong Zhang cols[0] = 0; cols[1] = 1; 469aad13602SShrirang Abhyankar vals[0] = 2*x[0]; vals[1] = 1.0; 470*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValues(JE,1,&zero,2,cols,vals,INSERT_VALUES)); 471*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(X,&x)); 472aad13602SShrirang Abhyankar } 473*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY)); 474*5f80ce2aSJacob Faibussowitsch CHKERRQ(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: 48409ee8bb0SRylee Sundermann args: -tao_converged_reason -tao_pdipm_kkt_shift_pd 485f560b561SHong Zhang requires: mumps 486aad13602SShrirang Abhyankar 487aad13602SShrirang Abhyankar test: 488aad13602SShrirang Abhyankar suffix: 2 489aad13602SShrirang Abhyankar nsize: 2 49009ee8bb0SRylee Sundermann args: -tao_converged_reason -tao_pdipm_kkt_shift_pd 491f560b561SHong Zhang requires: mumps 492aad13602SShrirang Abhyankar 4938e58fa1dSresundermann test: 4948e58fa1dSresundermann suffix: 3 4958e58fa1dSresundermann args: -tao_converged_reason -no_eq 496f560b561SHong Zhang requires: mumps 4978e58fa1dSresundermann 4988e58fa1dSresundermann test: 4998e58fa1dSresundermann suffix: 4 5008e58fa1dSresundermann nsize: 2 5018e58fa1dSresundermann args: -tao_converged_reason -no_eq 502f560b561SHong Zhang requires: mumps 5038e58fa1dSresundermann 504661095bbSAlp Dener test: 505661095bbSAlp Dener suffix: 5 506661095bbSAlp Dener args: -tao_cmonitor -tao_type almm 507f560b561SHong Zhang requires: mumps 508661095bbSAlp Dener 509661095bbSAlp Dener test: 510661095bbSAlp Dener suffix: 6 511661095bbSAlp Dener args: -tao_cmonitor -tao_type almm -tao_almm_type phr 512f560b561SHong Zhang requires: mumps 513661095bbSAlp Dener 514661095bbSAlp Dener test: 515661095bbSAlp Dener suffix: 7 516661095bbSAlp Dener nsize: 2 517f560b561SHong Zhang requires: mumps 518661095bbSAlp Dener args: -tao_cmonitor -tao_type almm 519661095bbSAlp Dener 520661095bbSAlp Dener test: 521661095bbSAlp Dener suffix: 8 522661095bbSAlp Dener nsize: 2 523f560b561SHong Zhang requires: cuda mumps 524661095bbSAlp Dener args: -tao_cmonitor -tao_type almm -vec_type cuda -mat_type aijcusparse 525661095bbSAlp Dener 526661095bbSAlp Dener test: 527661095bbSAlp Dener suffix: 9 528661095bbSAlp Dener nsize: 2 529661095bbSAlp Dener args: -tao_cmonitor -tao_type almm -no_eq 530f560b561SHong Zhang requires: mumps 531661095bbSAlp Dener 532661095bbSAlp Dener test: 533661095bbSAlp Dener suffix: 10 534661095bbSAlp Dener nsize: 2 535661095bbSAlp Dener args: -tao_cmonitor -tao_type almm -tao_almm_type phr -no_eq 536f560b561SHong Zhang requires: mumps 537661095bbSAlp Dener 538aad13602SShrirang Abhyankar TEST*/ 539