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