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"
11 int 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"
55 int 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;
59  DSDPTerminationReason reason;
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"
128 int 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"
149 int 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"
232 int 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"
255  DSDPFunctionBegin;
256  DSDPValid(dsdp);
257  *pdfeasible=dsdp->pdfeasible;
258  DSDPFunctionReturn(0);
259 }
276 #undef __FUNCT__
277 #define __FUNCT__ "DSDPGetTraceX"
278 int 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"
297 int 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"
343 int 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"
365 int 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"
386 int 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"
409 int 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"
434 int 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"
455 int 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"
483 int 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"
511 int 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 
int DSDPGetFinalErrors(DSDP dsdp, double err[6])
Copy six different error measurements into an array.
Definition: dsdpx.c:297
int DSDPGetDYMakeX(DSDP dsdp, double dy[], int m)
Copies the variables dy used to construct X into an array.
Definition: dsdpx.c:483
struct DSDPVec_C DSDPVec
This object hold m+2 variables: a scaling of C, the y variables, and r.
Definition: dsdpvec.h:25
Error handling, printing, and profiling.
int DSDPGetPTolerance(DSDP dsdp, double *inftol)
Copy the feasibility tolerance.
Definition: dsdpx.c:386
Internal structures for the DSDP solver.
Definition: dsdp.h:65
DSDPTerminationReason
There are many reasons to terminate the solver.
int DSDPComputeDualityGap(DSDP, double, double *)
Compute the current duality gap.
Definition: dualimpl.c:230
The API to DSDP for those applications using DSDP as a subroutine library.
int DSDPSetRTolerance(DSDP dsdp, double inftol)
Classify (D) as feasible only if the variable r is less than this tolerance.
Definition: dsdpx.c:409
int DSDPSetPTolerance(DSDP dsdp, double inftol)
Classify (P) as feasible only if the infeasibility is less than this tolerance.
Definition: dsdpx.c:365
int DSDPGetSolutionType(DSDP dsdp, DSDPSolutionType *pdfeasible)
Solutions can be bounded, infeasible, or unbounded.
Definition: dsdpx.c:254
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 DSDPGetPInfeasibility(DSDP dsdp, double *pperror)
Copy the infeasibility in (P).
Definition: dsdpx.c:343
int DSDPGetDObjective(DSDP dsdp, double *dobj)
Copy the objective value (D).
Definition: dsdpsetdata.c:502
int DSDPStopReason(DSDP dsdp, DSDPTerminationReason *reason)
Copy the reason why the solver terminated.
Definition: dsdpsetdata.c:582
int DSDPComputeDY(DSDP, double, DSDPVec, double *)
Compute the step direction.
Definition: dualimpl.c:45
int DSDPComputeXVariables(DSDP, double, DSDPVec, DSDPVec, DSDPVec, double *)
Compute the X variables in each cone.
Definition: dsdpcops.c:654
int DSDPGetPObjective(DSDP dsdp, double *pobj)
Copy the objective value (P).
Definition: dsdpx.c:232
int DSDPGetRR(DSDP, double *)
Get variable r.
Definition: dualimpl.c:361
int DSDPGetMuMakeX(DSDP dsdp, double *mu)
Copies the value of mu used to construct X.
Definition: dsdpx.c:511
int DSDPSaveYForX(DSDP dsdp, double mu, double pstep)
Save the current solution for later computation of X.
Definition: dsdpx.c:149
int DSDPGetDDObjective(DSDP dsdp, double *ddobj)
Copy the objective value (DD).
Definition: dsdpsetdata.c:523
int DSDPGetRTolerance(DSDP dsdp, double *inftol)
Copy the maximum infeasibility allowed (D).
Definition: dsdpx.c:434
int DSDPGetTraceX(DSDP dsdp, double *tracex)
Copy the trace of the variables X in (P).
Definition: dsdpx.c:278
int DSDPGetR(DSDP dsdp, double *res)
Copy the infeasibility in (D), or the variable r in (DD).
Definition: dsdpsetdata.c:601
int DSDPComputeX(DSDP dsdp)
Compute the X variables.
Definition: dsdpx.c:55
int DSDPGetScale(DSDP dsdp, double *scale)
Copy the internal scaling factor from the solver.
Definition: dsdpsetdata.c:128
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 DSDPSetConvergenceFlag(DSDP dsdp, DSDPTerminationReason reason)
Monitor each iteration of the solver.
Definition: dsdpsetdata.c:968
int DSDPGetYMakeX(DSDP dsdp, double y[], int m)
Copies the variables y used to construct X into an array.
Definition: dsdpx.c:455
DSDPSolutionType
Formulations (P) and (D) can be feasible and bounded, feasible and unbounded, or infeasible.
int DSDPGetMaxYElement(DSDP, double *)
Copy the the infinity norm of the variables y.
Definition: dsdpsetdata.c:645