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 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 68 PetscCallMPI(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 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"---- Constrained Problem -----\n")); 72 PetscCall(InitializeProblem(&user)); /* sets up problem, function below */ 73 PetscCall(TaoCreate(PETSC_COMM_WORLD,&tao)); 74 PetscCall(TaoSetType(tao,TAOPDIPM)); 75 PetscCall(TaoSetSolution(tao,user.x)); /* gets solution vector from problem */ 76 PetscCall(TaoSetVariableBounds(tao,user.xl,user.xu)); /* sets lower upper bounds from given solution */ 77 PetscCall(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void*)&user)); 78 79 if (!user.noeqflag) { 80 PetscCall(TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user)); 81 } 82 PetscCall(TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user)); 83 if (!user.noeqflag) { 84 PetscCall(TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user)); /* equality jacobian */ 85 } 86 PetscCall(TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user)); /* inequality jacobian */ 87 PetscCall(TaoSetTolerances(tao,1.e-6,1.e-6,1.e-6)); 88 PetscCall(TaoSetConstraintTolerances(tao,1.e-6,1.e-6)); 89 90 PetscCall(TaoGetKSP(tao,&ksp)); 91 PetscCall(KSPGetPC(ksp,&pc)); 92 PetscCall(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 PetscCall(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 PetscCall(KSPSetType(ksp,KSPPREONLY)); 103 PetscCall(KSPSetFromOptions(ksp)); 104 105 PetscCall(TaoSetFromOptions(tao)); 106 PetscCall(TaoGetType(tao,&type)); 107 PetscCall(PetscObjectTypeCompare((PetscObject)tao,TAOPDIPM,&pdipm)); 108 if (pdipm) { 109 PetscCall(TaoSetHessian(tao,user.H,user.H,FormHessian,(void*)&user)); 110 if (user.initview) { 111 PetscCall(TaoSetUp(tao)); 112 PetscCall(VecDuplicate(user.x, &G)); 113 PetscCall(FormFunctionGradient(tao, user.x, &f, G, (void*)&user)); 114 PetscCall(FormHessian(tao, user.x, user.H, user.H, (void*)&user)); 115 PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD)); 116 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nInitial point X:\n")); 117 PetscCall(VecView(user.x, PETSC_VIEWER_STDOUT_WORLD)); 118 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nInitial objective f(x) = %g\n",(double)f)); 119 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nInitial gradient and Hessian:\n")); 120 PetscCall(VecView(G, PETSC_VIEWER_STDOUT_WORLD)); 121 PetscCall(MatView(user.H, PETSC_VIEWER_STDOUT_WORLD)); 122 PetscCall(VecDestroy(&G)); 123 PetscCall(FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void*)&user)); 124 PetscCall(MatCreateVecs(user.Ai, NULL, &CI)); 125 PetscCall(FormInequalityConstraints(tao, user.x, CI, (void*)&user)); 126 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nInitial inequality constraints and Jacobian:\n")); 127 PetscCall(VecView(CI, PETSC_VIEWER_STDOUT_WORLD)); 128 PetscCall(MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD)); 129 PetscCall(VecDestroy(&CI)); 130 if (!user.noeqflag) { 131 PetscCall(FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void*)&user)); 132 PetscCall(MatCreateVecs(user.Ae, NULL, &CE)); 133 PetscCall(FormEqualityConstraints(tao, user.x, CE, (void*)&user)); 134 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nInitial equality constraints and Jacobian:\n")); 135 PetscCall(VecView(CE, PETSC_VIEWER_STDOUT_WORLD)); 136 PetscCall(MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD)); 137 PetscCall(VecDestroy(&CE)); 138 } 139 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n")); 140 PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD)); 141 } 142 } 143 144 PetscCall(TaoSolve(tao)); 145 PetscCall(TaoGetSolution(tao,&x)); 146 PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); 147 148 /* Free objects */ 149 PetscCall(DestroyProblem(&user)); 150 PetscCall(TaoDestroy(&tao)); 151 PetscCall(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 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 163 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 164 user->noeqflag = PETSC_FALSE; 165 PetscCall(PetscOptionsGetBool(NULL,NULL,"-no_eq",&user->noeqflag,NULL)); 166 user->initview = PETSC_FALSE; 167 PetscCall(PetscOptionsGetBool(NULL,NULL,"-init_view",&user->initview,NULL)); 168 169 if (!user->noeqflag) { 170 /* Tell user the correct solution, not an error checking */ 171 PetscCall(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 PetscCall(VecCreate(PETSC_COMM_WORLD,&user->x)); 178 PetscCall(VecSetSizes(user->x,nloc,user->n)); 179 PetscCall(VecSetFromOptions(user->x)); 180 PetscCall(VecSet(user->x,0)); 181 182 /* create and set lower and upper bound vectors */ 183 PetscCall(VecDuplicate(user->x,&user->xl)); 184 PetscCall(VecDuplicate(user->x,&user->xu)); 185 PetscCall(VecSet(user->xl,-1.0)); 186 PetscCall(VecSet(user->xu,2.0)); 187 188 /* create scater to zero */ 189 PetscCall(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 PetscCall(VecCreate(PETSC_COMM_WORLD,&user->ce)); /* a 1x1 vec for equality constraints */ 198 PetscCall(VecSetSizes(user->ce,neloc,user->ne)); 199 PetscCall(VecSetFromOptions(user->ce)); 200 PetscCall(VecSetUp(user->ce)); 201 } 202 203 PetscCall(VecCreate(PETSC_COMM_WORLD,&user->ci)); /* a 2x1 vec for inequality constraints */ 204 PetscCall(VecSetSizes(user->ci,niloc,user->ni)); 205 PetscCall(VecSetFromOptions(user->ci)); 206 PetscCall(VecSetUp(user->ci)); 207 208 /* nexn & nixn matricies for equally and inequalty constraints */ 209 if (!user->noeqflag) { 210 PetscCall(MatCreate(PETSC_COMM_WORLD,&user->Ae)); 211 PetscCall(MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n)); 212 PetscCall(MatSetFromOptions(user->Ae)); 213 PetscCall(MatSetUp(user->Ae)); 214 } 215 216 PetscCall(MatCreate(PETSC_COMM_WORLD,&user->Ai)); 217 PetscCall(MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n)); 218 PetscCall(MatSetFromOptions(user->Ai)); 219 PetscCall(MatSetUp(user->Ai)); 220 221 PetscCall(MatCreate(PETSC_COMM_WORLD,&user->H)); 222 PetscCall(MatSetSizes(user->H,nloc,nloc,user->n,user->n)); 223 PetscCall(MatSetFromOptions(user->H)); 224 PetscCall(MatSetUp(user->H)); 225 PetscFunctionReturn(0); 226 } 227 228 PetscErrorCode DestroyProblem(AppCtx *user) 229 { 230 PetscFunctionBegin; 231 if (!user->noeqflag) { 232 PetscCall(MatDestroy(&user->Ae)); 233 } 234 PetscCall(MatDestroy(&user->Ai)); 235 PetscCall(MatDestroy(&user->H)); 236 237 PetscCall(VecDestroy(&user->x)); 238 if (!user->noeqflag) { 239 PetscCall(VecDestroy(&user->ce)); 240 } 241 PetscCall(VecDestroy(&user->ci)); 242 PetscCall(VecDestroy(&user->xl)); 243 PetscCall(VecDestroy(&user->xu)); 244 PetscCall(VecDestroy(&user->Xseq)); 245 PetscCall(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 PetscCall(PetscObjectGetComm((PetscObject)tao,&comm)); 267 PetscCallMPI(MPI_Comm_rank(comm,&rank)); 268 269 PetscCall(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD)); 270 PetscCall(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD)); 271 272 fin = 0.0; 273 if (rank == 0) { 274 PetscCall(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 PetscCall(VecSetValue(G,0,g,INSERT_VALUES)); 278 g = 2.0*(x[1]-2.0) - 2.0; 279 PetscCall(VecSetValue(G,1,g,INSERT_VALUES)); 280 PetscCall(VecRestoreArrayRead(Xseq,&x)); 281 } 282 PetscCallMPI(MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm)); 283 PetscCall(VecAssemblyBegin(G)); 284 PetscCall(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 PetscCall(TaoGetDualVariables(tao,&DE,&DI)); 308 309 PetscCall(PetscObjectGetComm((PetscObject)tao,&comm)); 310 PetscCallMPI(MPI_Comm_rank(comm,&rank)); 311 312 if (!user->noeqflag) { 313 PetscCall(VecScatterCreateToZero(DE,&Descat,&Deseq)); 314 PetscCall(VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD)); 315 PetscCall(VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD)); 316 } 317 PetscCall(VecScatterCreateToZero(DI,&Discat,&Diseq)); 318 PetscCall(VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD)); 319 PetscCall(VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD)); 320 321 if (rank == 0) { 322 if (!user->noeqflag) { 323 PetscCall(VecGetArrayRead(Deseq,&de)); /* places equality constraint dual into array */ 324 } 325 PetscCall(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 PetscCall(VecRestoreArrayRead(Deseq,&de)); 330 PetscCall(VecRestoreArrayRead(Diseq,&di)); 331 } else { 332 val = 2.0 * (1 - di[0] + di[1]); 333 } 334 PetscCall(VecRestoreArrayRead(Diseq,&di)); 335 PetscCall(MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES)); 336 PetscCall(MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES)); 337 } 338 PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY)); 339 PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY)); 340 if (!user->noeqflag) { 341 PetscCall(VecScatterDestroy(&Descat)); 342 PetscCall(VecDestroy(&Deseq)); 343 } 344 PetscCall(VecScatterDestroy(&Discat)); 345 PetscCall(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 PetscCall(PetscObjectGetComm((PetscObject)tao,&comm)); 365 PetscCallMPI(MPI_Comm_rank(comm,&rank)); 366 367 PetscCall(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD)); 368 PetscCall(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD)); 369 370 if (rank == 0) { 371 PetscCall(VecGetArrayRead(Xseq,&x)); 372 ci = x[0]*x[0] - x[1]; 373 PetscCall(VecSetValue(CI,0,ci,INSERT_VALUES)); 374 ci = -x[0]*x[0] + x[1] + 1.0; 375 PetscCall(VecSetValue(CI,1,ci,INSERT_VALUES)); 376 PetscCall(VecRestoreArrayRead(Xseq,&x)); 377 } 378 PetscCall(VecAssemblyBegin(CI)); 379 PetscCall(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 PetscCall(PetscObjectGetComm((PetscObject)tao,&comm)); 398 PetscCallMPI(MPI_Comm_rank(comm,&rank)); 399 400 PetscCall(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD)); 401 PetscCall(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD)); 402 403 if (rank == 0) { 404 PetscCall(VecGetArrayRead(Xseq,&x)); 405 ce = x[0]*x[0] + x[1] - 2.0; 406 PetscCall(VecSetValue(CE,0,ce,INSERT_VALUES)); 407 PetscCall(VecRestoreArrayRead(Xseq,&x)); 408 } 409 PetscCall(VecAssemblyBegin(CE)); 410 PetscCall(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 PetscCall(PetscObjectGetComm((PetscObject)tao,&comm)); 431 PetscCallMPI(MPI_Comm_rank(comm,&rank)); 432 PetscCall(VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD)); 433 PetscCall(VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD)); 434 435 PetscCall(VecGetArrayRead(Xseq,&x)); 436 if (rank == 0) { 437 cols[0] = 0; cols[1] = 1; 438 vals[0] = 2*x[0]; vals[1] = -1.0; 439 PetscCall(MatSetValues(JI,1,&zero,2,cols,vals,INSERT_VALUES)); 440 vals[0] = -2*x[0]; vals[1] = 1.0; 441 PetscCall(MatSetValues(JI,1,&one,2,cols,vals,INSERT_VALUES)); 442 } 443 PetscCall(VecRestoreArrayRead(Xseq,&x)); 444 PetscCall(MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY)); 445 PetscCall(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 PetscCall(PetscObjectGetComm((PetscObject)tao,&comm)); 463 PetscCallMPI(MPI_Comm_rank(comm,&rank)); 464 465 if (rank == 0) { 466 PetscCall(VecGetArrayRead(X,&x)); 467 cols[0] = 0; cols[1] = 1; 468 vals[0] = 2*x[0]; vals[1] = 1.0; 469 PetscCall(MatSetValues(JE,1,&zero,2,cols,vals,INSERT_VALUES)); 470 PetscCall(VecRestoreArrayRead(X,&x)); 471 } 472 PetscCall(MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY)); 473 PetscCall(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