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 543 PetscFunctionBegin; 544 PetscCall(IPMComputeKKT(tao)); 545 *f = ipmP->phi; 546 PetscFunctionReturn(PETSC_SUCCESS); 547 } 548 */ 549 550 /* 551 f = d'x + 0.5 * x' * H * x 552 rd = H*x + d + Ae'*lame - Ai'*lami 553 Ai = jac_ineq 554 I (w/lb) 555 -I (w/ub) 556 557 rpe = ce 558 rpi = ci - s; 559 com = s.*lami 560 mu = yi' * lami/mi; 561 562 phi = ||rd|| + ||rpe|| + ||rpi|| + ||com|| 563 */ 564 static PetscErrorCode IPMComputeKKT(Tao tao) 565 { 566 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 567 PetscScalar norm; 568 569 PetscFunctionBegin; 570 PetscCall(VecCopy(tao->gradient, ipmP->rd)); 571 572 if (ipmP->me > 0) { 573 /* rd = gradient + Ae'*lambdae */ 574 PetscCall(MatMultTranspose(tao->jacobian_equality, ipmP->lambdae, ipmP->work)); 575 PetscCall(VecAXPY(ipmP->rd, 1.0, ipmP->work)); 576 577 /* rpe = ce(x) */ 578 PetscCall(VecCopy(tao->constraints_equality, ipmP->rpe)); 579 } 580 if (ipmP->nb > 0) { 581 /* rd = rd - Ai'*lambdai */ 582 PetscCall(MatMultTranspose(ipmP->Ai, ipmP->lambdai, ipmP->work)); 583 PetscCall(VecAXPY(ipmP->rd, -1.0, ipmP->work)); 584 585 /* rpi = cin - s */ 586 PetscCall(VecCopy(ipmP->ci, ipmP->rpi)); 587 PetscCall(VecAXPY(ipmP->rpi, -1.0, ipmP->s)); 588 589 /* com = s .* lami */ 590 PetscCall(VecPointwiseMult(ipmP->complementarity, ipmP->s, ipmP->lambdai)); 591 } 592 /* phi = ||rd; rpe; rpi; com|| */ 593 PetscCall(VecDot(ipmP->rd, ipmP->rd, &norm)); 594 ipmP->phi = norm; 595 if (ipmP->me > 0) { 596 PetscCall(VecDot(ipmP->rpe, ipmP->rpe, &norm)); 597 ipmP->phi += norm; 598 } 599 if (ipmP->nb > 0) { 600 PetscCall(VecDot(ipmP->rpi, ipmP->rpi, &norm)); 601 ipmP->phi += norm; 602 PetscCall(VecDot(ipmP->complementarity, ipmP->complementarity, &norm)); 603 ipmP->phi += norm; 604 /* mu = s'*lami/nb */ 605 PetscCall(VecDot(ipmP->s, ipmP->lambdai, &ipmP->mu)); 606 ipmP->mu /= ipmP->nb; 607 } else { 608 ipmP->mu = 1.0; 609 } 610 611 ipmP->phi = PetscSqrtScalar(ipmP->phi); 612 PetscFunctionReturn(PETSC_SUCCESS); 613 } 614 615 /* evaluate user info at current point */ 616 static PetscErrorCode IPMEvaluate(Tao tao) 617 { 618 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 619 620 PetscFunctionBegin; 621 PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &ipmP->kkt_f, tao->gradient)); 622 PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre)); 623 if (ipmP->me > 0) { 624 PetscCall(TaoComputeEqualityConstraints(tao, tao->solution, tao->constraints_equality)); 625 PetscCall(TaoComputeJacobianEquality(tao, tao->solution, tao->jacobian_equality, tao->jacobian_equality_pre)); 626 } 627 if (ipmP->mi > 0) { 628 PetscCall(TaoComputeInequalityConstraints(tao, tao->solution, tao->constraints_inequality)); 629 PetscCall(TaoComputeJacobianInequality(tao, tao->solution, tao->jacobian_inequality, tao->jacobian_inequality_pre)); 630 } 631 if (ipmP->nb > 0) { 632 /* Ai' = jac_ineq | I (w/lb) | -I (w/ub) */ 633 PetscCall(IPMUpdateAi(tao)); 634 } 635 PetscFunctionReturn(PETSC_SUCCESS); 636 } 637 638 /* Push initial point away from bounds */ 639 static PetscErrorCode IPMPushInitialPoint(Tao tao) 640 { 641 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 642 643 PetscFunctionBegin; 644 PetscCall(TaoComputeVariableBounds(tao)); 645 PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution)); 646 if (ipmP->nb > 0) { 647 PetscCall(VecSet(ipmP->s, ipmP->pushs)); 648 PetscCall(VecSet(ipmP->lambdai, ipmP->pushnu)); 649 if (ipmP->mi > 0) PetscCall(VecSet(tao->DI, ipmP->pushnu)); 650 } 651 if (ipmP->me > 0) { 652 PetscCall(VecSet(tao->DE, 1.0)); 653 PetscCall(VecSet(ipmP->lambdae, 1.0)); 654 } 655 PetscFunctionReturn(PETSC_SUCCESS); 656 } 657 658 static PetscErrorCode IPMUpdateAi(Tao tao) 659 { 660 /* Ai = Ji 661 I (w/lb) 662 -I (w/ub) */ 663 664 /* Ci = user->ci 665 Xi - lb (w/lb) 666 -Xi + ub (w/ub) */ 667 668 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 669 MPI_Comm comm; 670 PetscInt i; 671 PetscScalar newval; 672 PetscInt newrow, newcol, ncols; 673 const PetscScalar *vals; 674 const PetscInt *cols; 675 PetscInt astart, aend, jstart, jend; 676 PetscInt *nonzeros; 677 PetscInt r2, r3, r4; 678 PetscMPIInt size; 679 Vec solu; 680 PetscInt nloc; 681 682 PetscFunctionBegin; 683 r2 = ipmP->mi; 684 r3 = r2 + ipmP->nxlb; 685 r4 = r3 + ipmP->nxub; 686 687 if (!ipmP->nb) PetscFunctionReturn(PETSC_SUCCESS); 688 689 /* Create Ai matrix if it doesn't exist yet */ 690 if (!ipmP->Ai) { 691 comm = ((PetscObject)tao->solution)->comm; 692 PetscCallMPI(MPI_Comm_size(comm, &size)); 693 if (size == 1) { 694 PetscCall(PetscMalloc1(ipmP->nb, &nonzeros)); 695 for (i = 0; i < ipmP->mi; i++) { 696 PetscCall(MatGetRow(tao->jacobian_inequality, i, &ncols, NULL, NULL)); 697 nonzeros[i] = ncols; 698 PetscCall(MatRestoreRow(tao->jacobian_inequality, i, &ncols, NULL, NULL)); 699 } 700 for (i = r2; i < r4; i++) nonzeros[i] = 1; 701 } 702 PetscCall(MatCreate(comm, &ipmP->Ai)); 703 PetscCall(MatSetType(ipmP->Ai, MATAIJ)); 704 705 PetscCall(TaoGetSolution(tao, &solu)); 706 PetscCall(VecGetLocalSize(solu, &nloc)); 707 PetscCall(MatSetSizes(ipmP->Ai, PETSC_DECIDE, nloc, ipmP->nb, PETSC_DECIDE)); 708 PetscCall(MatSetFromOptions(ipmP->Ai)); 709 PetscCall(MatMPIAIJSetPreallocation(ipmP->Ai, ipmP->nb, NULL, ipmP->nb, NULL)); 710 PetscCall(MatSeqAIJSetPreallocation(ipmP->Ai, PETSC_DEFAULT, nonzeros)); 711 if (size == 1) PetscCall(PetscFree(nonzeros)); 712 } 713 714 /* Copy values from user jacobian to Ai */ 715 PetscCall(MatGetOwnershipRange(ipmP->Ai, &astart, &aend)); 716 717 /* Ai w/lb */ 718 if (ipmP->mi) { 719 PetscCall(MatZeroEntries(ipmP->Ai)); 720 PetscCall(MatGetOwnershipRange(tao->jacobian_inequality, &jstart, &jend)); 721 for (i = jstart; i < jend; i++) { 722 PetscCall(MatGetRow(tao->jacobian_inequality, i, &ncols, &cols, &vals)); 723 newrow = i; 724 PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, ncols, cols, vals, INSERT_VALUES)); 725 PetscCall(MatRestoreRow(tao->jacobian_inequality, i, &ncols, &cols, &vals)); 726 } 727 } 728 729 /* I w/ xlb */ 730 if (ipmP->nxlb) { 731 for (i = 0; i < ipmP->nxlb; i++) { 732 if (i >= astart && i < aend) { 733 newrow = i + r2; 734 newcol = i; 735 newval = 1.0; 736 PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES)); 737 } 738 } 739 } 740 if (ipmP->nxub) { 741 /* I w/ xub */ 742 for (i = 0; i < ipmP->nxub; i++) { 743 if (i >= astart && i < aend) { 744 newrow = i + r3; 745 newcol = i; 746 newval = -1.0; 747 PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES)); 748 } 749 } 750 } 751 752 PetscCall(MatAssemblyBegin(ipmP->Ai, MAT_FINAL_ASSEMBLY)); 753 PetscCall(MatAssemblyEnd(ipmP->Ai, MAT_FINAL_ASSEMBLY)); 754 CHKMEMQ; 755 756 PetscCall(VecSet(ipmP->ci, 0.0)); 757 758 /* user ci */ 759 if (ipmP->mi > 0) { 760 PetscCall(VecScatterBegin(ipmP->ci_scat, tao->constraints_inequality, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD)); 761 PetscCall(VecScatterEnd(ipmP->ci_scat, tao->constraints_inequality, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD)); 762 } 763 if (!ipmP->work) PetscCall(VecDuplicate(tao->solution, &ipmP->work)); 764 PetscCall(VecCopy(tao->solution, ipmP->work)); 765 if (tao->XL) { 766 PetscCall(VecAXPY(ipmP->work, -1.0, tao->XL)); 767 768 /* lower bounds on variables */ 769 if (ipmP->nxlb > 0) { 770 PetscCall(VecScatterBegin(ipmP->xl_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD)); 771 PetscCall(VecScatterEnd(ipmP->xl_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD)); 772 } 773 } 774 if (tao->XU) { 775 /* upper bounds on variables */ 776 PetscCall(VecCopy(tao->solution, ipmP->work)); 777 PetscCall(VecScale(ipmP->work, -1.0)); 778 PetscCall(VecAXPY(ipmP->work, 1.0, tao->XU)); 779 if (ipmP->nxub > 0) { 780 PetscCall(VecScatterBegin(ipmP->xu_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD)); 781 PetscCall(VecScatterEnd(ipmP->xu_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD)); 782 } 783 } 784 PetscFunctionReturn(PETSC_SUCCESS); 785 } 786 787 /* create K = [ Hlag , 0 , Ae', -Ai']; 788 [Ae , 0, 0 , 0]; 789 [Ai ,-I, 0 , 0]; 790 [ 0 , S , 0, Y ]; */ 791 static PetscErrorCode IPMUpdateK(Tao tao) 792 { 793 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 794 MPI_Comm comm; 795 PetscMPIInt size; 796 PetscInt i, j, row; 797 PetscInt ncols, newcol, newcols[2], newrow; 798 const PetscInt *cols; 799 const PetscReal *vals; 800 const PetscReal *l, *y; 801 PetscReal *newvals; 802 PetscReal newval; 803 PetscInt subsize; 804 const PetscInt *indices; 805 PetscInt *nonzeros, *d_nonzeros, *o_nonzeros; 806 PetscInt bigsize; 807 PetscInt r1, r2, r3; 808 PetscInt c1, c2, c3; 809 PetscInt klocalsize; 810 PetscInt hstart, hend, kstart, kend; 811 PetscInt aistart, aiend, aestart, aeend; 812 PetscInt sstart, send; 813 814 PetscFunctionBegin; 815 comm = ((PetscObject)tao->solution)->comm; 816 PetscCallMPI(MPI_Comm_size(comm, &size)); 817 PetscCall(IPMUpdateAi(tao)); 818 819 /* allocate workspace */ 820 subsize = PetscMax(ipmP->n, ipmP->nb); 821 subsize = PetscMax(ipmP->me, subsize); 822 subsize = PetscMax(2, subsize); 823 PetscCall(PetscMalloc1(subsize, &indices)); 824 PetscCall(PetscMalloc1(subsize, &newvals)); 825 826 r1 = c1 = ipmP->n; 827 r2 = r1 + ipmP->me; 828 c2 = c1 + ipmP->nb; 829 r3 = c3 = r2 + ipmP->nb; 830 831 bigsize = ipmP->n + 2 * ipmP->nb + ipmP->me; 832 PetscCall(VecGetOwnershipRange(ipmP->bigrhs, &kstart, &kend)); 833 PetscCall(MatGetOwnershipRange(tao->hessian, &hstart, &hend)); 834 klocalsize = kend - kstart; 835 if (!ipmP->K) { 836 if (size == 1) { 837 PetscCall(PetscMalloc1(kend - kstart, &nonzeros)); 838 for (i = 0; i < bigsize; i++) { 839 if (i < r1) { 840 PetscCall(MatGetRow(tao->hessian, i, &ncols, NULL, NULL)); 841 nonzeros[i] = ncols; 842 PetscCall(MatRestoreRow(tao->hessian, i, &ncols, NULL, NULL)); 843 nonzeros[i] += ipmP->me + ipmP->nb; 844 } else if (i < r2) { 845 nonzeros[i - kstart] = ipmP->n; 846 } else if (i < r3) { 847 nonzeros[i - kstart] = ipmP->n + 1; 848 } else if (i < bigsize) { 849 nonzeros[i - kstart] = 2; 850 } 851 } 852 PetscCall(MatCreate(comm, &ipmP->K)); 853 PetscCall(MatSetType(ipmP->K, MATSEQAIJ)); 854 PetscCall(MatSetSizes(ipmP->K, klocalsize, klocalsize, PETSC_DETERMINE, PETSC_DETERMINE)); 855 PetscCall(MatSeqAIJSetPreallocation(ipmP->K, 0, nonzeros)); 856 PetscCall(MatSetFromOptions(ipmP->K)); 857 PetscCall(PetscFree(nonzeros)); 858 } else { 859 PetscCall(PetscMalloc1(kend - kstart, &d_nonzeros)); 860 PetscCall(PetscMalloc1(kend - kstart, &o_nonzeros)); 861 for (i = kstart; i < kend; i++) { 862 if (i < r1) { 863 /* TODO fix preallocation for mpi mats */ 864 d_nonzeros[i - kstart] = PetscMin(ipmP->n + ipmP->me + ipmP->nb, kend - kstart); 865 o_nonzeros[i - kstart] = PetscMin(ipmP->n + ipmP->me + ipmP->nb, bigsize - (kend - kstart)); 866 } else if (i < r2) { 867 d_nonzeros[i - kstart] = PetscMin(ipmP->n, kend - kstart); 868 o_nonzeros[i - kstart] = PetscMin(ipmP->n, bigsize - (kend - kstart)); 869 } else if (i < r3) { 870 d_nonzeros[i - kstart] = PetscMin(ipmP->n + 2, kend - kstart); 871 o_nonzeros[i - kstart] = PetscMin(ipmP->n + 2, bigsize - (kend - kstart)); 872 } else { 873 d_nonzeros[i - kstart] = PetscMin(2, kend - kstart); 874 o_nonzeros[i - kstart] = PetscMin(2, bigsize - (kend - kstart)); 875 } 876 } 877 PetscCall(MatCreate(comm, &ipmP->K)); 878 PetscCall(MatSetType(ipmP->K, MATMPIAIJ)); 879 PetscCall(MatSetSizes(ipmP->K, klocalsize, klocalsize, PETSC_DETERMINE, PETSC_DETERMINE)); 880 PetscCall(MatMPIAIJSetPreallocation(ipmP->K, 0, d_nonzeros, 0, o_nonzeros)); 881 PetscCall(PetscFree(d_nonzeros)); 882 PetscCall(PetscFree(o_nonzeros)); 883 PetscCall(MatSetFromOptions(ipmP->K)); 884 } 885 } 886 887 PetscCall(MatZeroEntries(ipmP->K)); 888 /* Copy H */ 889 for (i = hstart; i < hend; i++) { 890 PetscCall(MatGetRow(tao->hessian, i, &ncols, &cols, &vals)); 891 if (ncols > 0) PetscCall(MatSetValues(ipmP->K, 1, &i, ncols, cols, vals, INSERT_VALUES)); 892 PetscCall(MatRestoreRow(tao->hessian, i, &ncols, &cols, &vals)); 893 } 894 895 /* Copy Ae and Ae' */ 896 if (ipmP->me > 0) { 897 PetscCall(MatGetOwnershipRange(tao->jacobian_equality, &aestart, &aeend)); 898 for (i = aestart; i < aeend; i++) { 899 PetscCall(MatGetRow(tao->jacobian_equality, i, &ncols, &cols, &vals)); 900 if (ncols > 0) { 901 /*Ae*/ 902 row = i + r1; 903 PetscCall(MatSetValues(ipmP->K, 1, &row, ncols, cols, vals, INSERT_VALUES)); 904 /*Ae'*/ 905 for (j = 0; j < ncols; j++) { 906 newcol = i + c2; 907 newrow = cols[j]; 908 newval = vals[j]; 909 PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES)); 910 } 911 } 912 PetscCall(MatRestoreRow(tao->jacobian_equality, i, &ncols, &cols, &vals)); 913 } 914 } 915 916 if (ipmP->nb > 0) { 917 PetscCall(MatGetOwnershipRange(ipmP->Ai, &aistart, &aiend)); 918 /* Copy Ai,and Ai' */ 919 for (i = aistart; i < aiend; i++) { 920 row = i + r2; 921 PetscCall(MatGetRow(ipmP->Ai, i, &ncols, &cols, &vals)); 922 if (ncols > 0) { 923 /*Ai*/ 924 PetscCall(MatSetValues(ipmP->K, 1, &row, ncols, cols, vals, INSERT_VALUES)); 925 /*-Ai'*/ 926 for (j = 0; j < ncols; j++) { 927 newcol = i + c3; 928 newrow = cols[j]; 929 newval = -vals[j]; 930 PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES)); 931 } 932 } 933 PetscCall(MatRestoreRow(ipmP->Ai, i, &ncols, &cols, &vals)); 934 } 935 936 /* -I */ 937 for (i = kstart; i < kend; i++) { 938 if (i >= r2 && i < r3) { 939 newrow = i; 940 newcol = i - r2 + c1; 941 newval = -1.0; 942 PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES)); 943 } 944 } 945 946 /* Copy L,Y */ 947 PetscCall(VecGetOwnershipRange(ipmP->s, &sstart, &send)); 948 PetscCall(VecGetArrayRead(ipmP->lambdai, &l)); 949 PetscCall(VecGetArrayRead(ipmP->s, &y)); 950 951 for (i = sstart; i < send; i++) { 952 newcols[0] = c1 + i; 953 newcols[1] = c3 + i; 954 newvals[0] = l[i - sstart]; 955 newvals[1] = y[i - sstart]; 956 newrow = r3 + i; 957 PetscCall(MatSetValues(ipmP->K, 1, &newrow, 2, newcols, newvals, INSERT_VALUES)); 958 } 959 960 PetscCall(VecRestoreArrayRead(ipmP->lambdai, &l)); 961 PetscCall(VecRestoreArrayRead(ipmP->s, &y)); 962 } 963 964 PetscCall(PetscFree(indices)); 965 PetscCall(PetscFree(newvals)); 966 PetscCall(MatAssemblyBegin(ipmP->K, MAT_FINAL_ASSEMBLY)); 967 PetscCall(MatAssemblyEnd(ipmP->K, MAT_FINAL_ASSEMBLY)); 968 PetscFunctionReturn(PETSC_SUCCESS); 969 } 970 971 static PetscErrorCode IPMGatherRHS(Tao tao, Vec RHS, Vec X1, Vec X2, Vec X3, Vec X4) 972 { 973 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 974 975 PetscFunctionBegin; 976 /* rhs = [x1 (n) 977 x2 (me) 978 x3 (nb) 979 x4 (nb)] */ 980 if (X1) { 981 PetscCall(VecScatterBegin(ipmP->rhs1, X1, RHS, INSERT_VALUES, SCATTER_FORWARD)); 982 PetscCall(VecScatterEnd(ipmP->rhs1, X1, RHS, INSERT_VALUES, SCATTER_FORWARD)); 983 } 984 if (ipmP->me > 0 && X2) { 985 PetscCall(VecScatterBegin(ipmP->rhs2, X2, RHS, INSERT_VALUES, SCATTER_FORWARD)); 986 PetscCall(VecScatterEnd(ipmP->rhs2, X2, RHS, INSERT_VALUES, SCATTER_FORWARD)); 987 } 988 if (ipmP->nb > 0) { 989 if (X3) { 990 PetscCall(VecScatterBegin(ipmP->rhs3, X3, RHS, INSERT_VALUES, SCATTER_FORWARD)); 991 PetscCall(VecScatterEnd(ipmP->rhs3, X3, RHS, INSERT_VALUES, SCATTER_FORWARD)); 992 } 993 if (X4) { 994 PetscCall(VecScatterBegin(ipmP->rhs4, X4, RHS, INSERT_VALUES, SCATTER_FORWARD)); 995 PetscCall(VecScatterEnd(ipmP->rhs4, X4, RHS, INSERT_VALUES, SCATTER_FORWARD)); 996 } 997 } 998 PetscFunctionReturn(PETSC_SUCCESS); 999 } 1000 1001 static PetscErrorCode IPMScatterStep(Tao tao, Vec STEP, Vec X1, Vec X2, Vec X3, Vec X4) 1002 { 1003 TAO_IPM *ipmP = (TAO_IPM *)tao->data; 1004 1005 PetscFunctionBegin; 1006 CHKMEMQ; 1007 /* [x1 (n) 1008 x2 (nb) may be 0 1009 x3 (me) may be 0 1010 x4 (nb) may be 0 */ 1011 if (X1) { 1012 PetscCall(VecScatterBegin(ipmP->step1, STEP, X1, INSERT_VALUES, SCATTER_FORWARD)); 1013 PetscCall(VecScatterEnd(ipmP->step1, STEP, X1, INSERT_VALUES, SCATTER_FORWARD)); 1014 } 1015 if (X2 && ipmP->nb > 0) { 1016 PetscCall(VecScatterBegin(ipmP->step2, STEP, X2, INSERT_VALUES, SCATTER_FORWARD)); 1017 PetscCall(VecScatterEnd(ipmP->step2, STEP, X2, INSERT_VALUES, SCATTER_FORWARD)); 1018 } 1019 if (X3 && ipmP->me > 0) { 1020 PetscCall(VecScatterBegin(ipmP->step3, STEP, X3, INSERT_VALUES, SCATTER_FORWARD)); 1021 PetscCall(VecScatterEnd(ipmP->step3, STEP, X3, INSERT_VALUES, SCATTER_FORWARD)); 1022 } 1023 if (X4 && ipmP->nb > 0) { 1024 PetscCall(VecScatterBegin(ipmP->step4, STEP, X4, INSERT_VALUES, SCATTER_FORWARD)); 1025 PetscCall(VecScatterEnd(ipmP->step4, STEP, X4, INSERT_VALUES, SCATTER_FORWARD)); 1026 } 1027 CHKMEMQ; 1028 PetscFunctionReturn(PETSC_SUCCESS); 1029 } 1030 1031 /*MC 1032 TAOIPM - Interior point algorithm for generally constrained optimization. 1033 1034 Options Database Keys: 1035 + -tao_ipm_pushnu - parameter to push initial dual variables away from bounds 1036 - -tao_ipm_pushs - parameter to push initial slack variables away from bounds 1037 1038 Notes: 1039 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. 1040 Level: beginner 1041 1042 .seealso: `Tao`, `TAOPDIPM`, `TaoType` 1043 M*/ 1044 1045 PETSC_EXTERN PetscErrorCode TaoCreate_IPM(Tao tao) 1046 { 1047 TAO_IPM *ipmP; 1048 1049 PetscFunctionBegin; 1050 tao->ops->setup = TaoSetup_IPM; 1051 tao->ops->solve = TaoSolve_IPM; 1052 tao->ops->view = TaoView_IPM; 1053 tao->ops->setfromoptions = TaoSetFromOptions_IPM; 1054 tao->ops->destroy = TaoDestroy_IPM; 1055 /* tao->ops->computedual = TaoComputeDual_IPM; */ 1056 1057 PetscCall(PetscNew(&ipmP)); 1058 tao->data = (void *)ipmP; 1059 1060 /* Override default settings (unless already changed) */ 1061 PetscCall(TaoParametersInitialize(tao)); 1062 PetscObjectParameterSetDefault(tao, max_it, 200); 1063 PetscObjectParameterSetDefault(tao, max_funcs, 500); 1064 1065 ipmP->dec = 10000; /* line search criteria */ 1066 ipmP->taumin = 0.995; 1067 ipmP->monitorkkt = PETSC_FALSE; 1068 ipmP->pushs = 100; 1069 ipmP->pushnu = 100; 1070 PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp)); 1071 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1)); 1072 PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix)); 1073 PetscFunctionReturn(PETSC_SUCCESS); 1074 } 1075