1 /* Program usage: mpiexec -n 2 ./ex1 [-help] [all TAO options] */ 2 3 /* ---------------------------------------------------------------------- 4 min f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1) 5 s.t. x0^2 + x1 - 2 = 0 6 0 <= x0^2 - x1 <= 1 7 -1 <= x0, x1 <= 2 8 --> 9 g(x) = 0 10 h(x) >= 0 11 -1 <= x0, x1 <= 2 12 where 13 g(x) = x0^2 + x1 - 2 14 h(x) = [x0^2 - x1 15 1 -(x0^2 - x1)] 16 ---------------------------------------------------------------------- */ 17 18 #include <petsctao.h> 19 20 static char help[]= "Solves constrained optimiztion problem using pdipm.\n\ 21 Input parameters include:\n\ 22 -tao_type pdipm : sets Tao solver\n\ 23 -no_eq : removes the equaility constraints from the problem\n\ 24 -init_view : view initial object setup\n\ 25 -snes_fd : snes with finite difference Jacobian (needed for pdipm)\n\ 26 -snes_compare_explicit : compare user Jacobian with finite difference Jacobian \n\ 27 -tao_cmonitor : convergence monitor with constraint norm \n\ 28 -tao_view_solution : view exact solution at each iteration\n\ 29 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"; 30 31 /* 32 User-defined application context - contains data needed by the application 33 */ 34 typedef struct { 35 PetscInt n; /* Global length of x */ 36 PetscInt ne; /* Global number of equality constraints */ 37 PetscInt ni; /* Global number of inequality constraints */ 38 PetscBool noeqflag, initview; 39 Vec x,xl,xu; 40 Vec ce,ci,bl,bu,Xseq; 41 Mat Ae,Ai,H; 42 VecScatter scat; 43 } AppCtx; 44 45 /* -------- User-defined Routines --------- */ 46 PetscErrorCode InitializeProblem(AppCtx *); 47 PetscErrorCode DestroyProblem(AppCtx *); 48 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *); 49 PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*); 50 PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*); 51 PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*); 52 PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*); 53 PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*); 54 55 PetscErrorCode main(int argc,char **argv) 56 { 57 PetscErrorCode ierr; /* used to check for functions returning nonzeros */ 58 Tao tao; 59 KSP ksp; 60 PC pc; 61 AppCtx user; /* application context */ 62 Vec x,G,CI,CE; 63 PetscMPIInt size; 64 TaoType type; 65 PetscReal f; 66 PetscBool pdipm; 67 68 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 69 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 70 if (size>2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,"More than 2 processors detected. Example written to use max of 2 processors.\n"); 71 72 ierr = PetscPrintf(PETSC_COMM_WORLD,"---- Constrained Problem -----\n");CHKERRQ(ierr); 73 ierr = InitializeProblem(&user);CHKERRQ(ierr); /* sets up problem, function below */ 74 ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr); 75 ierr = TaoSetType(tao,TAOPDIPM);CHKERRQ(ierr); 76 ierr = TaoSetInitialVector(tao,user.x);CHKERRQ(ierr); /* gets solution vector from problem */ 77 ierr = TaoSetVariableBounds(tao,user.xl,user.xu);CHKERRQ(ierr); /* sets lower upper bounds from given solution */ 78 ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);CHKERRQ(ierr); 79 80 if (!user.noeqflag) { 81 ierr = TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user);CHKERRQ(ierr); 82 } 83 ierr = TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user);CHKERRQ(ierr); 84 if (!user.noeqflag) { 85 ierr = TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user);CHKERRQ(ierr); /* equality jacobian */ 86 } 87 ierr = TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user);CHKERRQ(ierr); /* inequality jacobian */ 88 ierr = TaoSetTolerances(tao,1.e-6,1.e-6,1.e-6);CHKERRQ(ierr); 89 ierr = TaoSetConstraintTolerances(tao,1.e-6,1.e-6);CHKERRQ(ierr); 90 91 ierr = TaoGetKSP(tao,&ksp);CHKERRQ(ierr); 92 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 93 ierr = PCSetType(pc,PCCHOLESKY);CHKERRQ(ierr); 94 /* 95 This algorithm produces matrices with zeros along the diagonal therefore we use 96 MUMPS which provides solver for indefinite matrices 97 */ 98 #if defined(PETSC_HAVE_MUMPS) 99 ierr = PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);CHKERRQ(ierr); /* requires mumps to solve pdipm */ 100 #else 101 if (size > 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Requires an external package that supports parallel PCCHOLESKY, e.g., MUMPS."); 102 #endif 103 ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr); 104 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 105 106 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr); 107 ierr = TaoGetType(tao,&type);CHKERRQ(ierr); 108 ierr = PetscObjectTypeCompare((PetscObject)tao,TAOPDIPM,&pdipm);CHKERRQ(ierr); 109 if (pdipm) { 110 ierr = TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);CHKERRQ(ierr); 111 if (user.initview) { 112 ierr = TaoSetUp(tao);CHKERRQ(ierr); 113 ierr = VecDuplicate(user.x, &G);CHKERRQ(ierr); 114 ierr = FormFunctionGradient(tao, user.x, &f, G, (void*)&user);CHKERRQ(ierr); 115 ierr = FormHessian(tao, user.x, user.H, user.H, (void*)&user);CHKERRQ(ierr); 116 ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 117 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nInitial point X:\n",f);CHKERRQ(ierr); 118 ierr = VecView(user.x, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 119 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nInitial objective f(x) = %g\n",f);CHKERRQ(ierr); 120 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nInitial gradient and Hessian:\n",f);CHKERRQ(ierr); 121 ierr = VecView(G, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 122 ierr = MatView(user.H, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 123 ierr = VecDestroy(&G);CHKERRQ(ierr); 124 ierr = FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void*)&user);CHKERRQ(ierr); 125 ierr = MatCreateVecs(user.Ai, NULL, &CI);CHKERRQ(ierr); 126 ierr = FormInequalityConstraints(tao, user.x, CI, (void*)&user);CHKERRQ(ierr); 127 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nInitial inequality constraints and Jacobian:\n",f);CHKERRQ(ierr); 128 ierr = VecView(CI, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 129 ierr = MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 130 ierr = VecDestroy(&CI);CHKERRQ(ierr); 131 if (!user.noeqflag) { 132 ierr = FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void*)&user);CHKERRQ(ierr); 133 ierr = MatCreateVecs(user.Ae, NULL, &CE);CHKERRQ(ierr); 134 ierr = FormEqualityConstraints(tao, user.x, CE, (void*)&user);CHKERRQ(ierr); 135 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nInitial equality constraints and Jacobian:\n",f);CHKERRQ(ierr); 136 ierr = VecView(CE, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 137 ierr = MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 138 ierr = VecDestroy(&CE);CHKERRQ(ierr); 139 } 140 ierr = PetscPrintf(PETSC_COMM_WORLD, "\n");CHKERRQ(ierr); 141 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 142 } 143 } 144 145 ierr = TaoSolve(tao);CHKERRQ(ierr); 146 ierr = TaoGetSolutionVector(tao,&x);CHKERRQ(ierr); 147 ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 148 149 /* Free objects */ 150 ierr = DestroyProblem(&user);CHKERRQ(ierr); 151 ierr = TaoDestroy(&tao);CHKERRQ(ierr); 152 ierr = PetscFinalize(); 153 return ierr; 154 } 155 156 PetscErrorCode InitializeProblem(AppCtx *user) 157 { 158 PetscErrorCode ierr; 159 PetscMPIInt size; 160 PetscMPIInt rank; 161 PetscInt nloc,neloc,niloc; 162 163 PetscFunctionBegin; 164 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 165 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 166 user->noeqflag = PETSC_FALSE; 167 ierr = PetscOptionsGetBool(NULL,NULL,"-no_eq",&user->noeqflag,NULL);CHKERRQ(ierr); 168 user->initview = PETSC_FALSE; 169 ierr = PetscOptionsGetBool(NULL,NULL,"-init_view",&user->initview,NULL);CHKERRQ(ierr); 170 171 if (!user->noeqflag) { 172 /* Tell user the correct solution, not an error checking */ 173 ierr = PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n");CHKERRQ(ierr); 174 } 175 176 /* create vector x and set initial values */ 177 user->n = 2; /* global length */ 178 nloc = (size==1)?user->n:1; 179 ierr = VecCreate(PETSC_COMM_WORLD,&user->x);CHKERRQ(ierr); 180 ierr = VecSetSizes(user->x,nloc,user->n);CHKERRQ(ierr); 181 ierr = VecSetFromOptions(user->x);CHKERRQ(ierr); 182 ierr = VecSet(user->x,0);CHKERRQ(ierr); 183 184 /* create and set lower and upper bound vectors */ 185 ierr = VecDuplicate(user->x,&user->xl);CHKERRQ(ierr); 186 ierr = VecDuplicate(user->x,&user->xu);CHKERRQ(ierr); 187 ierr = VecSet(user->xl,-1.0);CHKERRQ(ierr); 188 ierr = VecSet(user->xu,2.0);CHKERRQ(ierr); 189 190 /* create scater to zero */ 191 ierr = VecScatterCreateToZero(user->x,&user->scat,&user->Xseq);CHKERRQ(ierr); 192 193 user->ne = 1; 194 user->ni = 2; 195 neloc = (rank==0)?user->ne:0; 196 niloc = (size==1)?user->ni:1; 197 198 if (!user->noeqflag) { 199 ierr = VecCreate(PETSC_COMM_WORLD,&user->ce);CHKERRQ(ierr); /* a 1x1 vec for equality constraints */ 200 ierr = VecSetSizes(user->ce,neloc,user->ne);CHKERRQ(ierr); 201 ierr = VecSetFromOptions(user->ce);CHKERRQ(ierr); 202 ierr = VecSetUp(user->ce);CHKERRQ(ierr); 203 } 204 205 ierr = VecCreate(PETSC_COMM_WORLD,&user->ci);CHKERRQ(ierr); /* a 2x1 vec for inequality constraints */ 206 ierr = VecSetSizes(user->ci,niloc,user->ni);CHKERRQ(ierr); 207 ierr = VecSetFromOptions(user->ci);CHKERRQ(ierr); 208 ierr = VecSetUp(user->ci);CHKERRQ(ierr); 209 210 /* nexn & nixn matricies for equally and inequalty constraints */ 211 if (!user->noeqflag) { 212 ierr = MatCreate(PETSC_COMM_WORLD,&user->Ae);CHKERRQ(ierr); 213 ierr = MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n);CHKERRQ(ierr); 214 ierr = MatSetFromOptions(user->Ae);CHKERRQ(ierr); 215 ierr = MatSetUp(user->Ae);CHKERRQ(ierr); 216 } 217 218 ierr = MatCreate(PETSC_COMM_WORLD,&user->Ai);CHKERRQ(ierr); 219 ierr = MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n);CHKERRQ(ierr); 220 ierr = MatSetFromOptions(user->Ai);CHKERRQ(ierr); 221 ierr = MatSetUp(user->Ai);CHKERRQ(ierr); 222 223 ierr = MatCreate(PETSC_COMM_WORLD,&user->H);CHKERRQ(ierr); 224 ierr = MatSetSizes(user->H,nloc,nloc,user->n,user->n);CHKERRQ(ierr); 225 ierr = MatSetFromOptions(user->H);CHKERRQ(ierr); 226 ierr = MatSetUp(user->H);CHKERRQ(ierr); 227 PetscFunctionReturn(0); 228 } 229 230 PetscErrorCode DestroyProblem(AppCtx *user) 231 { 232 PetscErrorCode ierr; 233 234 PetscFunctionBegin; 235 if (!user->noeqflag) { 236 ierr = MatDestroy(&user->Ae);CHKERRQ(ierr); 237 } 238 ierr = MatDestroy(&user->Ai);CHKERRQ(ierr); 239 ierr = MatDestroy(&user->H);CHKERRQ(ierr); 240 241 ierr = VecDestroy(&user->x);CHKERRQ(ierr); 242 if (!user->noeqflag) { 243 ierr = VecDestroy(&user->ce);CHKERRQ(ierr); 244 } 245 ierr = VecDestroy(&user->ci);CHKERRQ(ierr); 246 ierr = VecDestroy(&user->xl);CHKERRQ(ierr); 247 ierr = VecDestroy(&user->xu);CHKERRQ(ierr); 248 ierr = VecDestroy(&user->Xseq);CHKERRQ(ierr); 249 ierr = VecScatterDestroy(&user->scat);CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 253 /* Evaluate 254 f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1) 255 G = grad f = [2*(x0 - 2) - 2; 256 2*(x1 - 2) - 2] 257 */ 258 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx) 259 { 260 PetscScalar g; 261 const PetscScalar *x; 262 MPI_Comm comm; 263 PetscMPIInt rank; 264 PetscErrorCode ierr; 265 PetscReal fin; 266 AppCtx *user=(AppCtx*)ctx; 267 Vec Xseq=user->Xseq; 268 VecScatter scat=user->scat; 269 270 PetscFunctionBegin; 271 ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 272 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 273 274 ierr = VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 275 ierr = VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 276 277 fin = 0.0; 278 if (rank == 0) { 279 ierr = VecGetArrayRead(Xseq,&x);CHKERRQ(ierr); 280 fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]); 281 g = 2.0*(x[0]-2.0) - 2.0; 282 ierr = VecSetValue(G,0,g,INSERT_VALUES);CHKERRQ(ierr); 283 g = 2.0*(x[1]-2.0) - 2.0; 284 ierr = VecSetValue(G,1,g,INSERT_VALUES);CHKERRQ(ierr); 285 ierr = VecRestoreArrayRead(Xseq,&x);CHKERRQ(ierr); 286 } 287 ierr = MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm);CHKERRMPI(ierr); 288 ierr = VecAssemblyBegin(G);CHKERRQ(ierr); 289 ierr = VecAssemblyEnd(G);CHKERRQ(ierr); 290 PetscFunctionReturn(0); 291 } 292 293 /* Evaluate 294 H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)] 295 = [ 2*(1+de[0]-di[0]+di[1]), 0; 296 0, 2] 297 */ 298 PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx) 299 { 300 AppCtx *user=(AppCtx*)ctx; 301 Vec DE,DI; 302 const PetscScalar *de,*di; 303 PetscInt zero=0,one=1; 304 PetscScalar two=2.0; 305 PetscScalar val=0.0; 306 Vec Deseq,Diseq; 307 VecScatter Descat,Discat; 308 PetscMPIInt rank; 309 MPI_Comm comm; 310 PetscErrorCode ierr; 311 312 PetscFunctionBegin; 313 ierr = TaoGetDualVariables(tao,&DE,&DI);CHKERRQ(ierr); 314 315 ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 316 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 317 318 if (!user->noeqflag) { 319 ierr = VecScatterCreateToZero(DE,&Descat,&Deseq);CHKERRQ(ierr); 320 ierr = VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 321 ierr = VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 322 } 323 ierr = VecScatterCreateToZero(DI,&Discat,&Diseq);CHKERRQ(ierr); 324 ierr = VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 325 ierr = VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 326 327 if (rank == 0) { 328 if (!user->noeqflag) { 329 ierr = VecGetArrayRead(Deseq,&de);CHKERRQ(ierr); /* places equality constraint dual into array */ 330 } 331 ierr = VecGetArrayRead(Diseq,&di);CHKERRQ(ierr); /* places inequality constraint dual into array */ 332 333 if (!user->noeqflag) { 334 val = 2.0 * (1 + de[0] - di[0] + di[1]); 335 ierr = VecRestoreArrayRead(Deseq,&de);CHKERRQ(ierr); 336 ierr = VecRestoreArrayRead(Diseq,&di);CHKERRQ(ierr); 337 } else { 338 val = 2.0 * (1 - di[0] + di[1]); 339 } 340 ierr = VecRestoreArrayRead(Diseq,&di);CHKERRQ(ierr); 341 ierr = MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES);CHKERRQ(ierr); 342 ierr = MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES);CHKERRQ(ierr); 343 } 344 ierr = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 345 ierr = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 346 if (!user->noeqflag) { 347 ierr = VecScatterDestroy(&Descat);CHKERRQ(ierr); 348 ierr = VecDestroy(&Deseq);CHKERRQ(ierr); 349 } 350 ierr = VecScatterDestroy(&Discat);CHKERRQ(ierr); 351 ierr = VecDestroy(&Diseq);CHKERRQ(ierr); 352 PetscFunctionReturn(0); 353 } 354 355 /* Evaluate 356 h = [ x0^2 - x1; 357 1 -(x0^2 - x1)] 358 */ 359 PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx) 360 { 361 const PetscScalar *x; 362 PetscScalar ci; 363 PetscErrorCode ierr; 364 MPI_Comm comm; 365 PetscMPIInt rank; 366 AppCtx *user=(AppCtx*)ctx; 367 Vec Xseq=user->Xseq; 368 VecScatter scat=user->scat; 369 370 PetscFunctionBegin; 371 ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 372 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 373 374 ierr = VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 375 ierr = VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 376 377 if (rank == 0) { 378 ierr = VecGetArrayRead(Xseq,&x);CHKERRQ(ierr); 379 ci = x[0]*x[0] - x[1]; 380 ierr = VecSetValue(CI,0,ci,INSERT_VALUES);CHKERRQ(ierr); 381 ci = -x[0]*x[0] + x[1] + 1.0; 382 ierr = VecSetValue(CI,1,ci,INSERT_VALUES);CHKERRQ(ierr); 383 ierr = VecRestoreArrayRead(Xseq,&x);CHKERRQ(ierr); 384 } 385 ierr = VecAssemblyBegin(CI);CHKERRQ(ierr); 386 ierr = VecAssemblyEnd(CI);CHKERRQ(ierr); 387 PetscFunctionReturn(0); 388 } 389 390 /* Evaluate 391 g = [ x0^2 + x1 - 2] 392 */ 393 PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx) 394 { 395 const PetscScalar *x; 396 PetscScalar ce; 397 PetscErrorCode ierr; 398 MPI_Comm comm; 399 PetscMPIInt rank; 400 AppCtx *user=(AppCtx*)ctx; 401 Vec Xseq=user->Xseq; 402 VecScatter scat=user->scat; 403 404 PetscFunctionBegin; 405 ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 406 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 407 408 ierr = VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 409 ierr = VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 410 411 if (rank == 0) { 412 ierr = VecGetArrayRead(Xseq,&x);CHKERRQ(ierr); 413 ce = x[0]*x[0] + x[1] - 2.0; 414 ierr = VecSetValue(CE,0,ce,INSERT_VALUES);CHKERRQ(ierr); 415 ierr = VecRestoreArrayRead(Xseq,&x);CHKERRQ(ierr); 416 } 417 ierr = VecAssemblyBegin(CE);CHKERRQ(ierr); 418 ierr = VecAssemblyEnd(CE);CHKERRQ(ierr); 419 PetscFunctionReturn(0); 420 } 421 422 /* 423 grad h = [ 2*x0, -1; 424 -2*x0, 1] 425 */ 426 PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx) 427 { 428 AppCtx *user=(AppCtx*)ctx; 429 PetscInt zero=0,one=1,cols[2]; 430 PetscScalar vals[2]; 431 const PetscScalar *x; 432 PetscErrorCode ierr; 433 Vec Xseq=user->Xseq; 434 VecScatter scat=user->scat; 435 MPI_Comm comm; 436 PetscMPIInt rank; 437 438 PetscFunctionBegin; 439 ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 440 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 441 ierr = VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 442 ierr = VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 443 444 ierr = VecGetArrayRead(Xseq,&x);CHKERRQ(ierr); 445 if (!rank) { 446 cols[0] = 0; cols[1] = 1; 447 vals[0] = 2*x[0]; vals[1] = -1.0; 448 ierr = MatSetValues(JI,1,&zero,2,cols,vals,INSERT_VALUES);CHKERRQ(ierr); 449 vals[0] = -2*x[0]; vals[1] = 1.0; 450 ierr = MatSetValues(JI,1,&one,2,cols,vals,INSERT_VALUES);CHKERRQ(ierr); 451 } 452 ierr = VecRestoreArrayRead(Xseq,&x);CHKERRQ(ierr); 453 ierr = MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 454 ierr = MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 455 PetscFunctionReturn(0); 456 } 457 458 /* 459 grad g = [2*x0 460 1.0 ] 461 */ 462 PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx) 463 { 464 PetscInt zero=0,cols[2]; 465 PetscScalar vals[2]; 466 const PetscScalar *x; 467 PetscMPIInt rank; 468 MPI_Comm comm; 469 PetscErrorCode ierr; 470 471 PetscFunctionBegin; 472 ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 473 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 474 475 if (rank == 0) { 476 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 477 cols[0] = 0; cols[1] = 1; 478 vals[0] = 2*x[0]; vals[1] = 1.0; 479 ierr = MatSetValues(JE,1,&zero,2,cols,vals,INSERT_VALUES);CHKERRQ(ierr); 480 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 481 } 482 ierr = MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 483 ierr = MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 484 PetscFunctionReturn(0); 485 } 486 487 /*TEST 488 489 build: 490 requires: !complex !defined(PETSC_USE_CXX) 491 492 test: 493 args: -tao_converged_reason -tao_pdipm_kkt_shift_pd 494 requires: mumps 495 496 test: 497 suffix: 2 498 nsize: 2 499 args: -tao_converged_reason -tao_pdipm_kkt_shift_pd 500 requires: mumps 501 502 test: 503 suffix: 3 504 args: -tao_converged_reason -no_eq 505 requires: mumps 506 507 test: 508 suffix: 4 509 nsize: 2 510 args: -tao_converged_reason -no_eq 511 requires: mumps 512 513 test: 514 suffix: 5 515 args: -tao_cmonitor -tao_type almm 516 requires: mumps 517 518 test: 519 suffix: 6 520 args: -tao_cmonitor -tao_type almm -tao_almm_type phr 521 requires: mumps 522 523 test: 524 suffix: 7 525 nsize: 2 526 requires: mumps 527 args: -tao_cmonitor -tao_type almm 528 529 test: 530 suffix: 8 531 nsize: 2 532 requires: cuda mumps 533 args: -tao_cmonitor -tao_type almm -vec_type cuda -mat_type aijcusparse 534 535 test: 536 suffix: 9 537 nsize: 2 538 args: -tao_cmonitor -tao_type almm -no_eq 539 requires: mumps 540 541 test: 542 suffix: 10 543 nsize: 2 544 args: -tao_cmonitor -tao_type almm -tao_almm_type phr -no_eq 545 requires: mumps 546 547 TEST*/ 548