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);
DSDPTruth
Boolean variables.
int DSDPComputeLogSDeterminant(DSDP, double *)
Compute the logarithmic barrier function for the dual varialbe S.
struct DSDPVec_C DSDPVec
This object hold m+2 variables: a scaling of C, the y variables, and r.
int DSDPInvertS(DSDP)
Invert the S variables in each cone.
int DSDPComputeSS(DSDP, DSDPVec, DSDPDualFactorMatrix, DSDPTruth *)
Compute the dual variables S in each cone.
int DSDPSchurMatFactor(DSDPSchurMat M, DSDPTruth *successful)
Factor M.
Error handling, printing, and profiling.
int DSDPComputePDY1(DSDP, double, DSDPVec)
Compute an affine step direction dy1.
int DSDPSolveDynamicRho(DSDP dsdp)
Apply dual-scaling algorithm.
int DSDPYStepLineSearch2(DSDP, double, double, DSDPVec)
Used for centering steps, the merit function of this line search is the objective function plus the b...
int DSDPComputeNewY(DSDP, double, DSDPVec)
Update the Y variables.
int DSDPComputeAndFactorS(DSDP dsdp, DSDPTruth *psdefinite)
Compute and factor the dual matrix variables.
Internal structures for the DSDP solver.
DSDPTerminationReason
There are many reasons to terminate the solver.
int DSDPCheckConvergence(DSDP, DSDPTerminationReason *)
Check for convergence and monitor solution.
int DSDPSetRR(DSDP, double)
Set variable r.
The API to DSDP for those applications using DSDP as a subroutine library.
int DSDPComputePY(DSDP, double, DSDPVec)
Compute PY = Y - beta DY for use in computing X.
Internal data structure for the DSDP solver.
int DSDPSaveYForX(DSDP, double, double)
Save the current solution for later computation of X.
int DSDPSetY(DSDP, double, double, DSDPVec)
Update the solver with these y variables.
int DSDPSchurMatView(DSDPSchurMat M)
Print the matrix.
int DSDPComputeDY(DSDP, double, DSDPVec, double *)
Compute the step direction.
int DSDPComputeHessian(DSDP, DSDPSchurMat, DSDPVec, DSDPVec)
Compute the Schur complement, or Gram, matrix for each cone.
int DSDPComputePotential(DSDP, DSDPVec, double, double *)
Compute the potential of the given point.
int DSDPYStepLineSearch(DSDP, double, double, DSDPVec)
Used for Newton step, the merit function of this line search is the dual potential function.
int DSDPComputeObjective(DSDP, DSDPVec, double *)
Compute the objective function (DD).
int DSDPGetRR(DSDP, double *)
Get variable r.
int DSDPComputeG(DSDP, DSDPVec, DSDPVec, DSDPVec)
Compute the gradient of the barrier for each cone.
int DSDPComputeMaxStepLength(DSDP, DSDPVec, DSDPDualFactorMatrix, double *)
Compute the maximum step length for the given step direction.
int DSDPComputeDualStepDirections(DSDP dsdp)
Compute the step direction by computing a linear system and solving it.
int DSDPSchurMatShiftDiagonal(DSDPSchurMat M, double dd)
Add a scalar to each diagonal element of the matrix.
int DSDPComputePDY(DSDP, double, DSDPVec, double *)
Compute the step direction.
int DSDPSetConvergenceFlag(DSDP dsdp, DSDPTerminationReason reason)
Monitor each iteration of the solver.
int DSDPComputePotential2(DSDP, DSDPVec, double, double, double *)
Compute the objective function plus the barrier function.
int DSDPCGSolve(DSDP, DSDPSchurMat, DSDPVec, DSDPVec, double, DSDPTruth *)
Apply CG to solve for the step directions.
int DSDPInitializeVariables(DSDP dsdp)
Initialize variables and factor S.
int DSDPChooseBarrierParameter(DSDP, double, double *, double *)
DSDP barrier heuristic choses the smalles value of mu such that X>0.
int DSDPResetY0(DSDP)
After 1 iteration, consider increasing the variable r.
int DSDPGetMaxYElement(DSDP, double *)
Copy the the infinity norm of the variables y.