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"
21int 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"
45int 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"
77int 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"
105int 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"
125int 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"
150int 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"
177int 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"
200int 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"
230int 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"
261int 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"
287int 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"
309int 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"
345int 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"
361int 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"
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"
399int 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
Internal data structure for the DSDP solver.
int DSDPComputeSS(DSDP, DSDPVec, DSDPDualFactorMatrix, DSDPTruth *)
Compute the dual variables S in each cone.
Definition dsdpcops.c:272
int DSDPSaveYForX(DSDP, double, double)
Save the current solution for later computation of X.
Definition dsdpx.c:149
@ PRIMAL_FACTOR
DSDPTruth
Boolean variables.
@ DSDP_FALSE
@ DSDP_TRUE
int DSDPSchurMatSetR(DSDPSchurMat M, double rr)
Set up the data structure.
int DSDPSchurMatVariableCompute(DSDPSchurMat, int, double *)
Determine with the cone should compute this diagonal element of M and RHS.
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 DSDPComputePotential(DSDP dsdp, DSDPVec y, double logdet, double *potential)
Compute the potential of the given point.
Definition dualimpl.c:261
int DSDPComputePY(DSDP dsdp, double beta, DSDPVec PY)
Compute PY = Y - beta DY for use in computing X.
Definition dualimpl.c:150
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
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
int DSDPComputePDY1(DSDP dsdp, double mur, DSDPVec DY1)
Compute an affine step direction dy1.
Definition dualimpl.c:105
int DSDPComputePDY(DSDP dsdp, double mu, DSDPVec DY, double *pnorm)
Compute the step direction.
Definition dualimpl.c:77
int DSDPSetRR(DSDP dsdp, double res)
Set variable r.
Definition dualimpl.c:345
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 DSDPGetRR(DSDP dsdp, double *res)
Get variable r.
Definition dualimpl.c:361
int DSDPComputePNorm(DSDP dsdp, double mu, DSDPVec DY, double *pnorm)
Compute proximity to a point on the central path.
Definition dualimpl.c:200
int DSDPComputeObjective(DSDP dsdp, DSDPVec Y, double *ddobj)
Compute the objective function (DD).
Definition dualimpl.c:21
int DSDPSetY(DSDP dsdp, double beta, double logdet, DSDPVec ynew)
Update the solver with these y variables.
Definition dualimpl.c:309
int DSDPComputeDY(DSDP dsdp, double mu, DSDPVec DY, double *pnorm)
Compute the step direction.
Definition dualimpl.c:45
int DSDPObjectiveGH(DSDP dsdp, DSDPSchurMat M, DSDPVec vrhs1)
Compute gradient of dual objective.
Definition dualimpl.c:381
Schur complement matrix whose solution is the Newton direction.
Internal structures for the DSDP solver.
Definition dsdp.h:65