Main Page | Modules | Alphabetical List | Data Structures | Directories | File List | Globals | Related Pages

dualimpl.c

00001 #include "dsdp.h"
00002 #include "dsdpsys.h"
00019 #undef __FUNCT__
00020 #define __FUNCT__ "DSDPComputeObjective"
00021 int DSDPComputeObjective(DSDP dsdp, DSDPVec Y, double *ddobj){
00022   int info;
00023   DSDPFunctionBegin;
00024   info = DSDPVecDot(Y,dsdp->b,ddobj);DSDPCHKERR(info);
00025   DSDPFunctionReturn(0);
00026 }
00027 
00028 
00043 #undef __FUNCT__
00044 #define __FUNCT__ "DSDPComputeDY"
00045 int DSDPComputeDY(DSDP dsdp, double mu, DSDPVec DY, double *pnorm){
00046   int info;
00047   double ppnorm,ddy1=fabs(1.0/mu*dsdp->schurmu),ddy2=-1.0;
00048   DSDPFunctionBegin;
00049   info=DSDPComputeRHS(dsdp,mu,dsdp->rhs); DSDPCHKERR(info);
00050   info=DSDPVecWAXPBY(DY,ddy1,dsdp->dy1,ddy2,dsdp->dy2);DSDPCHKERR(info);
00051   info=DSDPComputePNorm(dsdp,mu,DY,&ppnorm);DSDPCHKERR(info);
00052   if (ppnorm<0){ /* If pnorm < 0 there are SMW numerical issues */
00053     DSDPLogInfo(0,2,"Problem with PNORM: %4.4e < 0 \n",ppnorm);
00054     /*    ppnorm=1.0; */
00055   }
00056   *pnorm=ppnorm;
00057   DSDPFunctionReturn(0);
00058 }
00059 
00075 #undef __FUNCT__
00076 #define __FUNCT__ "DSDPComputePDY"
00077 int DSDPComputePDY(DSDP dsdp, double mu, DSDPVec DY, double *pnorm){
00078   int info;
00079   double ppnorm,ddy1=-fabs(1.0/mu*dsdp->schurmu),ddy2=1.0;
00080   DSDPFunctionBegin;
00081   info=DSDPComputeRHS(dsdp,mu,dsdp->rhs); DSDPCHKERR(info);
00082   info=DSDPVecWAXPBY(DY,ddy1,dsdp->dy1,ddy2,dsdp->dy2);DSDPCHKERR(info);
00083   info=DSDPComputePNorm(dsdp,mu,DY,&ppnorm);DSDPCHKERR(info);
00084   if (ppnorm<0){ /* If pnorm < 0 there are SMW numerical issues */
00085     DSDPLogInfo(0,2,"Problem with PNORM: %4.4e < 0 \n",ppnorm);
00086     /*    ppnorm=1.0; */
00087   }
00088   *pnorm=ppnorm;
00089   DSDPFunctionReturn(0);
00090 }
00091 
00103 #undef __FUNCT__
00104 #define __FUNCT__ "DSDPComputePDY1"
00105 int DSDPComputePDY1(DSDP dsdp, double mur, DSDPVec DY1){
00106   int info;
00107   double ddy1=-fabs(mur*dsdp->schurmu);
00108   DSDPFunctionBegin;
00109   info=DSDPVecScaleCopy(dsdp->dy1,ddy1,DY1); DSDPCHKERR(info);
00110   DSDPFunctionReturn(0);
00111 }
00112 
00123 #undef __FUNCT__
00124 #define __FUNCT__ "DSDPComputeNewY"
00125 int DSDPComputeNewY(DSDP dsdp, double beta, DSDPVec Y){
00126   int info;
00127   double rtemp;
00128   DSDPFunctionBegin;
00129   info=DSDPVecWAXPY(Y,beta,dsdp->dy,dsdp->y);DSDPCHKERR(info);
00130   info=DSDPVecGetR(Y,&rtemp);DSDPCHKERR(info);
00131   rtemp=DSDPMin(0,rtemp);
00132   info=DSDPSchurMatSetR(dsdp->M,rtemp);DSDPCHKERR(info);
00133   info=DSDPVecSetR(Y,rtemp);DSDPCHKERR(info);
00134   info=DSDPApplyFixedVariables(dsdp->M,Y);DSDPCHKERR(info);
00135   DSDPFunctionReturn(0);
00136 }
00137 
00148 #undef __FUNCT__
00149 #define __FUNCT__ "DSDPComputePY"
00150 int DSDPComputePY(DSDP dsdp, double beta, DSDPVec PY){
00151   int info;
00152   DSDPFunctionBegin;
00153   info=DSDPVecWAXPY(PY,beta,dsdp->dy,dsdp->y);DSDPCHKERR(info);
00154   info=DSDPApplyFixedVariables(dsdp->M,PY);DSDPCHKERR(info);
00155   DSDPFunctionReturn(0);
00156 }
00157 
00175 #undef __FUNCT__
00176 #define __FUNCT__ "DSDPComputeRHS"
00177 int DSDPComputeRHS(DSDP dsdp, double mu, DSDPVec RHS){
00178   int info;
00179   double ddrhs1=1.0/mu*dsdp->schurmu,ddrhs2=-( mu/fabs(mu) );
00180   DSDPFunctionBegin;
00181   info=DSDPVecWAXPBY(RHS,ddrhs1,dsdp->rhs1,ddrhs2,dsdp->rhs2);DSDPCHKERR(info);
00182   DSDPFunctionReturn(0);
00183 }
00184 
00185 #undef __FUNCT__
00186 #define __FUNCT__ "DSDPComputePNorm"
00187 
00200 int DSDPComputePNorm(DSDP dsdp, double mu, DSDPVec DY, double *pnorm){
00201   int info;
00202   double ppnorm=0;
00203   DSDPFunctionBegin;
00204   info=DSDPComputeRHS(dsdp,mu,dsdp->rhs); DSDPCHKERR(info);
00205   info = DSDPVecDot(dsdp->rhs,DY,&ppnorm);DSDPCHKERR(info);
00206   ppnorm/=dsdp->schurmu;
00207   if (ppnorm >= 0){  /* Theoretically True */
00208     *pnorm=sqrt(ppnorm);
00209   } else {
00210     DSDPLogInfo(0,2,"Problem with PNORM: %4.4e is not positive.\n",ppnorm);
00211     *pnorm=ppnorm;
00212   }
00213   if (*pnorm!=*pnorm){DSDPSETERR1(1,"Problem with PNORM: %4.4e is not positive.\n",ppnorm);}
00214   DSDPFunctionReturn(0);
00215 }
00216 
00228 #undef __FUNCT__
00229 #define __FUNCT__ "DSDPComputeDualityGap"
00230 int DSDPComputeDualityGap(DSDP dsdp, double mu, double *gap){
00231   int info;
00232   double newgap=0,pnorm;
00233   double smu=1.0/dsdp->schurmu;
00234   DSDPFunctionBegin;
00235   info=DSDPComputeDY(dsdp,mu,dsdp->dy,&pnorm); DSDPCHKERR(info);
00236   info=DSDPVecDot(dsdp->dy,dsdp->rhs2,&newgap);DSDPCHKERR(info);
00237   newgap = (newgap*smu+dsdp->np)*mu;
00238   if (newgap<=0){
00239     DSDPLogInfo(0,2,"GAP :%4.4e<0: Problem\n",newgap);
00240   } else {
00241     DSDPLogInfo(0,2,"Duality Gap: %12.8e, Update primal objective: %12.8e\n",newgap,dsdp->ddobj+newgap);
00242   }
00243   newgap=DSDPMax(0,newgap);
00244   *gap=newgap;
00245   DSDPFunctionReturn(0);
00246 }
00247 
00259 #undef __FUNCT__
00260 #define __FUNCT__ "DSDPComputePotential"
00261 int DSDPComputePotential(DSDP dsdp, DSDPVec y, double  logdet, double *potential){
00262   int info;
00263   double dpotential,gap,ddobj;
00264   DSDPFunctionBegin;
00265   info=DSDPComputeObjective(dsdp,y,&ddobj);DSDPCHKERR(info);
00266   gap=dsdp->ppobj-ddobj;
00267   if (gap>0) dpotential=dsdp->rho*log(gap)-logdet;
00268   else {dpotential=dsdp->potential+1;}
00269   *potential=dpotential;
00270   DSDPLogInfo(0,9,"Gap: %4.4e, Log Determinant: %4.4e, Log Gap: %4.4e\n",gap,logdet,log(gap));
00271   DSDPFunctionReturn(0);
00272 }
00273 
00285 #undef __FUNCT__
00286 #define __FUNCT__ "DSDPComputePotential2"
00287 int DSDPComputePotential2(DSDP dsdp, DSDPVec y, double mu, double logdet, double *potential){
00288   int info;
00289   double ddobj;
00290   DSDPFunctionBegin;
00291   info=DSDPComputeObjective(dsdp,y,&ddobj);DSDPCHKERR(info);
00292   *potential=-(ddobj + mu*logdet)*dsdp->schurmu;
00293   *potential=-(ddobj/mu + logdet)*dsdp->schurmu;
00294   DSDPFunctionReturn(0);
00295 }
00296 
00297 
00307 #undef __FUNCT__
00308 #define __FUNCT__ "DSDPSetY"
00309 int DSDPSetY(DSDP dsdp, double beta, double logdet, DSDPVec ynew){
00310   int info;
00311   double r1,r2,rr,pp;
00312   DSDPFunctionBegin;
00313   info=DSDPVecGetR(dsdp->y,&r1);DSDPCHKERR(info);
00314   info=DSDPVecGetR(ynew,&r2);DSDPCHKERR(info);
00315   if (r2==0&&r1!=0){dsdp->rflag=1;} else {dsdp->rflag=0;};
00316   info=DSDPVecCopy(ynew,dsdp->y);DSDPCHKERR(info);
00317   info = DSDPComputeObjective(dsdp,dsdp->y,&dsdp->ddobj);DSDPCHKERR(info);
00318   /* Correct ppobj if ppobj < ddobj , which can happen when dual infeasibility is present */
00319   if (dsdp->ppobj<=dsdp->ddobj){
00320     dsdp->ppobj=dsdp->ddobj+2*dsdp->mu * dsdp->np;
00321     DSDPLogInfo(0,2,"Primal Objective Not Right.  Assigned: %8.8e\n",dsdp->ppobj);
00322   }
00323   info=DSDPVecGetR(ynew,&rr);DSDPCHKERR(info);
00324   info=DSDPVecGetR(dsdp->b,&pp);DSDPCHKERR(info);
00325   dsdp->dobj=dsdp->ddobj-rr*pp;
00326   DSDPLogInfo(0,2,"Duality Gap: %4.4e, Potential: %4.4e \n",dsdp->dualitygap,dsdp->potential);
00327   dsdp->dualitygap=dsdp->ppobj-dsdp->ddobj;
00328   dsdp->mu=(dsdp->dualitygap)/(dsdp->np);
00329   dsdp->dstep=beta;
00330   dsdp->logdet=logdet;
00331   info=DSDPComputePotential(dsdp,dsdp->y,dsdp->logdet,&dsdp->potential);DSDPCHKERR(info);
00332   DSDPLogInfo(0,2,"Duality Gap: %4.4e, Potential: %4.4e \n",dsdp->dualitygap,dsdp->potential);
00333   DSDPFunctionReturn(0);
00334 }
00335 
00336 #undef __FUNCT__
00337 #define __FUNCT__ "DSDPSaveBackupYForX"
00338 int DSDPSaveBackupYForX(DSDP dsdp, int count,double mu, double pstep){
00339   int info;
00340   double pnorm;
00341   DSDPFunctionBegin;
00342   info=DSDPVecCopy(dsdp->y,dsdp->xmaker[count].y);DSDPCHKERR(info);
00343   info=DSDPComputeDY(dsdp,mu,dsdp->xmaker[count].dy,&pnorm); DSDPCHKERR(info);
00344   dsdp->xmaker[count].pstep=pstep;
00345   dsdp->xmaker[count].mu = mu;
00346   DSDPFunctionReturn(0);
00347 }
00348 
00357 #undef __FUNCT__
00358 #define __FUNCT__ "DSDPSaveYForX"
00359 int DSDPSaveYForX(DSDP dsdp, double mu, double pstep){
00360   int info;
00361   double pnorm,newgap,ymax;
00362   DSDPFunctionBegin;
00363   info=DSDPGetMaxYElement(dsdp,&ymax);DSDPCHKERR(info);
00364   if (pstep==0){
00365     info=DSDPVecCopy(dsdp->y,dsdp->xmaker[0].y);DSDPCHKERR(info);
00366     dsdp->xmaker[0].pstep=pstep;
00367   } else if (dsdp->Mshift*ymax>dsdp->pinfeastol*10){
00368     info=DSDPComputeDualityGap(dsdp,mu,&newgap);DSDPCHKERR(info);
00369     if (pstep==1 && newgap>0){
00370       dsdp->ppobj = dsdp->ddobj + newgap; dsdp->mu=(newgap)/(dsdp->np);
00371       dsdp->dualitygap=newgap;
00372     }
00373     /* Not good enough to save */
00374   } else {
00375     info=DSDPVecCopy(dsdp->y,dsdp->xmaker[0].y);DSDPCHKERR(info);
00376     dsdp->xmaker[0].pstep=pstep;
00377     info=DSDPComputeRHS(dsdp,mu,dsdp->xmakerrhs); DSDPCHKERR(info);
00378     info=DSDPComputeDY(dsdp,mu,dsdp->xmaker[0].dy,&pnorm); DSDPCHKERR(info);
00379     dsdp->xmaker[0].mu = mu;
00380     info=DSDPComputeDualityGap(dsdp,mu,&newgap);DSDPCHKERR(info);
00381     if (pstep==1 && newgap>0){
00382       dsdp->ppobj = dsdp->ddobj + newgap; dsdp->mu=(newgap)/(dsdp->np);
00383       dsdp->dualitygap=newgap;
00384     }
00385   }
00386 
00387   if (pstep==1.0 && dsdp->rgap>5.0e-1) {
00388     info=DSDPSaveBackupYForX(dsdp,MAX_XMAKERS-1,mu,pstep);DSDPCHKERR(info);
00389   }
00390   if (pstep==1.0 && dsdp->rgap>1.0e-3) {
00391     info=DSDPSaveBackupYForX(dsdp,2,mu,pstep);DSDPCHKERR(info);
00392   }
00393   if (pstep==1.0 && dsdp->rgap>1.0e-5) {
00394     info=DSDPSaveBackupYForX(dsdp,1,mu,pstep);DSDPCHKERR(info);
00395   }
00396   
00397   DSDPFunctionReturn(0);
00398 }
00399 
00400 #undef __FUNCT__  
00401 #define __FUNCT__ "DSDPSetRR"
00402 
00408 int DSDPSetRR(DSDP dsdp,double res){ 
00409   int info;
00410   DSDPFunctionBegin;
00411   DSDPValid(dsdp);
00412   info=DSDPVecSetR(dsdp->y,-res);DSDPCHKERR(info);
00413   DSDPFunctionReturn(0);
00414 }
00415 
00416 #undef __FUNCT__  
00417 #define __FUNCT__ "DSDPGetRR"
00418 
00424 int DSDPGetRR(DSDP dsdp,double *res){ 
00425   int info;
00426   DSDPFunctionBegin;
00427   DSDPValid(dsdp);
00428   info=DSDPVecGetR(dsdp->y,res);DSDPCHKERR(info);
00429   *res=-*res;
00430   if (*res==0) *res=0;
00431   DSDPFunctionReturn(0);
00432 }
00433 
00434 
00435 #undef __FUNCT__  
00436 #define __FUNCT__ "DSDPObjectiveGH"
00437 
00444 int DSDPObjectiveGH( DSDP dsdp , DSDPSchurMat M, DSDPVec vrhs1){
00445   int i,info,m;
00446   double rtemp,dd;
00447 
00448   DSDPFunctionBegin;
00449   info=DSDPVecGetSize(vrhs1,&m); DSDPCHKERR(info);
00450   for (i=0;i<m;i++){
00451     info=DSDPSchurMatVariableCompute(M,i,&dd); DSDPCHKERR(info);
00452     if (dd){
00453       info=DSDPVecGetElement(dsdp->b,i,&rtemp);DSDPCHKERR(info);
00454       info=DSDPVecAddElement(vrhs1,i,rtemp);DSDPCHKERR(info);
00455     }
00456   }
00457   DSDPFunctionReturn(0);
00458 }
00459 
00460 #undef __FUNCT__
00461 #define __FUNCT__ "DSDPCheckForUnboundedObjective"
00462 int DSDPCheckForUnboundedObjective(DSDP dsdp, DSDPTruth *unbounded){
00463   int info;
00464   double dtemp;
00465   DSDPTruth psdefinite;
00466   DSDPFunctionBegin;
00467   *unbounded=DSDP_FALSE;
00468   info = DSDPVecDot(dsdp->b,dsdp->dy2,&dtemp);DSDPCHKERR(info);
00469   if ( dtemp < 0 /* && dsdp->r==0 && dsdp->ddobj > 0 */) {  
00470     info = DSDPVecScaleCopy(dsdp->dy2,-1.0,dsdp->ytemp); DSDPCHKERR(info);
00471     info = DSDPComputeSS(dsdp,dsdp->ytemp,PRIMAL_FACTOR,&psdefinite);DSDPCHKERR(info);
00472     if (psdefinite == DSDP_TRUE){
00473       psdefinite=DSDP_FALSE;
00474       while(psdefinite==DSDP_FALSE){ /* Dual point should be a certificate of dual unboundedness, and be dual feasible */
00475         info=DSDPComputeSS(dsdp,dsdp->ytemp,PRIMAL_FACTOR,&psdefinite);DSDPCHKERR(info);
00476         if (psdefinite == DSDP_TRUE) break;
00477         info=DSDPVecScale(2.0,dsdp->ytemp); DSDPCHKERR(info);
00478       }
00479       info = DSDPVecCopy(dsdp->ytemp,dsdp->y); DSDPCHKERR(info);
00480       info = DSDPSaveYForX(dsdp,1.0e-12,1.0);DSDPCHKERR(info);
00481       info = DSDPComputeObjective(dsdp,dsdp->y,&dsdp->ddobj);DSDPCHKERR(info);
00482       info = DSDPVecNormalize(dsdp->y); DSDPCHKERR(info);
00483       *unbounded=DSDP_TRUE;
00484     }
00485   }
00486   DSDPFunctionReturn(0);
00487 }
00488 

Generated on Sat Oct 15 11:05:40 2005 for DSDP by  doxygen 1.4.2