DSDP
dsdpsetup.c
Go to the documentation of this file.
1#include "dsdp.h"
2#include "dsdpsys.h"
3#include "dsdp5.h"
28#undef __FUNCT__
29#define __FUNCT__ "DSDPCreate"
30int DSDPCreate(int m,DSDP* dsdpnew){
31
32 DSDP dsdp;
33 int info;
34
35 DSDPFunctionBegin;
36
37 DSDPCALLOC1(&dsdp,PD_DSDP,&info);DSDPCHKERR(info);
38 *dsdpnew=dsdp;
39 dsdp->keyid=DSDPKEY;
40
41 /* Initialize some parameters */
42 DSDPEventLogInitialize();
43 dsdp->m=m;
44 dsdp->maxcones=0;
45 dsdp->ncones=0;
46 dsdp->K=0;
47 dsdp->setupcalled=DSDP_FALSE;
48 dsdp->ybcone=0;
49 dsdp->ndroutines=0;
50 /* info = DSDPSetStandardMonitor(dsdp);DSDPCHKERR(info); */
51 info = DSDPVecCreateSeq(m+2,&dsdp->b);DSDPCHKERR(info);
52 info = DSDPVecZero(dsdp->b);DSDPCHKERR(info);
53 info = DSDPVecDuplicate(dsdp->b,&dsdp->y);DSDPCHKERR(info);
54 info = DSDPVecDuplicate(dsdp->b,&dsdp->ytemp);DSDPCHKERR(info);
55 info = DSDPVecZero(dsdp->y);DSDPCHKERR(info);
56 info = DSDPVecSetC(dsdp->y,-1.0);DSDPCHKERR(info);
57
58 info = DSDPAddRCone(dsdp,&dsdp->rcone);DSDPCHKERR(info);
59 info = DSDPCreateLUBoundsCone(dsdp,&dsdp->ybcone);DSDPCHKERR(info);
60
61 info=DSDPSetDefaultStatistics(dsdp);DSDPCHKERR(info);
62 info=DSDPSetDefaultParameters(dsdp);DSDPCHKERR(info);
63 info=DSDPSetDefaultMonitors(dsdp);DSDPCHKERR(info);
64
65 /* info = DSDPMatInitialize(m,m,&dsdp->Q);DSDPCHKERR(info); */
66 info = DSDPSchurMatInitialize(&dsdp->M);DSDPCHKERR(info);
67 info = DSDPSetDefaultSchurMatrixStructure(dsdp); DSDPCHKERR(info);
68 info = DSDPCGInitialize(&dsdp->sles); DSDPCHKERR(info);
69
70 /* Set the one global variable
71 sdat=dsdp;
72 */
73 DSDPFunctionReturn(0);
74}
75
76
77#undef __FUNCT__
78#define __FUNCT__ "DSDPSetDefaultStatistics"
85
86 int i;
87 DSDPFunctionBegin;
88 DSDPValid(dsdp);
89 dsdp->reason=CONTINUE_ITERATING;
90 dsdp->pdfeasible=DSDP_PDUNKNOWN;
91 dsdp->itnow=0;
92 dsdp->pobj= 1.0e10;
93 dsdp->ppobj= 1.0e10;
94 dsdp->dobj= -1.0e+9;
95 dsdp->ddobj= -1.0e+9;
96 dsdp->dualitygap=dsdp->ppobj-dsdp->ddobj;
97 dsdp->pstep=1.0;
98 dsdp->dstep=0.0;
99 for (i=0;i<MAX_XMAKERS;i++){
100 dsdp->xmaker[i].mu=1.0e200;
101 dsdp->xmaker[i].pstep=0.0;
102 }
103 dsdp->pnorm=0.001;
104 dsdp->mu=1000.0;
105 dsdp->np=0;
106 dsdp->anorm=0;
107 dsdp->bnorm=0;
108 dsdp->cnorm=0;
109 dsdp->tracex=0;
110 dsdp->tracexs=0;
111 dsdp->Mshift=0;
112 dsdp->goty0=DSDP_FALSE;
113 DSDPFunctionReturn(0);
114}
115#undef __FUNCT__
116#define __FUNCT__ "DSDPSetDefaultParameters"
123
124 int info;
125 DSDPFunctionBegin;
126 DSDPValid(dsdp);
127
128 /* Stopping parameters */
129 info=DSDPSetMaxIts(dsdp,500);DSDPCHKERR(info);
130 info=DSDPSetGapTolerance(dsdp,1.0e-6);DSDPCHKERR(info);
131 info=DSDPSetPNormTolerance(dsdp,1.0e30);DSDPCHKERR(info);
132 if (dsdp->m<100){info=DSDPSetGapTolerance(dsdp,1.0e-7);DSDPCHKERR(info);}
133 if (dsdp->m>3000){info=DSDPSetGapTolerance(dsdp,5.0e-6);DSDPCHKERR(info);}
134 info=RConeSetType(dsdp->rcone,DSDPInfeasible);DSDPCHKERR(info);
135 info=DSDPSetDualBound(dsdp,1.0e20);DSDPCHKERR(info);
136 info=DSDPSetStepTolerance(dsdp,5.0e-2);DSDPCHKERR(info);
137 info=DSDPSetRTolerance(dsdp,1.0e-6);DSDPCHKERR(info);
138 info=DSDPSetPTolerance(dsdp,1.0e-4);DSDPCHKERR(info);
139 /* Solver options */
140 info=DSDPSetMaxTrustRadius(dsdp,1.0e10);DSDPCHKERR(info);
141 info=DSDPUsePenalty(dsdp,0);DSDPCHKERR(info);
142 info=DSDPSetInitialBarrierParameter(dsdp,-1.0);DSDPCHKERR(info);
143 info=DSDPSetPotentialParameter(dsdp,3.0);DSDPCHKERR(info);
144 info=DSDPUseDynamicRho(dsdp,1);DSDPCHKERR(info);
145 info=DSDPSetR0(dsdp,-1.0);DSDPCHKERR(info);
146 info=DSDPSetPenaltyParameter(dsdp,1.0e8);DSDPCHKERR(info);
147 info=DSDPReuseMatrix(dsdp,4);DSDPCHKERR(info);
148 if (dsdp->m>100){info=DSDPReuseMatrix(dsdp,7);DSDPCHKERR(info);}
149 if (dsdp->m>1000){info=DSDPReuseMatrix(dsdp,10);DSDPCHKERR(info);}
150 if (dsdp->m<=100){info=DSDPSetPotentialParameter(dsdp,5.0);DSDPCHKERR(info);DSDPCHKERR(info);}
151 dsdp->maxschurshift=1.0e-11;
152 dsdp->mu0=-1.0;
153 dsdp->slestype=2;
154 info = DSDPSetYBounds(dsdp,-1e7,1e7);DSDPCHKERR(info);
155 DSDPFunctionReturn(0);
156}
157
158#undef __FUNCT__
159#define __FUNCT__ "DSDPSetDefaultMonitors"
166
167 int info;
168
169 DSDPFunctionBegin;
170 DSDPValid(dsdp);
171 dsdp->nmonitors=0;
172 info=DSDPSetMonitor(dsdp,DSDPDefaultConvergence,(void*)&dsdp->conv); DSDPCHKERR(info);
173 DSDPFunctionReturn(0);
174}
175
191#undef __FUNCT__
192#define __FUNCT__ "DSDPSetUp"
193int DSDPSetup(DSDP dsdp){
194
195 int i,info;
196 DSDPFunctionBegin;
197 DSDPValid(dsdp);
198
199 /* Create the Work Vectors */
200 info = DSDPVecDuplicate(dsdp->y,&dsdp->rhs1);DSDPCHKERR(info);
201 info = DSDPVecDuplicate(dsdp->y,&dsdp->rhs2);DSDPCHKERR(info);
202 info = DSDPVecDuplicate(dsdp->y,&dsdp->rhs);DSDPCHKERR(info);
203 info = DSDPVecDuplicate(dsdp->y,&dsdp->rhstemp);DSDPCHKERR(info);
204 info = DSDPVecDuplicate(dsdp->y,&dsdp->dy1);DSDPCHKERR(info);
205 info = DSDPVecDuplicate(dsdp->y,&dsdp->dy2);DSDPCHKERR(info);
206 info = DSDPVecDuplicate(dsdp->y,&dsdp->dy);DSDPCHKERR(info);
207 info = DSDPVecDuplicate(dsdp->y,&dsdp->y0);DSDPCHKERR(info);
208 info = DSDPVecDuplicate(dsdp->y,&dsdp->xmakerrhs);DSDPCHKERR(info);
209 for (i=0;i<MAX_XMAKERS;i++){
210 info = DSDPVecDuplicate(dsdp->y,&dsdp->xmaker[i].y);DSDPCHKERR(info);
211 info = DSDPVecDuplicate(dsdp->y,&dsdp->xmaker[i].dy);DSDPCHKERR(info);
212 info = DSDPVecDuplicate(dsdp->y,&dsdp->xmaker[i].rhs);DSDPCHKERR(info);
213 }
214
215 /* Create M */
216 info = DSDPSetUpCones(dsdp);DSDPCHKERR(info);
217 info = DSDPSchurMatSetup(dsdp->M,dsdp->ytemp);DSDPCHKERR(info);
218
219 info = DSDPCGSetup(dsdp->sles,dsdp->ytemp); DSDPCHKERR(info);
220
221 info = DSDPSetUpCones2(dsdp,dsdp->y,dsdp->M);DSDPCHKERR(info);
222 info = DSDPGetConicDimension(dsdp,&dsdp->np);DSDPCHKERR(info);
223
224 info=DSDPComputeDataNorms(dsdp);DSDPCHKERR(info);
225 dsdp->pinfeas=dsdp->bnorm+1;
226 dsdp->perror=dsdp->bnorm+1;
227 info=DSDPScaleData(dsdp);DSDPCHKERR(info);
228
229 info=DSDPGetConicDimension(dsdp,&dsdp->np);DSDPCHKERR(info);
230 dsdp->solvetime=0;
231 dsdp->cgtime=0;
232 dsdp->ptime=0;
233 dsdp->dtime=0;
234 dsdp->ctime=0;
235 info=DSDPEventLogRegister("Primal Step",&dsdp->ptime);
236 info=DSDPEventLogRegister("Dual Step",&dsdp->dtime);
237 info=DSDPEventLogRegister("Corrector Step",&dsdp->ctime);
238 info=DSDPEventLogRegister("CG Solve",&dsdp->cgtime);
239 info=DSDPEventLogRegister("DSDP Solve",&dsdp->solvetime);
240 dsdp->setupcalled=DSDP_TRUE;
241 DSDPFunctionReturn(0);
242}
243
244
245
246#undef __FUNCT__
247#define __FUNCT__ "DSDPGetSchurMatrix"
248int DSDPGetSchurMatrix(DSDP dsdp, DSDPSchurMat *M){
249 DSDPFunctionBegin;
250 DSDPValid(dsdp);
251 *M=dsdp->M;
252 DSDPFunctionReturn(0);
253}
254
255#undef __FUNCT__
256#define __FUNCT__ "DSDPGetConvergenceMonitor"
268int DSDPGetConvergenceMonitor(DSDP dsdp, ConvergenceMonitor**ctx){
269 DSDPFunctionBegin;
270 DSDPValid(dsdp);
271 *ctx=&dsdp->conv;
272 DSDPFunctionReturn(0);
273}
274
275
276#undef __FUNCT__
277#define __FUNCT__ "DSDPComputeDataNorms"
284 int info;
285 DSDPVec ytemp=dsdp->ytemp;
286 DSDPFunctionBegin;
287 DSDPValid(dsdp);
288 info = DSDPComputeANorm2(dsdp,ytemp);DSDPCHKERR(info);
289 info = DSDPFixedVariablesNorm(dsdp->M,ytemp);DSDPCHKERR(info);
290 info = DSDPVecGetC(ytemp,&dsdp->cnorm);DSDPCHKERR(info);
291 dsdp->cnorm=sqrt(dsdp->cnorm);
292 info = DSDPVecSetR(ytemp,0);DSDPCHKERR(info);
293 info = DSDPVecSetC(ytemp,0);DSDPCHKERR(info);
294 info = DSDPVecNorm1(ytemp,&dsdp->anorm);DSDPCHKERR(info);
295 dsdp->anorm=sqrt(dsdp->anorm);
296 DSDPLogInfo(0,2,"Norm of data: %4.2e\n",dsdp->anorm);
297 info=DSDPVecCopy(dsdp->b,ytemp);DSDPCHKERR(info);
298 info = DSDPVecSetR(ytemp,0);DSDPCHKERR(info);
299 info = DSDPVecSetC(ytemp,0);DSDPCHKERR(info);
300 info = DSDPVecNorm2(ytemp,&dsdp->bnorm);DSDPCHKERR(info);
301 DSDPFunctionReturn(0);
302}
303
304#undef __FUNCT__
305#define __FUNCT__ "DSDPScaleData"
312 int info;
313 double scale;
314 DSDPFunctionBegin;
315 DSDPValid(dsdp);
316 scale=1.0*dsdp->anorm;
317 if (dsdp->bnorm){ scale/=dsdp->bnorm;}
318 if (dsdp->cnorm){ scale/=dsdp->cnorm;}
319 scale=DSDPMin(scale,1.0);
320 scale=DSDPMax(scale,1.0e-6);
321 if (dsdp->cnorm==0){ scale=1;}
322 info=DSDPSetScale(dsdp,scale);DSDPCHKERR(info);
323 DSDPFunctionReturn(0);
324}
325
341#undef __FUNCT__
342#define __FUNCT__ "DSDPSolve"
343int DSDPSolve(DSDP dsdp){
344 int info;
345 DSDPFunctionBegin;
346 info=DSDPEventLogBegin(dsdp->solvetime);
347 dsdp->pdfeasible=DSDP_PDUNKNOWN;
348 info=DSDPSetConvergenceFlag(dsdp,CONTINUE_ITERATING);DSDPCHKERR(info);
349 info=DSDPInitializeVariables(dsdp);DSDPCHKERR(info);
350 info=DSDPSolveDynamicRho(dsdp);DSDPCHKERR(info);
351 if (dsdp->pstep==1){info=DSDPRefineStepDirection(dsdp,dsdp->xmakerrhs,dsdp->xmaker[0].dy);DSDPCHKERR(info);}
352 if (dsdp->pdfeasible==DSDP_PDUNKNOWN) dsdp->pdfeasible=DSDP_PDFEASIBLE;
353 info=DSDPEventLogEnd(dsdp->solvetime);
354 DSDPFunctionReturn(0);
355}
356
357
358#undef __FUNCT__
359#define __FUNCT__ "DSDPCallMonitors"
367int DSDPCallMonitors(DSDP dsdp,DMonitor dmonitor[], int ndmonitors){
368 int i,info;
369 DSDPFunctionBegin;
370 for (i=0; i<ndmonitors;i++){
371 info=(dmonitor[i].monitor)(dsdp,dmonitor[i].monitorctx); DSDPCHKERR(info);
372 }
373 DSDPFunctionReturn(0);
374}
375/* ---------------------------------------------------------- */
376#undef __FUNCT__
377#define __FUNCT__ "DSDPCheckConvergence"
385 int info;
386 DSDPTruth unbounded;
387
388 DSDPFunctionBegin;
389 info = DSDPGetConicDimension(dsdp,&dsdp->np);DSDPCHKERR(info);
390 dsdp->rgap=(dsdp->ppobj-dsdp->ddobj)/(1.0+fabs(dsdp->ppobj)+fabs(dsdp->ddobj));
391 dsdp->pstepold=dsdp->pstep;
392 if (dsdp->reason==CONTINUE_ITERATING){
393 if (dsdp->itnow>0){
394 info=DSDPCheckForUnboundedObjective(dsdp,&unbounded);DSDPCHKERR(info);
395 if (unbounded==DSDP_TRUE){
396 dsdp->pdfeasible=DSDP_UNBOUNDED;
397 info=DSDPSetConvergenceFlag(dsdp,DSDP_CONVERGED); DSDPCHKERR(info);
398 }
399 }
400 if (dsdp->reason==CONTINUE_ITERATING){
401 if (dsdp->muold<dsdp->mutarget && dsdp->pstep==1 && dsdp->dstep==1 && dsdp->rgap<1e-5){
402 info=DSDPSetConvergenceFlag(dsdp,DSDP_NUMERICAL_ERROR); DSDPCHKERR(info);
403 DSDPLogInfo(0,2,"DSDP Finished: Numerical issues: Increase in Barrier function. \n");}
404 if (dsdp->itnow >= dsdp->maxiter){
405 info=DSDPSetConvergenceFlag(dsdp,DSDP_MAX_IT); DSDPCHKERR(info);}
406 if (dsdp->Mshift>dsdp->maxschurshift){
407 info = DSDPSetConvergenceFlag(dsdp,DSDP_INDEFINITE_SCHUR_MATRIX); DSDPCHKERR(info);
408 }
409 }
410 info=DSDPCallMonitors(dsdp,dsdp->dmonitor,dsdp->nmonitors);DSDPCHKERR(info);
411 info=DSDPMonitorCones(dsdp,0); DSDPCHKERR(info);
412 }
413 dsdp->muold=dsdp->mutarget;
414 info = DSDPStopReason(dsdp,reason); DSDPCHKERR(info);
415 DSDPFunctionReturn(0);
416}
417
418
419
420/* ---------------------------------------------------------- */
421#undef __FUNCT__
422#define __FUNCT__ "DSDPTakeDown"
429
430 int i,info;
431
432 DSDPFunctionBegin;
433 DSDPValid(dsdp);
434 info = DSDPVecDestroy(&dsdp->rhs);DSDPCHKERR(info);
435 info = DSDPVecDestroy(&dsdp->rhs1);DSDPCHKERR(info);
436 info = DSDPVecDestroy(&dsdp->rhs2);DSDPCHKERR(info);
437 info = DSDPVecDestroy(&dsdp->rhstemp);DSDPCHKERR(info);
438 info = DSDPVecDestroy(&dsdp->y);DSDPCHKERR(info);
439 info = DSDPVecDestroy(&dsdp->ytemp);DSDPCHKERR(info);
440 info = DSDPVecDestroy(&dsdp->dy1);DSDPCHKERR(info);
441 info = DSDPVecDestroy(&dsdp->dy2);DSDPCHKERR(info);
442 info = DSDPVecDestroy(&dsdp->dy);DSDPCHKERR(info);
443 for (i=0;i<MAX_XMAKERS;i++){
444 info = DSDPVecDestroy(&dsdp->xmaker[i].y);DSDPCHKERR(info);
445 info = DSDPVecDestroy(&dsdp->xmaker[i].dy);DSDPCHKERR(info);
446 info = DSDPVecDestroy(&dsdp->xmaker[i].rhs);DSDPCHKERR(info);
447 }
448 info = DSDPVecDestroy(&dsdp->xmakerrhs);DSDPCHKERR(info);
449 info = DSDPVecDestroy(&dsdp->y0);DSDPCHKERR(info);
450 info = DSDPVecDestroy(&dsdp->b);DSDPCHKERR(info);
451
452 info = DSDPCGDestroy(&dsdp->sles);DSDPCHKERR(info);
453 info = DSDPDestroyCones(dsdp);DSDPCHKERR(info);
454 info = DSDPSchurMatDestroy(&dsdp->M);DSDPCHKERR(info);
455 info = DSDPGetConicDimension(dsdp,&dsdp->np);DSDPCHKERR(info);
456 dsdp->setupcalled=DSDP_FALSE;
457 DSDPFunctionReturn(0);
458}
459
469int DSDPSetDestroyRoutine(DSDP dsdp, int (*fd)(void*), void* ctx){
470 int nd=dsdp->ndroutines;
471 if (nd<10){
472 dsdp->droutine[nd].f=fd;
473 dsdp->droutine[nd].ptr=ctx;
474 dsdp->ndroutines++;
475 } else {
476 printf("TOO MANY Destroy routines\n");
477 return 1;
478 }
479 return 0;
480}
481
482
494#undef __FUNCT__
495#define __FUNCT__ "DSDPDestroy"
496int DSDPDestroy(DSDP dsdp){
497 int i,info;
498 DSDPFunctionBegin;
499 DSDPValid(dsdp);
500 for (i=0;i<dsdp->ndroutines;i++){
501 info=(*dsdp->droutine[i].f)(dsdp->droutine[i].ptr);DSDPCHKERR(info);
502 }
503 info=DSDPTakeDown(dsdp);DSDPCHKERR(info);
504 DSDPFREE(&dsdp,&info);DSDPCHKERR(info);
505 DSDPFunctionReturn(0);
506}
int DSDPCreateLUBoundsCone(DSDP dsdp, LUBounds *dspcone)
Create bounds cone.
Definition allbounds.c:566
The API to DSDP for those applications using DSDP as a subroutine library.
Internal data structure for the DSDP solver.
int DSDPSetUpCones(DSDP)
Each cone should factor data or allocate internal data structures.
Definition dsdpcops.c:58
int DSDPComputeANorm2(DSDP, DSDPVec)
Compute norm of A and C.
Definition dsdpcops.c:246
int DSDPMonitorCones(DSDP, int)
This routine is called once per iteration.
Definition dsdpcops.c:450
int DSDPSolveDynamicRho(DSDP)
Apply dual-scaling algorithm.
Definition dualalg.c:121
int DSDPGetConicDimension(DSDP, double *)
Get the total dimension of the cones.
Definition dsdpcops.c:401
int DSDPDestroyCones(DSDP)
Each cone shoudl free its data structures.
Definition dsdpcops.c:107
int DSDPInitializeVariables(DSDP)
Initialize variables and factor S.
Definition dualalg.c:475
int DSDPSetUpCones2(DSDP, DSDPVec, DSDPSchurMat)
Each cone should allocate its data structures .
Definition dsdpcops.c:84
DSDPTerminationReason
There are many reasons to terminate the solver.
@ DSDP_INDEFINITE_SCHUR_MATRIX
@ CONTINUE_ITERATING
@ DSDP_MAX_IT
@ DSDP_CONVERGED
@ DSDP_NUMERICAL_ERROR
@ DSDP_UNBOUNDED
@ DSDP_PDFEASIBLE
@ DSDP_PDUNKNOWN
DSDPTruth
Boolean variables.
@ DSDP_FALSE
@ DSDP_TRUE
int DSDPSchurMatSetup(DSDPSchurMat M, DSDPVec Y)
Set up the data structure.
int DSDPSchurMatInitialize(DSDPSchurMat *M)
Initialize pointers to null.
int DSDPSchurMatDestroy(DSDPSchurMat *M)
Free the memory in the data structure.
int DSDPSetDefaultParameters(DSDP dsdp)
Set default parameters.
Definition dsdpsetup.c:122
int DSDPScaleData(DSDP dsdp)
Scale the matrix C.
Definition dsdpsetup.c:311
int DSDPTakeDown(DSDP dsdp)
Destroy internal data structures.
Definition dsdpsetup.c:428
int DSDPCheckConvergence(DSDP dsdp, DSDPTerminationReason *reason)
Check for convergence and monitor solution.
Definition dsdpsetup.c:384
int DSDPCallMonitors(DSDP dsdp, DMonitor dmonitor[], int ndmonitors)
Call the monitor routines.
Definition dsdpsetup.c:367
int DSDPSetDefaultMonitors(DSDP dsdp)
Set convergence monitor.
Definition dsdpsetup.c:165
int DSDPComputeDataNorms(DSDP dsdp)
Compute norms of A,C, and b.
Definition dsdpsetup.c:283
int DSDPSetDefaultStatistics(DSDP dsdp)
Set default statistics.
Definition dsdpsetup.c:84
int DSDPSetDestroyRoutine(DSDP dsdp, int(*fd)(void *), void *ctx)
Set a routine that will be called during DSDPDestroy().
Definition dsdpsetup.c:469
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
Schur complement matrix whose solution is the Newton direction.
Internal structures for the DSDP solver.
Definition dsdp.h:65