00001 #include "dsdp.h"
00002 #include "dsdpsys.h"
00003 #include "dsdp5.h"
00009 #undef __FUNCT__
00010 #define __FUNCT__ "DSDPInspectXY"
00011 int DSDPInspectXY(DSDP dsdp, double xmakermu, DSDPVec xmakery, DSDPVec xmakerdy, DSDPVec AX, double *tracexs2, double *pobj2, double *rpinfeas2){
00012 int info;
00013 DSDPFunctionBegin;
00014
00015 info=BoundYConeAddX(dsdp->ybcone,xmakermu,xmakery,xmakerdy,AX,tracexs2); DSDPCHKERR(info);
00016 info=DSDPVecGetC(AX,pobj2);DSDPCHKERR(info);
00017
00018 info=DSDPVecSetC(AX,0);DSDPCHKERR(info);
00019 info=DSDPVecSetR(AX,0);DSDPCHKERR(info);
00020 info=DSDPVecNorm1(AX,rpinfeas2); DSDPCHKERR(info);
00021 DSDPFunctionReturn(0);
00022 }
00023
00053 #undef __FUNCT__
00054 #define __FUNCT__ "DSDPComputeX"
00055 int DSDPComputeX(DSDP dsdp){
00056 int i,info;
00057 double pobj=0,ppobj2=0,ddobj,tracexs=0,tracexs2=0,rpinfeas=0,rpinfeas2=0,rpobjerr=0;
00058 double err1,cc,rrr,bigM,ymax,pfeastol=dsdp->pinfeastol;
00059 DSDPTerminationReason reason;
00060 DSDPVec AX=dsdp->ytemp;
00061
00062 DSDPFunctionBegin;
00063 info=DSDPStopReason(dsdp,&reason);DSDPCHKERR(info);
00064 info=DSDPGetDDObjective(dsdp,&ddobj);DSDPCHKERR(info);
00065 info=DSDPGetMaxYElement(dsdp,&ymax);DSDPCHKERR(info);
00066 info=DSDPGetR(dsdp,&rrr); DSDPCHKERR(info);
00067 info=DSDPGetPenalty(dsdp,&bigM);DSDPCHKERR(info);
00068 info=DSDPGetScale(dsdp,&cc);DSDPCHKERR(info);
00069
00070 dsdp->pdfeasible=DSDP_PDFEASIBLE;
00071 for (i=0;i<MAX_XMAKERS;i++){
00072 if (i>0 && dsdp->xmaker[i].pstep<1) continue;
00073 info=DSDPComputeXVariables(dsdp,dsdp->xmaker[i].mu,dsdp->xmaker[i].y,dsdp->xmaker[i].dy,AX,&tracexs);DSDPCHKERR(info);
00074 info=DSDPVecGetC(AX,&pobj); DSDPCHKERR(info);
00075 info=DSDPVecGetR(AX,&dsdp->tracex); DSDPCHKERR(info);
00076 info=DSDPVecSetC(AX,0);DSDPCHKERR(info);
00077 info=DSDPVecSetR(AX,0);DSDPCHKERR(info);
00078 info=DSDPVecNormInfinity(AX,&rpinfeas);DSDPCHKERR(info);
00079 rpinfeas=rpinfeas/(dsdp->tracex+1);
00080
00081 DSDPLogInfo(0,2,"POBJ: %4.4e, DOBJ: %4.4e\n",pobj,ddobj/cc);
00082
00083 info=DSDPVecNorm2(AX,&err1);DSDPCHKERR(info);
00084 dsdp->tracexs=tracexs;
00085 dsdp->perror=err1;
00086 dsdp->pobj=cc*pobj;
00087
00088 info=DSDPInspectXY(dsdp,dsdp->xmaker[i].mu,dsdp->xmaker[i].y,dsdp->xmaker[i].dy,AX,&tracexs2,&ppobj2,&rpinfeas2);DSDPCHKERR(info);
00089 rpinfeas2=rpinfeas2/(dsdp->tracex+1);
00090
00091
00092 DSDPLogInfo(0,2,"X P Infeas: %4.2e , PObj: %4.8e\n",rpinfeas,pobj*(cc));
00093 DSDPLogInfo(0,2,"TOTAL P Infeas: %4.2e PObj: %4.8e\n",rpinfeas2,ppobj2*(cc));
00094 rpobjerr= fabs(pobj-dsdp->ppobj)/(1+fabs(dsdp->ppobj));
00095
00096 if (rpinfeas2 < pfeastol){
00097
00098 if (dsdp->rgap<0.1){
00099 if (rpinfeas>pfeastol/100 && fabs(rrr)>dsdp->dinfeastol){
00100 dsdp->pdfeasible=DSDP_PDUNKNOWN;
00101 DSDPLogInfo(0,2,"Warning: Try Increasing penalty parameter\n");
00102 } else if (rpinfeas>pfeastol && ddobj>0 && ppobj2<0 && fabs(rrr)<dsdp->dinfeastol){
00103 dsdp->pdfeasible=DSDP_UNBOUNDED;
00104 DSDPLogInfo(0,2,"Warning: D probably unbounded\n");
00105
00106 } else if ( fabs(rrr)>dsdp->dinfeastol){
00107 dsdp->pdfeasible=DSDP_INFEASIBLE;
00108 DSDPLogInfo(0,2,"Warning: D probably infeasible \n");
00109 }
00110 }
00111 i=i+10;
00112 break;
00113
00114 } else {
00115
00116 DSDPLogInfo(0,2,"Try backup X\n");
00117 info=DSDPSetConvergenceFlag(dsdp,DSDP_NUMERICAL_ERROR); DSDPCHKERR(info);
00118 }
00119
00120 }
00121
00122 DSDPFunctionReturn(0);
00123 }
00124
00125
00126 #undef __FUNCT__
00127 #define __FUNCT__ "DSDPSaveBackupYForX"
00128 int DSDPSaveBackupYForX(DSDP dsdp, int count,double mu, double pstep){
00129 int info;
00130 double pnorm;
00131 DSDPFunctionBegin;
00132 info=DSDPVecCopy(dsdp->y,dsdp->xmaker[count].y);DSDPCHKERR(info);
00133 info=DSDPComputeDY(dsdp,mu,dsdp->xmaker[count].dy,&pnorm); DSDPCHKERR(info);
00134 dsdp->xmaker[count].pstep=pstep;
00135 dsdp->xmaker[count].mu = mu;
00136 DSDPFunctionReturn(0);
00137 }
00138
00147 #undef __FUNCT__
00148 #define __FUNCT__ "DSDPSaveYForX"
00149 int DSDPSaveYForX(DSDP dsdp, double mu, double pstep){
00150 int info,isgood=0;
00151 double pnorm,newgap,ymax,dd=0;
00152 DSDPFunctionBegin;
00153 dsdp->pstepold=dsdp->pstep;
00154 info=DSDPGetMaxYElement(dsdp,&ymax);DSDPCHKERR(info);
00155 if (pstep==0){
00156 info=DSDPVecCopy(dsdp->y,dsdp->xmaker[0].y);DSDPCHKERR(info);
00157 dsdp->xmaker[0].pstep=pstep;
00158 } else if (dsdp->Mshift*ymax>dsdp->pinfeastol*10){
00159 info=DSDPComputeDualityGap(dsdp,mu,&newgap);DSDPCHKERR(info);
00160 if (pstep==1 && newgap>0){
00161 dsdp->ppobj = dsdp->ddobj + newgap; dsdp->mu=(newgap)/(dsdp->np);
00162 dsdp->dualitygap=newgap;
00163 }
00164 info=DSDPVecZero(dsdp->rhstemp); DSDPCHKERR(info);
00165 info=BoundYConeAddX(dsdp->ybcone,dsdp->xmaker[0].mu,dsdp->xmaker[0].y,dsdp->xmaker[0].dy,dsdp->rhstemp,&dd); DSDPCHKERR(info);
00166 info=DSDPVecSetC(dsdp->rhstemp,0);
00167 info=DSDPVecSetR(dsdp->rhstemp,0);
00168 info=DSDPVecNormInfinity(dsdp->rhstemp,&dsdp->pinfeas); DSDPCHKERR(info);
00169 dsdp->pinfeas+=dsdp->Mshift*ymax;
00170 if (0==1){info=DSDPVecView(dsdp->rhstemp);}
00171
00172 } else {
00173 info=DSDPVecCopy(dsdp->y,dsdp->xmaker[0].y);DSDPCHKERR(info);
00174 dsdp->xmaker[0].pstep=pstep;
00175 info=DSDPComputeRHS(dsdp,mu,dsdp->xmakerrhs); DSDPCHKERR(info);
00176 info=DSDPComputeDY(dsdp,mu,dsdp->xmaker[0].dy,&pnorm); DSDPCHKERR(info);
00177 dsdp->xmaker[0].mu = mu;
00178 info=DSDPComputeDualityGap(dsdp,mu,&newgap);DSDPCHKERR(info);
00179 if (pstep==1 && newgap>0){
00180 dsdp->ppobj = dsdp->ddobj + newgap; dsdp->mu=(newgap)/(dsdp->np);
00181 dsdp->dualitygap=newgap;
00182
00183 info=DSDPVecZero(dsdp->rhstemp); DSDPCHKERR(info);
00184 info=BoundYConeAddX(dsdp->ybcone,dsdp->xmaker[0].mu,dsdp->xmaker[0].y,dsdp->xmaker[0].dy,dsdp->rhstemp,&dd); DSDPCHKERR(info);
00185 info=DSDPVecSetC(dsdp->rhstemp,0);
00186 info=DSDPVecSetR(dsdp->rhstemp,0);
00187 info=DSDPVecNormInfinity(dsdp->rhstemp,&dsdp->pinfeas); DSDPCHKERR(info);
00188 dsdp->pinfeas+=dsdp->Mshift*ymax;
00189 if (0==1){info=DSDPVecView(dsdp->rhstemp);}
00190
00191 }
00192 isgood=1;
00193 }
00194
00195 if (isgood==1){
00196 double penalty,r,rx;
00197 info=DSDPPassXVectors(dsdp,dsdp->xmaker[0].mu,dsdp->xmaker[0].y,dsdp->xmaker[0].dy); DSDPCHKERR(info);
00198 info=DSDPGetRR(dsdp,&r);DSDPCHKERR(info);
00199 if (r&& dsdp->rgap<0.1){
00200 info=RConeGetRX(dsdp->rcone,&rx);DSDPCHKERR(info);
00201 info=DSDPGetPenalty(dsdp,&penalty);DSDPCHKERR(info);
00202 dsdp->pinfeas = dsdp->pinfeas *(1+fabs(penalty-rx));
00203 }
00204 }
00205
00206 if (pstep==1.0 && dsdp->rgap>5.0e-1) {
00207 info=DSDPSaveBackupYForX(dsdp,MAX_XMAKERS-1,mu,pstep);DSDPCHKERR(info);
00208 }
00209 if (pstep==1.0 && dsdp->rgap>1.0e-3) {
00210 info=DSDPSaveBackupYForX(dsdp,2,mu,pstep);DSDPCHKERR(info);
00211 }
00212 if (pstep==1.0 && dsdp->rgap>1.0e-5) {
00213 info=DSDPSaveBackupYForX(dsdp,1,mu,pstep);DSDPCHKERR(info);
00214 }
00215
00216 DSDPFunctionReturn(0);
00217 }
00218
00230 #undef __FUNCT__
00231 #define __FUNCT__ "DSDPGetPObjective"
00232 int DSDPGetPObjective(DSDP dsdp,double *pobj){
00233 int info;
00234 double scale;
00235 DSDPFunctionBegin;
00236 DSDPValid(dsdp);
00237 info=DSDPGetScale(dsdp,&scale);DSDPCHKERR(info);
00238 *pobj=(dsdp->pobj)/scale;
00239 DSDPFunctionReturn(0);
00240 }
00241
00252 #undef __FUNCT__
00253 #define __FUNCT__ "DSDPGetSolutionType"
00254 int DSDPGetSolutionType(DSDP dsdp,DSDPSolutionType *pdfeasible){
00255 DSDPFunctionBegin;
00256 DSDPValid(dsdp);
00257 *pdfeasible=dsdp->pdfeasible;
00258 DSDPFunctionReturn(0);
00259 }
00276 #undef __FUNCT__
00277 #define __FUNCT__ "DSDPGetTraceX"
00278 int DSDPGetTraceX(DSDP dsdp, double *tracex){
00279 DSDPFunctionBegin;
00280 DSDPValid(dsdp);
00281 *tracex=dsdp->tracex;
00282 DSDPFunctionReturn(0);
00283 }
00284
00295 #undef __FUNCT__
00296 #define __FUNCT__ "DSDPGetFinalErrors"
00297 int DSDPGetFinalErrors(DSDP dsdp, double err[6]){
00298 int info;
00299 double scale,rr,bnorm,dobj=0,pobj=0;
00300 DSDPFunctionBegin;
00301 DSDPValid(dsdp);
00302 info=DSDPGetScale(dsdp,&scale);DSDPCHKERR(info);
00303 info=DSDPVecGetR(dsdp->y,&rr); DSDPCHKERR(info);
00304 info=DSDPGetPObjective(dsdp,&pobj);DSDPCHKERR(info);
00305 info=DSDPGetDObjective(dsdp,&dobj);DSDPCHKERR(info);
00306 err[0]=dsdp->perror;
00307 err[1]=0;
00308 err[2]=fabs(rr)/scale;
00309 err[3]=0;
00310 err[4]=pobj - dobj;
00311 err[5]=dsdp->tracexs/scale;
00312
00313 err[2] /= (1.0+dsdp->cnorm);
00314 info=DSDPVecCopy(dsdp->b,dsdp->ytemp);DSDPCHKERR(info);
00315 info=DSDPVecSetC(dsdp->ytemp,0);DSDPCHKERR(info);
00316 info=DSDPVecSetR(dsdp->ytemp,0);DSDPCHKERR(info);
00317 info=DSDPVecNormInfinity(dsdp->ytemp,&bnorm);
00318 err[0]=dsdp->perror/(1.0+bnorm);
00319
00320 err[4]=(err[4])/(1.0+fabs(pobj)+fabs(dobj));
00321 err[5]=(err[5])/(1.0+fabs(pobj)+fabs(dobj));
00322 DSDPFunctionReturn(0);
00323 }
00324
00341 #undef __FUNCT__
00342 #define __FUNCT__ "DSDPGetPInfeasibility"
00343 int DSDPGetPInfeasibility(DSDP dsdp, double *pperror){
00344 DSDPFunctionBegin;
00345 DSDPValid(dsdp);
00346 if (pperror) *pperror=dsdp->pinfeas;
00347 DSDPFunctionReturn(0);
00348 }
00349
00363 #undef __FUNCT__
00364 #define __FUNCT__ "DSDPSetPTolerance"
00365 int DSDPSetPTolerance(DSDP dsdp,double inftol){
00366 DSDPFunctionBegin;
00367 DSDPValid(dsdp);
00368 if (inftol > 0) dsdp->pinfeastol = inftol;
00369 DSDPLogInfo(0,2,"Set P Infeasibility Tolerance: %4.4e\n",inftol);
00370 DSDPFunctionReturn(0);
00371 }
00372
00384 #undef __FUNCT__
00385 #define __FUNCT__ "DSDPGetPTolerance"
00386 int DSDPGetPTolerance(DSDP dsdp,double *inftol){
00387 DSDPFunctionBegin;
00388 DSDPValid(dsdp);
00389 if (inftol) *inftol=dsdp->pinfeastol;
00390 DSDPFunctionReturn(0);
00391 }
00392
00393
00407 #undef __FUNCT__
00408 #define __FUNCT__ "DSDPSetRTolerance"
00409 int DSDPSetRTolerance(DSDP dsdp,double inftol){
00410 DSDPFunctionBegin;
00411 DSDPValid(dsdp);
00412 if (inftol > 0) dsdp->dinfeastol = inftol;
00413 DSDPLogInfo(0,2,"Set D Infeasibility Tolerance: %4.4e\n",inftol);
00414 DSDPFunctionReturn(0);
00415 }
00416
00432 #undef __FUNCT__
00433 #define __FUNCT__ "DSDPGetRTolerance"
00434 int DSDPGetRTolerance(DSDP dsdp,double *inftol){
00435 DSDPFunctionBegin;
00436 DSDPValid(dsdp);
00437 *inftol=dsdp->dinfeastol;
00438 DSDPFunctionReturn(0);
00439 }
00440
00453 #undef __FUNCT__
00454 #define __FUNCT__ "DSDPGetYMakeX"
00455 int DSDPGetYMakeX(DSDP dsdp,double y[], int m){
00456 int i,info;
00457 double scale,*yy;
00458 DSDPFunctionBegin;
00459 DSDPValid(dsdp);
00460 if (dsdp->m < m-1) DSDPFunctionReturn(1);
00461 if (dsdp->m > m) DSDPFunctionReturn(1);
00462 info=DSDPVecCopy(dsdp->xmaker[0].y,dsdp->ytemp); DSDPCHKERR(info);
00463 info=DSDPGetScale(dsdp,&scale);DSDPCHKERR(info);
00464 info=DSDPVecGetArray(dsdp->ytemp,&yy);DSDPCHKERR(info);
00465 for (i=0;i<m;i++) y[i]=yy[i+1]/scale;
00466 info=DSDPVecRestoreArray(dsdp->ytemp,&yy);DSDPCHKERR(info);
00467 DSDPFunctionReturn(0);
00468 }
00469
00481 #undef __FUNCT__
00482 #define __FUNCT__ "DSDPGetDYMakeX"
00483 int DSDPGetDYMakeX(DSDP dsdp,double dy[], int m){
00484 int i,info;
00485 double scale,*yy;
00486 DSDPFunctionBegin;
00487 DSDPValid(dsdp);
00488 if (dsdp->m < m-1) DSDPFunctionReturn(1);
00489 if (dsdp->m > m) DSDPFunctionReturn(1);
00490 info=DSDPVecCopy(dsdp->xmaker[0].dy,dsdp->ytemp); DSDPCHKERR(info);
00491 info=DSDPGetScale(dsdp,&scale);DSDPCHKERR(info);
00492 info=DSDPVecGetArray(dsdp->ytemp,&yy);DSDPCHKERR(info);
00493 for (i=0;i<m;i++) dy[i]=yy[i+1]/scale;
00494 info=DSDPVecRestoreArray(dsdp->ytemp,&yy);DSDPCHKERR(info);
00495 DSDPFunctionReturn(0);
00496 }
00497
00509 #undef __FUNCT__
00510 #define __FUNCT__ "DSDPGetMuMakeX"
00511 int DSDPGetMuMakeX(DSDP dsdp,double *mu){
00512 int info;
00513 double scale;
00514 DSDPFunctionBegin;
00515 DSDPValid(dsdp);
00516 info=DSDPGetScale(dsdp,&scale);DSDPCHKERR(info);
00517 *mu=dsdp->xmaker[0].mu/scale;
00518 DSDPFunctionReturn(0);
00519 }
00520