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