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