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 PetscReal f; 65 PetscBool pdipm; 66 67 PetscFunctionBeginUser; 68 PetscCall(PetscInitialize(&argc, &argv, NULL, 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 processes detected. Example written to use at most 2 processes."); 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 if (!user.noboundflag) 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(PetscObjectTypeCompare((PetscObject)tao, TAOPDIPM, &pdipm)); 93 if (pdipm) { 94 /* 95 PDIPM produces an indefinite KKT matrix with some zeros along the diagonal 96 Inverting this indefinite matrix requires PETSc to be configured with MUMPS 97 */ 98 PetscCheck(PetscDefined(HAVE_MUMPS), PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "TAOPDIPM requires PETSc to be configured with MUMPS (--download-mumps)"); 99 PetscCall(TaoGetKSP(tao, &ksp)); 100 PetscCall(KSPGetPC(ksp, &pc)); 101 PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS)); 102 PetscCall(KSPSetType(ksp, KSPPREONLY)); 103 PetscCall(PCSetType(pc, PCCHOLESKY)); 104 PetscCall(KSPSetFromOptions(ksp)); 105 PetscCall(TaoSetHessian(tao, user.H, user.H, FormPDIPMHessian, (void *)&user)); 106 } 107 108 /* Print out an initial view of the problem */ 109 if (user.initview) { 110 PetscCall(TaoSetUp(tao)); 111 PetscCall(VecDuplicate(user.x, &G)); 112 PetscCall(FormFunctionGradient(tao, user.x, &f, G, (void *)&user)); 113 if (pdipm) PetscCall(FormPDIPMHessian(tao, user.x, user.H, user.H, (void *)&user)); 114 PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD)); 115 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial point X:\n")); 116 PetscCall(VecView(user.x, PETSC_VIEWER_STDOUT_WORLD)); 117 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial objective f(x) = %g\n", (double)f)); 118 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial gradient and Hessian:\n")); 119 PetscCall(VecView(G, PETSC_VIEWER_STDOUT_WORLD)); 120 if (pdipm) PetscCall(MatView(user.H, PETSC_VIEWER_STDOUT_WORLD)); 121 PetscCall(VecDestroy(&G)); 122 PetscCall(FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void *)&user)); 123 PetscCall(MatCreateVecs(user.Ai, NULL, &CI)); 124 PetscCall(FormInequalityConstraints(tao, user.x, CI, (void *)&user)); 125 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial inequality constraints and Jacobian:\n")); 126 PetscCall(VecView(CI, PETSC_VIEWER_STDOUT_WORLD)); 127 PetscCall(MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD)); 128 PetscCall(VecDestroy(&CI)); 129 if (!user.noeqflag) { 130 PetscCall(FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void *)&user)); 131 PetscCall(MatCreateVecs(user.Ae, NULL, &CE)); 132 PetscCall(FormEqualityConstraints(tao, user.x, CE, (void *)&user)); 133 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial equality constraints and Jacobian:\n")); 134 PetscCall(VecView(CE, PETSC_VIEWER_STDOUT_WORLD)); 135 PetscCall(MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD)); 136 PetscCall(VecDestroy(&CE)); 137 } 138 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n")); 139 PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD)); 140 } 141 142 PetscCall(TaoSolve(tao)); 143 PetscCall(TaoGetSolution(tao, &x)); 144 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Found solution:\n")); 145 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD)); 146 147 /* Free objects */ 148 PetscCall(DestroyProblem(&user)); 149 PetscCall(TaoDestroy(&tao)); 150 PetscCall(PetscFinalize()); 151 return PETSC_SUCCESS; 152 } 153 154 PetscErrorCode InitializeProblem(AppCtx *user) 155 { 156 PetscMPIInt size; 157 PetscMPIInt rank; 158 PetscInt nloc, neloc, niloc; 159 160 PetscFunctionBegin; 161 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 162 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 163 user->noeqflag = PETSC_FALSE; 164 user->noboundflag = PETSC_FALSE; 165 PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_eq", &user->noeqflag, NULL)); 166 PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_bound", &user->noboundflag, 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 /* With equality constraint, bound or no bound makes no different to the solution */ 173 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution: f(1, 1) = -2\n")); 174 } else { 175 if (user->noboundflag) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution (-no_eq, -no_bound): f(2.05655, 3.22938) = -9.05728\n")); 176 else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution (-no_eq): f(1.73205, 2) = -7.3923\n")); 177 } 178 179 /* create vector x and set initial values */ 180 user->n = 2; /* global length */ 181 nloc = (size == 1) ? user->n : 1; 182 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->x)); 183 PetscCall(VecSetSizes(user->x, nloc, user->n)); 184 PetscCall(VecSetFromOptions(user->x)); 185 PetscCall(VecSet(user->x, 0.0)); 186 187 /* create and set lower and upper bound vectors */ 188 PetscCall(VecDuplicate(user->x, &user->xl)); 189 PetscCall(VecDuplicate(user->x, &user->xu)); 190 PetscCall(VecSet(user->xl, -1.0)); 191 PetscCall(VecSet(user->xu, 2.0)); 192 193 /* create scater to zero */ 194 PetscCall(VecScatterCreateToZero(user->x, &user->scat, &user->Xseq)); 195 196 user->ne = 1; 197 user->ni = 2; 198 neloc = (rank == 0) ? user->ne : 0; 199 niloc = (size == 1) ? user->ni : 1; 200 201 if (!user->noeqflag) { 202 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ce)); /* a 1x1 vec for equality constraints */ 203 PetscCall(VecSetSizes(user->ce, neloc, user->ne)); 204 PetscCall(VecSetFromOptions(user->ce)); 205 PetscCall(VecSetUp(user->ce)); 206 } 207 208 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ci)); /* a 2x1 vec for inequality constraints */ 209 PetscCall(VecSetSizes(user->ci, niloc, user->ni)); 210 PetscCall(VecSetFromOptions(user->ci)); 211 PetscCall(VecSetUp(user->ci)); 212 213 /* nexn & nixn matrices for equally and inequalty constraints */ 214 if (!user->noeqflag) { 215 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ae)); 216 PetscCall(MatSetSizes(user->Ae, neloc, nloc, user->ne, user->n)); 217 PetscCall(MatSetFromOptions(user->Ae)); 218 PetscCall(MatSetUp(user->Ae)); 219 } 220 221 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ai)); 222 PetscCall(MatSetSizes(user->Ai, niloc, nloc, user->ni, user->n)); 223 PetscCall(MatSetFromOptions(user->Ai)); 224 PetscCall(MatSetUp(user->Ai)); 225 226 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->H)); 227 PetscCall(MatSetSizes(user->H, nloc, nloc, user->n, user->n)); 228 PetscCall(MatSetFromOptions(user->H)); 229 PetscCall(MatSetUp(user->H)); 230 PetscFunctionReturn(PETSC_SUCCESS); 231 } 232 233 PetscErrorCode DestroyProblem(AppCtx *user) 234 { 235 PetscFunctionBegin; 236 if (!user->noeqflag) PetscCall(MatDestroy(&user->Ae)); 237 PetscCall(MatDestroy(&user->Ai)); 238 PetscCall(MatDestroy(&user->H)); 239 240 PetscCall(VecDestroy(&user->x)); 241 if (!user->noeqflag) PetscCall(VecDestroy(&user->ce)); 242 PetscCall(VecDestroy(&user->ci)); 243 PetscCall(VecDestroy(&user->xl)); 244 PetscCall(VecDestroy(&user->xu)); 245 PetscCall(VecDestroy(&user->Xseq)); 246 PetscCall(VecScatterDestroy(&user->scat)); 247 PetscFunctionReturn(PETSC_SUCCESS); 248 } 249 250 /* Evaluate 251 f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1) 252 G = grad f = [2*(x0 - 2) - 2, 2*(x1 - 2) - 2] = 2*X - 6 253 */ 254 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, PetscCtx ctx) 255 { 256 const PetscScalar *x; 257 MPI_Comm comm; 258 PetscMPIInt rank; 259 PetscReal fin; 260 AppCtx *user = (AppCtx *)ctx; 261 Vec Xseq = user->Xseq; 262 VecScatter scat = user->scat; 263 264 PetscFunctionBegin; 265 /* f = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1) */ 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 PetscCall(VecRestoreArrayRead(Xseq, &x)); 277 } 278 PetscCallMPI(MPIU_Allreduce(&fin, f, 1, MPIU_REAL, MPIU_SUM, comm)); 279 280 /* G = 2*X - 6 */ 281 PetscCall(VecSet(G, -6.0)); 282 PetscCall(VecAXPY(G, 2.0, X)); 283 PetscFunctionReturn(PETSC_SUCCESS); 284 } 285 286 /* Evaluate PDIPM Hessian, see Eqn(22) in http://doi.org/10.1049/gtd2.12708 287 H = Wxx = fxx + Jacobian(grad g^T*DE) - Jacobian(grad h^T*DI)] 288 = fxx + Jacobin([2*x0; 1]de[0]) - Jacobian([2*x0, -2*x0; -1, 1][di[0] di[1]]^T) 289 = [2 0; 0 2] + [2*de[0] 0; 0 0] - [2*di[0]-2*di[1], 0; 0, 0] 290 = [ 2*(1+de[0]-di[0]+di[1]), 0; 291 0, 2] 292 */ 293 PetscErrorCode FormPDIPMHessian(Tao tao, Vec x, Mat H, Mat Hpre, PetscCtx 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) PetscCall(VecGetArrayRead(Deseq, &de)); /* places equality constraint dual into array */ 323 PetscCall(VecGetArrayRead(Diseq, &di)); /* places inequality constraint dual into array */ 324 325 if (!user->noeqflag) { 326 val = 2.0 * (1 + de[0] - di[0] + di[1]); 327 PetscCall(VecRestoreArrayRead(Deseq, &de)); 328 PetscCall(VecRestoreArrayRead(Diseq, &di)); 329 } else { 330 val = 2.0 * (1 - di[0] + di[1]); 331 } 332 PetscCall(VecRestoreArrayRead(Diseq, &di)); 333 PetscCall(MatSetValues(H, 1, &zero, 1, &zero, &val, INSERT_VALUES)); 334 PetscCall(MatSetValues(H, 1, &one, 1, &one, &two, INSERT_VALUES)); 335 } 336 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY)); 337 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY)); 338 if (!user->noeqflag) { 339 PetscCall(VecScatterDestroy(&Descat)); 340 PetscCall(VecDestroy(&Deseq)); 341 } 342 PetscCall(VecScatterDestroy(&Discat)); 343 PetscCall(VecDestroy(&Diseq)); 344 PetscFunctionReturn(PETSC_SUCCESS); 345 } 346 347 /* Evaluate 348 h = [ x0^2 - x1; 349 1 -(x0^2 - x1)] 350 */ 351 PetscErrorCode FormInequalityConstraints(Tao tao, Vec X, Vec CI, PetscCtx ctx) 352 { 353 const PetscScalar *x; 354 PetscScalar ci; 355 MPI_Comm comm; 356 PetscMPIInt rank; 357 AppCtx *user = (AppCtx *)ctx; 358 Vec Xseq = user->Xseq; 359 VecScatter scat = user->scat; 360 361 PetscFunctionBegin; 362 PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 363 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 364 365 PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 366 PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 367 368 if (rank == 0) { 369 PetscCall(VecGetArrayRead(Xseq, &x)); 370 ci = x[0] * x[0] - x[1]; 371 PetscCall(VecSetValue(CI, 0, ci, INSERT_VALUES)); 372 ci = -x[0] * x[0] + x[1] + 1.0; 373 PetscCall(VecSetValue(CI, 1, ci, INSERT_VALUES)); 374 PetscCall(VecRestoreArrayRead(Xseq, &x)); 375 } 376 PetscCall(VecAssemblyBegin(CI)); 377 PetscCall(VecAssemblyEnd(CI)); 378 PetscFunctionReturn(PETSC_SUCCESS); 379 } 380 381 /* Evaluate 382 g = [ x0^2 + x1 - 2] 383 */ 384 PetscErrorCode FormEqualityConstraints(Tao tao, Vec X, Vec CE, PetscCtx ctx) 385 { 386 const PetscScalar *x; 387 PetscScalar ce; 388 MPI_Comm comm; 389 PetscMPIInt rank; 390 AppCtx *user = (AppCtx *)ctx; 391 Vec Xseq = user->Xseq; 392 VecScatter scat = user->scat; 393 394 PetscFunctionBegin; 395 PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 396 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 397 398 PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 399 PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 400 401 if (rank == 0) { 402 PetscCall(VecGetArrayRead(Xseq, &x)); 403 ce = x[0] * x[0] + x[1] - 2.0; 404 PetscCall(VecSetValue(CE, 0, ce, INSERT_VALUES)); 405 PetscCall(VecRestoreArrayRead(Xseq, &x)); 406 } 407 PetscCall(VecAssemblyBegin(CE)); 408 PetscCall(VecAssemblyEnd(CE)); 409 PetscFunctionReturn(PETSC_SUCCESS); 410 } 411 412 /* 413 grad h = [ 2*x0, -1; 414 -2*x0, 1] 415 */ 416 PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, PetscCtx ctx) 417 { 418 AppCtx *user = (AppCtx *)ctx; 419 PetscInt zero = 0, one = 1, cols[2]; 420 PetscScalar vals[2]; 421 const PetscScalar *x; 422 Vec Xseq = user->Xseq; 423 VecScatter scat = user->scat; 424 MPI_Comm comm; 425 PetscMPIInt rank; 426 427 PetscFunctionBegin; 428 PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 429 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 430 PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 431 PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD)); 432 433 PetscCall(VecGetArrayRead(Xseq, &x)); 434 if (rank == 0) { 435 cols[0] = 0; 436 cols[1] = 1; 437 vals[0] = 2 * x[0]; 438 vals[1] = -1.0; 439 PetscCall(MatSetValues(JI, 1, &zero, 2, cols, vals, INSERT_VALUES)); 440 vals[0] = -2 * x[0]; 441 vals[1] = 1.0; 442 PetscCall(MatSetValues(JI, 1, &one, 2, cols, vals, INSERT_VALUES)); 443 } 444 PetscCall(VecRestoreArrayRead(Xseq, &x)); 445 PetscCall(MatAssemblyBegin(JI, MAT_FINAL_ASSEMBLY)); 446 PetscCall(MatAssemblyEnd(JI, MAT_FINAL_ASSEMBLY)); 447 PetscFunctionReturn(PETSC_SUCCESS); 448 } 449 450 /* 451 grad g = [2*x0, 1.0] 452 */ 453 PetscErrorCode FormEqualityJacobian(Tao tao, Vec X, Mat JE, Mat JEpre, PetscCtx 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; 468 cols[1] = 1; 469 vals[0] = 2 * x[0]; 470 vals[1] = 1.0; 471 PetscCall(MatSetValues(JE, 1, &zero, 2, cols, vals, INSERT_VALUES)); 472 PetscCall(VecRestoreArrayRead(X, &x)); 473 } 474 PetscCall(MatAssemblyBegin(JE, MAT_FINAL_ASSEMBLY)); 475 PetscCall(MatAssemblyEnd(JE, MAT_FINAL_ASSEMBLY)); 476 PetscFunctionReturn(PETSC_SUCCESS); 477 } 478 479 /*TEST 480 481 build: 482 requires: !complex !defined(PETSC_USE_CXX) 483 484 test: 485 args: -tao_converged_reason -tao_gatol 1.e-6 -tao_type pdipm -tao_pdipm_kkt_shift_pd 486 requires: mumps !single 487 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 488 489 test: 490 suffix: pdipm_2 491 requires: mumps 492 nsize: 2 493 args: -tao_converged_reason -tao_gatol 1.e-6 -tao_type pdipm 494 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 495 496 test: 497 suffix: 2 498 args: -tao_converged_reason 499 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 500 501 test: 502 suffix: 3 503 args: -tao_converged_reason -no_eq 504 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 505 506 test: 507 suffix: 4 508 args: -tao_converged_reason -tao_almm_type classic 509 requires: !single 510 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 511 512 test: 513 suffix: 5 514 args: -tao_converged_reason -tao_almm_type classic -no_eq 515 requires: !single 516 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 517 518 test: 519 suffix: 6 520 args: -tao_converged_reason -tao_almm_subsolver_tao_type bqnktr 521 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 522 523 test: 524 suffix: 7 525 args: -tao_converged_reason -tao_almm_subsolver_tao_type bncg 526 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 527 528 test: 529 suffix: 8 530 nsize: 2 531 args: -tao_converged_reason 532 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 533 534 test: 535 suffix: 9 536 nsize: 2 537 args: -tao_converged_reason -vec_type cuda -mat_type aijcusparse 538 requires: cuda 539 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 540 541 test: 542 suffix: 10 543 args: -tao_converged_reason -tao_almm_type classic -no_eq -no_bound 544 requires: !single 545 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 546 547 test: 548 suffix: 11 549 args: -tao_converged_reason -tao_almm_type classic -no_bound 550 requires: !single 551 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 552 553 test: 554 suffix: 12 555 args: -tao_converged_reason -tao_almm_type phr -no_eq -no_bound 556 requires: !single 557 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 558 559 test: 560 suffix: 13 561 args: -tao_converged_reason -tao_almm_type phr -no_bound 562 requires: !single 563 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g" 564 565 TEST*/ 566