1 /* Program usage: mpiexec -n 2 ./ex1 [-help] [all TAO options] */ 2 3 /* ---------------------------------------------------------------------- 4 min f = (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 10 #include <petsctao.h> 11 12 static char help[]= "Solves constrained optimiztion problem using pdipm.\n\ 13 Input parameters include:\n\ 14 -tao_type pdipm : sets Tao solver\n\ 15 -no_eq : removes the equaility constraints from the problem\n\ 16 -snes_fd : snes with finite difference Jacobian (needed for pdipm)\n\ 17 -tao_cmonitor : convergence monitor with constraint norm \n\ 18 -tao_view_solution : view exact solution at each itteration\n\ 19 Note: external package superlu_dist is requried to run either for ipm or pdipm. Additionally This is designed for a maximum of 2 processors, the code will error if size > 2.\n\n"; 20 21 /* 22 User-defined application context - contains data needed by the 23 application-provided call-back routines, FormFunction(), 24 FormGradient(), and FormHessian(). 25 */ 26 27 /* 28 x,d in R^n 29 f in R 30 bin in R^mi 31 beq in R^me 32 Aeq in R^(me x n) 33 Ain in R^(mi x n) 34 H in R^(n x n) 35 min f=(1/2)*x'*H*x + d'*x 36 s.t. Aeq*x == beq 37 Ain*x >= bin 38 */ 39 typedef struct { 40 PetscInt n; /* Global length of x */ 41 PetscInt ne; /* Global number of equality constraints */ 42 PetscInt ni; /* Global number of inequality constraints */ 43 PetscBool noeqflag; 44 Vec x,xl,xu; 45 Vec ce,ci,bl,bu,Xseq; 46 Mat Ae,Ai,H; 47 VecScatter scat; 48 } AppCtx; 49 50 /* -------- User-defined Routines --------- */ 51 PetscErrorCode InitializeProblem(AppCtx *); 52 PetscErrorCode DestroyProblem(AppCtx *); 53 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *); 54 PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*); 55 PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*); 56 PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*); 57 PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*); 58 PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*); 59 60 PetscErrorCode main(int argc,char **argv) 61 { 62 PetscErrorCode ierr; /* used to check for functions returning nonzeros */ 63 Tao tao; 64 KSP ksp; 65 PC pc; 66 AppCtx user; /* application context */ 67 Vec x; 68 PetscMPIInt rank; 69 70 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 71 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 72 if (rank>1){ 73 SETERRQ(PETSC_COMM_SELF,1,"More than 2 processors detected. Example written to use max of 2 processors.\n"); 74 } 75 76 ierr = PetscPrintf(PETSC_COMM_WORLD,"---- Constrained Problem -----\n");CHKERRQ(ierr); 77 ierr = InitializeProblem(&user);CHKERRQ(ierr); /* sets up problem, function below */ 78 ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr); 79 ierr = TaoSetType(tao,TAOPDIPM);CHKERRQ(ierr); 80 ierr = TaoSetInitialVector(tao,user.x);CHKERRQ(ierr); /* gets solution vector from problem */ 81 ierr = TaoSetVariableBounds(tao,user.xl,user.xu);CHKERRQ(ierr); /* sets lower upper bounds from given solution */ 82 ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);CHKERRQ(ierr); 83 84 if (!user.noeqflag){ 85 ierr = TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user);CHKERRQ(ierr); 86 } 87 ierr = TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user);CHKERRQ(ierr); 88 89 if (!user.noeqflag){ 90 ierr = TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user);CHKERRQ(ierr); /* equality jacobian */ 91 } 92 ierr = TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user);CHKERRQ(ierr); /* inequality jacobian */ 93 ierr = TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);CHKERRQ(ierr); 94 ierr = TaoSetTolerances(tao,1.e-6,1.e-6,1.e-6);CHKERRQ(ierr); 95 ierr = TaoSetConstraintTolerances(tao,1.e-6,1.e-6);CHKERRQ(ierr); 96 97 ierr = TaoGetKSP(tao,&ksp);CHKERRQ(ierr); 98 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 99 ierr = PCSetType(pc,PCLU);CHKERRQ(ierr); 100 /* 101 This algorithm produces matrices with zeros along the diagonal therefore we use 102 SuperLU_DIST which provides shift to the diagonal 103 */ 104 ierr = PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU_DIST);CHKERRQ(ierr); /* requires superlu_dist to solve pdipm */ 105 ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr); 106 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 107 108 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr); 109 110 ierr = TaoSolve(tao);CHKERRQ(ierr); 111 ierr = TaoGetSolutionVector(tao,&x);CHKERRQ(ierr); 112 ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 113 114 /* Free objects */ 115 ierr = DestroyProblem(&user);CHKERRQ(ierr); 116 ierr = TaoDestroy(&tao);CHKERRQ(ierr); 117 ierr = PetscFinalize(); 118 return ierr; 119 } 120 121 PetscErrorCode InitializeProblem(AppCtx *user) 122 { 123 PetscErrorCode ierr; 124 PetscMPIInt size; 125 PetscMPIInt rank; 126 PetscInt nloc,neloc,niloc; 127 128 PetscFunctionBegin; 129 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 130 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 131 user->noeqflag = PETSC_FALSE; 132 ierr = PetscOptionsGetBool(NULL,NULL,"-no_eq",&user->noeqflag,NULL);CHKERRQ(ierr); 133 if (!user->noeqflag) { 134 ierr = PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n");CHKERRQ(ierr); 135 } 136 137 /* create vector x and set initial values */ 138 user->n = 2; /* global length */ 139 nloc = (rank==0)?user->n:0; 140 ierr = VecCreateMPI(PETSC_COMM_WORLD,nloc,user->n,&user->x);CHKERRQ(ierr); 141 ierr = VecSet(user->x,0);CHKERRQ(ierr); 142 143 /* create and set lower and upper bound vectors */ 144 ierr = VecDuplicate(user->x,&user->xl);CHKERRQ(ierr); 145 ierr = VecDuplicate(user->x,&user->xu);CHKERRQ(ierr); 146 ierr = VecSet(user->xl,-1.0);CHKERRQ(ierr); 147 ierr = VecSet(user->xu,2.0);CHKERRQ(ierr); 148 149 /* create scater to zero */ 150 ierr = VecScatterCreateToZero(user->x,&user->scat,&user->Xseq);CHKERRQ(ierr); 151 152 user->ne = 1; 153 neloc = (rank==0)?user->ne:0; 154 155 if (!user->noeqflag){ 156 ierr = VecCreateMPI(PETSC_COMM_WORLD,neloc,user->ne,&user->ce);CHKERRQ(ierr); /* a 1x1 vec for equality constraints */ 157 } 158 user->ni = 2; 159 niloc = (rank==0)?user->ni:0; 160 ierr = VecCreateMPI(PETSC_COMM_WORLD,niloc,user->ni,&user->ci);CHKERRQ(ierr); /* a 2x1 vec for inequality constraints */ 161 162 /* nexn & nixn matricies for equaly and inequalty constriants */ 163 if (!user->noeqflag){ 164 ierr = MatCreate(PETSC_COMM_WORLD,&user->Ae);CHKERRQ(ierr); 165 ierr = MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n);CHKERRQ(ierr); 166 ierr = MatSetUp(user->Ae);CHKERRQ(ierr); 167 ierr = MatSetFromOptions(user->Ae);CHKERRQ(ierr); 168 } 169 170 ierr = MatCreate(PETSC_COMM_WORLD,&user->Ai);CHKERRQ(ierr); 171 ierr = MatCreate(PETSC_COMM_WORLD,&user->H);CHKERRQ(ierr); 172 173 ierr = MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n);CHKERRQ(ierr); 174 ierr = MatSetSizes(user->H,nloc,nloc,user->n,user->n);CHKERRQ(ierr); 175 176 ierr = MatSetUp(user->Ai);CHKERRQ(ierr); 177 ierr = MatSetUp(user->H);CHKERRQ(ierr); 178 179 ierr = MatSetFromOptions(user->Ai);CHKERRQ(ierr); 180 ierr = MatSetFromOptions(user->H);CHKERRQ(ierr); 181 PetscFunctionReturn(0); 182 } 183 184 PetscErrorCode DestroyProblem(AppCtx *user) 185 { 186 PetscErrorCode ierr; 187 188 PetscFunctionBegin; 189 if (!user->noeqflag){ 190 ierr = MatDestroy(&user->Ae);CHKERRQ(ierr); 191 } 192 ierr = MatDestroy(&user->Ai);CHKERRQ(ierr); 193 ierr = MatDestroy(&user->H);CHKERRQ(ierr); 194 195 ierr = VecDestroy(&user->x);CHKERRQ(ierr); 196 if (!user->noeqflag){ 197 ierr = VecDestroy(&user->ce);CHKERRQ(ierr); 198 } 199 ierr = VecDestroy(&user->ci);CHKERRQ(ierr); 200 ierr = VecDestroy(&user->xl);CHKERRQ(ierr); 201 ierr = VecDestroy(&user->xu);CHKERRQ(ierr); 202 ierr = VecDestroy(&user->Xseq);CHKERRQ(ierr); 203 ierr = VecScatterDestroy(&user->scat);CHKERRQ(ierr); 204 PetscFunctionReturn(0); 205 } 206 207 /* 208 f(X) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1) 209 G(X) = fx = [2.0*(x[0]-2.0) - 2.0; 2.0*(x[1]-2.0) - 2.0] 210 */ 211 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx) 212 { 213 PetscScalar g; 214 const PetscScalar *x; 215 MPI_Comm comm; 216 PetscMPIInt rank; 217 PetscErrorCode ierr; 218 PetscReal fin; 219 AppCtx *user=(AppCtx*)ctx; 220 Vec Xseq=user->Xseq; 221 VecScatter scat=user->scat; 222 223 PetscFunctionBegin; 224 ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 225 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 226 227 ierr = VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 228 ierr = VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 229 230 fin = 0.0; 231 if (!rank) { 232 ierr = VecGetArrayRead(Xseq,&x);CHKERRQ(ierr); 233 fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]); 234 g = 2.0*(x[0]-2.0) - 2.0; 235 ierr = VecSetValue(G,0,g,INSERT_VALUES);CHKERRQ(ierr); 236 g = 2.0*(x[1]-2.0) - 2.0; 237 ierr = VecSetValue(G,1,g,INSERT_VALUES);CHKERRQ(ierr); 238 ierr = VecRestoreArrayRead(Xseq,&x);CHKERRQ(ierr); 239 } 240 ierr = MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm);CHKERRQ(ierr); 241 ierr = VecAssemblyBegin(G);CHKERRQ(ierr); 242 ierr = VecAssemblyEnd(G);CHKERRQ(ierr); 243 PetscFunctionReturn(0); 244 } 245 246 PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx) 247 { 248 AppCtx *user=(AppCtx*)ctx; 249 Vec DE,DI; 250 const PetscScalar *de, *di; 251 PetscInt zero=0,one=1; 252 PetscScalar two=2.0; 253 PetscScalar val=0.0; 254 Vec Deseq,Diseq=user->Xseq; 255 VecScatter Descat,Discat=user->scat; 256 PetscMPIInt rank; 257 MPI_Comm comm; 258 PetscErrorCode ierr; 259 260 PetscFunctionBegin; 261 ierr = TaoGetDualVariables(tao,&DE,&DI);CHKERRQ(ierr); 262 263 ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 264 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 265 266 if (!user->noeqflag){ 267 ierr = VecScatterCreateToZero(DE,&Descat,&Deseq);CHKERRQ(ierr); 268 ierr = VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 269 ierr = VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 270 } 271 272 ierr = VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 273 ierr = VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 274 275 276 if (!rank){ 277 if (!user->noeqflag){ 278 ierr = VecGetArrayRead(Deseq,&de);CHKERRQ(ierr); /* places equality constraint dual into array */ 279 } 280 281 ierr = VecGetArrayRead(Diseq,&di);CHKERRQ(ierr); /* places inequality constraint dual into array */ 282 283 if (!user->noeqflag){ 284 val = 2.0 * (1 + de[0] + di[0] - di[1]); 285 ierr = VecRestoreArrayRead(Deseq,&de);CHKERRQ(ierr); 286 }else{ 287 val = 2.0 * (1 + di[0] - di[1]); 288 } 289 ierr = VecRestoreArrayRead(Diseq,&di);CHKERRQ(ierr); 290 ierr = MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES);CHKERRQ(ierr); 291 ierr = MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES);CHKERRQ(ierr); 292 } 293 294 ierr = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 295 ierr = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 296 if (!user->noeqflag){ 297 ierr = VecScatterDestroy(&Descat);CHKERRQ(ierr); 298 ierr = VecDestroy(&Deseq);CHKERRQ(ierr); 299 } 300 PetscFunctionReturn(0); 301 } 302 303 PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx) 304 { 305 const PetscScalar *x; 306 PetscScalar ci; 307 PetscErrorCode ierr; 308 MPI_Comm comm; 309 PetscMPIInt rank; 310 AppCtx *user=(AppCtx*)ctx; 311 Vec Xseq=user->Xseq; 312 VecScatter scat=user->scat; 313 314 PetscFunctionBegin; 315 ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 316 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 317 318 ierr = VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 319 ierr = VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 320 321 if (!rank) { 322 ierr = VecGetArrayRead(Xseq,&x);CHKERRQ(ierr); 323 ci = x[0]*x[0] - x[1]; 324 ierr = VecSetValue(CI,0,ci,INSERT_VALUES);CHKERRQ(ierr); 325 ci = -x[0]*x[0] + x[1] + 1.0; 326 ierr = VecSetValue(CI,1,ci,INSERT_VALUES);CHKERRQ(ierr); 327 ierr = VecRestoreArrayRead(Xseq,&x);CHKERRQ(ierr); 328 } 329 ierr = VecAssemblyBegin(CI);CHKERRQ(ierr); 330 ierr = VecAssemblyEnd(CI);CHKERRQ(ierr); 331 PetscFunctionReturn(0); 332 } 333 334 PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx) 335 { 336 const PetscScalar *x; 337 PetscScalar ce; 338 PetscErrorCode ierr; 339 MPI_Comm comm; 340 PetscMPIInt rank; 341 AppCtx *user=(AppCtx*)ctx; 342 Vec Xseq=user->Xseq; 343 VecScatter scat=user->scat; 344 345 PetscFunctionBegin; 346 ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 347 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 348 349 ierr = VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 350 ierr = VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 351 352 if (!rank) { 353 ierr = VecGetArrayRead(Xseq,&x);CHKERRQ(ierr); 354 ce = x[0]*x[0] + x[1] - 2.0; 355 ierr = VecSetValue(CE,0,ce,INSERT_VALUES);CHKERRQ(ierr); 356 ierr = VecRestoreArrayRead(Xseq,&x);CHKERRQ(ierr); 357 } 358 ierr = VecAssemblyBegin(CE);CHKERRQ(ierr); 359 ierr = VecAssemblyEnd(CE);CHKERRQ(ierr); 360 PetscFunctionReturn(0); 361 } 362 363 PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx) 364 { 365 AppCtx *user=(AppCtx*)ctx; 366 PetscInt cols[2],min,max,i; 367 PetscScalar vals[2]; 368 const PetscScalar *x; 369 PetscErrorCode ierr; 370 Vec Xseq=user->Xseq; 371 VecScatter scat=user->scat; 372 MPI_Comm comm; 373 PetscMPIInt rank; 374 375 PetscFunctionBegin; 376 ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 377 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 378 ierr = VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 379 ierr = VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 380 381 ierr = VecGetArrayRead(Xseq,&x);CHKERRQ(ierr); 382 ierr = MatGetOwnershipRange(JI,&min,&max);CHKERRQ(ierr); 383 384 cols[0] = 0; cols[1] = 1; 385 for (i=min;i<max;i++) { 386 if (i==0){ 387 vals[0] = +2*x[0]; vals[1] = -1.0; 388 ierr = MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);CHKERRQ(ierr); 389 } 390 if (i==1) { 391 vals[0] = -2*x[0]; vals[1] = +1.0; 392 ierr = MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);CHKERRQ(ierr); 393 } 394 } 395 ierr = VecRestoreArrayRead(Xseq,&x);CHKERRQ(ierr); 396 397 ierr = MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 398 ierr = MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 399 PetscFunctionReturn(0); 400 } 401 402 PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx) 403 { 404 PetscInt rows[2]; 405 PetscScalar vals[2]; 406 const PetscScalar *x; 407 PetscMPIInt rank; 408 MPI_Comm comm; 409 PetscErrorCode ierr; 410 411 PetscFunctionBegin; 412 ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 413 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 414 415 if (!rank) { 416 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 417 rows[0] = 0; rows[1] = 1; 418 vals[0] = 2*x[0]; vals[1] = 1.0; 419 ierr = MatSetValues(JE,1,rows,2,rows,vals,INSERT_VALUES);CHKERRQ(ierr); 420 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 421 } 422 ierr = MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 423 ierr = MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 424 PetscFunctionReturn(0); 425 } 426 427 428 /*TEST 429 430 build: 431 requires: !complex !define(PETSC_USE_CXX) 432 433 test: 434 requires: superlu_dist 435 args: -tao_converged_reason 436 437 test: 438 suffix: 2 439 requires: superlu_dist 440 nsize: 2 441 args: -tao_converged_reason 442 443 test: 444 suffix: 3 445 requires: superlu_dist 446 args: -tao_converged_reason -no_eq 447 448 test: 449 suffix: 4 450 requires: superlu_dist 451 nsize: 2 452 args: -tao_converged_reason -no_eq 453 454 TEST*/ 455