14 #define __FUNCT__ "DSDPYStepLineSearch"
26 int info,attempt,maxattempts=30;
27 double dstep,newpotential,logdet;
28 double better=0.05, steptol=1e-8,maxmaxstep=0;
33 if (dsdp->pnorm<0.5) better=0.0;
34 dstep=DSDPMin(dstep0,0.95*maxmaxstep);
35 if (dstep * dsdp->pnorm > dsdp->maxtrustradius) dstep=dsdp->maxtrustradius/dsdp->pnorm;
36 DSDPLogInfo(0,8,
"Full Dual StepLength %4.4e, %4.4e\n",maxmaxstep,dstep);
44 if (newpotential>dsdp->potential-better && dstep > 0.001/dsdp->pnorm ){
45 DSDPLogInfo(0,2,
"Not sufficient reduction. Reduce stepsize. Trust Radius: %4.4e\n",dstep*dsdp->pnorm);
50 DSDPLogInfo(0,2,
"Dual Matrix not Positive Definite: Reduce step %4.4e",dstep);
52 if (dstep*dsdp->pnorm < steptol && dstep < steptol)
break;
55 info=
DSDPSetY(dsdp,dstep,logdet,dsdp->ytemp);DSDPCHKERR(info);
57 info=
DSDPSetY(dsdp,0,dsdp->logdet,dsdp->y);DSDPCHKERR(info);
59 DSDPFunctionReturn(0);
63 #define __FUNCT__ "DSDPYStepLineSearch2"
76 int info, attempt, maxattempts=10;
77 double dstep,newpotential,bdotdy,oldpotential,logdet;
78 double maxmaxstep=0,steptol=1e-6;
84 info=DSDPVecDot(dsdp->rhs,dy,&bdotdy);DSDPCHKERR(info);
85 dstep=DSDPMin(dstep0,0.95*maxmaxstep);
86 if (dstep * dsdp->pnorm > dsdp->maxtrustradius) dstep=dsdp->maxtrustradius/dsdp->pnorm;
87 DSDPLogInfo(0,8,
"Full Dual StepLength %4.4e, %4.4e\n",maxmaxstep,dstep);
89 if (dstep < steptol)
break;
95 b=bdotdy; a=2*(newpotential-oldpotential+bdotdy*dstep)/(dstep*dstep);
96 if (newpotential>oldpotential-0.1*dstep*bdotdy ){
97 DSDPLogInfo(0,2,
"Not sufficient reduction. Reduce stepsize. Step:: %4.4e\n",dstep);
99 if (b/a<dstep && b/a>0){ dstep=b/a;}
else { dstep=dstep/2; }
103 DSDPLogInfo(0,2,
"Dual Matrix not Positive Definite: Reduce step %4.4e",dstep);
106 if (psdefinite==
DSDP_TRUE && dstep>=steptol){
107 info=
DSDPSetY(dsdp,dstep,logdet,dsdp->ytemp);DSDPCHKERR(info);
109 info=
DSDPSetY(dsdp,0,dsdp->logdet,dsdp->y);DSDPCHKERR(info);
111 DSDPFunctionReturn(0);
115 #define __FUNCT__ "DSDPSolveDynmaicRho"
123 int info,attempt,maxattempts;
124 double dd1,dd2,mutarget,ppnorm;
131 info=DSDPVecCopy(dsdp->y,dsdp->y0);DSDPCHKERR(info);
132 for (dsdp->itnow=0; dsdp->itnow <= dsdp->maxiter+1 ; dsdp->itnow++){
137 if (dsdp->mu0>0){dsdp->mutarget=DSDPMin(dsdp->mutarget,dsdp->mu0);}
143 info=
DSDPComputePDY(dsdp,dsdp->mutarget,dsdp->dy,&dsdp->pnorm); DSDPCHKERR(info);
145 DSDPEventLogBegin(dsdp->ptime);
150 info=
DSDPSaveYForX(dsdp,dsdp->mutarget,dsdp->pstep);DSDPCHKERR(info);
156 dsdp->rho=dsdp->rhon*dsdp->np;
157 mutarget=(dsdp->ppobj-dsdp->ddobj)/(dsdp->rho);
161 dsdp->rho=(dsdp->ppobj-dsdp->ddobj)/(mutarget);
163 DSDPEventLogEnd(dsdp->ptime);
165 DSDPLogInfo(0,6,
"Current mu=%4.8e, Target X with mu=%4.8e, Rho: %8.4e, Primal Step Length: %4.8f, pnorm: %4.8e\n",dsdp->mu,mutarget,dsdp->rho/dsdp->np,dsdp->pstep, dsdp->pnorm);
170 DSDPEventLogBegin(dsdp->dtime);
171 info=
DSDPComputeDY(dsdp,mutarget,dsdp->dy,&dsdp->pnorm); DSDPCHKERR(info);
172 if (dsdp->pnorm<0.1){ mutarget/=10; info=
DSDPComputeDY(dsdp,mutarget,dsdp->dy,&dsdp->pnorm); DSDPCHKERR(info);}
175 DSDPEventLogEnd(dsdp->dtime);
177 maxattempts=dsdp->reuseM;
178 if (dsdp->dstep<1 && dsdp->rgap<1e-5) maxattempts=0;
179 if (dsdp->dstep<1e-13) maxattempts=0;
180 if (dsdp->rgap<1e-6) maxattempts=0;
181 if (maxattempts>dsdp->reuseM) maxattempts=dsdp->reuseM;
182 for (attempt=0;attempt<maxattempts;attempt++){
184 if (attempt>0 && dsdp->pnorm < 0.1)
break;
185 if (attempt > 0 && dsdp->dstep<1e-4)
break;
186 if (dsdp->rflag)
break;
187 DSDPEventLogBegin(dsdp->ctime);
188 DSDPLogInfo(0,2,
"Reuse Matrix %d: Ddobj: %12.8e, Pnorm: %4.2f, Step: %4.2f\n",attempt,dsdp->ddobj,dsdp->pnorm,dsdp->dstep);
190 info=
DSDPComputeG(dsdp,dsdp->rhstemp,dsdp->rhs1,dsdp->rhs2);DSDPCHKERR(info);
191 if (dsdp->slestype==2 || dsdp->slestype==3){
192 if (dsdp->rflag){info=
DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs1,dsdp->dy1,cgtol,&cg1);DSDPCHKERR(info);}
193 info=
DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs2,dsdp->dy2,cgtol,&cg1);DSDPCHKERR(info);
195 info=DSDPVecDot(dsdp->b,dsdp->dy1,&dd1);DSDPCHKERR(info);
196 info=DSDPVecDot(dsdp->b,dsdp->dy2,&dd2);DSDPCHKERR(info);
198 mutarget=DSDPMin(mutarget,dd1/dd2*dsdp->schurmu);
200 mutarget=mutarget*(dsdp->np/(dsdp->np+sqrt(dsdp->np)));
201 info=
DSDPComputeDY(dsdp,mutarget,dsdp->dy, &ppnorm);DSDPCHKERR(info);
202 if (ppnorm<=0){ DSDPEventLogEnd(dsdp->ctime);
break; }
205 DSDPEventLogEnd(dsdp->ctime);
207 if (attempt>0)dsdp->dstep=1.0;
209 dsdp->mutarget=DSDPMin(dsdp->mu,mutarget);
211 info=
DSDPGetRR(dsdp,&dd1);DSDPCHKERR(info);
212 if (dsdp->itnow==0 && dsdp->xmaker[0].pstep<1.0 && dd1> 0 && dsdp->pstep<1.0 && dsdp->goty0==
DSDP_FALSE){
220 DSDPFunctionReturn(0);
226 #define __FUNCT__ "DSDPChooseBarrierParameter"
241 int attempt,info,count=0;
242 double pnorm,pstep=*ppstep,pmumin,mur, dmury1, mutarget2;
247 *nextmutarget=mutarget;
255 info=
DSDPComputePDY(dsdp,mutarget,dsdp->dy,&pnorm); DSDPCHKERR(info);
258 if (pstep<1.0) {pstep=DSDPMin(0.97*pstep,1.0);}
else {pstep=DSDPMin(1.0*pstep,1.0);}
260 if (count > 2 && pstep <1e-8){pstep=0;
break;}
264 if (count>1) pstep=0.5*pstep;
else pstep=0.97*pstep;
265 DSDPLogInfo(0,2,
"Reducing pstep: %8.8e\n",pstep);
270 if (pstep > dsdp->xmaker[0].pstep || mutarget < dsdp->xmaker[0].mu * 1.0e-8){
274 DSDPFunctionReturn(0);
284 dmury1 = DSDPMin(1000,0.97*dmury1);
288 pmumin=mutarget / (1.0 + 1.0 * dmury1);
290 pmumin=mutarget / (1.0 + 1.0 * dmury1);
291 if (attempt>2){pmumin=mutarget;}
292 info=
DSDPComputePDY(dsdp,pmumin,dsdp->dy,&pnorm); DSDPCHKERR(info);
293 info=
DSDPComputePY(dsdp,pstep,dsdp->ytemp); DSDPCHKERR(info);
298 DSDPLogInfo(0,2,
"GOT X: Smallest Mu for feasible X: %4.4e.\n",pmumin);
301 DSDPLogInfo(0,6,
"GOT X: Smallest Mu for feasible X: %4.4e \n",pmumin);
308 mutarget2 = mutarget;
309 mutarget2 = (pstep)*mutarget + (1.0-pstep)*dsdp->mu;
310 mutarget2 = (pstep)*pmumin + (1.0-pstep)*dsdp->mu;
313 mutarget2=DSDPMax(dsdp->mu/dsdp->rhon,mutarget2);
314 if (dsdp->mu0>0){mutarget2=DSDPMin(mutarget2,dsdp->mu0);}
316 *nextmutarget=mutarget2;
317 DSDPFunctionReturn(0);
321 #define __FUNCT__ "DSDPResetY0"
333 info=
DSDPComputeDY(dsdp,dsdp->mutarget,dsdp->dy,&dsdp->pnorm); DSDPCHKERR(info);
334 info=DSDPVecCopy(dsdp->y0,dsdp->y);DSDPCHKERR(info);
335 info=
DSDPGetRR(dsdp,&r);DSDPCHKERR(info);
336 rr=DSDPMax(r*10000,1e12);
337 info=
DSDPSetRR(dsdp,rr);DSDPCHKERR(info);
340 info=
DSDPSetY(dsdp,1.0,dsdp->logdet,dsdp->y);DSDPCHKERR(info);
341 info=DSDPVecGetR(dsdp->b,&dd);DSDPCHKERR(info);
343 dsdp->mutarget=fabs(rr*dd);
344 dsdp->mu=fabs(rr*dd);
348 DSDPLogInfo(0,2,
"Restart algorithm\n");
349 DSDPFunctionReturn(0);
369 #define __FUNCT__ "DSDPComputeDualStepDirections"
372 double madd,ymax,cgtol=1e-7;
376 if (dsdp->itnow>30) dsdp->slestype=3;
377 if (dsdp->rgap<1e-3) dsdp->slestype=3;
378 if (dsdp->m<40) dsdp->slestype=3;
379 if (0 && dsdp->itnow>20 && dsdp->m<500) dsdp->slestype=3;
381 if (dsdp->slestype==1){
384 info=
DSDPComputeG(dsdp,dsdp->rhstemp,dsdp->rhs1,dsdp->rhs2);DSDPCHKERR(info);
385 info=
DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs1,dsdp->dy1,cgtol,&cg1);DSDPCHKERR(info);
386 if (cg1==
DSDP_TRUE){info=
DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs2,dsdp->dy2,cgtol,&cg2);DSDPCHKERR(info);}
389 if (dsdp->slestype==2){
391 DSDPLogInfo(0,9,
"Compute Hessian\n");
395 DSDPLogInfo(0,9,
"Apply CG\n");
396 info=
DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs1,dsdp->dy1,cgtol,&cg1);DSDPCHKERR(info);
397 if (cg1==
DSDP_TRUE){info=
DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs2,dsdp->dy2,cgtol,&cg2);DSDPCHKERR(info);}
401 if (dsdp->slestype==3){
402 DSDPLogInfo(0,9,
"Factor Hessian\n");
404 if (dsdp->Mshift < 1e-12 || dsdp->rgap<0.1 || dsdp->Mshift > 1e-6){
413 if (0==1 && dsdp->Mshift>dsdp->maxschurshift){
417 if (0 && dsdp->Mshift*ymax>dsdp->pinfeastol/10){
421 if (madd*ymax>dsdp->pinfeastol*1000){
432 madd=madd*4 + 1.0e-13;
437 info=
DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs1,dsdp->dy1,cgtol,&cg1);DSDPCHKERR(info);
438 info=
DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs2,dsdp->dy2,cgtol,&cg2);DSDPCHKERR(info);
442 DSDPFunctionReturn(0);
445 #define __FUNCT__ "DSDPComputeDualStepDirections"
454 }
else if (dsdp->pstep<1){
457 info=
DSDPCGSolve(dsdp,dsdp->M,rhs,dy,cgtol,&cg1);DSDPCHKERR(info);
461 DSDPFunctionReturn(0);
468 #define __FUNCT__ "DSDPInitializeVariables"
478 double r0,mutarget=dsdp->mutarget,penalty;
482 info=
DSDPGetRR(dsdp,&r0);DSDPCHKERR(info);
483 dsdp->rho=dsdp->np*dsdp->rhon;
487 if (mutarget<0) mutarget=(dsdp->ppobj-dsdp->ddobj)/(dsdp->rho);
489 info=DSDPGetPenalty(dsdp,&penalty);DSDPCHKERR(info);
490 r0=0.1/(1.0+dsdp->cnorm);
493 DSDPLogInfo(0,9,
"Set Initial R0 %4.2e\n",r0);
494 info=
DSDPSetRR(dsdp,r0);DSDPCHKERR(info);
498 if (dsdp->cnorm>0 && dsdp->anorm>0 && dsdp->cnorm/dsdp->anorm<1){ r0=r0/(dsdp->cnorm/dsdp->anorm);}
501 mutarget=(dsdp->ppobj-dsdp->ddobj)/(dsdp->rho);
502 mutarget=(dsdp->ppobj-dsdp->ddobj)/(dsdp->np);
505 DSDPLogInfo(0,9,
"Set Initial R0 %4.2e\n",r0);
506 info=
DSDPSetRR(dsdp,r0);DSDPCHKERR(info);
518 info=
DSDPSetY(dsdp,0,dsdp->logdet,dsdp->y);DSDPCHKERR(info);
519 info=
DSDPSaveYForX(dsdp,dsdp->xmaker[0].mu,0);DSDPCHKERR(info);
520 dsdp->mutarget=mutarget;
524 DSDPFunctionReturn(0);
539 #define __FUNCT__ "DSDPComputeAndFactorS"
544 DSDPFunctionReturn(0);