10#define __FUNCT__ "DSDPInspectXY"
11int DSDPInspectXY(
DSDP dsdp,
double xmakermu,
DSDPVec xmakery,
DSDPVec xmakerdy,
DSDPVec AX,
double *tracexs2,
double *pobj2,
double *rpinfeas2){
15 info=BoundYConeAddX(dsdp->ybcone,xmakermu,xmakery,xmakerdy,AX,tracexs2); DSDPCHKERR(info);
16 info=DSDPVecGetC(AX,pobj2);DSDPCHKERR(info);
18 info=DSDPVecSetC(AX,0);DSDPCHKERR(info);
19 info=DSDPVecSetR(AX,0);DSDPCHKERR(info);
20 info=DSDPVecNorm1(AX,rpinfeas2); DSDPCHKERR(info);
21 DSDPFunctionReturn(0);
54#define __FUNCT__ "DSDPComputeX"
55int DSDPComputeX(
DSDP dsdp){
57 double pobj=0,ppobj2=0,ddobj,tracexs=0,tracexs2=0,rpinfeas=0,rpinfeas2=0,rpobjerr=0;
58 double err1,cc,rrr,bigM,ymax,pfeastol=dsdp->pinfeastol;
63 info=DSDPStopReason(dsdp,&reason);DSDPCHKERR(info);
64 info=DSDPGetDDObjective(dsdp,&ddobj);DSDPCHKERR(info);
65 info=DSDPGetMaxYElement(dsdp,&ymax);DSDPCHKERR(info);
66 info=DSDPGetR(dsdp,&rrr); DSDPCHKERR(info);
67 info=DSDPGetPenalty(dsdp,&bigM);DSDPCHKERR(info);
68 info=DSDPGetScale(dsdp,&cc);DSDPCHKERR(info);
71 for (i=0;i<MAX_XMAKERS;i++){
72 if (i>0 && dsdp->xmaker[i].pstep<1)
continue;
73 info=
DSDPComputeXVariables(dsdp,dsdp->xmaker[i].mu,dsdp->xmaker[i].y,dsdp->xmaker[i].dy,AX,&tracexs);DSDPCHKERR(info);
74 info=DSDPVecGetC(AX,&pobj); DSDPCHKERR(info);
75 info=DSDPVecGetR(AX,&dsdp->tracex); DSDPCHKERR(info);
76 info=DSDPVecSetC(AX,0);DSDPCHKERR(info);
77 info=DSDPVecSetR(AX,0);DSDPCHKERR(info);
78 info=DSDPVecNormInfinity(AX,&rpinfeas);DSDPCHKERR(info);
79 rpinfeas=rpinfeas/(dsdp->tracex+1);
81 DSDPLogInfo(0,2,
"POBJ: %4.4e, DOBJ: %4.4e\n",pobj,ddobj/cc);
83 info=DSDPVecNorm2(AX,&err1);DSDPCHKERR(info);
84 dsdp->tracexs=tracexs;
88 info=DSDPInspectXY(dsdp,dsdp->xmaker[i].mu,dsdp->xmaker[i].y,dsdp->xmaker[i].dy,AX,&tracexs2,&ppobj2,&rpinfeas2);DSDPCHKERR(info);
89 rpinfeas2=rpinfeas2/(dsdp->tracex+1);
92 DSDPLogInfo(0,2,
"X P Infeas: %4.2e , PObj: %4.8e\n",rpinfeas,pobj*(cc));
93 DSDPLogInfo(0,2,
"TOTAL P Infeas: %4.2e PObj: %4.8e\n",rpinfeas2,ppobj2*(cc));
94 rpobjerr= fabs(pobj-dsdp->ppobj)/(1+fabs(dsdp->ppobj));
96 if (rpinfeas2 < pfeastol){
99 if (rpinfeas>pfeastol/100 && fabs(rrr)>dsdp->dinfeastol){
101 DSDPLogInfo(0,2,
"Warning: Try Increasing penalty parameter\n");
102 }
else if (rpinfeas>pfeastol && ddobj>0 && ppobj2<0 && fabs(rrr)<dsdp->dinfeastol){
104 DSDPLogInfo(0,2,
"Warning: D probably unbounded\n");
106 }
else if ( fabs(rrr)>dsdp->dinfeastol){
108 DSDPLogInfo(0,2,
"Warning: D probably infeasible \n");
116 DSDPLogInfo(0,2,
"Try backup X\n");
122 DSDPFunctionReturn(0);
127#define __FUNCT__ "DSDPSaveBackupYForX"
128int DSDPSaveBackupYForX(
DSDP dsdp,
int count,
double mu,
double pstep){
132 info=DSDPVecCopy(dsdp->y,dsdp->xmaker[count].y);DSDPCHKERR(info);
133 info=
DSDPComputeDY(dsdp,mu,dsdp->xmaker[count].dy,&pnorm); DSDPCHKERR(info);
134 dsdp->xmaker[count].pstep=pstep;
135 dsdp->xmaker[count].mu = mu;
136 DSDPFunctionReturn(0);
148#define __FUNCT__ "DSDPSaveYForX"
151 double pnorm,newgap,ymax,dd=0;
153 dsdp->pstepold=dsdp->pstep;
154 info=DSDPGetMaxYElement(dsdp,&ymax);DSDPCHKERR(info);
156 info=DSDPVecCopy(dsdp->y,dsdp->xmaker[0].y);DSDPCHKERR(info);
157 dsdp->xmaker[0].pstep=pstep;
158 }
else if (dsdp->Mshift*ymax>dsdp->pinfeastol*10){
160 if (pstep==1 && newgap>0){
161 dsdp->ppobj = dsdp->ddobj + newgap; dsdp->mu=(newgap)/(dsdp->np);
162 dsdp->dualitygap=newgap;
164 info=DSDPVecZero(dsdp->rhstemp); DSDPCHKERR(info);
165 info=BoundYConeAddX(dsdp->ybcone,dsdp->xmaker[0].mu,dsdp->xmaker[0].y,dsdp->xmaker[0].dy,dsdp->rhstemp,&dd); DSDPCHKERR(info);
166 info=DSDPVecSetC(dsdp->rhstemp,0);
167 info=DSDPVecSetR(dsdp->rhstemp,0);
168 info=DSDPVecNormInfinity(dsdp->rhstemp,&dsdp->pinfeas); DSDPCHKERR(info);
169 dsdp->pinfeas+=dsdp->Mshift*ymax;
170 if (0==1){info=DSDPVecView(dsdp->rhstemp);}
173 info=DSDPVecCopy(dsdp->y,dsdp->xmaker[0].y);DSDPCHKERR(info);
174 dsdp->xmaker[0].pstep=pstep;
176 info=
DSDPComputeDY(dsdp,mu,dsdp->xmaker[0].dy,&pnorm); DSDPCHKERR(info);
177 dsdp->xmaker[0].mu = mu;
179 if (pstep==1 && newgap>0){
180 dsdp->ppobj = dsdp->ddobj + newgap; dsdp->mu=(newgap)/(dsdp->np);
181 dsdp->dualitygap=newgap;
183 info=DSDPVecZero(dsdp->rhstemp); DSDPCHKERR(info);
184 info=BoundYConeAddX(dsdp->ybcone,dsdp->xmaker[0].mu,dsdp->xmaker[0].y,dsdp->xmaker[0].dy,dsdp->rhstemp,&dd); DSDPCHKERR(info);
185 info=DSDPVecSetC(dsdp->rhstemp,0);
186 info=DSDPVecSetR(dsdp->rhstemp,0);
187 info=DSDPVecNormInfinity(dsdp->rhstemp,&dsdp->pinfeas); DSDPCHKERR(info);
188 dsdp->pinfeas+=dsdp->Mshift*ymax;
189 if (0==1){info=DSDPVecView(dsdp->rhstemp);}
197 info=
DSDPPassXVectors(dsdp,dsdp->xmaker[0].mu,dsdp->xmaker[0].y,dsdp->xmaker[0].dy); DSDPCHKERR(info);
198 info=
DSDPGetRR(dsdp,&r);DSDPCHKERR(info);
199 if (r&& dsdp->rgap<0.1){
200 info=RConeGetRX(dsdp->rcone,&rx);DSDPCHKERR(info);
201 info=DSDPGetPenalty(dsdp,&penalty);DSDPCHKERR(info);
202 dsdp->pinfeas = dsdp->pinfeas *(1+fabs(penalty-rx));
206 if (pstep==1.0 && dsdp->rgap>5.0e-1) {
207 info=DSDPSaveBackupYForX(dsdp,MAX_XMAKERS-1,mu,pstep);DSDPCHKERR(info);
209 if (pstep==1.0 && dsdp->rgap>1.0e-3) {
210 info=DSDPSaveBackupYForX(dsdp,2,mu,pstep);DSDPCHKERR(info);
212 if (pstep==1.0 && dsdp->rgap>1.0e-5) {
213 info=DSDPSaveBackupYForX(dsdp,1,mu,pstep);DSDPCHKERR(info);
216 DSDPFunctionReturn(0);
231#define __FUNCT__ "DSDPGetPObjective"
232int DSDPGetPObjective(
DSDP dsdp,
double *pobj){
237 info=DSDPGetScale(dsdp,&scale);DSDPCHKERR(info);
238 *pobj=(dsdp->pobj)/scale;
239 DSDPFunctionReturn(0);
253#define __FUNCT__ "DSDPGetSolutionType"
257 *pdfeasible=dsdp->pdfeasible;
258 DSDPFunctionReturn(0);
277#define __FUNCT__ "DSDPGetTraceX"
278int DSDPGetTraceX(
DSDP dsdp,
double *tracex){
281 *tracex=dsdp->tracex;
282 DSDPFunctionReturn(0);
296#define __FUNCT__ "DSDPGetFinalErrors"
297int DSDPGetFinalErrors(
DSDP dsdp,
double err[6]){
299 double scale,rr,bnorm,dobj=0,pobj=0;
302 info=DSDPGetScale(dsdp,&scale);DSDPCHKERR(info);
303 info=DSDPVecGetR(dsdp->y,&rr); DSDPCHKERR(info);
304 info=DSDPGetPObjective(dsdp,&pobj);DSDPCHKERR(info);
305 info=DSDPGetDObjective(dsdp,&dobj);DSDPCHKERR(info);
308 err[2]=fabs(rr)/scale;
311 err[5]=dsdp->tracexs/scale;
313 err[2] /= (1.0+dsdp->cnorm);
314 info=DSDPVecCopy(dsdp->b,dsdp->ytemp);DSDPCHKERR(info);
315 info=DSDPVecSetC(dsdp->ytemp,0);DSDPCHKERR(info);
316 info=DSDPVecSetR(dsdp->ytemp,0);DSDPCHKERR(info);
317 info=DSDPVecNormInfinity(dsdp->ytemp,&bnorm);
318 err[0]=dsdp->perror/(1.0+bnorm);
320 err[4]=(err[4])/(1.0+fabs(pobj)+fabs(dobj));
321 err[5]=(err[5])/(1.0+fabs(pobj)+fabs(dobj));
322 DSDPFunctionReturn(0);
342#define __FUNCT__ "DSDPGetPInfeasibility"
343int DSDPGetPInfeasibility(
DSDP dsdp,
double *pperror){
346 if (pperror) *pperror=dsdp->pinfeas;
347 DSDPFunctionReturn(0);
364#define __FUNCT__ "DSDPSetPTolerance"
365int DSDPSetPTolerance(
DSDP dsdp,
double inftol){
368 if (inftol > 0) dsdp->pinfeastol = inftol;
369 DSDPLogInfo(0,2,
"Set P Infeasibility Tolerance: %4.4e\n",inftol);
370 DSDPFunctionReturn(0);
385#define __FUNCT__ "DSDPGetPTolerance"
386int DSDPGetPTolerance(
DSDP dsdp,
double *inftol){
389 if (inftol) *inftol=dsdp->pinfeastol;
390 DSDPFunctionReturn(0);
408#define __FUNCT__ "DSDPSetRTolerance"
409int DSDPSetRTolerance(
DSDP dsdp,
double inftol){
412 if (inftol > 0) dsdp->dinfeastol = inftol;
413 DSDPLogInfo(0,2,
"Set D Infeasibility Tolerance: %4.4e\n",inftol);
414 DSDPFunctionReturn(0);
433#define __FUNCT__ "DSDPGetRTolerance"
434int DSDPGetRTolerance(
DSDP dsdp,
double *inftol){
437 *inftol=dsdp->dinfeastol;
438 DSDPFunctionReturn(0);
454#define __FUNCT__ "DSDPGetYMakeX"
455int DSDPGetYMakeX(
DSDP dsdp,
double y[],
int m){
460 if (dsdp->m < m-1) DSDPFunctionReturn(1);
461 if (dsdp->m > m) DSDPFunctionReturn(1);
462 info=DSDPVecCopy(dsdp->xmaker[0].y,dsdp->ytemp); DSDPCHKERR(info);
463 info=DSDPGetScale(dsdp,&scale);DSDPCHKERR(info);
464 info=DSDPVecGetArray(dsdp->ytemp,&yy);DSDPCHKERR(info);
465 for (i=0;i<m;i++) y[i]=yy[i+1]/scale;
466 info=DSDPVecRestoreArray(dsdp->ytemp,&yy);DSDPCHKERR(info);
467 DSDPFunctionReturn(0);
482#define __FUNCT__ "DSDPGetDYMakeX"
483int DSDPGetDYMakeX(
DSDP dsdp,
double dy[],
int m){
488 if (dsdp->m < m-1) DSDPFunctionReturn(1);
489 if (dsdp->m > m) DSDPFunctionReturn(1);
490 info=DSDPVecCopy(dsdp->xmaker[0].dy,dsdp->ytemp); DSDPCHKERR(info);
491 info=DSDPGetScale(dsdp,&scale);DSDPCHKERR(info);
492 info=DSDPVecGetArray(dsdp->ytemp,&yy);DSDPCHKERR(info);
493 for (i=0;i<m;i++) dy[i]=yy[i+1]/scale;
494 info=DSDPVecRestoreArray(dsdp->ytemp,&yy);DSDPCHKERR(info);
495 DSDPFunctionReturn(0);
510#define __FUNCT__ "DSDPGetMuMakeX"
511int DSDPGetMuMakeX(
DSDP dsdp,
double *mu){
516 info=DSDPGetScale(dsdp,&scale);DSDPCHKERR(info);
517 *mu=dsdp->xmaker[0].mu/scale;
518 DSDPFunctionReturn(0);
The API to DSDP for those applications using DSDP as a subroutine library.
Internal data structure for the DSDP solver.
int DSDPComputeRHS(DSDP, double, DSDPVec)
Compute the right-hand side of the linear system that determines the step direction.
int DSDPPassXVectors(DSDP, double, DSDPVec, DSDPVec)
Pass the information needed to compute the variables X in each cone but do not compute X.
int DSDPComputeDualityGap(DSDP, double, double *)
Compute the current duality gap.
int DSDPComputeXVariables(DSDP, double, DSDPVec, DSDPVec, DSDPVec, double *)
Compute the X variables in each cone.
int DSDPComputeDY(DSDP, double, DSDPVec, double *)
Compute the step direction.
int DSDPGetRR(DSDP, double *)
Get variable r.
DSDPTerminationReason
There are many reasons to terminate the solver.
DSDPSolutionType
Formulations (P) and (D) can be feasible and bounded, feasible and unbounded, or infeasible.
Error handling, printing, and profiling.
struct DSDPVec_C DSDPVec
This object hold m+2 variables: a scaling of C, the y variables, and r.
int DSDPSaveYForX(DSDP dsdp, double mu, double pstep)
Save the current solution for later computation of X.
Internal structures for the DSDP solver.