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
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 }
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
00075
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 }
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
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
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
00169
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 }
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
00249
00250 if (*ppstep >=1){
00251 pstep=1.0;
00252
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
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
00279 mur=pstep/mutarget;
00280 info=DSDPComputePDY1(dsdp,mur,dsdp->rhstemp); DSDPCHKERR(info);
00281
00282
00283 info=DSDPComputeMaxStepLength(dsdp,dsdp->rhstemp,PRIMAL_FACTOR,&dmury1);DSDPCHKERR(info);
00284 dmury1 = DSDPMin(1000,0.97*dmury1);
00285
00286
00287 attempt=0;psdefinite=DSDP_FALSE;
00288 pmumin=mutarget / (1.0 + 1.0 * dmury1);
00289 while ( 0 && psdefinite==DSDP_FALSE){
00290 pmumin=mutarget / (1.0 + 1.0 * dmury1);
00291 if (attempt>2){pmumin=mutarget;}
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; }
00296 else { }
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
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
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) {
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
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
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 }