DSDP
|
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 /* rpinfeas is infeasibility of (P) while rpinfeas2 is infeasibility of (PP) */ 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){ /* (PP) must be close to feasible */ 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(bigM)-dsdp->tracex < fabs(rrr) && rpinfeas<pfeastol */ 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 /* Step direction was not accurate enough to compute X from Schur complement */ 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 /* Not good enough to save */ 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){ /* Estimate error in X */ 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