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