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