dsdpx.c

Go to the documentation of this file.
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 

Generated on Sun Mar 23 07:30:49 2008 for DSDP by  doxygen 1.5.5