DSDP
|
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 00337 #undef __FUNCT__ 00338 #define __FUNCT__ "DSDPSetRR" 00339 00345 int DSDPSetRR(DSDP dsdp,double res){ 00346 int info; 00347 DSDPFunctionBegin; 00348 DSDPValid(dsdp); 00349 info=DSDPVecSetR(dsdp->y,-res);DSDPCHKERR(info); 00350 DSDPFunctionReturn(0); 00351 } 00352 00353 #undef __FUNCT__ 00354 #define __FUNCT__ "DSDPGetRR" 00355 00361 int DSDPGetRR(DSDP dsdp,double *res){ 00362 int info; 00363 DSDPFunctionBegin; 00364 DSDPValid(dsdp); 00365 info=DSDPVecGetR(dsdp->y,res);DSDPCHKERR(info); 00366 *res=-*res; 00367 if (*res==0) *res=0; 00368 DSDPFunctionReturn(0); 00369 } 00370 00371 00372 #undef __FUNCT__ 00373 #define __FUNCT__ "DSDPObjectiveGH" 00374 00381 int DSDPObjectiveGH( DSDP dsdp , DSDPSchurMat M, DSDPVec vrhs1){ 00382 int i,info,m; 00383 double rtemp,dd; 00384 00385 DSDPFunctionBegin; 00386 info=DSDPVecGetSize(vrhs1,&m); DSDPCHKERR(info); 00387 for (i=0;i<m;i++){ 00388 info=DSDPSchurMatVariableCompute(M,i,&dd); DSDPCHKERR(info); 00389 if (dd){ 00390 info=DSDPVecGetElement(dsdp->b,i,&rtemp);DSDPCHKERR(info); 00391 info=DSDPVecAddElement(vrhs1,i,rtemp);DSDPCHKERR(info); 00392 } 00393 } 00394 DSDPFunctionReturn(0); 00395 } 00396 00397 #undef __FUNCT__ 00398 #define __FUNCT__ "DSDPCheckForUnboundedObjective" 00399 int DSDPCheckForUnboundedObjective(DSDP dsdp, DSDPTruth *unbounded){ 00400 int info; 00401 double dtemp; 00402 DSDPTruth psdefinite; 00403 DSDPFunctionBegin; 00404 *unbounded=DSDP_FALSE; 00405 info = DSDPVecDot(dsdp->b,dsdp->dy2,&dtemp);DSDPCHKERR(info); 00406 if ( dtemp < 0 /* && dsdp->r==0 && dsdp->ddobj > 0 */) { 00407 info = DSDPVecScaleCopy(dsdp->dy2,-1.0,dsdp->ytemp); DSDPCHKERR(info); 00408 info = DSDPComputeSS(dsdp,dsdp->ytemp,PRIMAL_FACTOR,&psdefinite);DSDPCHKERR(info); 00409 if (psdefinite == DSDP_TRUE){ 00410 psdefinite=DSDP_FALSE; 00411 while(psdefinite==DSDP_FALSE){ /* Dual point should be a certificate of dual unboundedness, and be dual feasible */ 00412 info=DSDPComputeSS(dsdp,dsdp->ytemp,PRIMAL_FACTOR,&psdefinite);DSDPCHKERR(info); 00413 if (psdefinite == DSDP_TRUE) break; 00414 info=DSDPVecScale(2.0,dsdp->ytemp); DSDPCHKERR(info); 00415 } 00416 info = DSDPVecCopy(dsdp->ytemp,dsdp->y); DSDPCHKERR(info); 00417 info = DSDPSaveYForX(dsdp,1.0e-12,1.0);DSDPCHKERR(info); 00418 info = DSDPComputeObjective(dsdp,dsdp->y,&dsdp->ddobj);DSDPCHKERR(info); 00419 info = DSDPVecNormalize(dsdp->y); DSDPCHKERR(info); 00420 *unbounded=DSDP_TRUE; 00421 } 00422 } 00423 DSDPFunctionReturn(0); 00424 } 00425