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