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