1 #include <petsctaolinesearch.h> 2 #include <../src/tao/constrained/impls/ipm/ipm.h> /*I "ipm.h" I*/ 3 4 /* 5 x,d in R^n 6 f in R 7 nb = mi + nlb+nub 8 s in R^nb is slack vector CI(x) / x-XL / -x+XU 9 bin in R^mi (tao->constraints_inequality) 10 beq in R^me (tao->constraints_equality) 11 lambdai in R^nb (ipmP->lambdai) 12 lambdae in R^me (ipmP->lambdae) 13 Jeq in R^(me x n) (tao->jacobian_equality) 14 Jin in R^(mi x n) (tao->jacobian_inequality) 15 Ai in R^(nb x n) (ipmP->Ai) 16 H in R^(n x n) (tao->hessian) 17 min f=(1/2)*x'*H*x + d'*x 18 s.t. CE(x) == 0 19 CI(x) >= 0 20 x >= tao->XL 21 -x >= -tao->XU 22 */ 23 24 static PetscErrorCode IPMComputeKKT(Tao tao); 25 static PetscErrorCode IPMPushInitialPoint(Tao tao); 26 static PetscErrorCode IPMEvaluate(Tao tao); 27 static PetscErrorCode IPMUpdateK(Tao tao); 28 static PetscErrorCode IPMUpdateAi(Tao tao); 29 static PetscErrorCode IPMGatherRHS(Tao tao, Vec, Vec, Vec, Vec, Vec); 30 static PetscErrorCode IPMScatterStep(Tao tao, Vec, Vec, Vec, Vec, Vec); 31 static PetscErrorCode IPMInitializeBounds(Tao tao); 32 33 static PetscErrorCode TaoSolve_IPM(Tao tao) 34 { 35 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 36 PetscInt its, i; 37 PetscScalar stepsize = 1.0; 38 PetscScalar step_s, step_l, alpha, tau, sigma, phi_target; 39 40 PetscFunctionBegin; 41 /* Push initial point away from bounds */ 42 PetscCall(IPMInitializeBounds(tao)); 43 PetscCall(IPMPushInitialPoint(tao)); 44 PetscCall(VecCopy(tao->solution, ipmP->rhs_x)); 45 PetscCall(IPMEvaluate(tao)); 46 PetscCall(IPMComputeKKT(tao)); 47 48 tao->reason = TAO_CONTINUE_ITERATING; 49 PetscCall(TaoLogConvergenceHistory(tao, ipmP->kkt_f, ipmP->phi, 0.0, tao->ksp_its)); 50 PetscCall(TaoMonitor(tao, tao->niter, ipmP->kkt_f, ipmP->phi, 0.0, 1.0)); 51 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 52 53 while (tao->reason == TAO_CONTINUE_ITERATING) { 54 /* Call general purpose update function */ 55 PetscTryTypeMethod(tao, update, tao->niter, tao->user_update); 56 57 tao->ksp_its = 0; 58 PetscCall(IPMUpdateK(tao)); 59 /* 60 rhs.x = -rd 61 rhs.lame = -rpe 62 rhs.lami = -rpi 63 rhs.com = -com 64 */ 65 66 PetscCall(VecCopy(ipmP->rd, ipmP->rhs_x)); 67 if (ipmP->me > 0) PetscCall(VecCopy(ipmP->rpe, ipmP->rhs_lambdae)); 68 if (ipmP->nb > 0) { 69 PetscCall(VecCopy(ipmP->rpi, ipmP->rhs_lambdai)); 70 PetscCall(VecCopy(ipmP->complementarity, ipmP->rhs_s)); 71 } 72 PetscCall(IPMGatherRHS(tao, ipmP->bigrhs, ipmP->rhs_x, ipmP->rhs_lambdae, ipmP->rhs_lambdai, ipmP->rhs_s)); 73 PetscCall(VecScale(ipmP->bigrhs, -1.0)); 74 75 /* solve K * step = rhs */ 76 PetscCall(KSPSetOperators(tao->ksp, ipmP->K, ipmP->K)); 77 PetscCall(KSPSolve(tao->ksp, ipmP->bigrhs, ipmP->bigstep)); 78 79 PetscCall(IPMScatterStep(tao, ipmP->bigstep, tao->stepdirection, ipmP->ds, ipmP->dlambdae, ipmP->dlambdai)); 80 PetscCall(KSPGetIterationNumber(tao->ksp, &its)); 81 tao->ksp_its += its; 82 tao->ksp_tot_its += its; 83 /* Find distance along step direction to closest bound */ 84 if (ipmP->nb > 0) { 85 PetscCall(VecStepBoundInfo(ipmP->s, ipmP->ds, ipmP->Zero_nb, ipmP->Inf_nb, &step_s, NULL, NULL)); 86 PetscCall(VecStepBoundInfo(ipmP->lambdai, ipmP->dlambdai, ipmP->Zero_nb, ipmP->Inf_nb, &step_l, NULL, NULL)); 87 alpha = PetscMin(step_s, step_l); 88 alpha = PetscMin(alpha, 1.0); 89 ipmP->alpha1 = alpha; 90 } else { 91 ipmP->alpha1 = alpha = 1.0; 92 } 93 94 /* x_aff = x + alpha*d */ 95 PetscCall(VecCopy(tao->solution, ipmP->save_x)); 96 if (ipmP->me > 0) PetscCall(VecCopy(ipmP->lambdae, ipmP->save_lambdae)); 97 if (ipmP->nb > 0) { 98 PetscCall(VecCopy(ipmP->lambdai, ipmP->save_lambdai)); 99 PetscCall(VecCopy(ipmP->s, ipmP->save_s)); 100 } 101 102 PetscCall(VecAXPY(tao->solution, alpha, tao->stepdirection)); 103 if (ipmP->me > 0) PetscCall(VecAXPY(ipmP->lambdae, alpha, ipmP->dlambdae)); 104 if (ipmP->nb > 0) { 105 PetscCall(VecAXPY(ipmP->lambdai, alpha, ipmP->dlambdai)); 106 PetscCall(VecAXPY(ipmP->s, alpha, ipmP->ds)); 107 } 108 109 /* Recompute kkt to find centering parameter sigma = (new_mu/old_mu)^3 */ 110 if (ipmP->mu == 0.0) { 111 sigma = 0.0; 112 } else { 113 sigma = 1.0 / ipmP->mu; 114 } 115 PetscCall(IPMComputeKKT(tao)); 116 sigma *= ipmP->mu; 117 sigma *= sigma * sigma; 118 119 /* revert kkt info */ 120 PetscCall(VecCopy(ipmP->save_x, tao->solution)); 121 if (ipmP->me > 0) PetscCall(VecCopy(ipmP->save_lambdae, ipmP->lambdae)); 122 if (ipmP->nb > 0) { 123 PetscCall(VecCopy(ipmP->save_lambdai, ipmP->lambdai)); 124 PetscCall(VecCopy(ipmP->save_s, ipmP->s)); 125 } 126 PetscCall(IPMComputeKKT(tao)); 127 128 /* update rhs with new complementarity vector */ 129 if (ipmP->nb > 0) { 130 PetscCall(VecCopy(ipmP->complementarity, ipmP->rhs_s)); 131 PetscCall(VecScale(ipmP->rhs_s, -1.0)); 132 PetscCall(VecShift(ipmP->rhs_s, sigma * ipmP->mu)); 133 } 134 PetscCall(IPMGatherRHS(tao, ipmP->bigrhs, NULL, NULL, NULL, ipmP->rhs_s)); 135 136 /* solve K * step = rhs */ 137 PetscCall(KSPSetOperators(tao->ksp, ipmP->K, ipmP->K)); 138 PetscCall(KSPSolve(tao->ksp, ipmP->bigrhs, ipmP->bigstep)); 139 140 PetscCall(IPMScatterStep(tao, ipmP->bigstep, tao->stepdirection, ipmP->ds, ipmP->dlambdae, ipmP->dlambdai)); 141 PetscCall(KSPGetIterationNumber(tao->ksp, &its)); 142 tao->ksp_its += its; 143 tao->ksp_tot_its += its; 144 if (ipmP->nb > 0) { 145 /* Get max step size and apply frac-to-boundary */ 146 tau = PetscMax(ipmP->taumin, 1.0 - ipmP->mu); 147 tau = PetscMin(tau, 1.0); 148 if (tau != 1.0) { 149 PetscCall(VecScale(ipmP->s, tau)); 150 PetscCall(VecScale(ipmP->lambdai, tau)); 151 } 152 PetscCall(VecStepBoundInfo(ipmP->s, ipmP->ds, ipmP->Zero_nb, ipmP->Inf_nb, &step_s, NULL, NULL)); 153 PetscCall(VecStepBoundInfo(ipmP->lambdai, ipmP->dlambdai, ipmP->Zero_nb, ipmP->Inf_nb, &step_l, NULL, NULL)); 154 if (tau != 1.0) { 155 PetscCall(VecCopy(ipmP->save_s, ipmP->s)); 156 PetscCall(VecCopy(ipmP->save_lambdai, ipmP->lambdai)); 157 } 158 alpha = PetscMin(step_s, step_l); 159 alpha = PetscMin(alpha, 1.0); 160 } else { 161 alpha = 1.0; 162 } 163 ipmP->alpha2 = alpha; 164 /* TODO make phi_target meaningful */ 165 phi_target = ipmP->dec * ipmP->phi; 166 for (i = 0; i < 11; i++) { 167 PetscCall(VecAXPY(tao->solution, alpha, tao->stepdirection)); 168 if (ipmP->nb > 0) { 169 PetscCall(VecAXPY(ipmP->s, alpha, ipmP->ds)); 170 PetscCall(VecAXPY(ipmP->lambdai, alpha, ipmP->dlambdai)); 171 } 172 if (ipmP->me > 0) PetscCall(VecAXPY(ipmP->lambdae, alpha, ipmP->dlambdae)); 173 174 /* update dual variables */ 175 if (ipmP->me > 0) PetscCall(VecCopy(ipmP->lambdae, tao->DE)); 176 177 PetscCall(IPMEvaluate(tao)); 178 PetscCall(IPMComputeKKT(tao)); 179 if (ipmP->phi <= phi_target) break; 180 alpha /= 2.0; 181 } 182 183 PetscCall(TaoLogConvergenceHistory(tao, ipmP->kkt_f, ipmP->phi, 0.0, tao->ksp_its)); 184 PetscCall(TaoMonitor(tao, tao->niter, ipmP->kkt_f, ipmP->phi, 0.0, stepsize)); 185 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 186 tao->niter++; 187 } 188 PetscFunctionReturn(PETSC_SUCCESS); 189 } 190 191 static PetscErrorCode TaoSetup_IPM(Tao tao) 192 { 193 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 194 195 PetscFunctionBegin; 196 ipmP->nb = ipmP->mi = ipmP->me = 0; 197 ipmP->K = NULL; 198 PetscCall(VecGetSize(tao->solution, &ipmP->n)); 199 if (!tao->gradient) { 200 PetscCall(VecDuplicate(tao->solution, &tao->gradient)); 201 PetscCall(VecDuplicate(tao->solution, &tao->stepdirection)); 202 PetscCall(VecDuplicate(tao->solution, &ipmP->rd)); 203 PetscCall(VecDuplicate(tao->solution, &ipmP->rhs_x)); 204 PetscCall(VecDuplicate(tao->solution, &ipmP->work)); 205 PetscCall(VecDuplicate(tao->solution, &ipmP->save_x)); 206 } 207 if (tao->constraints_equality) { 208 PetscCall(VecGetSize(tao->constraints_equality, &ipmP->me)); 209 PetscCall(VecDuplicate(tao->constraints_equality, &ipmP->lambdae)); 210 PetscCall(VecDuplicate(tao->constraints_equality, &ipmP->dlambdae)); 211 PetscCall(VecDuplicate(tao->constraints_equality, &ipmP->rhs_lambdae)); 212 PetscCall(VecDuplicate(tao->constraints_equality, &ipmP->save_lambdae)); 213 PetscCall(VecDuplicate(tao->constraints_equality, &ipmP->rpe)); 214 PetscCall(VecDuplicate(tao->constraints_equality, &tao->DE)); 215 } 216 if (tao->constraints_inequality) PetscCall(VecDuplicate(tao->constraints_inequality, &tao->DI)); 217 PetscFunctionReturn(PETSC_SUCCESS); 218 } 219 220 static PetscErrorCode IPMInitializeBounds(Tao tao) 221 { 222 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 223 Vec xtmp; 224 PetscInt xstart, xend; 225 PetscInt ucstart, ucend; /* user ci */ 226 PetscInt ucestart, uceend; /* user ce */ 227 PetscInt sstart = 0, send = 0; 228 PetscInt bigsize; 229 PetscInt i, counter, nloc; 230 PetscInt *cind, *xind, *ucind, *uceind, *stepind; 231 VecType vtype; 232 const PetscInt *xli, *xui; 233 PetscInt xl_offset, xu_offset; 234 IS bigxl, bigxu, isuc, isc, isx, sis, is1; 235 MPI_Comm comm; 236 237 PetscFunctionBegin; 238 cind = xind = ucind = uceind = stepind = NULL; 239 ipmP->mi = 0; 240 ipmP->nxlb = 0; 241 ipmP->nxub = 0; 242 ipmP->nb = 0; 243 ipmP->nslack = 0; 244 245 PetscCall(VecDuplicate(tao->solution, &xtmp)); 246 PetscCall(TaoComputeVariableBounds(tao)); 247 if (tao->XL) { 248 PetscCall(VecSet(xtmp, PETSC_NINFINITY)); 249 PetscCall(VecWhichGreaterThan(tao->XL, xtmp, &ipmP->isxl)); 250 PetscCall(ISGetSize(ipmP->isxl, &ipmP->nxlb)); 251 } else { 252 ipmP->nxlb = 0; 253 } 254 if (tao->XU) { 255 PetscCall(VecSet(xtmp, PETSC_INFINITY)); 256 PetscCall(VecWhichLessThan(tao->XU, xtmp, &ipmP->isxu)); 257 PetscCall(ISGetSize(ipmP->isxu, &ipmP->nxub)); 258 } else { 259 ipmP->nxub = 0; 260 } 261 PetscCall(VecDestroy(&xtmp)); 262 if (tao->constraints_inequality) { 263 PetscCall(VecGetSize(tao->constraints_inequality, &ipmP->mi)); 264 } else { 265 ipmP->mi = 0; 266 } 267 ipmP->nb = ipmP->nxlb + ipmP->nxub + ipmP->mi; 268 269 PetscCall(PetscObjectGetComm((PetscObject)tao->solution, &comm)); 270 271 bigsize = ipmP->n + 2 * ipmP->nb + ipmP->me; 272 PetscCall(PetscMalloc1(bigsize, &stepind)); 273 PetscCall(PetscMalloc1(ipmP->n, &xind)); 274 PetscCall(PetscMalloc1(ipmP->me, &uceind)); 275 PetscCall(VecGetOwnershipRange(tao->solution, &xstart, &xend)); 276 277 if (ipmP->nb > 0) { 278 PetscCall(VecCreate(comm, &ipmP->s)); 279 PetscCall(VecSetSizes(ipmP->s, PETSC_DECIDE, ipmP->nb)); 280 PetscCall(VecSetFromOptions(ipmP->s)); 281 PetscCall(VecDuplicate(ipmP->s, &ipmP->ds)); 282 PetscCall(VecDuplicate(ipmP->s, &ipmP->rhs_s)); 283 PetscCall(VecDuplicate(ipmP->s, &ipmP->complementarity)); 284 PetscCall(VecDuplicate(ipmP->s, &ipmP->ci)); 285 286 PetscCall(VecDuplicate(ipmP->s, &ipmP->lambdai)); 287 PetscCall(VecDuplicate(ipmP->s, &ipmP->dlambdai)); 288 PetscCall(VecDuplicate(ipmP->s, &ipmP->rhs_lambdai)); 289 PetscCall(VecDuplicate(ipmP->s, &ipmP->save_lambdai)); 290 291 PetscCall(VecDuplicate(ipmP->s, &ipmP->save_s)); 292 PetscCall(VecDuplicate(ipmP->s, &ipmP->rpi)); 293 PetscCall(VecDuplicate(ipmP->s, &ipmP->Zero_nb)); 294 PetscCall(VecSet(ipmP->Zero_nb, 0.0)); 295 PetscCall(VecDuplicate(ipmP->s, &ipmP->One_nb)); 296 PetscCall(VecSet(ipmP->One_nb, 1.0)); 297 PetscCall(VecDuplicate(ipmP->s, &ipmP->Inf_nb)); 298 PetscCall(VecSet(ipmP->Inf_nb, PETSC_INFINITY)); 299 300 PetscCall(PetscMalloc1(ipmP->nb, &cind)); 301 PetscCall(PetscMalloc1(ipmP->mi, &ucind)); 302 PetscCall(VecGetOwnershipRange(ipmP->s, &sstart, &send)); 303 304 if (ipmP->mi > 0) { 305 PetscCall(VecGetOwnershipRange(tao->constraints_inequality, &ucstart, &ucend)); 306 counter = 0; 307 for (i = ucstart; i < ucend; i++) cind[counter++] = i; 308 PetscCall(ISCreateGeneral(comm, counter, cind, PETSC_COPY_VALUES, &isuc)); 309 PetscCall(ISCreateGeneral(comm, counter, cind, PETSC_COPY_VALUES, &isc)); 310 PetscCall(VecScatterCreate(tao->constraints_inequality, isuc, ipmP->ci, isc, &ipmP->ci_scat)); 311 312 PetscCall(ISDestroy(&isuc)); 313 PetscCall(ISDestroy(&isc)); 314 } 315 /* need to know how may xbound indices are on each process */ 316 /* TODO better way */ 317 if (ipmP->nxlb) { 318 PetscCall(ISAllGather(ipmP->isxl, &bigxl)); 319 PetscCall(ISGetIndices(bigxl, &xli)); 320 /* find offsets for this processor */ 321 xl_offset = ipmP->mi; 322 for (i = 0; i < ipmP->nxlb; i++) { 323 if (xli[i] < xstart) { 324 xl_offset++; 325 } else break; 326 } 327 PetscCall(ISRestoreIndices(bigxl, &xli)); 328 329 PetscCall(ISGetIndices(ipmP->isxl, &xli)); 330 PetscCall(ISGetLocalSize(ipmP->isxl, &nloc)); 331 for (i = 0; i < nloc; i++) { 332 xind[i] = xli[i]; 333 cind[i] = xl_offset + i; 334 } 335 336 PetscCall(ISCreateGeneral(comm, nloc, xind, PETSC_COPY_VALUES, &isx)); 337 PetscCall(ISCreateGeneral(comm, nloc, cind, PETSC_COPY_VALUES, &isc)); 338 PetscCall(VecScatterCreate(tao->XL, isx, ipmP->ci, isc, &ipmP->xl_scat)); 339 PetscCall(ISDestroy(&isx)); 340 PetscCall(ISDestroy(&isc)); 341 PetscCall(ISDestroy(&bigxl)); 342 } 343 344 if (ipmP->nxub) { 345 PetscCall(ISAllGather(ipmP->isxu, &bigxu)); 346 PetscCall(ISGetIndices(bigxu, &xui)); 347 /* find offsets for this processor */ 348 xu_offset = ipmP->mi + ipmP->nxlb; 349 for (i = 0; i < ipmP->nxub; i++) { 350 if (xui[i] < xstart) { 351 xu_offset++; 352 } else break; 353 } 354 PetscCall(ISRestoreIndices(bigxu, &xui)); 355 356 PetscCall(ISGetIndices(ipmP->isxu, &xui)); 357 PetscCall(ISGetLocalSize(ipmP->isxu, &nloc)); 358 for (i = 0; i < nloc; i++) { 359 xind[i] = xui[i]; 360 cind[i] = xu_offset + i; 361 } 362 363 PetscCall(ISCreateGeneral(comm, nloc, xind, PETSC_COPY_VALUES, &isx)); 364 PetscCall(ISCreateGeneral(comm, nloc, cind, PETSC_COPY_VALUES, &isc)); 365 PetscCall(VecScatterCreate(tao->XU, isx, ipmP->ci, isc, &ipmP->xu_scat)); 366 PetscCall(ISDestroy(&isx)); 367 PetscCall(ISDestroy(&isc)); 368 PetscCall(ISDestroy(&bigxu)); 369 } 370 } 371 PetscCall(VecCreate(comm, &ipmP->bigrhs)); 372 PetscCall(VecGetType(tao->solution, &vtype)); 373 PetscCall(VecSetType(ipmP->bigrhs, vtype)); 374 PetscCall(VecSetSizes(ipmP->bigrhs, PETSC_DECIDE, bigsize)); 375 PetscCall(VecSetFromOptions(ipmP->bigrhs)); 376 PetscCall(VecDuplicate(ipmP->bigrhs, &ipmP->bigstep)); 377 378 /* create scatters for step->x and x->rhs */ 379 for (i = xstart; i < xend; i++) { 380 stepind[i - xstart] = i; 381 xind[i - xstart] = i; 382 } 383 PetscCall(ISCreateGeneral(comm, xend - xstart, stepind, PETSC_COPY_VALUES, &sis)); 384 PetscCall(ISCreateGeneral(comm, xend - xstart, xind, PETSC_COPY_VALUES, &is1)); 385 PetscCall(VecScatterCreate(ipmP->bigstep, sis, tao->solution, is1, &ipmP->step1)); 386 PetscCall(VecScatterCreate(tao->solution, is1, ipmP->bigrhs, sis, &ipmP->rhs1)); 387 PetscCall(ISDestroy(&sis)); 388 PetscCall(ISDestroy(&is1)); 389 390 if (ipmP->nb > 0) { 391 for (i = sstart; i < send; i++) { 392 stepind[i - sstart] = i + ipmP->n; 393 cind[i - sstart] = i; 394 } 395 PetscCall(ISCreateGeneral(comm, send - sstart, stepind, PETSC_COPY_VALUES, &sis)); 396 PetscCall(ISCreateGeneral(comm, send - sstart, cind, PETSC_COPY_VALUES, &is1)); 397 PetscCall(VecScatterCreate(ipmP->bigstep, sis, ipmP->s, is1, &ipmP->step2)); 398 PetscCall(ISDestroy(&sis)); 399 400 for (i = sstart; i < send; i++) { 401 stepind[i - sstart] = i + ipmP->n + ipmP->me; 402 cind[i - sstart] = i; 403 } 404 PetscCall(ISCreateGeneral(comm, send - sstart, stepind, PETSC_COPY_VALUES, &sis)); 405 PetscCall(VecScatterCreate(ipmP->s, is1, ipmP->bigrhs, sis, &ipmP->rhs3)); 406 PetscCall(ISDestroy(&sis)); 407 PetscCall(ISDestroy(&is1)); 408 } 409 410 if (ipmP->me > 0) { 411 PetscCall(VecGetOwnershipRange(tao->constraints_equality, &ucestart, &uceend)); 412 for (i = ucestart; i < uceend; i++) { 413 stepind[i - ucestart] = i + ipmP->n + ipmP->nb; 414 uceind[i - ucestart] = i; 415 } 416 417 PetscCall(ISCreateGeneral(comm, uceend - ucestart, stepind, PETSC_COPY_VALUES, &sis)); 418 PetscCall(ISCreateGeneral(comm, uceend - ucestart, uceind, PETSC_COPY_VALUES, &is1)); 419 PetscCall(VecScatterCreate(ipmP->bigstep, sis, tao->constraints_equality, is1, &ipmP->step3)); 420 PetscCall(ISDestroy(&sis)); 421 422 for (i = ucestart; i < uceend; i++) stepind[i - ucestart] = i + ipmP->n; 423 424 PetscCall(ISCreateGeneral(comm, uceend - ucestart, stepind, PETSC_COPY_VALUES, &sis)); 425 PetscCall(VecScatterCreate(tao->constraints_equality, is1, ipmP->bigrhs, sis, &ipmP->rhs2)); 426 PetscCall(ISDestroy(&sis)); 427 PetscCall(ISDestroy(&is1)); 428 } 429 430 if (ipmP->nb > 0) { 431 for (i = sstart; i < send; i++) { 432 stepind[i - sstart] = i + ipmP->n + ipmP->nb + ipmP->me; 433 cind[i - sstart] = i; 434 } 435 PetscCall(ISCreateGeneral(comm, send - sstart, cind, PETSC_COPY_VALUES, &is1)); 436 PetscCall(ISCreateGeneral(comm, send - sstart, stepind, PETSC_COPY_VALUES, &sis)); 437 PetscCall(VecScatterCreate(ipmP->bigstep, sis, ipmP->s, is1, &ipmP->step4)); 438 PetscCall(VecScatterCreate(ipmP->s, is1, ipmP->bigrhs, sis, &ipmP->rhs4)); 439 PetscCall(ISDestroy(&sis)); 440 PetscCall(ISDestroy(&is1)); 441 } 442 443 PetscCall(PetscFree(stepind)); 444 PetscCall(PetscFree(cind)); 445 PetscCall(PetscFree(ucind)); 446 PetscCall(PetscFree(uceind)); 447 PetscCall(PetscFree(xind)); 448 PetscFunctionReturn(PETSC_SUCCESS); 449 } 450 451 static PetscErrorCode TaoDestroy_IPM(Tao tao) 452 { 453 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 454 455 PetscFunctionBegin; 456 PetscCall(VecDestroy(&ipmP->rd)); 457 PetscCall(VecDestroy(&ipmP->rpe)); 458 PetscCall(VecDestroy(&ipmP->rpi)); 459 PetscCall(VecDestroy(&ipmP->work)); 460 PetscCall(VecDestroy(&ipmP->lambdae)); 461 PetscCall(VecDestroy(&ipmP->lambdai)); 462 PetscCall(VecDestroy(&ipmP->s)); 463 PetscCall(VecDestroy(&ipmP->ds)); 464 PetscCall(VecDestroy(&ipmP->ci)); 465 466 PetscCall(VecDestroy(&ipmP->rhs_x)); 467 PetscCall(VecDestroy(&ipmP->rhs_lambdae)); 468 PetscCall(VecDestroy(&ipmP->rhs_lambdai)); 469 PetscCall(VecDestroy(&ipmP->rhs_s)); 470 471 PetscCall(VecDestroy(&ipmP->save_x)); 472 PetscCall(VecDestroy(&ipmP->save_lambdae)); 473 PetscCall(VecDestroy(&ipmP->save_lambdai)); 474 PetscCall(VecDestroy(&ipmP->save_s)); 475 476 PetscCall(VecScatterDestroy(&ipmP->step1)); 477 PetscCall(VecScatterDestroy(&ipmP->step2)); 478 PetscCall(VecScatterDestroy(&ipmP->step3)); 479 PetscCall(VecScatterDestroy(&ipmP->step4)); 480 481 PetscCall(VecScatterDestroy(&ipmP->rhs1)); 482 PetscCall(VecScatterDestroy(&ipmP->rhs2)); 483 PetscCall(VecScatterDestroy(&ipmP->rhs3)); 484 PetscCall(VecScatterDestroy(&ipmP->rhs4)); 485 486 PetscCall(VecScatterDestroy(&ipmP->ci_scat)); 487 PetscCall(VecScatterDestroy(&ipmP->xl_scat)); 488 PetscCall(VecScatterDestroy(&ipmP->xu_scat)); 489 490 PetscCall(VecDestroy(&ipmP->dlambdai)); 491 PetscCall(VecDestroy(&ipmP->dlambdae)); 492 PetscCall(VecDestroy(&ipmP->Zero_nb)); 493 PetscCall(VecDestroy(&ipmP->One_nb)); 494 PetscCall(VecDestroy(&ipmP->Inf_nb)); 495 PetscCall(VecDestroy(&ipmP->complementarity)); 496 497 PetscCall(VecDestroy(&ipmP->bigrhs)); 498 PetscCall(VecDestroy(&ipmP->bigstep)); 499 PetscCall(MatDestroy(&ipmP->Ai)); 500 PetscCall(MatDestroy(&ipmP->K)); 501 PetscCall(ISDestroy(&ipmP->isxu)); 502 PetscCall(ISDestroy(&ipmP->isxl)); 503 PetscCall(KSPDestroy(&tao->ksp)); 504 PetscCall(PetscFree(tao->data)); 505 PetscFunctionReturn(PETSC_SUCCESS); 506 } 507 508 static PetscErrorCode TaoSetFromOptions_IPM(Tao tao, PetscOptionItems *PetscOptionsObject) 509 { 510 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 511 512 PetscFunctionBegin; 513 PetscOptionsHeadBegin(PetscOptionsObject, "IPM method for constrained optimization"); 514 PetscCall(PetscOptionsBool("-tao_ipm_monitorkkt", "monitor kkt status", NULL, ipmP->monitorkkt, &ipmP->monitorkkt, NULL)); 515 PetscCall(PetscOptionsReal("-tao_ipm_pushs", "parameter to push initial slack variables away from bounds", NULL, ipmP->pushs, &ipmP->pushs, NULL)); 516 PetscCall(PetscOptionsReal("-tao_ipm_pushnu", "parameter to push initial (inequality) dual variables away from bounds", NULL, ipmP->pushnu, &ipmP->pushnu, NULL)); 517 PetscOptionsHeadEnd(); 518 PetscCall(KSPSetFromOptions(tao->ksp)); 519 PetscFunctionReturn(PETSC_SUCCESS); 520 } 521 522 static PetscErrorCode TaoView_IPM(Tao tao, PetscViewer viewer) 523 { 524 return PETSC_SUCCESS; 525 } 526 527 /* IPMObjectiveAndGradient() 528 f = d'x + 0.5 * x' * H * x 529 rd = H*x + d + Ae'*lame - Ai'*lami 530 rpe = Ae*x - be 531 rpi = Ai*x - yi - bi 532 mu = yi' * lami/mi; 533 com = yi.*lami 534 535 phi = ||rd|| + ||rpe|| + ||rpi|| + ||com|| 536 */ 537 /* 538 static PetscErrorCode IPMObjective(TaoLineSearch ls, Vec X, PetscReal *f, void *tptr) 539 { 540 Tao tao = (Tao)tptr; 541 TAO_IPM *ipmP = (TAO_IPM*)tao->data; 542 PetscFunctionBegin; 543 PetscCall(IPMComputeKKT(tao)); 544 *f = ipmP->phi; 545 PetscFunctionReturn(PETSC_SUCCESS); 546 } 547 */ 548 549 /* 550 f = d'x + 0.5 * x' * H * x 551 rd = H*x + d + Ae'*lame - Ai'*lami 552 Ai = jac_ineq 553 I (w/lb) 554 -I (w/ub) 555 556 rpe = ce 557 rpi = ci - s; 558 com = s.*lami 559 mu = yi' * lami/mi; 560 561 phi = ||rd|| + ||rpe|| + ||rpi|| + ||com|| 562 */ 563 static PetscErrorCode IPMComputeKKT(Tao tao) 564 { 565 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 566 PetscScalar norm; 567 568 PetscFunctionBegin; 569 PetscCall(VecCopy(tao->gradient, ipmP->rd)); 570 571 if (ipmP->me > 0) { 572 /* rd = gradient + Ae'*lambdae */ 573 PetscCall(MatMultTranspose(tao->jacobian_equality, ipmP->lambdae, ipmP->work)); 574 PetscCall(VecAXPY(ipmP->rd, 1.0, ipmP->work)); 575 576 /* rpe = ce(x) */ 577 PetscCall(VecCopy(tao->constraints_equality, ipmP->rpe)); 578 } 579 if (ipmP->nb > 0) { 580 /* rd = rd - Ai'*lambdai */ 581 PetscCall(MatMultTranspose(ipmP->Ai, ipmP->lambdai, ipmP->work)); 582 PetscCall(VecAXPY(ipmP->rd, -1.0, ipmP->work)); 583 584 /* rpi = cin - s */ 585 PetscCall(VecCopy(ipmP->ci, ipmP->rpi)); 586 PetscCall(VecAXPY(ipmP->rpi, -1.0, ipmP->s)); 587 588 /* com = s .* lami */ 589 PetscCall(VecPointwiseMult(ipmP->complementarity, ipmP->s, ipmP->lambdai)); 590 } 591 /* phi = ||rd; rpe; rpi; com|| */ 592 PetscCall(VecDot(ipmP->rd, ipmP->rd, &norm)); 593 ipmP->phi = norm; 594 if (ipmP->me > 0) { 595 PetscCall(VecDot(ipmP->rpe, ipmP->rpe, &norm)); 596 ipmP->phi += norm; 597 } 598 if (ipmP->nb > 0) { 599 PetscCall(VecDot(ipmP->rpi, ipmP->rpi, &norm)); 600 ipmP->phi += norm; 601 PetscCall(VecDot(ipmP->complementarity, ipmP->complementarity, &norm)); 602 ipmP->phi += norm; 603 /* mu = s'*lami/nb */ 604 PetscCall(VecDot(ipmP->s, ipmP->lambdai, &ipmP->mu)); 605 ipmP->mu /= ipmP->nb; 606 } else { 607 ipmP->mu = 1.0; 608 } 609 610 ipmP->phi = PetscSqrtScalar(ipmP->phi); 611 PetscFunctionReturn(PETSC_SUCCESS); 612 } 613 614 /* evaluate user info at current point */ 615 static PetscErrorCode IPMEvaluate(Tao tao) 616 { 617 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 618 619 PetscFunctionBegin; 620 PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &ipmP->kkt_f, tao->gradient)); 621 PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre)); 622 if (ipmP->me > 0) { 623 PetscCall(TaoComputeEqualityConstraints(tao, tao->solution, tao->constraints_equality)); 624 PetscCall(TaoComputeJacobianEquality(tao, tao->solution, tao->jacobian_equality, tao->jacobian_equality_pre)); 625 } 626 if (ipmP->mi > 0) { 627 PetscCall(TaoComputeInequalityConstraints(tao, tao->solution, tao->constraints_inequality)); 628 PetscCall(TaoComputeJacobianInequality(tao, tao->solution, tao->jacobian_inequality, tao->jacobian_inequality_pre)); 629 } 630 if (ipmP->nb > 0) { 631 /* Ai' = jac_ineq | I (w/lb) | -I (w/ub) */ 632 PetscCall(IPMUpdateAi(tao)); 633 } 634 PetscFunctionReturn(PETSC_SUCCESS); 635 } 636 637 /* Push initial point away from bounds */ 638 static PetscErrorCode IPMPushInitialPoint(Tao tao) 639 { 640 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 641 642 PetscFunctionBegin; 643 PetscCall(TaoComputeVariableBounds(tao)); 644 PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution)); 645 if (ipmP->nb > 0) { 646 PetscCall(VecSet(ipmP->s, ipmP->pushs)); 647 PetscCall(VecSet(ipmP->lambdai, ipmP->pushnu)); 648 if (ipmP->mi > 0) PetscCall(VecSet(tao->DI, ipmP->pushnu)); 649 } 650 if (ipmP->me > 0) { 651 PetscCall(VecSet(tao->DE, 1.0)); 652 PetscCall(VecSet(ipmP->lambdae, 1.0)); 653 } 654 PetscFunctionReturn(PETSC_SUCCESS); 655 } 656 657 static PetscErrorCode IPMUpdateAi(Tao tao) 658 { 659 /* Ai = Ji 660 I (w/lb) 661 -I (w/ub) */ 662 663 /* Ci = user->ci 664 Xi - lb (w/lb) 665 -Xi + ub (w/ub) */ 666 667 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 668 MPI_Comm comm; 669 PetscInt i; 670 PetscScalar newval; 671 PetscInt newrow, newcol, ncols; 672 const PetscScalar *vals; 673 const PetscInt *cols; 674 PetscInt astart, aend, jstart, jend; 675 PetscInt *nonzeros; 676 PetscInt r2, r3, r4; 677 PetscMPIInt size; 678 Vec solu; 679 PetscInt nloc; 680 681 PetscFunctionBegin; 682 r2 = ipmP->mi; 683 r3 = r2 + ipmP->nxlb; 684 r4 = r3 + ipmP->nxub; 685 686 if (!ipmP->nb) PetscFunctionReturn(PETSC_SUCCESS); 687 688 /* Create Ai matrix if it doesn't exist yet */ 689 if (!ipmP->Ai) { 690 comm = ((PetscObject)(tao->solution))->comm; 691 PetscCallMPI(MPI_Comm_size(comm, &size)); 692 if (size == 1) { 693 PetscCall(PetscMalloc1(ipmP->nb, &nonzeros)); 694 for (i = 0; i < ipmP->mi; i++) { 695 PetscCall(MatGetRow(tao->jacobian_inequality, i, &ncols, NULL, NULL)); 696 nonzeros[i] = ncols; 697 PetscCall(MatRestoreRow(tao->jacobian_inequality, i, &ncols, NULL, NULL)); 698 } 699 for (i = r2; i < r4; i++) nonzeros[i] = 1; 700 } 701 PetscCall(MatCreate(comm, &ipmP->Ai)); 702 PetscCall(MatSetType(ipmP->Ai, MATAIJ)); 703 704 PetscCall(TaoGetSolution(tao, &solu)); 705 PetscCall(VecGetLocalSize(solu, &nloc)); 706 PetscCall(MatSetSizes(ipmP->Ai, PETSC_DECIDE, nloc, ipmP->nb, PETSC_DECIDE)); 707 PetscCall(MatSetFromOptions(ipmP->Ai)); 708 PetscCall(MatMPIAIJSetPreallocation(ipmP->Ai, ipmP->nb, NULL, ipmP->nb, NULL)); 709 PetscCall(MatSeqAIJSetPreallocation(ipmP->Ai, PETSC_DEFAULT, nonzeros)); 710 if (size == 1) PetscCall(PetscFree(nonzeros)); 711 } 712 713 /* Copy values from user jacobian to Ai */ 714 PetscCall(MatGetOwnershipRange(ipmP->Ai, &astart, &aend)); 715 716 /* Ai w/lb */ 717 if (ipmP->mi) { 718 PetscCall(MatZeroEntries(ipmP->Ai)); 719 PetscCall(MatGetOwnershipRange(tao->jacobian_inequality, &jstart, &jend)); 720 for (i = jstart; i < jend; i++) { 721 PetscCall(MatGetRow(tao->jacobian_inequality, i, &ncols, &cols, &vals)); 722 newrow = i; 723 PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, ncols, cols, vals, INSERT_VALUES)); 724 PetscCall(MatRestoreRow(tao->jacobian_inequality, i, &ncols, &cols, &vals)); 725 } 726 } 727 728 /* I w/ xlb */ 729 if (ipmP->nxlb) { 730 for (i = 0; i < ipmP->nxlb; i++) { 731 if (i >= astart && i < aend) { 732 newrow = i + r2; 733 newcol = i; 734 newval = 1.0; 735 PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES)); 736 } 737 } 738 } 739 if (ipmP->nxub) { 740 /* I w/ xub */ 741 for (i = 0; i < ipmP->nxub; i++) { 742 if (i >= astart && i < aend) { 743 newrow = i + r3; 744 newcol = i; 745 newval = -1.0; 746 PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES)); 747 } 748 } 749 } 750 751 PetscCall(MatAssemblyBegin(ipmP->Ai, MAT_FINAL_ASSEMBLY)); 752 PetscCall(MatAssemblyEnd(ipmP->Ai, MAT_FINAL_ASSEMBLY)); 753 CHKMEMQ; 754 755 PetscCall(VecSet(ipmP->ci, 0.0)); 756 757 /* user ci */ 758 if (ipmP->mi > 0) { 759 PetscCall(VecScatterBegin(ipmP->ci_scat, tao->constraints_inequality, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD)); 760 PetscCall(VecScatterEnd(ipmP->ci_scat, tao->constraints_inequality, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD)); 761 } 762 if (!ipmP->work) PetscCall(VecDuplicate(tao->solution, &ipmP->work)); 763 PetscCall(VecCopy(tao->solution, ipmP->work)); 764 if (tao->XL) { 765 PetscCall(VecAXPY(ipmP->work, -1.0, tao->XL)); 766 767 /* lower bounds on variables */ 768 if (ipmP->nxlb > 0) { 769 PetscCall(VecScatterBegin(ipmP->xl_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD)); 770 PetscCall(VecScatterEnd(ipmP->xl_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD)); 771 } 772 } 773 if (tao->XU) { 774 /* upper bounds on variables */ 775 PetscCall(VecCopy(tao->solution, ipmP->work)); 776 PetscCall(VecScale(ipmP->work, -1.0)); 777 PetscCall(VecAXPY(ipmP->work, 1.0, tao->XU)); 778 if (ipmP->nxub > 0) { 779 PetscCall(VecScatterBegin(ipmP->xu_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD)); 780 PetscCall(VecScatterEnd(ipmP->xu_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD)); 781 } 782 } 783 PetscFunctionReturn(PETSC_SUCCESS); 784 } 785 786 /* create K = [ Hlag , 0 , Ae', -Ai']; 787 [Ae , 0, 0 , 0]; 788 [Ai ,-I, 0 , 0]; 789 [ 0 , S , 0, Y ]; */ 790 static PetscErrorCode IPMUpdateK(Tao tao) 791 { 792 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 793 MPI_Comm comm; 794 PetscMPIInt size; 795 PetscInt i, j, row; 796 PetscInt ncols, newcol, newcols[2], newrow; 797 const PetscInt *cols; 798 const PetscReal *vals; 799 const PetscReal *l, *y; 800 PetscReal *newvals; 801 PetscReal newval; 802 PetscInt subsize; 803 const PetscInt *indices; 804 PetscInt *nonzeros, *d_nonzeros, *o_nonzeros; 805 PetscInt bigsize; 806 PetscInt r1, r2, r3; 807 PetscInt c1, c2, c3; 808 PetscInt klocalsize; 809 PetscInt hstart, hend, kstart, kend; 810 PetscInt aistart, aiend, aestart, aeend; 811 PetscInt sstart, send; 812 813 PetscFunctionBegin; 814 comm = ((PetscObject)(tao->solution))->comm; 815 PetscCallMPI(MPI_Comm_size(comm, &size)); 816 PetscCall(IPMUpdateAi(tao)); 817 818 /* allocate workspace */ 819 subsize = PetscMax(ipmP->n, ipmP->nb); 820 subsize = PetscMax(ipmP->me, subsize); 821 subsize = PetscMax(2, subsize); 822 PetscCall(PetscMalloc1(subsize, (PetscInt **)&indices)); 823 PetscCall(PetscMalloc1(subsize, &newvals)); 824 825 r1 = c1 = ipmP->n; 826 r2 = r1 + ipmP->me; 827 c2 = c1 + ipmP->nb; 828 r3 = c3 = r2 + ipmP->nb; 829 830 bigsize = ipmP->n + 2 * ipmP->nb + ipmP->me; 831 PetscCall(VecGetOwnershipRange(ipmP->bigrhs, &kstart, &kend)); 832 PetscCall(MatGetOwnershipRange(tao->hessian, &hstart, &hend)); 833 klocalsize = kend - kstart; 834 if (!ipmP->K) { 835 if (size == 1) { 836 PetscCall(PetscMalloc1(kend - kstart, &nonzeros)); 837 for (i = 0; i < bigsize; i++) { 838 if (i < r1) { 839 PetscCall(MatGetRow(tao->hessian, i, &ncols, NULL, NULL)); 840 nonzeros[i] = ncols; 841 PetscCall(MatRestoreRow(tao->hessian, i, &ncols, NULL, NULL)); 842 nonzeros[i] += ipmP->me + ipmP->nb; 843 } else if (i < r2) { 844 nonzeros[i - kstart] = ipmP->n; 845 } else if (i < r3) { 846 nonzeros[i - kstart] = ipmP->n + 1; 847 } else if (i < bigsize) { 848 nonzeros[i - kstart] = 2; 849 } 850 } 851 PetscCall(MatCreate(comm, &ipmP->K)); 852 PetscCall(MatSetType(ipmP->K, MATSEQAIJ)); 853 PetscCall(MatSetSizes(ipmP->K, klocalsize, klocalsize, PETSC_DETERMINE, PETSC_DETERMINE)); 854 PetscCall(MatSeqAIJSetPreallocation(ipmP->K, 0, nonzeros)); 855 PetscCall(MatSetFromOptions(ipmP->K)); 856 PetscCall(PetscFree(nonzeros)); 857 } else { 858 PetscCall(PetscMalloc1(kend - kstart, &d_nonzeros)); 859 PetscCall(PetscMalloc1(kend - kstart, &o_nonzeros)); 860 for (i = kstart; i < kend; i++) { 861 if (i < r1) { 862 /* TODO fix preallocation for mpi mats */ 863 d_nonzeros[i - kstart] = PetscMin(ipmP->n + ipmP->me + ipmP->nb, kend - kstart); 864 o_nonzeros[i - kstart] = PetscMin(ipmP->n + ipmP->me + ipmP->nb, bigsize - (kend - kstart)); 865 } else if (i < r2) { 866 d_nonzeros[i - kstart] = PetscMin(ipmP->n, kend - kstart); 867 o_nonzeros[i - kstart] = PetscMin(ipmP->n, bigsize - (kend - kstart)); 868 } else if (i < r3) { 869 d_nonzeros[i - kstart] = PetscMin(ipmP->n + 2, kend - kstart); 870 o_nonzeros[i - kstart] = PetscMin(ipmP->n + 2, bigsize - (kend - kstart)); 871 } else { 872 d_nonzeros[i - kstart] = PetscMin(2, kend - kstart); 873 o_nonzeros[i - kstart] = PetscMin(2, bigsize - (kend - kstart)); 874 } 875 } 876 PetscCall(MatCreate(comm, &ipmP->K)); 877 PetscCall(MatSetType(ipmP->K, MATMPIAIJ)); 878 PetscCall(MatSetSizes(ipmP->K, klocalsize, klocalsize, PETSC_DETERMINE, PETSC_DETERMINE)); 879 PetscCall(MatMPIAIJSetPreallocation(ipmP->K, 0, d_nonzeros, 0, o_nonzeros)); 880 PetscCall(PetscFree(d_nonzeros)); 881 PetscCall(PetscFree(o_nonzeros)); 882 PetscCall(MatSetFromOptions(ipmP->K)); 883 } 884 } 885 886 PetscCall(MatZeroEntries(ipmP->K)); 887 /* Copy H */ 888 for (i = hstart; i < hend; i++) { 889 PetscCall(MatGetRow(tao->hessian, i, &ncols, &cols, &vals)); 890 if (ncols > 0) PetscCall(MatSetValues(ipmP->K, 1, &i, ncols, cols, vals, INSERT_VALUES)); 891 PetscCall(MatRestoreRow(tao->hessian, i, &ncols, &cols, &vals)); 892 } 893 894 /* Copy Ae and Ae' */ 895 if (ipmP->me > 0) { 896 PetscCall(MatGetOwnershipRange(tao->jacobian_equality, &aestart, &aeend)); 897 for (i = aestart; i < aeend; i++) { 898 PetscCall(MatGetRow(tao->jacobian_equality, i, &ncols, &cols, &vals)); 899 if (ncols > 0) { 900 /*Ae*/ 901 row = i + r1; 902 PetscCall(MatSetValues(ipmP->K, 1, &row, ncols, cols, vals, INSERT_VALUES)); 903 /*Ae'*/ 904 for (j = 0; j < ncols; j++) { 905 newcol = i + c2; 906 newrow = cols[j]; 907 newval = vals[j]; 908 PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES)); 909 } 910 } 911 PetscCall(MatRestoreRow(tao->jacobian_equality, i, &ncols, &cols, &vals)); 912 } 913 } 914 915 if (ipmP->nb > 0) { 916 PetscCall(MatGetOwnershipRange(ipmP->Ai, &aistart, &aiend)); 917 /* Copy Ai,and Ai' */ 918 for (i = aistart; i < aiend; i++) { 919 row = i + r2; 920 PetscCall(MatGetRow(ipmP->Ai, i, &ncols, &cols, &vals)); 921 if (ncols > 0) { 922 /*Ai*/ 923 PetscCall(MatSetValues(ipmP->K, 1, &row, ncols, cols, vals, INSERT_VALUES)); 924 /*-Ai'*/ 925 for (j = 0; j < ncols; j++) { 926 newcol = i + c3; 927 newrow = cols[j]; 928 newval = -vals[j]; 929 PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES)); 930 } 931 } 932 PetscCall(MatRestoreRow(ipmP->Ai, i, &ncols, &cols, &vals)); 933 } 934 935 /* -I */ 936 for (i = kstart; i < kend; i++) { 937 if (i >= r2 && i < r3) { 938 newrow = i; 939 newcol = i - r2 + c1; 940 newval = -1.0; 941 PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES)); 942 } 943 } 944 945 /* Copy L,Y */ 946 PetscCall(VecGetOwnershipRange(ipmP->s, &sstart, &send)); 947 PetscCall(VecGetArrayRead(ipmP->lambdai, &l)); 948 PetscCall(VecGetArrayRead(ipmP->s, &y)); 949 950 for (i = sstart; i < send; i++) { 951 newcols[0] = c1 + i; 952 newcols[1] = c3 + i; 953 newvals[0] = l[i - sstart]; 954 newvals[1] = y[i - sstart]; 955 newrow = r3 + i; 956 PetscCall(MatSetValues(ipmP->K, 1, &newrow, 2, newcols, newvals, INSERT_VALUES)); 957 } 958 959 PetscCall(VecRestoreArrayRead(ipmP->lambdai, &l)); 960 PetscCall(VecRestoreArrayRead(ipmP->s, &y)); 961 } 962 963 PetscCall(PetscFree(indices)); 964 PetscCall(PetscFree(newvals)); 965 PetscCall(MatAssemblyBegin(ipmP->K, MAT_FINAL_ASSEMBLY)); 966 PetscCall(MatAssemblyEnd(ipmP->K, MAT_FINAL_ASSEMBLY)); 967 PetscFunctionReturn(PETSC_SUCCESS); 968 } 969 970 static PetscErrorCode IPMGatherRHS(Tao tao, Vec RHS, Vec X1, Vec X2, Vec X3, Vec X4) 971 { 972 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 973 974 PetscFunctionBegin; 975 /* rhs = [x1 (n) 976 x2 (me) 977 x3 (nb) 978 x4 (nb)] */ 979 if (X1) { 980 PetscCall(VecScatterBegin(ipmP->rhs1, X1, RHS, INSERT_VALUES, SCATTER_FORWARD)); 981 PetscCall(VecScatterEnd(ipmP->rhs1, X1, RHS, INSERT_VALUES, SCATTER_FORWARD)); 982 } 983 if (ipmP->me > 0 && X2) { 984 PetscCall(VecScatterBegin(ipmP->rhs2, X2, RHS, INSERT_VALUES, SCATTER_FORWARD)); 985 PetscCall(VecScatterEnd(ipmP->rhs2, X2, RHS, INSERT_VALUES, SCATTER_FORWARD)); 986 } 987 if (ipmP->nb > 0) { 988 if (X3) { 989 PetscCall(VecScatterBegin(ipmP->rhs3, X3, RHS, INSERT_VALUES, SCATTER_FORWARD)); 990 PetscCall(VecScatterEnd(ipmP->rhs3, X3, RHS, INSERT_VALUES, SCATTER_FORWARD)); 991 } 992 if (X4) { 993 PetscCall(VecScatterBegin(ipmP->rhs4, X4, RHS, INSERT_VALUES, SCATTER_FORWARD)); 994 PetscCall(VecScatterEnd(ipmP->rhs4, X4, RHS, INSERT_VALUES, SCATTER_FORWARD)); 995 } 996 } 997 PetscFunctionReturn(PETSC_SUCCESS); 998 } 999 1000 static PetscErrorCode IPMScatterStep(Tao tao, Vec STEP, Vec X1, Vec X2, Vec X3, Vec X4) 1001 { 1002 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 1003 1004 PetscFunctionBegin; 1005 CHKMEMQ; 1006 /* [x1 (n) 1007 x2 (nb) may be 0 1008 x3 (me) may be 0 1009 x4 (nb) may be 0 */ 1010 if (X1) { 1011 PetscCall(VecScatterBegin(ipmP->step1, STEP, X1, INSERT_VALUES, SCATTER_FORWARD)); 1012 PetscCall(VecScatterEnd(ipmP->step1, STEP, X1, INSERT_VALUES, SCATTER_FORWARD)); 1013 } 1014 if (X2 && ipmP->nb > 0) { 1015 PetscCall(VecScatterBegin(ipmP->step2, STEP, X2, INSERT_VALUES, SCATTER_FORWARD)); 1016 PetscCall(VecScatterEnd(ipmP->step2, STEP, X2, INSERT_VALUES, SCATTER_FORWARD)); 1017 } 1018 if (X3 && ipmP->me > 0) { 1019 PetscCall(VecScatterBegin(ipmP->step3, STEP, X3, INSERT_VALUES, SCATTER_FORWARD)); 1020 PetscCall(VecScatterEnd(ipmP->step3, STEP, X3, INSERT_VALUES, SCATTER_FORWARD)); 1021 } 1022 if (X4 && ipmP->nb > 0) { 1023 PetscCall(VecScatterBegin(ipmP->step4, STEP, X4, INSERT_VALUES, SCATTER_FORWARD)); 1024 PetscCall(VecScatterEnd(ipmP->step4, STEP, X4, INSERT_VALUES, SCATTER_FORWARD)); 1025 } 1026 CHKMEMQ; 1027 PetscFunctionReturn(PETSC_SUCCESS); 1028 } 1029 1030 /*MC 1031 TAOIPM - Interior point algorithm for generally constrained optimization. 1032 1033 Options Database Keys: 1034 + -tao_ipm_pushnu - parameter to push initial dual variables away from bounds 1035 - -tao_ipm_pushs - parameter to push initial slack variables away from bounds 1036 1037 Notes: 1038 This algorithm is more of a place-holder for future constrained optimization algorithms and should not yet be used for large problems or production code. 1039 Level: beginner 1040 1041 .seealso: `Tao`, `TAOPDIPM`, `TaoType` 1042 M*/ 1043 1044 PETSC_EXTERN PetscErrorCode TaoCreate_IPM(Tao tao) 1045 { 1046 TAO_IPM *ipmP; 1047 1048 PetscFunctionBegin; 1049 tao->ops->setup = TaoSetup_IPM; 1050 tao->ops->solve = TaoSolve_IPM; 1051 tao->ops->view = TaoView_IPM; 1052 tao->ops->setfromoptions = TaoSetFromOptions_IPM; 1053 tao->ops->destroy = TaoDestroy_IPM; 1054 /* tao->ops->computedual = TaoComputeDual_IPM; */ 1055 1056 PetscCall(PetscNew(&ipmP)); 1057 tao->data = (void *)ipmP; 1058 1059 /* Override default settings (unless already changed) */ 1060 if (!tao->max_it_changed) tao->max_it = 200; 1061 if (!tao->max_funcs_changed) tao->max_funcs = 500; 1062 1063 ipmP->dec = 10000; /* line search criteria */ 1064 ipmP->taumin = 0.995; 1065 ipmP->monitorkkt = PETSC_FALSE; 1066 ipmP->pushs = 100; 1067 ipmP->pushnu = 100; 1068 PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp)); 1069 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1)); 1070 PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix)); 1071 PetscFunctionReturn(PETSC_SUCCESS); 1072 } 1073