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