DSDP
dualalg.c
Go to the documentation of this file.
1#include "dsdp.h"
2#include "dsdpsys.h"
8int DSDPChooseBarrierParameter(DSDP,double,double*,double*);
9int DSDPYStepLineSearch(DSDP,double,double,DSDPVec);
10int DSDPYStepLineSearch2(DSDP,double,double,DSDPVec);
11int DSDPResetY0(DSDP);
12
13#undef __FUNCT__
14#define __FUNCT__ "DSDPYStepLineSearch"
24int DSDPYStepLineSearch(DSDP dsdp, double mutarget, double dstep0, DSDPVec dy){
25 /* The merit function is the dual potential function */
26 int info,attempt,maxattempts=30;
27 double dstep,newpotential,logdet;
28 double better=0.05, steptol=1e-8,maxmaxstep=0;
29 DSDPTruth psdefinite;
30 DSDPFunctionBegin;
31 info=DSDPComputeMaxStepLength(dsdp,dy,DUAL_FACTOR,&maxmaxstep);DSDPCHKERR(info);
32 info=DSDPComputePotential(dsdp,dsdp->y,dsdp->logdet,&dsdp->potential);DSDPCHKERR(info);
33 if (dsdp->pnorm<0.5) better=0.0;
34 dstep=DSDPMin(dstep0,0.95*maxmaxstep);
35 if (dstep * dsdp->pnorm > dsdp->maxtrustradius) dstep=dsdp->maxtrustradius/dsdp->pnorm;
36 DSDPLogInfo(0,8,"Full Dual StepLength %4.4e, %4.4e\n",maxmaxstep,dstep);
37 psdefinite=DSDP_FALSE;
38 for (psdefinite=DSDP_FALSE,attempt=0; attempt<maxattempts && psdefinite==DSDP_FALSE; attempt++){
39 info=DSDPComputeNewY(dsdp,dstep,dsdp->ytemp);DSDPCHKERR(info);
40 info=DSDPComputeSS(dsdp,dsdp->ytemp,DUAL_FACTOR,&psdefinite);DSDPCHKERR(info);
41 if (psdefinite==DSDP_TRUE){
42 info=DSDPComputeLogSDeterminant(dsdp,&logdet);DSDPCHKERR(info);
43 info=DSDPComputePotential(dsdp,dsdp->ytemp,logdet,&newpotential);DSDPCHKERR(info);
44 if (newpotential>dsdp->potential-better && dstep > 0.001/dsdp->pnorm ){
45 DSDPLogInfo(0,2,"Not sufficient reduction. Reduce stepsize. Trust Radius: %4.4e\n",dstep*dsdp->pnorm);
46 psdefinite=DSDP_FALSE; dstep=0.3*dstep;
47 }
48 } else {
49 dstep=dstep/3.0;
50 DSDPLogInfo(0,2,"Dual Matrix not Positive Definite: Reduce step %4.4e",dstep);
51 }
52 if (dstep*dsdp->pnorm < steptol && dstep < steptol) break;
53 } /* Hopefully, the initial step size works and only go through loop once */
54 if (psdefinite==DSDP_TRUE){
55 info=DSDPSetY(dsdp,dstep,logdet,dsdp->ytemp);DSDPCHKERR(info);
56 } else {
57 info=DSDPSetY(dsdp,0,dsdp->logdet,dsdp->y);DSDPCHKERR(info);
58 }
59 DSDPFunctionReturn(0);
60}
61
62#undef __FUNCT__
63#define __FUNCT__ "DSDPYStepLineSearch2"
73int DSDPYStepLineSearch2(DSDP dsdp, double mutarget, double dstep0, DSDPVec dy){
74 /* The merit function is the objective (DD) plus the barrier function */
75 /* This line search is used in the corrector steps */
76 int info, attempt, maxattempts=10;
77 double dstep,newpotential,bdotdy,oldpotential,logdet;
78 double maxmaxstep=0,steptol=1e-6;
79 double a,b;
80 DSDPTruth psdefinite;
81 DSDPFunctionBegin;
82 info=DSDPComputeMaxStepLength(dsdp,dy,DUAL_FACTOR,&maxmaxstep);DSDPCHKERR(info);
83 info=DSDPComputePotential2(dsdp,dsdp->y,mutarget, dsdp->logdet,&oldpotential);DSDPCHKERR(info);
84 info=DSDPVecDot(dsdp->rhs,dy,&bdotdy);DSDPCHKERR(info);
85 dstep=DSDPMin(dstep0,0.95*maxmaxstep);
86 if (dstep * dsdp->pnorm > dsdp->maxtrustradius) dstep=dsdp->maxtrustradius/dsdp->pnorm;
87 DSDPLogInfo(0,8,"Full Dual StepLength %4.4e, %4.4e\n",maxmaxstep,dstep);
88 for (psdefinite=DSDP_FALSE,attempt=0; attempt<maxattempts && psdefinite==DSDP_FALSE; attempt++){
89 if (dstep < steptol) break;
90 info=DSDPComputeNewY(dsdp,dstep,dsdp->ytemp);DSDPCHKERR(info);
91 info=DSDPComputeSS(dsdp,dsdp->ytemp,DUAL_FACTOR,&psdefinite);DSDPCHKERR(info);
92 if (psdefinite==DSDP_TRUE){
93 info=DSDPComputeLogSDeterminant(dsdp,&logdet);DSDPCHKERR(info);
94 info=DSDPComputePotential2(dsdp,dsdp->ytemp,mutarget,logdet,&newpotential);DSDPCHKERR(info);
95 b=bdotdy; a=2*(newpotential-oldpotential+bdotdy*dstep)/(dstep*dstep);
96 if (newpotential>oldpotential-0.1*dstep*bdotdy ){
97 DSDPLogInfo(0,2,"Not sufficient reduction. Reduce stepsize. Step:: %4.4e\n",dstep);
98 psdefinite=DSDP_FALSE;
99 if (b/a<dstep && b/a>0){ dstep=b/a;} else { dstep=dstep/2; }
100 }
101 } else {
102 dstep=dstep/2.0;
103 DSDPLogInfo(0,2,"Dual Matrix not Positive Definite: Reduce step %4.4e",dstep);
104 }
105 } /* Hopefully, the initial step size works and only go through loop once */
106 if (psdefinite==DSDP_TRUE && dstep>=steptol){
107 info=DSDPSetY(dsdp,dstep,logdet,dsdp->ytemp);DSDPCHKERR(info);
108 } else {
109 info=DSDPSetY(dsdp,0,dsdp->logdet,dsdp->y);DSDPCHKERR(info);
110 }
111 DSDPFunctionReturn(0);
112}
113
114#undef __FUNCT__
115#define __FUNCT__ "DSDPSolveDynmaicRho"
122
123 int info,attempt,maxattempts;
124 double dd1,dd2,mutarget,ppnorm;
126 DSDPTruth cg1;
127 DSDPTruth psdefinite;
128
129 DSDPFunctionBegin;
130
131 info=DSDPVecCopy(dsdp->y,dsdp->y0);DSDPCHKERR(info);
132 for (dsdp->itnow=0; dsdp->itnow <= dsdp->maxiter+1 ; dsdp->itnow++){
133
134 /* Check Convergence, and print information if desired */
135 info=DSDPCheckConvergence(dsdp,&reason);DSDPCHKERR(info);
136 if (reason != CONTINUE_ITERATING){break;}
137 if (dsdp->mu0>0){dsdp->mutarget=DSDPMin(dsdp->mutarget,dsdp->mu0);}
138
139 /* Compute the Gram matrix M and rhs */
140 info=DSDPComputeDualStepDirections(dsdp); DSDPCHKERR(info);
141 if (dsdp->reason==DSDP_INDEFINITE_SCHUR_MATRIX){continue;}
142
143 info=DSDPComputePDY(dsdp,dsdp->mutarget,dsdp->dy,&dsdp->pnorm); DSDPCHKERR(info);
144
145 DSDPEventLogBegin(dsdp->ptime);
146 info=DSDPComputePY(dsdp,1.0,dsdp->ytemp);DSDPCHKERR(info);
147 info=DSDPComputeSS(dsdp,dsdp->ytemp,PRIMAL_FACTOR,&psdefinite);DSDPCHKERR(info);
148 if (psdefinite==DSDP_TRUE){
149 dsdp->pstep=1.0;
150 info=DSDPSaveYForX(dsdp,dsdp->mutarget,dsdp->pstep);DSDPCHKERR(info);
151 } else {
152 dsdp->pstep=0.0;
153 }
154
155 if (dsdp->usefixedrho==DSDP_TRUE){
156 dsdp->rho=dsdp->rhon*dsdp->np;
157 mutarget=(dsdp->ppobj-dsdp->ddobj)/(dsdp->rho);
158 dsdp->pstep=0.5;
159 } else {
160 info = DSDPChooseBarrierParameter(dsdp,dsdp->mutarget,&dsdp->pstep,&mutarget);DSDPCHKERR(info);
161 dsdp->rho=(dsdp->ppobj-dsdp->ddobj)/(mutarget);
162 }
163 DSDPEventLogEnd(dsdp->ptime);
164
165 DSDPLogInfo(0,6,"Current mu=%4.8e, Target X with mu=%4.8e, Rho: %8.4e, Primal Step Length: %4.8f, pnorm: %4.8e\n",dsdp->mu,mutarget,dsdp->rho/dsdp->np,dsdp->pstep, dsdp->pnorm);
166
167
168 /* Take Dual Step */
169 /* Compute distance from chosen point on central path Pnorm */
170 DSDPEventLogBegin(dsdp->dtime);
171 info=DSDPComputeDY(dsdp,mutarget,dsdp->dy,&dsdp->pnorm); DSDPCHKERR(info);
172 if (dsdp->pnorm<0.1){ mutarget/=10; info=DSDPComputeDY(dsdp,mutarget,dsdp->dy,&dsdp->pnorm); DSDPCHKERR(info);}
173
174 info=DSDPYStepLineSearch(dsdp, mutarget, 1.0, dsdp->dy);DSDPCHKERR(info);
175 DSDPEventLogEnd(dsdp->dtime);
176
177 maxattempts=dsdp->reuseM;
178 if (dsdp->dstep<1 && dsdp->rgap<1e-5) maxattempts=0;
179 if (dsdp->dstep<1e-13) maxattempts=0;
180 if (dsdp->rgap<1e-6) maxattempts=0;
181 if (maxattempts>dsdp->reuseM) maxattempts=dsdp->reuseM;
182 for (attempt=0;attempt<maxattempts;attempt++){
183 double cgtol=1e-6;
184 if (attempt>0 && dsdp->pnorm < 0.1) break;
185 if (attempt > 0 && dsdp->dstep<1e-4) break;
186 if (dsdp->rflag) break;
187 DSDPEventLogBegin(dsdp->ctime);
188 DSDPLogInfo(0,2,"Reuse Matrix %d: Ddobj: %12.8e, Pnorm: %4.2f, Step: %4.2f\n",attempt,dsdp->ddobj,dsdp->pnorm,dsdp->dstep);
189 info=DSDPInvertS(dsdp);DSDPCHKERR(info);
190 info=DSDPComputeG(dsdp,dsdp->rhstemp,dsdp->rhs1,dsdp->rhs2);DSDPCHKERR(info);
191 if (dsdp->slestype==2 || dsdp->slestype==3){
192 if (dsdp->rflag){info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs1,dsdp->dy1,cgtol,&cg1);DSDPCHKERR(info);}
193 info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs2,dsdp->dy2,cgtol,&cg1);DSDPCHKERR(info);
194 }
195 info=DSDPVecDot(dsdp->b,dsdp->dy1,&dd1);DSDPCHKERR(info);
196 info=DSDPVecDot(dsdp->b,dsdp->dy2,&dd2);DSDPCHKERR(info);
197 if (dd1>0 && dd2>0){
198 mutarget=DSDPMin(mutarget,dd1/dd2*dsdp->schurmu);
199 }
200 mutarget=mutarget*(dsdp->np/(dsdp->np+sqrt(dsdp->np)));
201 info=DSDPComputeDY(dsdp,mutarget,dsdp->dy, &ppnorm);DSDPCHKERR(info);
202 if (ppnorm<=0){ DSDPEventLogEnd(dsdp->ctime); break; }
203 dsdp->pnorm=ppnorm;
204 info=DSDPYStepLineSearch2(dsdp, mutarget, dsdp->dstep, dsdp->dy);DSDPCHKERR(info);
205 DSDPEventLogEnd(dsdp->ctime);
206 }
207 if (attempt>0)dsdp->dstep=1.0;
208
209 dsdp->mutarget=DSDPMin(dsdp->mu,mutarget);
210
211 info=DSDPGetRR(dsdp,&dd1);DSDPCHKERR(info);
212 if (dsdp->itnow==0 && dsdp->xmaker[0].pstep<1.0 && dd1> 0 && dsdp->pstep<1.0 && dsdp->goty0==DSDP_FALSE){
213 info=DSDPResetY0(dsdp);DSDPCHKERR(info); continue;
214 dsdp->goty0=DSDP_FALSE;
215 }
216
217 } /* End of Dual Scaling Algorithm */
218
219
220 DSDPFunctionReturn(0);
221}
222
223
224
225#undef __FUNCT__
226#define __FUNCT__ "DSDPChooseBarrierParameter"
240int DSDPChooseBarrierParameter(DSDP dsdp, double mutarget, double *ppstep, double *nextmutarget){
241 int attempt,info,count=0;
242 double pnorm,pstep=*ppstep,pmumin,mur, dmury1, mutarget2;
243 DSDPTruth psdefinite=DSDP_FALSE;
244
245 DSDPFunctionBegin;
246
247 *nextmutarget=mutarget;
248 /* Compute a feasible primal matrix with the smallest possible value of \mu */
249 /* Start by finding a psd feasible matrix */
250 if (*ppstep >=1){
251 pstep=1.0;
252 /* PS should alread be formed and factored */
253 } else {
254 mur=-1.0/mutarget;
255 info=DSDPComputePDY(dsdp,mutarget,dsdp->dy,&pnorm); DSDPCHKERR(info);
256 info=DSDPComputeMaxStepLength(dsdp,dsdp->dy,DUAL_FACTOR,&pstep);DSDPCHKERR(info);
257 /* pstep=DSDPMin(0.97*pstep,1.0); */
258 if (pstep<1.0) {pstep=DSDPMin(0.97*pstep,1.0);} else {pstep=DSDPMin(1.0*pstep,1.0);}
259 while (psdefinite==DSDP_FALSE){
260 if (count > 2 && pstep <1e-8){pstep=0;break;}
261 info=DSDPComputePY(dsdp,pstep,dsdp->ytemp);DSDPCHKERR(info);
262 info=DSDPComputeSS(dsdp,dsdp->ytemp,PRIMAL_FACTOR,&psdefinite);DSDPCHKERR(info);
263 if (psdefinite==DSDP_FALSE){
264 if (count>1) pstep=0.5*pstep; else pstep=0.97*pstep;
265 DSDPLogInfo(0,2,"Reducing pstep: %8.8e\n",pstep);
266 count++;
267 }
268 }
269 *ppstep=pstep;
270 if (pstep > dsdp->xmaker[0].pstep || mutarget < dsdp->xmaker[0].mu * 1.0e-8){
271 info=DSDPSaveYForX(dsdp,mutarget,pstep);DSDPCHKERR(info);
272 }
273 if (pstep==0){
274 DSDPFunctionReturn(0);
275 }
276 }
277
278 /* Now determine how much smaller of mu can be used */
279 mur=pstep/mutarget;
280 info=DSDPComputePDY1(dsdp,mur,dsdp->rhstemp); DSDPCHKERR(info);
281
282 /* Smallest value of mu that gives a positive definite matrix */
283 info=DSDPComputeMaxStepLength(dsdp,dsdp->rhstemp,PRIMAL_FACTOR,&dmury1);DSDPCHKERR(info);
284 dmury1 = DSDPMin(1000,0.97*dmury1);
285
286 /* We could test the point, My tests say its not necessary -- its good! */
287 attempt=0;psdefinite=DSDP_FALSE;
288 pmumin=mutarget / (1.0 + 1.0 * dmury1); /* This should be positive definite */
289 while ( 0 && psdefinite==DSDP_FALSE){
290 pmumin=mutarget / (1.0 + 1.0 * dmury1); /* This should be positive definite */
291 if (attempt>2){pmumin=mutarget;}/* We have actually factored this one. It is PSD. */
292 info=DSDPComputePDY(dsdp,pmumin,dsdp->dy,&pnorm); DSDPCHKERR(info);
293 info=DSDPComputePY(dsdp,pstep,dsdp->ytemp); DSDPCHKERR(info);
294 info=DSDPComputeSS(dsdp,dsdp->ytemp,PRIMAL_FACTOR,&psdefinite);DSDPCHKERR(info);
295 if (psdefinite==DSDP_FALSE){ dmury1*=0.9; /* printf("NO GOOD \n"); */ }
296 else { /* printf("ITS GOOD \n"); */}
297 attempt++;
298 DSDPLogInfo(0,2,"GOT X: Smallest Mu for feasible X: %4.4e.\n",pmumin);
299 }
300
301 DSDPLogInfo(0,6,"GOT X: Smallest Mu for feasible X: %4.4e \n",pmumin);
302
303 mutarget2=mutarget;
304 if (dsdp->pstep==1){
305 mutarget2 = pmumin;
306 } else {
307 /* printf("PMUMIN: %4.4e MUTARGET: %4.4e \n",pmumin,dsdp->mutarget); */
308 mutarget2 = mutarget;
309 mutarget2 = (pstep)*mutarget + (1.0-pstep)*dsdp->mu;
310 mutarget2 = (pstep)*pmumin + (1.0-pstep)*dsdp->mu;
311 }
312
313 mutarget2=DSDPMax(dsdp->mu/dsdp->rhon,mutarget2);
314 if (dsdp->mu0>0){mutarget2=DSDPMin(mutarget2,dsdp->mu0);}
315
316 *nextmutarget=mutarget2;
317 DSDPFunctionReturn(0);
318}
319
320#undef __FUNCT__
321#define __FUNCT__ "DSDPResetY0"
329 int info;
330 double r,rr,dd;
331 DSDPTruth psdefinite;
332 DSDPFunctionBegin;
333 info=DSDPComputeDY(dsdp,dsdp->mutarget,dsdp->dy,&dsdp->pnorm); DSDPCHKERR(info);
334 info=DSDPVecCopy(dsdp->y0,dsdp->y);DSDPCHKERR(info);
335 info=DSDPGetRR(dsdp,&r);DSDPCHKERR(info);
336 rr=DSDPMax(r*10000,1e12);
337 info=DSDPSetRR(dsdp,rr);DSDPCHKERR(info);
338 info=DSDPComputeSS(dsdp,dsdp->y,DUAL_FACTOR,&psdefinite);DSDPCHKERR(info);
339 info=DSDPComputeLogSDeterminant(dsdp,&dsdp->logdet);DSDPCHKERR(info);
340 info=DSDPSetY(dsdp,1.0,dsdp->logdet,dsdp->y);DSDPCHKERR(info);
341 info=DSDPVecGetR(dsdp->b,&dd);DSDPCHKERR(info);
342
343 dsdp->mutarget=fabs(rr*dd);
344 dsdp->mu=fabs(rr*dd);
345
346 /* dsdp->par.mu0=mutarget; */
347 dsdp->goty0=DSDP_TRUE;
348 DSDPLogInfo(0,2,"Restart algorithm\n");
349 DSDPFunctionReturn(0);
350}
351
352
368#undef __FUNCT__
369#define __FUNCT__ "DSDPComputeDualStepDirections"
371 int info,computem=1;
372 double madd,ymax,cgtol=1e-7;
373 DSDPTruth cg1,cg2,psdefinite;
374 DSDPFunctionBegin;
375
376 if (dsdp->itnow>30) dsdp->slestype=3;
377 if (dsdp->rgap<1e-3) dsdp->slestype=3;
378 if (dsdp->m<40) dsdp->slestype=3;
379 if (0 && dsdp->itnow>20 && dsdp->m<500) dsdp->slestype=3;
380 info=DSDPGetMaxYElement(dsdp,&ymax);DSDPCHKERR(info);
381 if (dsdp->slestype==1){
382 cg1=DSDP_TRUE; cg2=DSDP_TRUE;
383 info=DSDPInvertS(dsdp);DSDPCHKERR(info);
384 info=DSDPComputeG(dsdp,dsdp->rhstemp,dsdp->rhs1,dsdp->rhs2);DSDPCHKERR(info);
385 info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs1,dsdp->dy1,cgtol,&cg1);DSDPCHKERR(info);
386 if (cg1==DSDP_TRUE){info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs2,dsdp->dy2,cgtol,&cg2);DSDPCHKERR(info);}
387 if (cg1==DSDP_FALSE || cg2==DSDP_FALSE) dsdp->slestype=2;
388 }
389 if (dsdp->slestype==2){
390 cg1=DSDP_TRUE; cg2=DSDP_TRUE;
391 DSDPLogInfo(0,9,"Compute Hessian\n");
392 info=DSDPInvertS(dsdp);DSDPCHKERR(info);
393 info=DSDPComputeHessian(dsdp,dsdp->M,dsdp->rhs1,dsdp->rhs2);DSDPCHKERR(info);
394 computem=0;
395 DSDPLogInfo(0,9,"Apply CG\n");
396 info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs1,dsdp->dy1,cgtol,&cg1);DSDPCHKERR(info);
397 if (cg1==DSDP_TRUE){info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs2,dsdp->dy2,cgtol,&cg2);DSDPCHKERR(info);}
398 if (cg1==DSDP_FALSE || cg2==DSDP_FALSE) dsdp->slestype=3;
399
400 }
401 if (dsdp->slestype==3){
402 DSDPLogInfo(0,9,"Factor Hessian\n");
403 psdefinite=DSDP_FALSE;
404 if (dsdp->Mshift < 1e-12 || dsdp->rgap<0.1 || dsdp->Mshift > 1e-6){
405 madd=dsdp->Mshift;
406 } else {
407 madd=1e-13;
408 }
409 if (computem){
410 info=DSDPInvertS(dsdp);DSDPCHKERR(info);
411 }
412 while (psdefinite==DSDP_FALSE){
413 if (0==1 && dsdp->Mshift>dsdp->maxschurshift){
414 info = DSDPSetConvergenceFlag(dsdp,DSDP_INDEFINITE_SCHUR_MATRIX); DSDPCHKERR(info);
415 break;
416 }
417 if (0 && dsdp->Mshift*ymax>dsdp->pinfeastol/10){
418 info = DSDPSetConvergenceFlag(dsdp,DSDP_INDEFINITE_SCHUR_MATRIX); DSDPCHKERR(info);
419 break;
420 }
421 if (madd*ymax>dsdp->pinfeastol*1000){
422 info = DSDPSetConvergenceFlag(dsdp,DSDP_INDEFINITE_SCHUR_MATRIX); DSDPCHKERR(info);
423 break;
424 }
425 if (computem){
426 info=DSDPComputeHessian(dsdp,dsdp->M,dsdp->rhs1,dsdp->rhs2);DSDPCHKERR(info);}
427 if (0==1){info=DSDPSchurMatView(dsdp->M);DSDPCHKERR(info);}
428 info = DSDPSchurMatShiftDiagonal(dsdp->M,madd);DSDPCHKERR(info);
429 info = DSDPSchurMatFactor(dsdp->M,&psdefinite); DSDPCHKERR(info);
430 computem=1;
431 if (psdefinite==DSDP_FALSE){
432 madd=madd*4 + 1.0e-13;
433 }
434 }
435 dsdp->Mshift=madd;
436 if (psdefinite==DSDP_TRUE){
437 info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs1,dsdp->dy1,cgtol,&cg1);DSDPCHKERR(info);
438 info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs2,dsdp->dy2,cgtol,&cg2);DSDPCHKERR(info);
439 }
440 }
441
442 DSDPFunctionReturn(0);
443}
444#undef __FUNCT__
445#define __FUNCT__ "DSDPComputeDualStepDirections"
446int DSDPRefineStepDirection(DSDP dsdp,DSDPVec rhs, DSDPVec dy){
447 int info;
448 double cgtol=1e-20;
449 DSDPTruth cg1;
450 DSDPFunctionBegin;
451
452 if (dsdp->reason==DSDP_INDEFINITE_SCHUR_MATRIX){
453 } else if (dsdp->reason==DSDP_SMALL_STEPS){
454 } else if (dsdp->pstep<1){
455 } else {
456 dsdp->slestype=4;
457 info=DSDPCGSolve(dsdp,dsdp->M,rhs,dy,cgtol,&cg1);DSDPCHKERR(info);
458 dsdp->slestype=3;
459 }
460
461 DSDPFunctionReturn(0);
462}
463
464#include "dsdp5.h"
465
466
467#undef __FUNCT__
468#define __FUNCT__ "DSDPInitializeVariables"
476
477 int info;
478 double r0,mutarget=dsdp->mutarget,penalty;
479 DSDPTruth psdefinite=DSDP_FALSE;
480
481 DSDPFunctionBegin;
482 info=DSDPGetRR(dsdp,&r0);DSDPCHKERR(info);
483 dsdp->rho=dsdp->np*dsdp->rhon;
484 if (r0>=0) { /* Use the specified starting point */
485 info=DSDPComputeObjective(dsdp,dsdp->y,&dsdp->ddobj);DSDPCHKERR(info);
486 info=DSDPComputeSS(dsdp,dsdp->y,DUAL_FACTOR,&psdefinite);DSDPCHKERR(info);
487 if (mutarget<0) mutarget=(dsdp->ppobj-dsdp->ddobj)/(dsdp->rho);
488 } else {
489 info=DSDPGetPenalty(dsdp,&penalty);DSDPCHKERR(info);
490 r0=0.1/(1.0+dsdp->cnorm);
491 while (psdefinite==DSDP_FALSE ){
492 r0=r0*100;
493 DSDPLogInfo(0,9,"Set Initial R0 %4.2e\n",r0);
494 info=DSDPSetRR(dsdp,r0);DSDPCHKERR(info);
495 info=DSDPComputeSS(dsdp,dsdp->y,DUAL_FACTOR,&psdefinite);DSDPCHKERR(info);
496 }
497 r0=r0*dsdp->np;
498 if (dsdp->cnorm>0 && dsdp->anorm>0 && dsdp->cnorm/dsdp->anorm<1){ r0=r0/(dsdp->cnorm/dsdp->anorm);}
499 dsdp->mu=r0*penalty;
500 if (mutarget<0){
501 mutarget=(dsdp->ppobj-dsdp->ddobj)/(dsdp->rho);
502 mutarget=(dsdp->ppobj-dsdp->ddobj)/(dsdp->np);
503 mutarget=r0*penalty;
504 }
505 DSDPLogInfo(0,9,"Set Initial R0 %4.2e\n",r0);
506 info=DSDPSetRR(dsdp,r0);DSDPCHKERR(info);
507 info=DSDPComputeObjective(dsdp,dsdp->y,&dsdp->ddobj);DSDPCHKERR(info);
508 info=DSDPComputeSS(dsdp,dsdp->y,DUAL_FACTOR,&psdefinite);DSDPCHKERR(info);
509 }
510 info=DSDPComputeObjective(dsdp,dsdp->y,&dsdp->ddobj);DSDPCHKERR(info);
511 if (psdefinite==DSDP_FALSE){
512 info=DSDPSetConvergenceFlag(dsdp,DSDP_INFEASIBLE_START);DSDPCHKERR(info);
513 } else {
514 info=DSDPComputeLogSDeterminant(dsdp,&dsdp->logdet);DSDPCHKERR(info);
515 info=DSDPComputePotential(dsdp,dsdp->y,dsdp->logdet,&dsdp->potential);DSDPCHKERR(info);
516 }
517 /* Tough guess, as all these rules suggest */
518 info=DSDPSetY(dsdp,0,dsdp->logdet,dsdp->y);DSDPCHKERR(info);
519 info=DSDPSaveYForX(dsdp,dsdp->xmaker[0].mu,0);DSDPCHKERR(info);
520 dsdp->mutarget=mutarget;
521 dsdp->pstep=0.0;
522 dsdp->pnorm=0;
523 /* dsdp->par.mu0=mutarget; */
524 DSDPFunctionReturn(0);
525}
526
538#undef __FUNCT__
539#define __FUNCT__ "DSDPComputeAndFactorS"
540int DSDPComputeAndFactorS(DSDP dsdp,DSDPTruth *psdefinite){
541 int info;
542 DSDPFunctionBegin;
543 info=DSDPComputeSS(dsdp,dsdp->y,DUAL_FACTOR,psdefinite);DSDPCHKERR(info);
544 DSDPFunctionReturn(0);
545}
The API to DSDP for those applications using DSDP as a subroutine library.
Internal data structure for the DSDP solver.
int DSDPComputePotential2(DSDP, DSDPVec, double, double, double *)
Compute the objective function plus the barrier function.
Definition dualimpl.c:287
int DSDPComputeSS(DSDP, DSDPVec, DSDPDualFactorMatrix, DSDPTruth *)
Compute the dual variables S in each cone.
Definition dsdpcops.c:272
int DSDPSetRR(DSDP, double)
Set variable r.
Definition dualimpl.c:345
int DSDPComputeG(DSDP, DSDPVec, DSDPVec, DSDPVec)
Compute the gradient of the barrier for each cone.
Definition dsdpcops.c:215
int DSDPComputeHessian(DSDP, DSDPSchurMat, DSDPVec, DSDPVec)
Compute the Schur complement, or Gram, matrix for each cone.
Definition dsdpcops.c:142
int DSDPComputeLogSDeterminant(DSDP, double *)
Compute the logarithmic barrier function for the dual varialbe S.
Definition dsdpcops.c:495
int DSDPInvertS(DSDP)
Invert the S variables in each cone.
Definition dsdpcops.c:307
int DSDPCGSolve(DSDP, DSDPSchurMat, DSDPVec, DSDPVec, double, DSDPTruth *)
Apply CG to solve for the step directions.
Definition dsdpcg.c:239
int DSDPComputeNewY(DSDP, double, DSDPVec)
Update the Y variables.
Definition dualimpl.c:125
int DSDPComputeObjective(DSDP, DSDPVec, double *)
Compute the objective function (DD).
Definition dualimpl.c:21
int DSDPComputePDY(DSDP, double, DSDPVec, double *)
Compute the step direction.
Definition dualimpl.c:77
int DSDPCheckConvergence(DSDP, DSDPTerminationReason *)
Check for convergence and monitor solution.
Definition dsdpsetup.c:384
int DSDPComputeDY(DSDP, double, DSDPVec, double *)
Compute the step direction.
Definition dualimpl.c:45
int DSDPComputePDY1(DSDP, double, DSDPVec)
Compute an affine step direction dy1.
Definition dualimpl.c:105
int DSDPComputePY(DSDP, double, DSDPVec)
Compute PY = Y - beta DY for use in computing X.
Definition dualimpl.c:150
int DSDPSaveYForX(DSDP, double, double)
Save the current solution for later computation of X.
Definition dsdpx.c:149
int DSDPComputeMaxStepLength(DSDP, DSDPVec, DSDPDualFactorMatrix, double *)
Compute the maximum step length for the given step direction.
Definition dsdpcops.c:336
int DSDPComputePotential(DSDP, DSDPVec, double, double *)
Compute the potential of the given point.
Definition dualimpl.c:261
int DSDPGetRR(DSDP, double *)
Get variable r.
Definition dualimpl.c:361
int DSDPSetY(DSDP, double, double, DSDPVec)
Update the solver with these y variables.
Definition dualimpl.c:309
DSDPTerminationReason
There are many reasons to terminate the solver.
@ DSDP_INDEFINITE_SCHUR_MATRIX
@ CONTINUE_ITERATING
@ DSDP_INFEASIBLE_START
@ DSDP_SMALL_STEPS
@ PRIMAL_FACTOR
@ DUAL_FACTOR
DSDPTruth
Boolean variables.
@ DSDP_FALSE
@ DSDP_TRUE
int DSDPSchurMatFactor(DSDPSchurMat M, DSDPTruth *successful)
Factor M.
int DSDPSchurMatView(DSDPSchurMat M)
Print the matrix.
int DSDPSchurMatShiftDiagonal(DSDPSchurMat M, double dd)
Add a scalar to each diagonal element of the matrix.
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 DSDPInitializeVariables(DSDP dsdp)
Initialize variables and factor S.
Definition dualalg.c:475
int DSDPComputeDualStepDirections(DSDP dsdp)
Compute the step direction by computing a linear system and solving it.
Definition dualalg.c:370
int DSDPYStepLineSearch(DSDP, double, double, DSDPVec)
Used for Newton step, the merit function of this line search is the dual potential function.
Definition dualalg.c:24
int DSDPYStepLineSearch2(DSDP, double, double, DSDPVec)
Used for centering steps, the merit function of this line search is the objective function plus the b...
Definition dualalg.c:73
int DSDPResetY0(DSDP)
After 1 iteration, consider increasing the variable r.
Definition dualalg.c:328
int DSDPSolveDynamicRho(DSDP dsdp)
Apply dual-scaling algorithm.
Definition dualalg.c:121
int DSDPChooseBarrierParameter(DSDP, double, double *, double *)
DSDP barrier heuristic choses the smalles value of mu such that X>0.
Definition dualalg.c:240
Internal structures for the DSDP solver.
Definition dsdp.h:65