1aad13602SShrirang Abhyankar /* Program usage: mpiexec -n 2 ./ex1 [-help] [all TAO options] */ 2aad13602SShrirang Abhyankar 3aad13602SShrirang Abhyankar /* ---------------------------------------------------------------------- 4*628da978Sresundermann 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 8*628da978Sresundermann --> 9*628da978Sresundermann g(x) = 0 10*628da978Sresundermann h(x) >= 0 11*628da978Sresundermann -1 <= x0, x1 <= 2 12*628da978Sresundermann where 13*628da978Sresundermann g(x) = x0^2 + x1 - 2 14*628da978Sresundermann h(x) = [x0^2 - x1 15*628da978Sresundermann 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\ 24aad13602SShrirang Abhyankar -snes_fd : snes with finite difference Jacobian (needed for pdipm)\n\ 25*628da978Sresundermann -snes_compare_explicit : compare user Jacobian with finite difference Jacobian \n\ 26aad13602SShrirang Abhyankar -tao_cmonitor : convergence monitor with constraint norm \n\ 27aad13602SShrirang Abhyankar -tao_view_solution : view exact solution at each itteration\n\ 28*628da978Sresundermann Note: external package MUMPS is required to run pdipm. This is designed for a maximum of 2 processors, the code will error if size > 2.\n"; 29aad13602SShrirang Abhyankar 30aad13602SShrirang Abhyankar /* 31*628da978Sresundermann User-defined application context - contains data needed by the application 32aad13602SShrirang Abhyankar */ 33aad13602SShrirang Abhyankar typedef struct { 34aad13602SShrirang Abhyankar PetscInt n; /* Global length of x */ 35aad13602SShrirang Abhyankar PetscInt ne; /* Global number of equality constraints */ 36aad13602SShrirang Abhyankar PetscInt ni; /* Global number of inequality constraints */ 378e58fa1dSresundermann PetscBool noeqflag; 38aad13602SShrirang Abhyankar Vec x,xl,xu; 39aad13602SShrirang Abhyankar Vec ce,ci,bl,bu,Xseq; 40aad13602SShrirang Abhyankar Mat Ae,Ai,H; 41aad13602SShrirang Abhyankar VecScatter scat; 42aad13602SShrirang Abhyankar } AppCtx; 43aad13602SShrirang Abhyankar 44aad13602SShrirang Abhyankar /* -------- User-defined Routines --------- */ 45aad13602SShrirang Abhyankar PetscErrorCode InitializeProblem(AppCtx *); 46aad13602SShrirang Abhyankar PetscErrorCode DestroyProblem(AppCtx *); 47aad13602SShrirang Abhyankar PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *); 48aad13602SShrirang Abhyankar PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*); 49aad13602SShrirang Abhyankar PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*); 50aad13602SShrirang Abhyankar PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*); 51aad13602SShrirang Abhyankar PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*); 52aad13602SShrirang Abhyankar PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*); 53aad13602SShrirang Abhyankar 54aad13602SShrirang Abhyankar PetscErrorCode main(int argc,char **argv) 55aad13602SShrirang Abhyankar { 56aad13602SShrirang Abhyankar PetscErrorCode ierr; /* used to check for functions returning nonzeros */ 57aad13602SShrirang Abhyankar Tao tao; 58aad13602SShrirang Abhyankar KSP ksp; 59aad13602SShrirang Abhyankar PC pc; 60aad13602SShrirang Abhyankar AppCtx user; /* application context */ 61aad13602SShrirang Abhyankar Vec x; 62aad13602SShrirang Abhyankar PetscMPIInt rank; 63661095bbSAlp Dener TaoType type; 64661095bbSAlp Dener PetscBool pdipm; 65aad13602SShrirang Abhyankar 66aad13602SShrirang Abhyankar ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 67ffc4695bSBarry Smith ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 68aad13602SShrirang Abhyankar if (rank>1){ 69aad13602SShrirang Abhyankar SETERRQ(PETSC_COMM_SELF,1,"More than 2 processors detected. Example written to use max of 2 processors.\n"); 70aad13602SShrirang Abhyankar } 71aad13602SShrirang Abhyankar 72aad13602SShrirang Abhyankar ierr = PetscPrintf(PETSC_COMM_WORLD,"---- Constrained Problem -----\n");CHKERRQ(ierr); 73aad13602SShrirang Abhyankar ierr = InitializeProblem(&user);CHKERRQ(ierr); /* sets up problem, function below */ 74aad13602SShrirang Abhyankar ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr); 75aad13602SShrirang Abhyankar ierr = TaoSetType(tao,TAOPDIPM);CHKERRQ(ierr); 76aad13602SShrirang Abhyankar ierr = TaoSetInitialVector(tao,user.x);CHKERRQ(ierr); /* gets solution vector from problem */ 77aad13602SShrirang Abhyankar ierr = TaoSetVariableBounds(tao,user.xl,user.xu);CHKERRQ(ierr); /* sets lower upper bounds from given solution */ 78aad13602SShrirang Abhyankar ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);CHKERRQ(ierr); 79aad13602SShrirang Abhyankar 808e58fa1dSresundermann if (!user.noeqflag){ 81aad13602SShrirang Abhyankar ierr = TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user);CHKERRQ(ierr); 828e58fa1dSresundermann } 83aad13602SShrirang Abhyankar ierr = TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user);CHKERRQ(ierr); 848e58fa1dSresundermann if (!user.noeqflag){ 85aad13602SShrirang Abhyankar ierr = TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user);CHKERRQ(ierr); /* equality jacobian */ 868e58fa1dSresundermann } 87aad13602SShrirang Abhyankar ierr = TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user);CHKERRQ(ierr); /* inequality jacobian */ 88aad13602SShrirang Abhyankar ierr = TaoSetTolerances(tao,1.e-6,1.e-6,1.e-6);CHKERRQ(ierr); 89aad13602SShrirang Abhyankar ierr = TaoSetConstraintTolerances(tao,1.e-6,1.e-6);CHKERRQ(ierr); 90aad13602SShrirang Abhyankar 91aad13602SShrirang Abhyankar ierr = TaoGetKSP(tao,&ksp);CHKERRQ(ierr); 92aad13602SShrirang Abhyankar ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 9312d688e0SRylee Sundermann ierr = PCSetType(pc,PCCHOLESKY);CHKERRQ(ierr); 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 */ 9812d688e0SRylee Sundermann ierr = PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);CHKERRQ(ierr); /* requires mumps to solve pdipm */ 99aad13602SShrirang Abhyankar ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr); 100aad13602SShrirang Abhyankar ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 101aad13602SShrirang Abhyankar 102aad13602SShrirang Abhyankar ierr = TaoSetFromOptions(tao);CHKERRQ(ierr); 103661095bbSAlp Dener ierr = TaoGetType(tao,&type);CHKERRQ(ierr); 104661095bbSAlp Dener ierr = PetscObjectTypeCompare((PetscObject)tao,TAOPDIPM,&pdipm);CHKERRQ(ierr); 105661095bbSAlp Dener if (pdipm) { 106661095bbSAlp Dener ierr = TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);CHKERRQ(ierr); 107661095bbSAlp Dener } 108aad13602SShrirang Abhyankar 109aad13602SShrirang Abhyankar ierr = TaoSolve(tao);CHKERRQ(ierr); 110aad13602SShrirang Abhyankar ierr = TaoGetSolutionVector(tao,&x);CHKERRQ(ierr); 111aad13602SShrirang Abhyankar ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 112aad13602SShrirang Abhyankar 113aad13602SShrirang Abhyankar /* Free objects */ 114aad13602SShrirang Abhyankar ierr = DestroyProblem(&user);CHKERRQ(ierr); 115aad13602SShrirang Abhyankar ierr = TaoDestroy(&tao);CHKERRQ(ierr); 116aad13602SShrirang Abhyankar ierr = PetscFinalize(); 117aad13602SShrirang Abhyankar return ierr; 118aad13602SShrirang Abhyankar } 119aad13602SShrirang Abhyankar 120aad13602SShrirang Abhyankar PetscErrorCode InitializeProblem(AppCtx *user) 121aad13602SShrirang Abhyankar { 122aad13602SShrirang Abhyankar PetscErrorCode ierr; 123aad13602SShrirang Abhyankar PetscMPIInt size; 124aad13602SShrirang Abhyankar PetscMPIInt rank; 125aad13602SShrirang Abhyankar PetscInt nloc,neloc,niloc; 126aad13602SShrirang Abhyankar 127aad13602SShrirang Abhyankar PetscFunctionBegin; 128ffc4695bSBarry Smith ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 129ffc4695bSBarry Smith ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 1308e58fa1dSresundermann user->noeqflag = PETSC_FALSE; 1318e58fa1dSresundermann ierr = PetscOptionsGetBool(NULL,NULL,"-no_eq",&user->noeqflag,NULL);CHKERRQ(ierr); 132*628da978Sresundermann 1338e58fa1dSresundermann if (!user->noeqflag) { 1348e58fa1dSresundermann ierr = PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n");CHKERRQ(ierr); 1358e58fa1dSresundermann } 136aad13602SShrirang Abhyankar 1372d4ee042Sprj- /* create vector x and set initial values */ 138aad13602SShrirang Abhyankar user->n = 2; /* global length */ 139aad13602SShrirang Abhyankar nloc = (rank==0)?user->n:0; 140661095bbSAlp Dener ierr = VecCreate(PETSC_COMM_WORLD,&user->x);CHKERRQ(ierr); 141661095bbSAlp Dener ierr = VecSetSizes(user->x,nloc,user->n);CHKERRQ(ierr); 142661095bbSAlp Dener ierr = VecSetFromOptions(user->x);CHKERRQ(ierr); 143aad13602SShrirang Abhyankar ierr = VecSet(user->x,0);CHKERRQ(ierr); 144aad13602SShrirang Abhyankar 145aad13602SShrirang Abhyankar /* create and set lower and upper bound vectors */ 146aad13602SShrirang Abhyankar ierr = VecDuplicate(user->x,&user->xl);CHKERRQ(ierr); 147aad13602SShrirang Abhyankar ierr = VecDuplicate(user->x,&user->xu);CHKERRQ(ierr); 148aad13602SShrirang Abhyankar ierr = VecSet(user->xl,-1.0);CHKERRQ(ierr); 149aad13602SShrirang Abhyankar ierr = VecSet(user->xu,2.0);CHKERRQ(ierr); 150aad13602SShrirang Abhyankar 151aad13602SShrirang Abhyankar /* create scater to zero */ 152aad13602SShrirang Abhyankar ierr = VecScatterCreateToZero(user->x,&user->scat,&user->Xseq);CHKERRQ(ierr); 153aad13602SShrirang Abhyankar 154aad13602SShrirang Abhyankar user->ne = 1; 155*628da978Sresundermann user->ni = 2; 156aad13602SShrirang Abhyankar neloc = (rank==0)?user->ne:0; 157*628da978Sresundermann niloc = (rank==0)?user->ni:0; 158aad13602SShrirang Abhyankar 1598e58fa1dSresundermann if (!user->noeqflag){ 160661095bbSAlp Dener ierr = VecCreate(PETSC_COMM_WORLD,&user->ce);CHKERRQ(ierr); /* a 1x1 vec for equality constraints */ 161661095bbSAlp Dener ierr = VecSetSizes(user->ce,neloc,user->ne);CHKERRQ(ierr); 162661095bbSAlp Dener ierr = VecSetFromOptions(user->ce);CHKERRQ(ierr); 163661095bbSAlp Dener ierr = VecSetUp(user->ce);CHKERRQ(ierr); 1648e58fa1dSresundermann } 165*628da978Sresundermann 166661095bbSAlp Dener ierr = VecCreate(PETSC_COMM_WORLD,&user->ci);CHKERRQ(ierr); /* a 2x1 vec for inequality constraints */ 167661095bbSAlp Dener ierr = VecSetSizes(user->ci,niloc,user->ni);CHKERRQ(ierr); 168661095bbSAlp Dener ierr = VecSetFromOptions(user->ci);CHKERRQ(ierr); 169661095bbSAlp Dener ierr = VecSetUp(user->ci);CHKERRQ(ierr); 170aad13602SShrirang Abhyankar 171aad13602SShrirang Abhyankar /* nexn & nixn matricies for equaly and inequalty constriants */ 1728e58fa1dSresundermann if (!user->noeqflag){ 173aad13602SShrirang Abhyankar ierr = MatCreate(PETSC_COMM_WORLD,&user->Ae);CHKERRQ(ierr); 1748e58fa1dSresundermann ierr = MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n);CHKERRQ(ierr); 1758e58fa1dSresundermann ierr = MatSetFromOptions(user->Ae);CHKERRQ(ierr); 176661095bbSAlp Dener ierr = MatSetUp(user->Ae);CHKERRQ(ierr); 1778e58fa1dSresundermann } 1788e58fa1dSresundermann 179aad13602SShrirang Abhyankar ierr = MatCreate(PETSC_COMM_WORLD,&user->Ai);CHKERRQ(ierr); 180aad13602SShrirang Abhyankar ierr = MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n);CHKERRQ(ierr); 181aad13602SShrirang Abhyankar ierr = MatSetFromOptions(user->Ai);CHKERRQ(ierr); 182661095bbSAlp Dener ierr = MatSetUp(user->Ai);CHKERRQ(ierr); 183*628da978Sresundermann 184*628da978Sresundermann ierr = MatCreate(PETSC_COMM_WORLD,&user->H);CHKERRQ(ierr); 185*628da978Sresundermann ierr = MatSetSizes(user->H,nloc,nloc,user->n,user->n);CHKERRQ(ierr); 186*628da978Sresundermann ierr = MatSetFromOptions(user->H);CHKERRQ(ierr); 187661095bbSAlp Dener ierr = MatSetUp(user->H);CHKERRQ(ierr); 188aad13602SShrirang Abhyankar PetscFunctionReturn(0); 189aad13602SShrirang Abhyankar } 190aad13602SShrirang Abhyankar 191aad13602SShrirang Abhyankar PetscErrorCode DestroyProblem(AppCtx *user) 192aad13602SShrirang Abhyankar { 193aad13602SShrirang Abhyankar PetscErrorCode ierr; 194aad13602SShrirang Abhyankar 195aad13602SShrirang Abhyankar PetscFunctionBegin; 1968e58fa1dSresundermann if (!user->noeqflag){ 197aad13602SShrirang Abhyankar ierr = MatDestroy(&user->Ae);CHKERRQ(ierr); 1988e58fa1dSresundermann } 199aad13602SShrirang Abhyankar ierr = MatDestroy(&user->Ai);CHKERRQ(ierr); 200aad13602SShrirang Abhyankar ierr = MatDestroy(&user->H);CHKERRQ(ierr); 201aad13602SShrirang Abhyankar 202aad13602SShrirang Abhyankar ierr = VecDestroy(&user->x);CHKERRQ(ierr); 2038e58fa1dSresundermann if (!user->noeqflag){ 204aad13602SShrirang Abhyankar ierr = VecDestroy(&user->ce);CHKERRQ(ierr); 2058e58fa1dSresundermann } 206aad13602SShrirang Abhyankar ierr = VecDestroy(&user->ci);CHKERRQ(ierr); 207aad13602SShrirang Abhyankar ierr = VecDestroy(&user->xl);CHKERRQ(ierr); 208aad13602SShrirang Abhyankar ierr = VecDestroy(&user->xu);CHKERRQ(ierr); 209aad13602SShrirang Abhyankar ierr = VecDestroy(&user->Xseq);CHKERRQ(ierr); 210aad13602SShrirang Abhyankar ierr = VecScatterDestroy(&user->scat);CHKERRQ(ierr); 211aad13602SShrirang Abhyankar PetscFunctionReturn(0); 212aad13602SShrirang Abhyankar } 213aad13602SShrirang Abhyankar 214*628da978Sresundermann /* Evaluate 215*628da978Sresundermann f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1) 216*628da978Sresundermann G = grad f = [2*(x0 - 2) - 2; 217*628da978Sresundermann 2*(x1 - 2) - 2] 218aad13602SShrirang Abhyankar */ 219aad13602SShrirang Abhyankar PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx) 220aad13602SShrirang Abhyankar { 221aad13602SShrirang Abhyankar PetscScalar g; 222aad13602SShrirang Abhyankar const PetscScalar *x; 223aad13602SShrirang Abhyankar MPI_Comm comm; 224aad13602SShrirang Abhyankar PetscMPIInt rank; 225aad13602SShrirang Abhyankar PetscErrorCode ierr; 226aad13602SShrirang Abhyankar PetscReal fin; 227aad13602SShrirang Abhyankar AppCtx *user=(AppCtx*)ctx; 228aad13602SShrirang Abhyankar Vec Xseq=user->Xseq; 229aad13602SShrirang Abhyankar VecScatter scat=user->scat; 230aad13602SShrirang Abhyankar 231aad13602SShrirang Abhyankar PetscFunctionBegin; 232aad13602SShrirang Abhyankar ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 233ffc4695bSBarry Smith ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 234aad13602SShrirang Abhyankar 235aad13602SShrirang Abhyankar ierr = VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 236aad13602SShrirang Abhyankar ierr = VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 237aad13602SShrirang Abhyankar 238aad13602SShrirang Abhyankar fin = 0.0; 239aad13602SShrirang Abhyankar if (!rank) { 240aad13602SShrirang Abhyankar ierr = VecGetArrayRead(Xseq,&x);CHKERRQ(ierr); 241aad13602SShrirang Abhyankar fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]); 242aad13602SShrirang Abhyankar g = 2.0*(x[0]-2.0) - 2.0; 243aad13602SShrirang Abhyankar ierr = VecSetValue(G,0,g,INSERT_VALUES);CHKERRQ(ierr); 244aad13602SShrirang Abhyankar g = 2.0*(x[1]-2.0) - 2.0; 245aad13602SShrirang Abhyankar ierr = VecSetValue(G,1,g,INSERT_VALUES);CHKERRQ(ierr); 246aad13602SShrirang Abhyankar ierr = VecRestoreArrayRead(Xseq,&x);CHKERRQ(ierr); 247aad13602SShrirang Abhyankar } 248ffc4695bSBarry Smith ierr = MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm);CHKERRMPI(ierr); 249aad13602SShrirang Abhyankar ierr = VecAssemblyBegin(G);CHKERRQ(ierr); 250aad13602SShrirang Abhyankar ierr = VecAssemblyEnd(G);CHKERRQ(ierr); 251aad13602SShrirang Abhyankar PetscFunctionReturn(0); 252aad13602SShrirang Abhyankar } 253aad13602SShrirang Abhyankar 254*628da978Sresundermann /* Evaluate 255*628da978Sresundermann H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)] 256*628da978Sresundermann = [ 2*(1+de[0]-di[0]+di[1]), 0; 257*628da978Sresundermann 0, 2] 258*628da978Sresundermann */ 259aad13602SShrirang Abhyankar PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx) 260aad13602SShrirang Abhyankar { 261aad13602SShrirang Abhyankar AppCtx *user=(AppCtx*)ctx; 262aad13602SShrirang Abhyankar Vec DE,DI; 263aad13602SShrirang Abhyankar const PetscScalar *de,*di; 264aad13602SShrirang Abhyankar PetscInt zero=0,one=1; 265aad13602SShrirang Abhyankar PetscScalar two=2.0; 266aad13602SShrirang Abhyankar PetscScalar val=0.0; 267aad13602SShrirang Abhyankar Vec Deseq,Diseq=user->Xseq; 268aad13602SShrirang Abhyankar VecScatter Descat,Discat=user->scat; 269aad13602SShrirang Abhyankar PetscMPIInt rank; 270aad13602SShrirang Abhyankar MPI_Comm comm; 271aad13602SShrirang Abhyankar PetscErrorCode ierr; 272aad13602SShrirang Abhyankar 273aad13602SShrirang Abhyankar PetscFunctionBegin; 274aad13602SShrirang Abhyankar ierr = TaoGetDualVariables(tao,&DE,&DI);CHKERRQ(ierr); 275aad13602SShrirang Abhyankar 276aad13602SShrirang Abhyankar ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 277ffc4695bSBarry Smith ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 278aad13602SShrirang Abhyankar 2798e58fa1dSresundermann if (!user->noeqflag){ 280aad13602SShrirang Abhyankar ierr = VecScatterCreateToZero(DE,&Descat,&Deseq);CHKERRQ(ierr); 281aad13602SShrirang Abhyankar ierr = VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 282aad13602SShrirang Abhyankar ierr = VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2838e58fa1dSresundermann } 2848e58fa1dSresundermann ierr = VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 285aad13602SShrirang Abhyankar ierr = VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 286aad13602SShrirang Abhyankar 287aad13602SShrirang Abhyankar if (!rank){ 2888e58fa1dSresundermann if (!user->noeqflag){ 289aad13602SShrirang Abhyankar ierr = VecGetArrayRead(Deseq,&de);CHKERRQ(ierr); /* places equality constraint dual into array */ 2908e58fa1dSresundermann } 291aad13602SShrirang Abhyankar ierr = VecGetArrayRead(Diseq,&di);CHKERRQ(ierr); /* places inequality constraint dual into array */ 292*628da978Sresundermann 2938e58fa1dSresundermann if (!user->noeqflag){ 294*628da978Sresundermann val = 2.0 * (1 + de[0] - di[0] + di[1]); 295aad13602SShrirang Abhyankar ierr = VecRestoreArrayRead(Deseq,&de);CHKERRQ(ierr); 296*628da978Sresundermann ierr = VecRestoreArrayRead(Diseq,&di);CHKERRQ(ierr); 2978e58fa1dSresundermann } else { 298*628da978Sresundermann val = 2.0 * (1 - di[0] + di[1]); 2998e58fa1dSresundermann } 300aad13602SShrirang Abhyankar ierr = VecRestoreArrayRead(Diseq,&di);CHKERRQ(ierr); 301aad13602SShrirang Abhyankar ierr = MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES);CHKERRQ(ierr); 302aad13602SShrirang Abhyankar ierr = MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES);CHKERRQ(ierr); 303aad13602SShrirang Abhyankar } 304aad13602SShrirang Abhyankar ierr = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 305aad13602SShrirang Abhyankar ierr = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3068e58fa1dSresundermann if (!user->noeqflag){ 307aad13602SShrirang Abhyankar ierr = VecScatterDestroy(&Descat);CHKERRQ(ierr); 308aad13602SShrirang Abhyankar ierr = VecDestroy(&Deseq);CHKERRQ(ierr); 3098e58fa1dSresundermann } 310aad13602SShrirang Abhyankar PetscFunctionReturn(0); 311aad13602SShrirang Abhyankar } 312aad13602SShrirang Abhyankar 313*628da978Sresundermann /* Evaluate 314*628da978Sresundermann h = [ x0^2 - x1; 315*628da978Sresundermann 1 -(x0^2 - x1)] 316*628da978Sresundermann */ 317aad13602SShrirang Abhyankar PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx) 318aad13602SShrirang Abhyankar { 319aad13602SShrirang Abhyankar const PetscScalar *x; 320aad13602SShrirang Abhyankar PetscScalar ci; 321aad13602SShrirang Abhyankar PetscErrorCode ierr; 322aad13602SShrirang Abhyankar MPI_Comm comm; 323aad13602SShrirang Abhyankar PetscMPIInt rank; 324aad13602SShrirang Abhyankar AppCtx *user=(AppCtx*)ctx; 325aad13602SShrirang Abhyankar Vec Xseq=user->Xseq; 326aad13602SShrirang Abhyankar VecScatter scat=user->scat; 327aad13602SShrirang Abhyankar 328aad13602SShrirang Abhyankar PetscFunctionBegin; 329aad13602SShrirang Abhyankar ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 330ffc4695bSBarry Smith ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 331aad13602SShrirang Abhyankar 332aad13602SShrirang Abhyankar ierr = VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 333aad13602SShrirang Abhyankar ierr = VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 334aad13602SShrirang Abhyankar 335aad13602SShrirang Abhyankar if (!rank) { 336aad13602SShrirang Abhyankar ierr = VecGetArrayRead(Xseq,&x);CHKERRQ(ierr); 337aad13602SShrirang Abhyankar ci = x[0]*x[0] - x[1]; 338aad13602SShrirang Abhyankar ierr = VecSetValue(CI,0,ci,INSERT_VALUES);CHKERRQ(ierr); 339aad13602SShrirang Abhyankar ci = -x[0]*x[0] + x[1] + 1.0; 340aad13602SShrirang Abhyankar ierr = VecSetValue(CI,1,ci,INSERT_VALUES);CHKERRQ(ierr); 341aad13602SShrirang Abhyankar ierr = VecRestoreArrayRead(Xseq,&x);CHKERRQ(ierr); 342aad13602SShrirang Abhyankar } 343aad13602SShrirang Abhyankar ierr = VecAssemblyBegin(CI);CHKERRQ(ierr); 344aad13602SShrirang Abhyankar ierr = VecAssemblyEnd(CI);CHKERRQ(ierr); 345aad13602SShrirang Abhyankar PetscFunctionReturn(0); 346aad13602SShrirang Abhyankar } 347aad13602SShrirang Abhyankar 348*628da978Sresundermann /* Evaluate 349*628da978Sresundermann g = [ x0^2 + x1 - 2] 350*628da978Sresundermann */ 351aad13602SShrirang Abhyankar PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx) 352aad13602SShrirang Abhyankar { 353aad13602SShrirang Abhyankar const PetscScalar *x; 354aad13602SShrirang Abhyankar PetscScalar ce; 355aad13602SShrirang Abhyankar PetscErrorCode ierr; 356aad13602SShrirang Abhyankar MPI_Comm comm; 357aad13602SShrirang Abhyankar PetscMPIInt rank; 358aad13602SShrirang Abhyankar AppCtx *user=(AppCtx*)ctx; 359aad13602SShrirang Abhyankar Vec Xseq=user->Xseq; 360aad13602SShrirang Abhyankar VecScatter scat=user->scat; 361aad13602SShrirang Abhyankar 362aad13602SShrirang Abhyankar PetscFunctionBegin; 363aad13602SShrirang Abhyankar ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 364ffc4695bSBarry Smith ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 365aad13602SShrirang Abhyankar 366aad13602SShrirang Abhyankar ierr = VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 367aad13602SShrirang Abhyankar ierr = VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 368aad13602SShrirang Abhyankar 369aad13602SShrirang Abhyankar if (!rank) { 370aad13602SShrirang Abhyankar ierr = VecGetArrayRead(Xseq,&x);CHKERRQ(ierr); 371aad13602SShrirang Abhyankar ce = x[0]*x[0] + x[1] - 2.0; 372aad13602SShrirang Abhyankar ierr = VecSetValue(CE,0,ce,INSERT_VALUES);CHKERRQ(ierr); 373aad13602SShrirang Abhyankar ierr = VecRestoreArrayRead(Xseq,&x);CHKERRQ(ierr); 374aad13602SShrirang Abhyankar } 375aad13602SShrirang Abhyankar ierr = VecAssemblyBegin(CE);CHKERRQ(ierr); 376aad13602SShrirang Abhyankar ierr = VecAssemblyEnd(CE);CHKERRQ(ierr); 377aad13602SShrirang Abhyankar PetscFunctionReturn(0); 378aad13602SShrirang Abhyankar } 379aad13602SShrirang Abhyankar 380*628da978Sresundermann /* 381*628da978Sresundermann grad h = [ 2*x0, -1; 382*628da978Sresundermann -2*x0, 1] 383*628da978Sresundermann */ 384aad13602SShrirang Abhyankar PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx) 385aad13602SShrirang Abhyankar { 386aad13602SShrirang Abhyankar AppCtx *user=(AppCtx*)ctx; 387aad13602SShrirang Abhyankar PetscInt cols[2],min,max,i; 388aad13602SShrirang Abhyankar PetscScalar vals[2]; 389aad13602SShrirang Abhyankar const PetscScalar *x; 390aad13602SShrirang Abhyankar PetscErrorCode ierr; 391aad13602SShrirang Abhyankar Vec Xseq=user->Xseq; 392aad13602SShrirang Abhyankar VecScatter scat=user->scat; 393aad13602SShrirang Abhyankar MPI_Comm comm; 394aad13602SShrirang Abhyankar PetscMPIInt rank; 395aad13602SShrirang Abhyankar 396aad13602SShrirang Abhyankar PetscFunctionBegin; 397aad13602SShrirang Abhyankar ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 398ffc4695bSBarry Smith ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 399aad13602SShrirang Abhyankar ierr = VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 400aad13602SShrirang Abhyankar ierr = VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 401aad13602SShrirang Abhyankar 402aad13602SShrirang Abhyankar ierr = VecGetArrayRead(Xseq,&x);CHKERRQ(ierr); 403aad13602SShrirang Abhyankar ierr = MatGetOwnershipRange(JI,&min,&max);CHKERRQ(ierr); 404aad13602SShrirang Abhyankar 405aad13602SShrirang Abhyankar cols[0] = 0; cols[1] = 1; 406aad13602SShrirang Abhyankar for (i=min;i<max;i++) { 407aad13602SShrirang Abhyankar if (i==0){ 408*628da978Sresundermann vals[0] = 2*x[0]; vals[1] = -1.0; 409aad13602SShrirang Abhyankar ierr = MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);CHKERRQ(ierr); 410aad13602SShrirang Abhyankar } 411aad13602SShrirang Abhyankar if (i==1) { 412*628da978Sresundermann vals[0] = -2*x[0]; vals[1] = 1.0; 413aad13602SShrirang Abhyankar ierr = MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);CHKERRQ(ierr); 414aad13602SShrirang Abhyankar } 415aad13602SShrirang Abhyankar } 416aad13602SShrirang Abhyankar ierr = VecRestoreArrayRead(Xseq,&x);CHKERRQ(ierr); 417aad13602SShrirang Abhyankar 418aad13602SShrirang Abhyankar ierr = MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 419aad13602SShrirang Abhyankar ierr = MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 420aad13602SShrirang Abhyankar PetscFunctionReturn(0); 421aad13602SShrirang Abhyankar } 422aad13602SShrirang Abhyankar 423*628da978Sresundermann /* 424*628da978Sresundermann grad g = [2*x0 425*628da978Sresundermann 1.0 ] 426*628da978Sresundermann */ 427aad13602SShrirang Abhyankar PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx) 428aad13602SShrirang Abhyankar { 429aad13602SShrirang Abhyankar PetscInt rows[2]; 430aad13602SShrirang Abhyankar PetscScalar vals[2]; 431aad13602SShrirang Abhyankar const PetscScalar *x; 432aad13602SShrirang Abhyankar PetscMPIInt rank; 433aad13602SShrirang Abhyankar MPI_Comm comm; 434aad13602SShrirang Abhyankar PetscErrorCode ierr; 435aad13602SShrirang Abhyankar 436aad13602SShrirang Abhyankar PetscFunctionBegin; 437aad13602SShrirang Abhyankar ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 438ffc4695bSBarry Smith ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 439aad13602SShrirang Abhyankar 440aad13602SShrirang Abhyankar if (!rank) { 441aad13602SShrirang Abhyankar ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 442aad13602SShrirang Abhyankar rows[0] = 0; rows[1] = 1; 443aad13602SShrirang Abhyankar vals[0] = 2*x[0]; vals[1] = 1.0; 444aad13602SShrirang Abhyankar ierr = MatSetValues(JE,1,rows,2,rows,vals,INSERT_VALUES);CHKERRQ(ierr); 445aad13602SShrirang Abhyankar ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 446aad13602SShrirang Abhyankar } 447aad13602SShrirang Abhyankar ierr = MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 448aad13602SShrirang Abhyankar ierr = MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 449aad13602SShrirang Abhyankar PetscFunctionReturn(0); 450aad13602SShrirang Abhyankar } 451aad13602SShrirang Abhyankar 452aad13602SShrirang Abhyankar 453aad13602SShrirang Abhyankar /*TEST 454aad13602SShrirang Abhyankar 455aad13602SShrirang Abhyankar build: 45612d688e0SRylee Sundermann requires: !complex !define(PETSC_USE_CXX) mumps 457aad13602SShrirang Abhyankar 458aad13602SShrirang Abhyankar test: 45909ee8bb0SRylee Sundermann args: -tao_converged_reason -tao_pdipm_kkt_shift_pd 460aad13602SShrirang Abhyankar 461aad13602SShrirang Abhyankar test: 462aad13602SShrirang Abhyankar suffix: 2 463aad13602SShrirang Abhyankar nsize: 2 46409ee8bb0SRylee Sundermann args: -tao_converged_reason -tao_pdipm_kkt_shift_pd 465aad13602SShrirang Abhyankar 4668e58fa1dSresundermann test: 4678e58fa1dSresundermann suffix: 3 4688e58fa1dSresundermann args: -tao_converged_reason -no_eq 4698e58fa1dSresundermann 4708e58fa1dSresundermann test: 4718e58fa1dSresundermann suffix: 4 4728e58fa1dSresundermann nsize: 2 4738e58fa1dSresundermann args: -tao_converged_reason -no_eq 4748e58fa1dSresundermann 475661095bbSAlp Dener test: 476661095bbSAlp Dener suffix: 5 477661095bbSAlp Dener args: -tao_cmonitor -tao_type almm 478661095bbSAlp Dener 479661095bbSAlp Dener test: 480661095bbSAlp Dener suffix: 6 481661095bbSAlp Dener args: -tao_cmonitor -tao_type almm -tao_almm_type phr 482661095bbSAlp Dener 483661095bbSAlp Dener test: 484661095bbSAlp Dener suffix: 7 485661095bbSAlp Dener nsize: 2 486661095bbSAlp Dener args: -tao_cmonitor -tao_type almm 487661095bbSAlp Dener 488661095bbSAlp Dener test: 489661095bbSAlp Dener suffix: 8 490661095bbSAlp Dener nsize: 2 491661095bbSAlp Dener requires: cuda 492661095bbSAlp Dener args: -tao_cmonitor -tao_type almm -vec_type cuda -mat_type aijcusparse 493661095bbSAlp Dener 494661095bbSAlp Dener test: 495661095bbSAlp Dener suffix: 9 496661095bbSAlp Dener nsize: 2 497661095bbSAlp Dener args: -tao_cmonitor -tao_type almm -no_eq 498661095bbSAlp Dener 499661095bbSAlp Dener test: 500661095bbSAlp Dener suffix: 10 501661095bbSAlp Dener nsize: 2 502661095bbSAlp Dener args: -tao_cmonitor -tao_type almm -tao_almm_type phr -no_eq 503661095bbSAlp Dener 504aad13602SShrirang Abhyankar TEST*/ 505