DSDP
dualimpl.c
Go to the documentation of this file.
1 #include "dsdp.h"
2 #include "dsdpsys.h"
19 #undef __FUNCT__
20 #define __FUNCT__ "DSDPComputeObjective"
21 int DSDPComputeObjective(DSDP dsdp, DSDPVec Y, double *ddobj){
22  int info;
23  DSDPFunctionBegin;
24  info = DSDPVecDot(Y,dsdp->b,ddobj);DSDPCHKERR(info);
25  DSDPFunctionReturn(0);
26 }
27 
28 
43 #undef __FUNCT__
44 #define __FUNCT__ "DSDPComputeDY"
45 int DSDPComputeDY(DSDP dsdp, double mu, DSDPVec DY, double *pnorm){
46  int info;
47  double ppnorm,ddy1=fabs(1.0/mu*dsdp->schurmu),ddy2=-1.0;
48  DSDPFunctionBegin;
49  info=DSDPComputeRHS(dsdp,mu,dsdp->rhs); DSDPCHKERR(info);
50  info=DSDPVecWAXPBY(DY,ddy1,dsdp->dy1,ddy2,dsdp->dy2);DSDPCHKERR(info);
51  info=DSDPComputePNorm(dsdp,mu,DY,&ppnorm);DSDPCHKERR(info);
52  if (ppnorm<0){ /* If pnorm < 0 there are SMW numerical issues */
53  DSDPLogInfo(0,2,"Problem with PNORM: %4.4e < 0 \n",ppnorm);
54  /* ppnorm=1.0; */
55  }
56  *pnorm=ppnorm;
57  DSDPFunctionReturn(0);
58 }
59 
75 #undef __FUNCT__
76 #define __FUNCT__ "DSDPComputePDY"
77 int DSDPComputePDY(DSDP dsdp, double mu, DSDPVec DY, double *pnorm){
78  int info;
79  double ppnorm,ddy1=-fabs(1.0/mu*dsdp->schurmu),ddy2=1.0;
80  DSDPFunctionBegin;
81  info=DSDPComputeRHS(dsdp,mu,dsdp->rhs); DSDPCHKERR(info);
82  info=DSDPVecWAXPBY(DY,ddy1,dsdp->dy1,ddy2,dsdp->dy2);DSDPCHKERR(info);
83  info=DSDPComputePNorm(dsdp,mu,DY,&ppnorm);DSDPCHKERR(info);
84  if (ppnorm<0){ /* If pnorm < 0 there are SMW numerical issues */
85  DSDPLogInfo(0,2,"Problem with PNORM: %4.4e < 0 \n",ppnorm);
86  /* ppnorm=1.0; */
87  }
88  *pnorm=ppnorm;
89  DSDPFunctionReturn(0);
90 }
91 
103 #undef __FUNCT__
104 #define __FUNCT__ "DSDPComputePDY1"
105 int DSDPComputePDY1(DSDP dsdp, double mur, DSDPVec DY1){
106  int info;
107  double ddy1=-fabs(mur*dsdp->schurmu);
108  DSDPFunctionBegin;
109  info=DSDPVecScaleCopy(dsdp->dy1,ddy1,DY1); DSDPCHKERR(info);
110  DSDPFunctionReturn(0);
111 }
112 
123 #undef __FUNCT__
124 #define __FUNCT__ "DSDPComputeNewY"
125 int DSDPComputeNewY(DSDP dsdp, double beta, DSDPVec Y){
126  int info;
127  double rtemp;
128  DSDPFunctionBegin;
129  info=DSDPVecWAXPY(Y,beta,dsdp->dy,dsdp->y);DSDPCHKERR(info);
130  info=DSDPVecGetR(Y,&rtemp);DSDPCHKERR(info);
131  rtemp=DSDPMin(0,rtemp);
132  info=DSDPSchurMatSetR(dsdp->M,rtemp);DSDPCHKERR(info);
133  info=DSDPVecSetR(Y,rtemp);DSDPCHKERR(info);
134  info=DSDPApplyFixedVariables(dsdp->M,Y);DSDPCHKERR(info);
135  DSDPFunctionReturn(0);
136 }
137 
148 #undef __FUNCT__
149 #define __FUNCT__ "DSDPComputePY"
150 int DSDPComputePY(DSDP dsdp, double beta, DSDPVec PY){
151  int info;
152  DSDPFunctionBegin;
153  info=DSDPVecWAXPY(PY,beta,dsdp->dy,dsdp->y);DSDPCHKERR(info);
154  info=DSDPApplyFixedVariables(dsdp->M,PY);DSDPCHKERR(info);
155  DSDPFunctionReturn(0);
156 }
157 
175 #undef __FUNCT__
176 #define __FUNCT__ "DSDPComputeRHS"
177 int DSDPComputeRHS(DSDP dsdp, double mu, DSDPVec RHS){
178  int info;
179  double ddrhs1=1.0/mu*dsdp->schurmu,ddrhs2=-( mu/fabs(mu) );
180  DSDPFunctionBegin;
181  info=DSDPVecWAXPBY(RHS,ddrhs1,dsdp->rhs1,ddrhs2,dsdp->rhs2);DSDPCHKERR(info);
182  DSDPFunctionReturn(0);
183 }
184 
185 #undef __FUNCT__
186 #define __FUNCT__ "DSDPComputePNorm"
187 
200 int DSDPComputePNorm(DSDP dsdp, double mu, DSDPVec DY, double *pnorm){
201  int info;
202  double ppnorm=0;
203  DSDPFunctionBegin;
204  info=DSDPComputeRHS(dsdp,mu,dsdp->rhs); DSDPCHKERR(info);
205  info = DSDPVecDot(dsdp->rhs,DY,&ppnorm);DSDPCHKERR(info);
206  ppnorm/=dsdp->schurmu;
207  if (ppnorm >= 0){ /* Theoretically True */
208  *pnorm=sqrt(ppnorm);
209  } else {
210  DSDPLogInfo(0,2,"Problem with PNORM: %4.4e is not positive.\n",ppnorm);
211  *pnorm=ppnorm;
212  }
213  if (*pnorm!=*pnorm){DSDPSETERR1(1,"Problem with PNORM: %4.4e is not positive.\n",ppnorm);}
214  DSDPFunctionReturn(0);
215 }
216 
228 #undef __FUNCT__
229 #define __FUNCT__ "DSDPComputeDualityGap"
230 int DSDPComputeDualityGap(DSDP dsdp, double mu, double *gap){
231  int info;
232  double newgap=0,pnorm;
233  double smu=1.0/dsdp->schurmu;
234  DSDPFunctionBegin;
235  info=DSDPComputeDY(dsdp,mu,dsdp->dy,&pnorm); DSDPCHKERR(info);
236  info=DSDPVecDot(dsdp->dy,dsdp->rhs2,&newgap);DSDPCHKERR(info);
237  newgap = (newgap*smu+dsdp->np)*mu;
238  if (newgap<=0){
239  DSDPLogInfo(0,2,"GAP :%4.4e<0: Problem\n",newgap);
240  } else {
241  DSDPLogInfo(0,2,"Duality Gap: %12.8e, Update primal objective: %12.8e\n",newgap,dsdp->ddobj+newgap);
242  }
243  newgap=DSDPMax(0,newgap);
244  *gap=newgap;
245  DSDPFunctionReturn(0);
246 }
247 
259 #undef __FUNCT__
260 #define __FUNCT__ "DSDPComputePotential"
261 int DSDPComputePotential(DSDP dsdp, DSDPVec y, double logdet, double *potential){
262  int info;
263  double dpotential,gap,ddobj;
264  DSDPFunctionBegin;
265  info=DSDPComputeObjective(dsdp,y,&ddobj);DSDPCHKERR(info);
266  gap=dsdp->ppobj-ddobj;
267  if (gap>0) dpotential=dsdp->rho*log(gap)-logdet;
268  else {dpotential=dsdp->potential+1;}
269  *potential=dpotential;
270  DSDPLogInfo(0,9,"Gap: %4.4e, Log Determinant: %4.4e, Log Gap: %4.4e\n",gap,logdet,log(gap));
271  DSDPFunctionReturn(0);
272 }
273 
285 #undef __FUNCT__
286 #define __FUNCT__ "DSDPComputePotential2"
287 int DSDPComputePotential2(DSDP dsdp, DSDPVec y, double mu, double logdet, double *potential){
288  int info;
289  double ddobj;
290  DSDPFunctionBegin;
291  info=DSDPComputeObjective(dsdp,y,&ddobj);DSDPCHKERR(info);
292  *potential=-(ddobj + mu*logdet)*dsdp->schurmu;
293  *potential=-(ddobj/mu + logdet)*dsdp->schurmu;
294  DSDPFunctionReturn(0);
295 }
296 
297 
307 #undef __FUNCT__
308 #define __FUNCT__ "DSDPSetY"
309 int DSDPSetY(DSDP dsdp, double beta, double logdet, DSDPVec ynew){
310  int info;
311  double r1,r2,rr,pp;
312  DSDPFunctionBegin;
313  info=DSDPVecGetR(dsdp->y,&r1);DSDPCHKERR(info);
314  info=DSDPVecGetR(ynew,&r2);DSDPCHKERR(info);
315  if (r2==0&&r1!=0){dsdp->rflag=1;} else {dsdp->rflag=0;};
316  info=DSDPVecCopy(ynew,dsdp->y);DSDPCHKERR(info);
317  info = DSDPComputeObjective(dsdp,dsdp->y,&dsdp->ddobj);DSDPCHKERR(info);
318  /* Correct ppobj if ppobj < ddobj , which can happen when dual infeasibility is present */
319  if (dsdp->ppobj<=dsdp->ddobj){
320  dsdp->ppobj=dsdp->ddobj+2*dsdp->mu * dsdp->np;
321  DSDPLogInfo(0,2,"Primal Objective Not Right. Assigned: %8.8e\n",dsdp->ppobj);
322  }
323  info=DSDPVecGetR(ynew,&rr);DSDPCHKERR(info);
324  info=DSDPVecGetR(dsdp->b,&pp);DSDPCHKERR(info);
325  dsdp->dobj=dsdp->ddobj-rr*pp;
326  DSDPLogInfo(0,2,"Duality Gap: %4.4e, Potential: %4.4e \n",dsdp->dualitygap,dsdp->potential);
327  dsdp->dualitygap=dsdp->ppobj-dsdp->ddobj;
328  dsdp->mu=(dsdp->dualitygap)/(dsdp->np);
329  dsdp->dstep=beta;
330  dsdp->logdet=logdet;
331  info=DSDPComputePotential(dsdp,dsdp->y,dsdp->logdet,&dsdp->potential);DSDPCHKERR(info);
332  DSDPLogInfo(0,2,"Duality Gap: %4.4e, Potential: %4.4e \n",dsdp->dualitygap,dsdp->potential);
333  DSDPFunctionReturn(0);
334 }
335 
336 
337 #undef __FUNCT__
338 #define __FUNCT__ "DSDPSetRR"
339 
345 int DSDPSetRR(DSDP dsdp,double res){
346  int info;
347  DSDPFunctionBegin;
348  DSDPValid(dsdp);
349  info=DSDPVecSetR(dsdp->y,-res);DSDPCHKERR(info);
350  DSDPFunctionReturn(0);
351 }
352 
353 #undef __FUNCT__
354 #define __FUNCT__ "DSDPGetRR"
355 
361 int DSDPGetRR(DSDP dsdp,double *res){
362  int info;
363  DSDPFunctionBegin;
364  DSDPValid(dsdp);
365  info=DSDPVecGetR(dsdp->y,res);DSDPCHKERR(info);
366  *res=-*res;
367  if (*res==0) *res=0;
368  DSDPFunctionReturn(0);
369 }
370 
371 
372 #undef __FUNCT__
373 #define __FUNCT__ "DSDPObjectiveGH"
374 
381 int DSDPObjectiveGH( DSDP dsdp , DSDPSchurMat M, DSDPVec vrhs1){
382  int i,info,m;
383  double rtemp,dd;
384 
385  DSDPFunctionBegin;
386  info=DSDPVecGetSize(vrhs1,&m); DSDPCHKERR(info);
387  for (i=0;i<m;i++){
388  info=DSDPSchurMatVariableCompute(M,i,&dd); DSDPCHKERR(info);
389  if (dd){
390  info=DSDPVecGetElement(dsdp->b,i,&rtemp);DSDPCHKERR(info);
391  info=DSDPVecAddElement(vrhs1,i,rtemp);DSDPCHKERR(info);
392  }
393  }
394  DSDPFunctionReturn(0);
395 }
396 
397 #undef __FUNCT__
398 #define __FUNCT__ "DSDPCheckForUnboundedObjective"
399 int DSDPCheckForUnboundedObjective(DSDP dsdp, DSDPTruth *unbounded){
400  int info;
401  double dtemp;
402  DSDPTruth psdefinite;
403  DSDPFunctionBegin;
404  *unbounded=DSDP_FALSE;
405  info = DSDPVecDot(dsdp->b,dsdp->dy2,&dtemp);DSDPCHKERR(info);
406  if ( dtemp < 0 /* && dsdp->r==0 && dsdp->ddobj > 0 */) {
407  info = DSDPVecScaleCopy(dsdp->dy2,-1.0,dsdp->ytemp); DSDPCHKERR(info);
408  info = DSDPComputeSS(dsdp,dsdp->ytemp,PRIMAL_FACTOR,&psdefinite);DSDPCHKERR(info);
409  if (psdefinite == DSDP_TRUE){
410  psdefinite=DSDP_FALSE;
411  while(psdefinite==DSDP_FALSE){ /* Dual point should be a certificate of dual unboundedness, and be dual feasible */
412  info=DSDPComputeSS(dsdp,dsdp->ytemp,PRIMAL_FACTOR,&psdefinite);DSDPCHKERR(info);
413  if (psdefinite == DSDP_TRUE) break;
414  info=DSDPVecScale(2.0,dsdp->ytemp); DSDPCHKERR(info);
415  }
416  info = DSDPVecCopy(dsdp->ytemp,dsdp->y); DSDPCHKERR(info);
417  info = DSDPSaveYForX(dsdp,1.0e-12,1.0);DSDPCHKERR(info);
418  info = DSDPComputeObjective(dsdp,dsdp->y,&dsdp->ddobj);DSDPCHKERR(info);
419  info = DSDPVecNormalize(dsdp->y); DSDPCHKERR(info);
420  *unbounded=DSDP_TRUE;
421  }
422  }
423  DSDPFunctionReturn(0);
424 }
425 
int DSDPObjectiveGH(DSDP dsdp, DSDPSchurMat M, DSDPVec vrhs1)
Compute gradient of dual objective.
Definition: dualimpl.c:381
DSDPTruth
Boolean variables.
int DSDPComputePY(DSDP dsdp, double beta, DSDPVec PY)
Compute PY = Y - beta DY for use in computing X.
Definition: dualimpl.c:150
int DSDPComputePDY(DSDP dsdp, double mu, DSDPVec DY, double *pnorm)
Compute the step direction.
Definition: dualimpl.c:77
struct DSDPVec_C DSDPVec
This object hold m+2 variables: a scaling of C, the y variables, and r.
Definition: dsdpvec.h:25
Schur complement matrix whose solution is the Newton direction.
Definition: dsdpschurmat.h:35
int DSDPComputeSS(DSDP, DSDPVec, DSDPDualFactorMatrix, DSDPTruth *)
Compute the dual variables S in each cone.
Definition: dsdpcops.c:272
int DSDPComputePotential(DSDP dsdp, DSDPVec y, double logdet, double *potential)
Compute the potential of the given point.
Definition: dualimpl.c:261
Error handling, printing, and profiling.
int DSDPComputePDY1(DSDP dsdp, double mur, DSDPVec DY1)
Compute an affine step direction dy1.
Definition: dualimpl.c:105
int DSDPSetY(DSDP dsdp, double beta, double logdet, DSDPVec ynew)
Update the solver with these y variables.
Definition: dualimpl.c:309
int DSDPGetRR(DSDP dsdp, double *res)
Get variable r.
Definition: dualimpl.c:361
Internal structures for the DSDP solver.
Definition: dsdp.h:65
int DSDPComputeObjective(DSDP dsdp, DSDPVec Y, double *ddobj)
Compute the objective function (DD).
Definition: dualimpl.c:21
int DSDPComputeRHS(DSDP dsdp, double mu, DSDPVec RHS)
Compute the right-hand side of the linear system that determines the step direction.
Definition: dualimpl.c:177
Internal data structure for the DSDP solver.
int DSDPSaveYForX(DSDP, double, double)
Save the current solution for later computation of X.
Definition: dsdpx.c:149
int DSDPComputePotential2(DSDP dsdp, DSDPVec y, double mu, double logdet, double *potential)
Compute the objective function plus the barrier function.
Definition: dualimpl.c:287
int DSDPComputeDY(DSDP dsdp, double mu, DSDPVec DY, double *pnorm)
Compute the step direction.
Definition: dualimpl.c:45
int DSDPComputePNorm(DSDP dsdp, double mu, DSDPVec DY, double *pnorm)
Compute proximity to a point on the central path.
Definition: dualimpl.c:200
int DSDPSchurMatSetR(DSDPSchurMat M, double rr)
Set up the data structure.
Definition: dsdpschurmat.c:338
int DSDPSetRR(DSDP dsdp, double res)
Set variable r.
Definition: dualimpl.c:345
int DSDPSchurMatVariableCompute(DSDPSchurMat, int, double *)
Determine with the cone should compute this diagonal element of M and RHS.
int DSDPComputeDualityGap(DSDP dsdp, double mu, double *gap)
Compute the current duality gap.
Definition: dualimpl.c:230
int DSDPComputeNewY(DSDP dsdp, double beta, DSDPVec Y)
Update the Y variables.
Definition: dualimpl.c:125