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