DSDP
dualalg.c
Go to the documentation of this file.
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 }