DSDP
dsdpx.c
Go to the documentation of this file.
1#include "dsdp.h"
2#include "dsdpsys.h"
3#include "dsdp5.h"
9#undef __FUNCT__
10#define __FUNCT__ "DSDPInspectXY"
11int DSDPInspectXY(DSDP dsdp, double xmakermu, DSDPVec xmakery, DSDPVec xmakerdy, DSDPVec AX, double *tracexs2, double *pobj2, double *rpinfeas2){
12 int info;
13 DSDPFunctionBegin;
14
15 info=BoundYConeAddX(dsdp->ybcone,xmakermu,xmakery,xmakerdy,AX,tracexs2); DSDPCHKERR(info);
16 info=DSDPVecGetC(AX,pobj2);DSDPCHKERR(info);
17
18 info=DSDPVecSetC(AX,0);DSDPCHKERR(info);
19 info=DSDPVecSetR(AX,0);DSDPCHKERR(info);
20 info=DSDPVecNorm1(AX,rpinfeas2); DSDPCHKERR(info);
21 DSDPFunctionReturn(0);
22}
23
53#undef __FUNCT__
54#define __FUNCT__ "DSDPComputeX"
55int DSDPComputeX(DSDP dsdp){
56 int i,info;
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;
60 DSDPVec AX=dsdp->ytemp;
61
62 DSDPFunctionBegin;
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);
69
70 dsdp->pdfeasible=DSDP_PDFEASIBLE;
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);
80
81 DSDPLogInfo(0,2,"POBJ: %4.4e, DOBJ: %4.4e\n",pobj,ddobj/cc);
82
83 info=DSDPVecNorm2(AX,&err1);DSDPCHKERR(info);
84 dsdp->tracexs=tracexs;
85 dsdp->perror=err1;
86 dsdp->pobj=cc*pobj;
87
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);
90 /* rpinfeas is infeasibility of (P) while rpinfeas2 is infeasibility of (PP) */
91
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));
95
96 if (rpinfeas2 < pfeastol){ /* (PP) must be close to feasible */
97
98 if (dsdp->rgap<0.1){
99 if (rpinfeas>pfeastol/100 && fabs(rrr)>dsdp->dinfeastol){
100 dsdp->pdfeasible=DSDP_PDUNKNOWN;
101 DSDPLogInfo(0,2,"Warning: Try Increasing penalty parameter\n");
102 } else if (rpinfeas>pfeastol && ddobj>0 && ppobj2<0 && fabs(rrr)<dsdp->dinfeastol){
103 dsdp->pdfeasible=DSDP_UNBOUNDED;
104 DSDPLogInfo(0,2,"Warning: D probably unbounded\n");
105
106 } else if (/* fabs(bigM)-dsdp->tracex < fabs(rrr) && rpinfeas<pfeastol */ fabs(rrr)>dsdp->dinfeastol){
107 dsdp->pdfeasible=DSDP_INFEASIBLE;
108 DSDPLogInfo(0,2,"Warning: D probably infeasible \n");
109 }
110 }
111 i=i+10;
112 break;
113
114 } else {
115 /* Step direction was not accurate enough to compute X from Schur complement */
116 DSDPLogInfo(0,2,"Try backup X\n");
117 info=DSDPSetConvergenceFlag(dsdp,DSDP_NUMERICAL_ERROR); DSDPCHKERR(info);
118 }
119
120 }
121
122 DSDPFunctionReturn(0);
123}
124
125
126#undef __FUNCT__
127#define __FUNCT__ "DSDPSaveBackupYForX"
128int DSDPSaveBackupYForX(DSDP dsdp, int count,double mu, double pstep){
129 int info;
130 double pnorm;
131 DSDPFunctionBegin;
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);
137}
138
147#undef __FUNCT__
148#define __FUNCT__ "DSDPSaveYForX"
149int DSDPSaveYForX(DSDP dsdp, double mu, double pstep){
150 int info,isgood=0;
151 double pnorm,newgap,ymax,dd=0;
152 DSDPFunctionBegin;
153 dsdp->pstepold=dsdp->pstep;
154 info=DSDPGetMaxYElement(dsdp,&ymax);DSDPCHKERR(info);
155 if (pstep==0){
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){
159 info=DSDPComputeDualityGap(dsdp,mu,&newgap);DSDPCHKERR(info);
160 if (pstep==1 && newgap>0){
161 dsdp->ppobj = dsdp->ddobj + newgap; dsdp->mu=(newgap)/(dsdp->np);
162 dsdp->dualitygap=newgap;
163 }
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);}
171 /* Not good enough to save */
172 } else {
173 info=DSDPVecCopy(dsdp->y,dsdp->xmaker[0].y);DSDPCHKERR(info);
174 dsdp->xmaker[0].pstep=pstep;
175 info=DSDPComputeRHS(dsdp,mu,dsdp->xmakerrhs); DSDPCHKERR(info);
176 info=DSDPComputeDY(dsdp,mu,dsdp->xmaker[0].dy,&pnorm); DSDPCHKERR(info);
177 dsdp->xmaker[0].mu = mu;
178 info=DSDPComputeDualityGap(dsdp,mu,&newgap);DSDPCHKERR(info);
179 if (pstep==1 && newgap>0){
180 dsdp->ppobj = dsdp->ddobj + newgap; dsdp->mu=(newgap)/(dsdp->np);
181 dsdp->dualitygap=newgap;
182
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);}
190
191 }
192 isgood=1;
193 }
194
195 if (isgood==1){
196 double penalty,r,rx;
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){ /* Estimate error in X */
200 info=RConeGetRX(dsdp->rcone,&rx);DSDPCHKERR(info);
201 info=DSDPGetPenalty(dsdp,&penalty);DSDPCHKERR(info);
202 dsdp->pinfeas = dsdp->pinfeas *(1+fabs(penalty-rx));
203 }
204 }
205
206 if (pstep==1.0 && dsdp->rgap>5.0e-1) {
207 info=DSDPSaveBackupYForX(dsdp,MAX_XMAKERS-1,mu,pstep);DSDPCHKERR(info);
208 }
209 if (pstep==1.0 && dsdp->rgap>1.0e-3) {
210 info=DSDPSaveBackupYForX(dsdp,2,mu,pstep);DSDPCHKERR(info);
211 }
212 if (pstep==1.0 && dsdp->rgap>1.0e-5) {
213 info=DSDPSaveBackupYForX(dsdp,1,mu,pstep);DSDPCHKERR(info);
214 }
215
216 DSDPFunctionReturn(0);
217}
218
230#undef __FUNCT__
231#define __FUNCT__ "DSDPGetPObjective"
232int DSDPGetPObjective(DSDP dsdp,double *pobj){
233 int info;
234 double scale;
235 DSDPFunctionBegin;
236 DSDPValid(dsdp);
237 info=DSDPGetScale(dsdp,&scale);DSDPCHKERR(info);
238 *pobj=(dsdp->pobj)/scale;
239 DSDPFunctionReturn(0);
240}
241
252#undef __FUNCT__
253#define __FUNCT__ "DSDPGetSolutionType"
254int DSDPGetSolutionType(DSDP dsdp,DSDPSolutionType *pdfeasible){
255 DSDPFunctionBegin;
256 DSDPValid(dsdp);
257 *pdfeasible=dsdp->pdfeasible;
258 DSDPFunctionReturn(0);
259}
276#undef __FUNCT__
277#define __FUNCT__ "DSDPGetTraceX"
278int DSDPGetTraceX(DSDP dsdp, double *tracex){
279 DSDPFunctionBegin;
280 DSDPValid(dsdp);
281 *tracex=dsdp->tracex;
282 DSDPFunctionReturn(0);
283}
284
295#undef __FUNCT__
296#define __FUNCT__ "DSDPGetFinalErrors"
297int DSDPGetFinalErrors(DSDP dsdp, double err[6]){
298 int info;
299 double scale,rr,bnorm,dobj=0,pobj=0;
300 DSDPFunctionBegin;
301 DSDPValid(dsdp);
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);
306 err[0]=dsdp->perror;
307 err[1]=0;
308 err[2]=fabs(rr)/scale;
309 err[3]=0;
310 err[4]=pobj - dobj;
311 err[5]=dsdp->tracexs/scale;
312
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);
319
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);
323}
324
341#undef __FUNCT__
342#define __FUNCT__ "DSDPGetPInfeasibility"
343int DSDPGetPInfeasibility(DSDP dsdp, double *pperror){
344 DSDPFunctionBegin;
345 DSDPValid(dsdp);
346 if (pperror) *pperror=dsdp->pinfeas;
347 DSDPFunctionReturn(0);
348}
349
363#undef __FUNCT__
364#define __FUNCT__ "DSDPSetPTolerance"
365int DSDPSetPTolerance(DSDP dsdp,double inftol){
366 DSDPFunctionBegin;
367 DSDPValid(dsdp);
368 if (inftol > 0) dsdp->pinfeastol = inftol;
369 DSDPLogInfo(0,2,"Set P Infeasibility Tolerance: %4.4e\n",inftol);
370 DSDPFunctionReturn(0);
371}
372
384#undef __FUNCT__
385#define __FUNCT__ "DSDPGetPTolerance"
386int DSDPGetPTolerance(DSDP dsdp,double *inftol){
387 DSDPFunctionBegin;
388 DSDPValid(dsdp);
389 if (inftol) *inftol=dsdp->pinfeastol;
390 DSDPFunctionReturn(0);
391}
392
393
407#undef __FUNCT__
408#define __FUNCT__ "DSDPSetRTolerance"
409int DSDPSetRTolerance(DSDP dsdp,double inftol){
410 DSDPFunctionBegin;
411 DSDPValid(dsdp);
412 if (inftol > 0) dsdp->dinfeastol = inftol;
413 DSDPLogInfo(0,2,"Set D Infeasibility Tolerance: %4.4e\n",inftol);
414 DSDPFunctionReturn(0);
415}
416
432#undef __FUNCT__
433#define __FUNCT__ "DSDPGetRTolerance"
434int DSDPGetRTolerance(DSDP dsdp,double *inftol){
435 DSDPFunctionBegin;
436 DSDPValid(dsdp);
437 *inftol=dsdp->dinfeastol;
438 DSDPFunctionReturn(0);
439}
440
453#undef __FUNCT__
454#define __FUNCT__ "DSDPGetYMakeX"
455int DSDPGetYMakeX(DSDP dsdp,double y[], int m){
456 int i,info;
457 double scale,*yy;
458 DSDPFunctionBegin;
459 DSDPValid(dsdp);
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);
468}
469
481#undef __FUNCT__
482#define __FUNCT__ "DSDPGetDYMakeX"
483int DSDPGetDYMakeX(DSDP dsdp,double dy[], int m){
484 int i,info;
485 double scale,*yy;
486 DSDPFunctionBegin;
487 DSDPValid(dsdp);
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);
496}
497
509#undef __FUNCT__
510#define __FUNCT__ "DSDPGetMuMakeX"
511int DSDPGetMuMakeX(DSDP dsdp,double *mu){
512 int info;
513 double scale;
514 DSDPFunctionBegin;
515 DSDPValid(dsdp);
516 info=DSDPGetScale(dsdp,&scale);DSDPCHKERR(info);
517 *mu=dsdp->xmaker[0].mu/scale;
518 DSDPFunctionReturn(0);
519}
520
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.
Definition dualimpl.c:177
int DSDPPassXVectors(DSDP, double, DSDPVec, DSDPVec)
Pass the information needed to compute the variables X in each cone but do not compute X.
Definition dsdpcops.c:378
int DSDPComputeDualityGap(DSDP, double, double *)
Compute the current duality gap.
Definition dualimpl.c:230
int DSDPComputeXVariables(DSDP, double, DSDPVec, DSDPVec, DSDPVec, double *)
Compute the X variables in each cone.
Definition dsdpcops.c:654
int DSDPComputeDY(DSDP, double, DSDPVec, double *)
Compute the step direction.
Definition dualimpl.c:45
int DSDPGetRR(DSDP, double *)
Get variable r.
Definition dualimpl.c:361
DSDPTerminationReason
There are many reasons to terminate the solver.
@ DSDP_NUMERICAL_ERROR
DSDPSolutionType
Formulations (P) and (D) can be feasible and bounded, feasible and unbounded, or infeasible.
@ DSDP_UNBOUNDED
@ DSDP_PDFEASIBLE
@ DSDP_PDUNKNOWN
@ DSDP_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.
Definition dsdpvec.h:25
int DSDPSaveYForX(DSDP dsdp, double mu, double pstep)
Save the current solution for later computation of X.
Definition dsdpx.c:149
Internal structures for the DSDP solver.
Definition dsdp.h:65