DSDP
dualalg.c
Go to the documentation of this file.
1 #include "dsdp.h"
2 #include "dsdpsys.h"
8 int DSDPChooseBarrierParameter(DSDP,double,double*,double*);
9 int DSDPYStepLineSearch(DSDP,double,double,DSDPVec);
10 int DSDPYStepLineSearch2(DSDP,double,double,DSDPVec);
11 int DSDPResetY0(DSDP);
12 
13 #undef __FUNCT__
14 #define __FUNCT__ "DSDPYStepLineSearch"
15 
24 int 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"
64 
73 int 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"
116 
122 
123  int info,attempt,maxattempts;
124  double dd1,dd2,mutarget,ppnorm;
125  DSDPTerminationReason reason;
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"
227 
240 int 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"
322 
328 int DSDPResetY0(DSDP dsdp){
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"
446 int 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"
469 
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"
540 int 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 }
DSDPTruth
Boolean variables.
int DSDPComputeLogSDeterminant(DSDP, double *)
Compute the logarithmic barrier function for the dual varialbe S.
Definition: dsdpcops.c:495
struct DSDPVec_C DSDPVec
This object hold m+2 variables: a scaling of C, the y variables, and r.
Definition: dsdpvec.h:25
int DSDPInvertS(DSDP)
Invert the S variables in each cone.
Definition: dsdpcops.c:307
int DSDPComputeSS(DSDP, DSDPVec, DSDPDualFactorMatrix, DSDPTruth *)
Compute the dual variables S in each cone.
Definition: dsdpcops.c:272
int DSDPSchurMatFactor(DSDPSchurMat M, DSDPTruth *successful)
Factor M.
Definition: dsdpschurmat.c:196
Error handling, printing, and profiling.
int DSDPComputePDY1(DSDP, double, DSDPVec)
Compute an affine step direction dy1.
Definition: dualimpl.c:105
int DSDPSolveDynamicRho(DSDP dsdp)
Apply dual-scaling algorithm.
Definition: dualalg.c:121
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 DSDPComputeNewY(DSDP, double, DSDPVec)
Update the Y variables.
Definition: dualimpl.c:125
int DSDPComputeAndFactorS(DSDP dsdp, DSDPTruth *psdefinite)
Compute and factor the dual matrix variables.
Definition: dualalg.c:540
Internal structures for the DSDP solver.
Definition: dsdp.h:65
DSDPTerminationReason
There are many reasons to terminate the solver.
int DSDPCheckConvergence(DSDP, DSDPTerminationReason *)
Check for convergence and monitor solution.
Definition: dsdpsetup.c:384
int DSDPSetRR(DSDP, double)
Set variable r.
Definition: dualimpl.c:345
The API to DSDP for those applications using DSDP as a subroutine library.
int DSDPComputePY(DSDP, double, DSDPVec)
Compute PY = Y - beta DY for use in computing X.
Definition: dualimpl.c:150
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 DSDPSetY(DSDP, double, double, DSDPVec)
Update the solver with these y variables.
Definition: dualimpl.c:309
int DSDPSchurMatView(DSDPSchurMat M)
Print the matrix.
Definition: dsdpschurmat.c:376
int DSDPComputeDY(DSDP, double, DSDPVec, double *)
Compute the step direction.
Definition: dualimpl.c:45
int DSDPComputeHessian(DSDP, DSDPSchurMat, DSDPVec, DSDPVec)
Compute the Schur complement, or Gram, matrix for each cone.
Definition: dsdpcops.c:142
int DSDPComputePotential(DSDP, DSDPVec, double, double *)
Compute the potential of the given point.
Definition: dualimpl.c:261
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 DSDPComputeObjective(DSDP, DSDPVec, double *)
Compute the objective function (DD).
Definition: dualimpl.c:21
int DSDPGetRR(DSDP, double *)
Get variable r.
Definition: dualimpl.c:361
int DSDPComputeG(DSDP, DSDPVec, DSDPVec, DSDPVec)
Compute the gradient of the barrier for each cone.
Definition: dsdpcops.c:215
int DSDPComputeMaxStepLength(DSDP, DSDPVec, DSDPDualFactorMatrix, double *)
Compute the maximum step length for the given step direction.
Definition: dsdpcops.c:336
int DSDPComputeDualStepDirections(DSDP dsdp)
Compute the step direction by computing a linear system and solving it.
Definition: dualalg.c:370
int DSDPSchurMatShiftDiagonal(DSDPSchurMat M, double dd)
Add a scalar to each diagonal element of the matrix.
Definition: dsdpschurmat.c:120
int DSDPComputePDY(DSDP, double, DSDPVec, double *)
Compute the step direction.
Definition: dualimpl.c:77
int DSDPSetConvergenceFlag(DSDP dsdp, DSDPTerminationReason reason)
Monitor each iteration of the solver.
Definition: dsdpsetdata.c:968
int DSDPComputePotential2(DSDP, DSDPVec, double, double, double *)
Compute the objective function plus the barrier function.
Definition: dualimpl.c:287
int DSDPCGSolve(DSDP, DSDPSchurMat, DSDPVec, DSDPVec, double, DSDPTruth *)
Apply CG to solve for the step directions.
Definition: dsdpcg.c:239
int DSDPInitializeVariables(DSDP dsdp)
Initialize variables and factor S.
Definition: dualalg.c:475
int DSDPChooseBarrierParameter(DSDP, double, double *, double *)
DSDP barrier heuristic choses the smalles value of mu such that X>0.
Definition: dualalg.c:240
int DSDPResetY0(DSDP)
After 1 iteration, consider increasing the variable r.
Definition: dualalg.c:328
int DSDPGetMaxYElement(DSDP, double *)
Copy the the infinity norm of the variables y.
Definition: dsdpsetdata.c:645