DSDP
|
00001 #include "dsdp.h" 00002 #include "dsdpsys.h" 00008 int DSDPChooseBarrierParameter(DSDP,double,double*,double*); 00009 int DSDPYStepLineSearch(DSDP,double,double,DSDPVec); 00010 int DSDPYStepLineSearch2(DSDP,double,double,DSDPVec); 00011 int DSDPResetY0(DSDP); 00012 00013 #undef __FUNCT__ 00014 #define __FUNCT__ "DSDPYStepLineSearch" 00015 00024 int DSDPYStepLineSearch(DSDP dsdp, double mutarget, double dstep0, DSDPVec dy){ 00025 /* The merit function is the dual potential function */ 00026 int info,attempt,maxattempts=30; 00027 double dstep,newpotential,logdet; 00028 double better=0.05, steptol=1e-8,maxmaxstep=0; 00029 DSDPTruth psdefinite; 00030 DSDPFunctionBegin; 00031 info=DSDPComputeMaxStepLength(dsdp,dy,DUAL_FACTOR,&maxmaxstep);DSDPCHKERR(info); 00032 info=DSDPComputePotential(dsdp,dsdp->y,dsdp->logdet,&dsdp->potential);DSDPCHKERR(info); 00033 if (dsdp->pnorm<0.5) better=0.0; 00034 dstep=DSDPMin(dstep0,0.95*maxmaxstep); 00035 if (dstep * dsdp->pnorm > dsdp->maxtrustradius) dstep=dsdp->maxtrustradius/dsdp->pnorm; 00036 DSDPLogInfo(0,8,"Full Dual StepLength %4.4e, %4.4e\n",maxmaxstep,dstep); 00037 psdefinite=DSDP_FALSE; 00038 for (psdefinite=DSDP_FALSE,attempt=0; attempt<maxattempts && psdefinite==DSDP_FALSE; attempt++){ 00039 info=DSDPComputeNewY(dsdp,dstep,dsdp->ytemp);DSDPCHKERR(info); 00040 info=DSDPComputeSS(dsdp,dsdp->ytemp,DUAL_FACTOR,&psdefinite);DSDPCHKERR(info); 00041 if (psdefinite==DSDP_TRUE){ 00042 info=DSDPComputeLogSDeterminant(dsdp,&logdet);DSDPCHKERR(info); 00043 info=DSDPComputePotential(dsdp,dsdp->ytemp,logdet,&newpotential);DSDPCHKERR(info); 00044 if (newpotential>dsdp->potential-better && dstep > 0.001/dsdp->pnorm ){ 00045 DSDPLogInfo(0,2,"Not sufficient reduction. Reduce stepsize. Trust Radius: %4.4e\n",dstep*dsdp->pnorm); 00046 psdefinite=DSDP_FALSE; dstep=0.3*dstep; 00047 } 00048 } else { 00049 dstep=dstep/3.0; 00050 DSDPLogInfo(0,2,"Dual Matrix not Positive Definite: Reduce step %4.4e",dstep); 00051 } 00052 if (dstep*dsdp->pnorm < steptol && dstep < steptol) break; 00053 } /* Hopefully, the initial step size works and only go through loop once */ 00054 if (psdefinite==DSDP_TRUE){ 00055 info=DSDPSetY(dsdp,dstep,logdet,dsdp->ytemp);DSDPCHKERR(info); 00056 } else { 00057 info=DSDPSetY(dsdp,0,dsdp->logdet,dsdp->y);DSDPCHKERR(info); 00058 } 00059 DSDPFunctionReturn(0); 00060 } 00061 00062 #undef __FUNCT__ 00063 #define __FUNCT__ "DSDPYStepLineSearch2" 00064 00073 int DSDPYStepLineSearch2(DSDP dsdp, double mutarget, double dstep0, DSDPVec dy){ 00074 /* The merit function is the objective (DD) plus the barrier function */ 00075 /* This line search is used in the corrector steps */ 00076 int info, attempt, maxattempts=10; 00077 double dstep,newpotential,bdotdy,oldpotential,logdet; 00078 double maxmaxstep=0,steptol=1e-6; 00079 double a,b; 00080 DSDPTruth psdefinite; 00081 DSDPFunctionBegin; 00082 info=DSDPComputeMaxStepLength(dsdp,dy,DUAL_FACTOR,&maxmaxstep);DSDPCHKERR(info); 00083 info=DSDPComputePotential2(dsdp,dsdp->y,mutarget, dsdp->logdet,&oldpotential);DSDPCHKERR(info); 00084 info=DSDPVecDot(dsdp->rhs,dy,&bdotdy);DSDPCHKERR(info); 00085 dstep=DSDPMin(dstep0,0.95*maxmaxstep); 00086 if (dstep * dsdp->pnorm > dsdp->maxtrustradius) dstep=dsdp->maxtrustradius/dsdp->pnorm; 00087 DSDPLogInfo(0,8,"Full Dual StepLength %4.4e, %4.4e\n",maxmaxstep,dstep); 00088 for (psdefinite=DSDP_FALSE,attempt=0; attempt<maxattempts && psdefinite==DSDP_FALSE; attempt++){ 00089 if (dstep < steptol) break; 00090 info=DSDPComputeNewY(dsdp,dstep,dsdp->ytemp);DSDPCHKERR(info); 00091 info=DSDPComputeSS(dsdp,dsdp->ytemp,DUAL_FACTOR,&psdefinite);DSDPCHKERR(info); 00092 if (psdefinite==DSDP_TRUE){ 00093 info=DSDPComputeLogSDeterminant(dsdp,&logdet);DSDPCHKERR(info); 00094 info=DSDPComputePotential2(dsdp,dsdp->ytemp,mutarget,logdet,&newpotential);DSDPCHKERR(info); 00095 b=bdotdy; a=2*(newpotential-oldpotential+bdotdy*dstep)/(dstep*dstep); 00096 if (newpotential>oldpotential-0.1*dstep*bdotdy ){ 00097 DSDPLogInfo(0,2,"Not sufficient reduction. Reduce stepsize. Step:: %4.4e\n",dstep); 00098 psdefinite=DSDP_FALSE; 00099 if (b/a<dstep && b/a>0){ dstep=b/a;} else { dstep=dstep/2; } 00100 } 00101 } else { 00102 dstep=dstep/2.0; 00103 DSDPLogInfo(0,2,"Dual Matrix not Positive Definite: Reduce step %4.4e",dstep); 00104 } 00105 } /* Hopefully, the initial step size works and only go through loop once */ 00106 if (psdefinite==DSDP_TRUE && dstep>=steptol){ 00107 info=DSDPSetY(dsdp,dstep,logdet,dsdp->ytemp);DSDPCHKERR(info); 00108 } else { 00109 info=DSDPSetY(dsdp,0,dsdp->logdet,dsdp->y);DSDPCHKERR(info); 00110 } 00111 DSDPFunctionReturn(0); 00112 } 00113 00114 #undef __FUNCT__ 00115 #define __FUNCT__ "DSDPSolveDynmaicRho" 00116 00121 int DSDPSolveDynamicRho(DSDP dsdp){ 00122 00123 int info,attempt,maxattempts; 00124 double dd1,dd2,mutarget,ppnorm; 00125 DSDPTerminationReason reason; 00126 DSDPTruth cg1; 00127 DSDPTruth psdefinite; 00128 00129 DSDPFunctionBegin; 00130 00131 info=DSDPVecCopy(dsdp->y,dsdp->y0);DSDPCHKERR(info); 00132 for (dsdp->itnow=0; dsdp->itnow <= dsdp->maxiter+1 ; dsdp->itnow++){ 00133 00134 /* Check Convergence, and print information if desired */ 00135 info=DSDPCheckConvergence(dsdp,&reason);DSDPCHKERR(info); 00136 if (reason != CONTINUE_ITERATING){break;} 00137 if (dsdp->mu0>0){dsdp->mutarget=DSDPMin(dsdp->mutarget,dsdp->mu0);} 00138 00139 /* Compute the Gram matrix M and rhs */ 00140 info=DSDPComputeDualStepDirections(dsdp); DSDPCHKERR(info); 00141 if (dsdp->reason==DSDP_INDEFINITE_SCHUR_MATRIX){continue;} 00142 00143 info=DSDPComputePDY(dsdp,dsdp->mutarget,dsdp->dy,&dsdp->pnorm); DSDPCHKERR(info); 00144 00145 DSDPEventLogBegin(dsdp->ptime); 00146 info=DSDPComputePY(dsdp,1.0,dsdp->ytemp);DSDPCHKERR(info); 00147 info=DSDPComputeSS(dsdp,dsdp->ytemp,PRIMAL_FACTOR,&psdefinite);DSDPCHKERR(info); 00148 if (psdefinite==DSDP_TRUE){ 00149 dsdp->pstep=1.0; 00150 info=DSDPSaveYForX(dsdp,dsdp->mutarget,dsdp->pstep);DSDPCHKERR(info); 00151 } else { 00152 dsdp->pstep=0.0; 00153 } 00154 00155 if (dsdp->usefixedrho==DSDP_TRUE){ 00156 dsdp->rho=dsdp->rhon*dsdp->np; 00157 mutarget=(dsdp->ppobj-dsdp->ddobj)/(dsdp->rho); 00158 dsdp->pstep=0.5; 00159 } else { 00160 info = DSDPChooseBarrierParameter(dsdp,dsdp->mutarget,&dsdp->pstep,&mutarget);DSDPCHKERR(info); 00161 dsdp->rho=(dsdp->ppobj-dsdp->ddobj)/(mutarget); 00162 } 00163 DSDPEventLogEnd(dsdp->ptime); 00164 00165 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); 00166 00167 00168 /* Take Dual Step */ 00169 /* Compute distance from chosen point on central path Pnorm */ 00170 DSDPEventLogBegin(dsdp->dtime); 00171 info=DSDPComputeDY(dsdp,mutarget,dsdp->dy,&dsdp->pnorm); DSDPCHKERR(info); 00172 if (dsdp->pnorm<0.1){ mutarget/=10; info=DSDPComputeDY(dsdp,mutarget,dsdp->dy,&dsdp->pnorm); DSDPCHKERR(info);} 00173 00174 info=DSDPYStepLineSearch(dsdp, mutarget, 1.0, dsdp->dy);DSDPCHKERR(info); 00175 DSDPEventLogEnd(dsdp->dtime); 00176 00177 maxattempts=dsdp->reuseM; 00178 if (dsdp->dstep<1 && dsdp->rgap<1e-5) maxattempts=0; 00179 if (dsdp->dstep<1e-13) maxattempts=0; 00180 if (dsdp->rgap<1e-6) maxattempts=0; 00181 if (maxattempts>dsdp->reuseM) maxattempts=dsdp->reuseM; 00182 for (attempt=0;attempt<maxattempts;attempt++){ 00183 double cgtol=1e-6; 00184 if (attempt>0 && dsdp->pnorm < 0.1) break; 00185 if (attempt > 0 && dsdp->dstep<1e-4) break; 00186 if (dsdp->rflag) break; 00187 DSDPEventLogBegin(dsdp->ctime); 00188 DSDPLogInfo(0,2,"Reuse Matrix %d: Ddobj: %12.8e, Pnorm: %4.2f, Step: %4.2f\n",attempt,dsdp->ddobj,dsdp->pnorm,dsdp->dstep); 00189 info=DSDPInvertS(dsdp);DSDPCHKERR(info); 00190 info=DSDPComputeG(dsdp,dsdp->rhstemp,dsdp->rhs1,dsdp->rhs2);DSDPCHKERR(info); 00191 if (dsdp->slestype==2 || dsdp->slestype==3){ 00192 if (dsdp->rflag){info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs1,dsdp->dy1,cgtol,&cg1);DSDPCHKERR(info);} 00193 info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs2,dsdp->dy2,cgtol,&cg1);DSDPCHKERR(info); 00194 } 00195 info=DSDPVecDot(dsdp->b,dsdp->dy1,&dd1);DSDPCHKERR(info); 00196 info=DSDPVecDot(dsdp->b,dsdp->dy2,&dd2);DSDPCHKERR(info); 00197 if (dd1>0 && dd2>0){ 00198 mutarget=DSDPMin(mutarget,dd1/dd2*dsdp->schurmu); 00199 } 00200 mutarget=mutarget*(dsdp->np/(dsdp->np+sqrt(dsdp->np))); 00201 info=DSDPComputeDY(dsdp,mutarget,dsdp->dy, &ppnorm);DSDPCHKERR(info); 00202 if (ppnorm<=0){ DSDPEventLogEnd(dsdp->ctime); break; } 00203 dsdp->pnorm=ppnorm; 00204 info=DSDPYStepLineSearch2(dsdp, mutarget, dsdp->dstep, dsdp->dy);DSDPCHKERR(info); 00205 DSDPEventLogEnd(dsdp->ctime); 00206 } 00207 if (attempt>0)dsdp->dstep=1.0; 00208 00209 dsdp->mutarget=DSDPMin(dsdp->mu,mutarget); 00210 00211 info=DSDPGetRR(dsdp,&dd1);DSDPCHKERR(info); 00212 if (dsdp->itnow==0 && dsdp->xmaker[0].pstep<1.0 && dd1> 0 && dsdp->pstep<1.0 && dsdp->goty0==DSDP_FALSE){ 00213 info=DSDPResetY0(dsdp);DSDPCHKERR(info); continue; 00214 dsdp->goty0=DSDP_FALSE; 00215 } 00216 00217 } /* End of Dual Scaling Algorithm */ 00218 00219 00220 DSDPFunctionReturn(0); 00221 } 00222 00223 00224 00225 #undef __FUNCT__ 00226 #define __FUNCT__ "DSDPChooseBarrierParameter" 00227 00240 int DSDPChooseBarrierParameter(DSDP dsdp, double mutarget, double *ppstep, double *nextmutarget){ 00241 int attempt,info,count=0; 00242 double pnorm,pstep=*ppstep,pmumin,mur, dmury1, mutarget2; 00243 DSDPTruth psdefinite=DSDP_FALSE; 00244 00245 DSDPFunctionBegin; 00246 00247 *nextmutarget=mutarget; 00248 /* Compute a feasible primal matrix with the smallest possible value of \mu */ 00249 /* Start by finding a psd feasible matrix */ 00250 if (*ppstep >=1){ 00251 pstep=1.0; 00252 /* PS should alread be formed and factored */ 00253 } else { 00254 mur=-1.0/mutarget; 00255 info=DSDPComputePDY(dsdp,mutarget,dsdp->dy,&pnorm); DSDPCHKERR(info); 00256 info=DSDPComputeMaxStepLength(dsdp,dsdp->dy,DUAL_FACTOR,&pstep);DSDPCHKERR(info); 00257 /* pstep=DSDPMin(0.97*pstep,1.0); */ 00258 if (pstep<1.0) {pstep=DSDPMin(0.97*pstep,1.0);} else {pstep=DSDPMin(1.0*pstep,1.0);} 00259 while (psdefinite==DSDP_FALSE){ 00260 if (count > 2 && pstep <1e-8){pstep=0;break;} 00261 info=DSDPComputePY(dsdp,pstep,dsdp->ytemp);DSDPCHKERR(info); 00262 info=DSDPComputeSS(dsdp,dsdp->ytemp,PRIMAL_FACTOR,&psdefinite);DSDPCHKERR(info); 00263 if (psdefinite==DSDP_FALSE){ 00264 if (count>1) pstep=0.5*pstep; else pstep=0.97*pstep; 00265 DSDPLogInfo(0,2,"Reducing pstep: %8.8e\n",pstep); 00266 count++; 00267 } 00268 } 00269 *ppstep=pstep; 00270 if (pstep > dsdp->xmaker[0].pstep || mutarget < dsdp->xmaker[0].mu * 1.0e-8){ 00271 info=DSDPSaveYForX(dsdp,mutarget,pstep);DSDPCHKERR(info); 00272 } 00273 if (pstep==0){ 00274 DSDPFunctionReturn(0); 00275 } 00276 } 00277 00278 /* Now determine how much smaller of mu can be used */ 00279 mur=pstep/mutarget; 00280 info=DSDPComputePDY1(dsdp,mur,dsdp->rhstemp); DSDPCHKERR(info); 00281 00282 /* Smallest value of mu that gives a positive definite matrix */ 00283 info=DSDPComputeMaxStepLength(dsdp,dsdp->rhstemp,PRIMAL_FACTOR,&dmury1);DSDPCHKERR(info); 00284 dmury1 = DSDPMin(1000,0.97*dmury1); 00285 00286 /* We could test the point, My tests say its not necessary -- its good! */ 00287 attempt=0;psdefinite=DSDP_FALSE; 00288 pmumin=mutarget / (1.0 + 1.0 * dmury1); /* This should be positive definite */ 00289 while ( 0 && psdefinite==DSDP_FALSE){ 00290 pmumin=mutarget / (1.0 + 1.0 * dmury1); /* This should be positive definite */ 00291 if (attempt>2){pmumin=mutarget;}/* We have actually factored this one. It is PSD. */ 00292 info=DSDPComputePDY(dsdp,pmumin,dsdp->dy,&pnorm); DSDPCHKERR(info); 00293 info=DSDPComputePY(dsdp,pstep,dsdp->ytemp); DSDPCHKERR(info); 00294 info=DSDPComputeSS(dsdp,dsdp->ytemp,PRIMAL_FACTOR,&psdefinite);DSDPCHKERR(info); 00295 if (psdefinite==DSDP_FALSE){ dmury1*=0.9; /* printf("NO GOOD \n"); */ } 00296 else { /* printf("ITS GOOD \n"); */} 00297 attempt++; 00298 DSDPLogInfo(0,2,"GOT X: Smallest Mu for feasible X: %4.4e.\n",pmumin); 00299 } 00300 00301 DSDPLogInfo(0,6,"GOT X: Smallest Mu for feasible X: %4.4e \n",pmumin); 00302 00303 mutarget2=mutarget; 00304 if (dsdp->pstep==1){ 00305 mutarget2 = pmumin; 00306 } else { 00307 /* printf("PMUMIN: %4.4e MUTARGET: %4.4e \n",pmumin,dsdp->mutarget); */ 00308 mutarget2 = mutarget; 00309 mutarget2 = (pstep)*mutarget + (1.0-pstep)*dsdp->mu; 00310 mutarget2 = (pstep)*pmumin + (1.0-pstep)*dsdp->mu; 00311 } 00312 00313 mutarget2=DSDPMax(dsdp->mu/dsdp->rhon,mutarget2); 00314 if (dsdp->mu0>0){mutarget2=DSDPMin(mutarget2,dsdp->mu0);} 00315 00316 *nextmutarget=mutarget2; 00317 DSDPFunctionReturn(0); 00318 } 00319 00320 #undef __FUNCT__ 00321 #define __FUNCT__ "DSDPResetY0" 00322 00328 int DSDPResetY0(DSDP dsdp){ 00329 int info; 00330 double r,rr,dd; 00331 DSDPTruth psdefinite; 00332 DSDPFunctionBegin; 00333 info=DSDPComputeDY(dsdp,dsdp->mutarget,dsdp->dy,&dsdp->pnorm); DSDPCHKERR(info); 00334 info=DSDPVecCopy(dsdp->y0,dsdp->y);DSDPCHKERR(info); 00335 info=DSDPGetRR(dsdp,&r);DSDPCHKERR(info); 00336 rr=DSDPMax(r*10000,1e12); 00337 info=DSDPSetRR(dsdp,rr);DSDPCHKERR(info); 00338 info=DSDPComputeSS(dsdp,dsdp->y,DUAL_FACTOR,&psdefinite);DSDPCHKERR(info); 00339 info=DSDPComputeLogSDeterminant(dsdp,&dsdp->logdet);DSDPCHKERR(info); 00340 info=DSDPSetY(dsdp,1.0,dsdp->logdet,dsdp->y);DSDPCHKERR(info); 00341 info=DSDPVecGetR(dsdp->b,&dd);DSDPCHKERR(info); 00342 00343 dsdp->mutarget=fabs(rr*dd); 00344 dsdp->mu=fabs(rr*dd); 00345 00346 /* dsdp->par.mu0=mutarget; */ 00347 dsdp->goty0=DSDP_TRUE; 00348 DSDPLogInfo(0,2,"Restart algorithm\n"); 00349 DSDPFunctionReturn(0); 00350 } 00351 00352 00368 #undef __FUNCT__ 00369 #define __FUNCT__ "DSDPComputeDualStepDirections" 00370 int DSDPComputeDualStepDirections(DSDP dsdp){ 00371 int info,computem=1; 00372 double madd,ymax,cgtol=1e-7; 00373 DSDPTruth cg1,cg2,psdefinite; 00374 DSDPFunctionBegin; 00375 00376 if (dsdp->itnow>30) dsdp->slestype=3; 00377 if (dsdp->rgap<1e-3) dsdp->slestype=3; 00378 if (dsdp->m<40) dsdp->slestype=3; 00379 if (0 && dsdp->itnow>20 && dsdp->m<500) dsdp->slestype=3; 00380 info=DSDPGetMaxYElement(dsdp,&ymax);DSDPCHKERR(info); 00381 if (dsdp->slestype==1){ 00382 cg1=DSDP_TRUE; cg2=DSDP_TRUE; 00383 info=DSDPInvertS(dsdp);DSDPCHKERR(info); 00384 info=DSDPComputeG(dsdp,dsdp->rhstemp,dsdp->rhs1,dsdp->rhs2);DSDPCHKERR(info); 00385 info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs1,dsdp->dy1,cgtol,&cg1);DSDPCHKERR(info); 00386 if (cg1==DSDP_TRUE){info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs2,dsdp->dy2,cgtol,&cg2);DSDPCHKERR(info);} 00387 if (cg1==DSDP_FALSE || cg2==DSDP_FALSE) dsdp->slestype=2; 00388 } 00389 if (dsdp->slestype==2){ 00390 cg1=DSDP_TRUE; cg2=DSDP_TRUE; 00391 DSDPLogInfo(0,9,"Compute Hessian\n"); 00392 info=DSDPInvertS(dsdp);DSDPCHKERR(info); 00393 info=DSDPComputeHessian(dsdp,dsdp->M,dsdp->rhs1,dsdp->rhs2);DSDPCHKERR(info); 00394 computem=0; 00395 DSDPLogInfo(0,9,"Apply CG\n"); 00396 info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs1,dsdp->dy1,cgtol,&cg1);DSDPCHKERR(info); 00397 if (cg1==DSDP_TRUE){info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs2,dsdp->dy2,cgtol,&cg2);DSDPCHKERR(info);} 00398 if (cg1==DSDP_FALSE || cg2==DSDP_FALSE) dsdp->slestype=3; 00399 00400 } 00401 if (dsdp->slestype==3){ 00402 DSDPLogInfo(0,9,"Factor Hessian\n"); 00403 psdefinite=DSDP_FALSE; 00404 if (dsdp->Mshift < 1e-12 || dsdp->rgap<0.1 || dsdp->Mshift > 1e-6){ 00405 madd=dsdp->Mshift; 00406 } else { 00407 madd=1e-13; 00408 } 00409 if (computem){ 00410 info=DSDPInvertS(dsdp);DSDPCHKERR(info); 00411 } 00412 while (psdefinite==DSDP_FALSE){ 00413 if (0==1 && dsdp->Mshift>dsdp->maxschurshift){ 00414 info = DSDPSetConvergenceFlag(dsdp,DSDP_INDEFINITE_SCHUR_MATRIX); DSDPCHKERR(info); 00415 break; 00416 } 00417 if (0 && dsdp->Mshift*ymax>dsdp->pinfeastol/10){ 00418 info = DSDPSetConvergenceFlag(dsdp,DSDP_INDEFINITE_SCHUR_MATRIX); DSDPCHKERR(info); 00419 break; 00420 } 00421 if (madd*ymax>dsdp->pinfeastol*1000){ 00422 info = DSDPSetConvergenceFlag(dsdp,DSDP_INDEFINITE_SCHUR_MATRIX); DSDPCHKERR(info); 00423 break; 00424 } 00425 if (computem){ 00426 info=DSDPComputeHessian(dsdp,dsdp->M,dsdp->rhs1,dsdp->rhs2);DSDPCHKERR(info);} 00427 if (0==1){info=DSDPSchurMatView(dsdp->M);DSDPCHKERR(info);} 00428 info = DSDPSchurMatShiftDiagonal(dsdp->M,madd);DSDPCHKERR(info); 00429 info = DSDPSchurMatFactor(dsdp->M,&psdefinite); DSDPCHKERR(info); 00430 computem=1; 00431 if (psdefinite==DSDP_FALSE){ 00432 madd=madd*4 + 1.0e-13; 00433 } 00434 } 00435 dsdp->Mshift=madd; 00436 if (psdefinite==DSDP_TRUE){ 00437 info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs1,dsdp->dy1,cgtol,&cg1);DSDPCHKERR(info); 00438 info=DSDPCGSolve(dsdp,dsdp->M,dsdp->rhs2,dsdp->dy2,cgtol,&cg2);DSDPCHKERR(info); 00439 } 00440 } 00441 00442 DSDPFunctionReturn(0); 00443 } 00444 #undef __FUNCT__ 00445 #define __FUNCT__ "DSDPComputeDualStepDirections" 00446 int DSDPRefineStepDirection(DSDP dsdp,DSDPVec rhs, DSDPVec dy){ 00447 int info; 00448 double cgtol=1e-20; 00449 DSDPTruth cg1; 00450 DSDPFunctionBegin; 00451 00452 if (dsdp->reason==DSDP_INDEFINITE_SCHUR_MATRIX){ 00453 } else if (dsdp->reason==DSDP_SMALL_STEPS){ 00454 } else if (dsdp->pstep<1){ 00455 } else { 00456 dsdp->slestype=4; 00457 info=DSDPCGSolve(dsdp,dsdp->M,rhs,dy,cgtol,&cg1);DSDPCHKERR(info); 00458 dsdp->slestype=3; 00459 } 00460 00461 DSDPFunctionReturn(0); 00462 } 00463 00464 #include "dsdp5.h" 00465 00466 00467 #undef __FUNCT__ 00468 #define __FUNCT__ "DSDPInitializeVariables" 00469 00475 int DSDPInitializeVariables( DSDP dsdp ){ 00476 00477 int info; 00478 double r0,mutarget=dsdp->mutarget,penalty; 00479 DSDPTruth psdefinite=DSDP_FALSE; 00480 00481 DSDPFunctionBegin; 00482 info=DSDPGetRR(dsdp,&r0);DSDPCHKERR(info); 00483 dsdp->rho=dsdp->np*dsdp->rhon; 00484 if (r0>=0) { /* Use the specified starting point */ 00485 info=DSDPComputeObjective(dsdp,dsdp->y,&dsdp->ddobj);DSDPCHKERR(info); 00486 info=DSDPComputeSS(dsdp,dsdp->y,DUAL_FACTOR,&psdefinite);DSDPCHKERR(info); 00487 if (mutarget<0) mutarget=(dsdp->ppobj-dsdp->ddobj)/(dsdp->rho); 00488 } else { 00489 info=DSDPGetPenalty(dsdp,&penalty);DSDPCHKERR(info); 00490 r0=0.1/(1.0+dsdp->cnorm); 00491 while (psdefinite==DSDP_FALSE ){ 00492 r0=r0*100; 00493 DSDPLogInfo(0,9,"Set Initial R0 %4.2e\n",r0); 00494 info=DSDPSetRR(dsdp,r0);DSDPCHKERR(info); 00495 info=DSDPComputeSS(dsdp,dsdp->y,DUAL_FACTOR,&psdefinite);DSDPCHKERR(info); 00496 } 00497 r0=r0*dsdp->np; 00498 if (dsdp->cnorm>0 && dsdp->anorm>0 && dsdp->cnorm/dsdp->anorm<1){ r0=r0/(dsdp->cnorm/dsdp->anorm);} 00499 dsdp->mu=r0*penalty; 00500 if (mutarget<0){ 00501 mutarget=(dsdp->ppobj-dsdp->ddobj)/(dsdp->rho); 00502 mutarget=(dsdp->ppobj-dsdp->ddobj)/(dsdp->np); 00503 mutarget=r0*penalty; 00504 } 00505 DSDPLogInfo(0,9,"Set Initial R0 %4.2e\n",r0); 00506 info=DSDPSetRR(dsdp,r0);DSDPCHKERR(info); 00507 info=DSDPComputeObjective(dsdp,dsdp->y,&dsdp->ddobj);DSDPCHKERR(info); 00508 info=DSDPComputeSS(dsdp,dsdp->y,DUAL_FACTOR,&psdefinite);DSDPCHKERR(info); 00509 } 00510 info=DSDPComputeObjective(dsdp,dsdp->y,&dsdp->ddobj);DSDPCHKERR(info); 00511 if (psdefinite==DSDP_FALSE){ 00512 info=DSDPSetConvergenceFlag(dsdp,DSDP_INFEASIBLE_START);DSDPCHKERR(info); 00513 } else { 00514 info=DSDPComputeLogSDeterminant(dsdp,&dsdp->logdet);DSDPCHKERR(info); 00515 info=DSDPComputePotential(dsdp,dsdp->y,dsdp->logdet,&dsdp->potential);DSDPCHKERR(info); 00516 } 00517 /* Tough guess, as all these rules suggest */ 00518 info=DSDPSetY(dsdp,0,dsdp->logdet,dsdp->y);DSDPCHKERR(info); 00519 info=DSDPSaveYForX(dsdp,dsdp->xmaker[0].mu,0);DSDPCHKERR(info); 00520 dsdp->mutarget=mutarget; 00521 dsdp->pstep=0.0; 00522 dsdp->pnorm=0; 00523 /* dsdp->par.mu0=mutarget; */ 00524 DSDPFunctionReturn(0); 00525 } 00526 00538 #undef __FUNCT__ 00539 #define __FUNCT__ "DSDPComputeAndFactorS" 00540 int DSDPComputeAndFactorS(DSDP dsdp,DSDPTruth *psdefinite){ 00541 int info; 00542 DSDPFunctionBegin; 00543 info=DSDPComputeSS(dsdp,dsdp->y,DUAL_FACTOR,psdefinite);DSDPCHKERR(info); 00544 DSDPFunctionReturn(0); 00545 }