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){
00053 DSDPLogInfo(0,2,"Problem with PNORM: %4.4e < 0 \n",ppnorm);
00054
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){
00085 DSDPLogInfo(0,2,"Problem with PNORM: %4.4e < 0 \n",ppnorm);
00086
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){
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
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){
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