DSDP
|
00001 #include "dsdp.h" 00002 #include "dsdpsys.h" 00003 #include "dsdp5.h" 00028 #undef __FUNCT__ 00029 #define __FUNCT__ "DSDPCreate" 00030 int DSDPCreate(int m,DSDP* dsdpnew){ 00031 00032 DSDP dsdp; 00033 int info; 00034 00035 DSDPFunctionBegin; 00036 00037 DSDPCALLOC1(&dsdp,PD_DSDP,&info);DSDPCHKERR(info); 00038 *dsdpnew=dsdp; 00039 dsdp->keyid=DSDPKEY; 00040 00041 /* Initialize some parameters */ 00042 DSDPEventLogInitialize(); 00043 dsdp->m=m; 00044 dsdp->maxcones=0; 00045 dsdp->ncones=0; 00046 dsdp->K=0; 00047 dsdp->setupcalled=DSDP_FALSE; 00048 dsdp->ybcone=0; 00049 dsdp->ndroutines=0; 00050 /* info = DSDPSetStandardMonitor(dsdp);DSDPCHKERR(info); */ 00051 info = DSDPVecCreateSeq(m+2,&dsdp->b);DSDPCHKERR(info); 00052 info = DSDPVecZero(dsdp->b);DSDPCHKERR(info); 00053 info = DSDPVecDuplicate(dsdp->b,&dsdp->y);DSDPCHKERR(info); 00054 info = DSDPVecDuplicate(dsdp->b,&dsdp->ytemp);DSDPCHKERR(info); 00055 info = DSDPVecZero(dsdp->y);DSDPCHKERR(info); 00056 info = DSDPVecSetC(dsdp->y,-1.0);DSDPCHKERR(info); 00057 00058 info = DSDPAddRCone(dsdp,&dsdp->rcone);DSDPCHKERR(info); 00059 info = DSDPCreateLUBoundsCone(dsdp,&dsdp->ybcone);DSDPCHKERR(info); 00060 00061 info=DSDPSetDefaultStatistics(dsdp);DSDPCHKERR(info); 00062 info=DSDPSetDefaultParameters(dsdp);DSDPCHKERR(info); 00063 info=DSDPSetDefaultMonitors(dsdp);DSDPCHKERR(info); 00064 00065 /* info = DSDPMatInitialize(m,m,&dsdp->Q);DSDPCHKERR(info); */ 00066 info = DSDPSchurMatInitialize(&dsdp->M);DSDPCHKERR(info); 00067 info = DSDPSetDefaultSchurMatrixStructure(dsdp); DSDPCHKERR(info); 00068 info = DSDPCGInitialize(&dsdp->sles); DSDPCHKERR(info); 00069 00070 /* Set the one global variable 00071 sdat=dsdp; 00072 */ 00073 DSDPFunctionReturn(0); 00074 } 00075 00076 00077 #undef __FUNCT__ 00078 #define __FUNCT__ "DSDPSetDefaultStatistics" 00079 00084 int DSDPSetDefaultStatistics(DSDP dsdp){ 00085 00086 int i; 00087 DSDPFunctionBegin; 00088 DSDPValid(dsdp); 00089 dsdp->reason=CONTINUE_ITERATING; 00090 dsdp->pdfeasible=DSDP_PDUNKNOWN; 00091 dsdp->itnow=0; 00092 dsdp->pobj= 1.0e10; 00093 dsdp->ppobj= 1.0e10; 00094 dsdp->dobj= -1.0e+9; 00095 dsdp->ddobj= -1.0e+9; 00096 dsdp->dualitygap=dsdp->ppobj-dsdp->ddobj; 00097 dsdp->pstep=1.0; 00098 dsdp->dstep=0.0; 00099 for (i=0;i<MAX_XMAKERS;i++){ 00100 dsdp->xmaker[i].mu=1.0e200; 00101 dsdp->xmaker[i].pstep=0.0; 00102 } 00103 dsdp->pnorm=0.001; 00104 dsdp->mu=1000.0; 00105 dsdp->np=0; 00106 dsdp->anorm=0; 00107 dsdp->bnorm=0; 00108 dsdp->cnorm=0; 00109 dsdp->tracex=0; 00110 dsdp->tracexs=0; 00111 dsdp->Mshift=0; 00112 dsdp->goty0=DSDP_FALSE; 00113 DSDPFunctionReturn(0); 00114 } 00115 #undef __FUNCT__ 00116 #define __FUNCT__ "DSDPSetDefaultParameters" 00117 00122 int DSDPSetDefaultParameters(DSDP dsdp){ 00123 00124 int info; 00125 DSDPFunctionBegin; 00126 DSDPValid(dsdp); 00127 00128 /* Stopping parameters */ 00129 info=DSDPSetMaxIts(dsdp,500);DSDPCHKERR(info); 00130 info=DSDPSetGapTolerance(dsdp,1.0e-6);DSDPCHKERR(info); 00131 info=DSDPSetPNormTolerance(dsdp,1.0e30);DSDPCHKERR(info); 00132 if (dsdp->m<100){info=DSDPSetGapTolerance(dsdp,1.0e-7);DSDPCHKERR(info);} 00133 if (dsdp->m>3000){info=DSDPSetGapTolerance(dsdp,5.0e-6);DSDPCHKERR(info);} 00134 info=RConeSetType(dsdp->rcone,DSDPInfeasible);DSDPCHKERR(info); 00135 info=DSDPSetDualBound(dsdp,1.0e20);DSDPCHKERR(info); 00136 info=DSDPSetStepTolerance(dsdp,5.0e-2);DSDPCHKERR(info); 00137 info=DSDPSetRTolerance(dsdp,1.0e-6);DSDPCHKERR(info); 00138 info=DSDPSetPTolerance(dsdp,1.0e-4);DSDPCHKERR(info); 00139 /* Solver options */ 00140 info=DSDPSetMaxTrustRadius(dsdp,1.0e10);DSDPCHKERR(info); 00141 info=DSDPUsePenalty(dsdp,0);DSDPCHKERR(info); 00142 info=DSDPSetInitialBarrierParameter(dsdp,-1.0);DSDPCHKERR(info); 00143 info=DSDPSetPotentialParameter(dsdp,3.0);DSDPCHKERR(info); 00144 info=DSDPUseDynamicRho(dsdp,1);DSDPCHKERR(info); 00145 info=DSDPSetR0(dsdp,-1.0);DSDPCHKERR(info); 00146 info=DSDPSetPenaltyParameter(dsdp,1.0e8);DSDPCHKERR(info); 00147 info=DSDPReuseMatrix(dsdp,4);DSDPCHKERR(info); 00148 if (dsdp->m>100){info=DSDPReuseMatrix(dsdp,7);DSDPCHKERR(info);} 00149 if (dsdp->m>1000){info=DSDPReuseMatrix(dsdp,10);DSDPCHKERR(info);} 00150 if (dsdp->m<=100){info=DSDPSetPotentialParameter(dsdp,5.0);DSDPCHKERR(info);DSDPCHKERR(info);} 00151 dsdp->maxschurshift=1.0e-11; 00152 dsdp->mu0=-1.0; 00153 dsdp->slestype=2; 00154 info = DSDPSetYBounds(dsdp,-1e7,1e7);DSDPCHKERR(info); 00155 DSDPFunctionReturn(0); 00156 } 00157 00158 #undef __FUNCT__ 00159 #define __FUNCT__ "DSDPSetDefaultMonitors" 00160 00165 int DSDPSetDefaultMonitors(DSDP dsdp){ 00166 00167 int info; 00168 00169 DSDPFunctionBegin; 00170 DSDPValid(dsdp); 00171 dsdp->nmonitors=0; 00172 info=DSDPSetMonitor(dsdp,DSDPDefaultConvergence,(void*)&dsdp->conv); DSDPCHKERR(info); 00173 DSDPFunctionReturn(0); 00174 } 00175 00191 #undef __FUNCT__ 00192 #define __FUNCT__ "DSDPSetUp" 00193 int DSDPSetup(DSDP dsdp){ 00194 00195 int i,info; 00196 DSDPFunctionBegin; 00197 DSDPValid(dsdp); 00198 00199 /* Create the Work Vectors */ 00200 info = DSDPVecDuplicate(dsdp->y,&dsdp->rhs1);DSDPCHKERR(info); 00201 info = DSDPVecDuplicate(dsdp->y,&dsdp->rhs2);DSDPCHKERR(info); 00202 info = DSDPVecDuplicate(dsdp->y,&dsdp->rhs);DSDPCHKERR(info); 00203 info = DSDPVecDuplicate(dsdp->y,&dsdp->rhstemp);DSDPCHKERR(info); 00204 info = DSDPVecDuplicate(dsdp->y,&dsdp->dy1);DSDPCHKERR(info); 00205 info = DSDPVecDuplicate(dsdp->y,&dsdp->dy2);DSDPCHKERR(info); 00206 info = DSDPVecDuplicate(dsdp->y,&dsdp->dy);DSDPCHKERR(info); 00207 info = DSDPVecDuplicate(dsdp->y,&dsdp->y0);DSDPCHKERR(info); 00208 info = DSDPVecDuplicate(dsdp->y,&dsdp->xmakerrhs);DSDPCHKERR(info); 00209 for (i=0;i<MAX_XMAKERS;i++){ 00210 info = DSDPVecDuplicate(dsdp->y,&dsdp->xmaker[i].y);DSDPCHKERR(info); 00211 info = DSDPVecDuplicate(dsdp->y,&dsdp->xmaker[i].dy);DSDPCHKERR(info); 00212 info = DSDPVecDuplicate(dsdp->y,&dsdp->xmaker[i].rhs);DSDPCHKERR(info); 00213 } 00214 00215 /* Create M */ 00216 info = DSDPSetUpCones(dsdp);DSDPCHKERR(info); 00217 info = DSDPSchurMatSetup(dsdp->M,dsdp->ytemp);DSDPCHKERR(info); 00218 00219 info = DSDPCGSetup(dsdp->sles,dsdp->ytemp); DSDPCHKERR(info); 00220 00221 info = DSDPSetUpCones2(dsdp,dsdp->y,dsdp->M);DSDPCHKERR(info); 00222 info = DSDPGetConicDimension(dsdp,&dsdp->np);DSDPCHKERR(info); 00223 00224 info=DSDPComputeDataNorms(dsdp);DSDPCHKERR(info); 00225 dsdp->pinfeas=dsdp->bnorm+1; 00226 dsdp->perror=dsdp->bnorm+1; 00227 info=DSDPScaleData(dsdp);DSDPCHKERR(info); 00228 00229 info=DSDPGetConicDimension(dsdp,&dsdp->np);DSDPCHKERR(info); 00230 dsdp->solvetime=0; 00231 dsdp->cgtime=0; 00232 dsdp->ptime=0; 00233 dsdp->dtime=0; 00234 dsdp->ctime=0; 00235 info=DSDPEventLogRegister("Primal Step",&dsdp->ptime); 00236 info=DSDPEventLogRegister("Dual Step",&dsdp->dtime); 00237 info=DSDPEventLogRegister("Corrector Step",&dsdp->ctime); 00238 info=DSDPEventLogRegister("CG Solve",&dsdp->cgtime); 00239 info=DSDPEventLogRegister("DSDP Solve",&dsdp->solvetime); 00240 dsdp->setupcalled=DSDP_TRUE; 00241 DSDPFunctionReturn(0); 00242 } 00243 00244 00245 00246 #undef __FUNCT__ 00247 #define __FUNCT__ "DSDPGetSchurMatrix" 00248 int DSDPGetSchurMatrix(DSDP dsdp, DSDPSchurMat *M){ 00249 DSDPFunctionBegin; 00250 DSDPValid(dsdp); 00251 *M=dsdp->M; 00252 DSDPFunctionReturn(0); 00253 } 00254 00255 #undef __FUNCT__ 00256 #define __FUNCT__ "DSDPGetConvergenceMonitor" 00257 00268 int DSDPGetConvergenceMonitor(DSDP dsdp, ConvergenceMonitor**ctx){ 00269 DSDPFunctionBegin; 00270 DSDPValid(dsdp); 00271 *ctx=&dsdp->conv; 00272 DSDPFunctionReturn(0); 00273 } 00274 00275 00276 #undef __FUNCT__ 00277 #define __FUNCT__ "DSDPComputeDataNorms" 00278 00283 int DSDPComputeDataNorms(DSDP dsdp){ 00284 int info; 00285 DSDPVec ytemp=dsdp->ytemp; 00286 DSDPFunctionBegin; 00287 DSDPValid(dsdp); 00288 info = DSDPComputeANorm2(dsdp,ytemp);DSDPCHKERR(info); 00289 info = DSDPFixedVariablesNorm(dsdp->M,ytemp);DSDPCHKERR(info); 00290 info = DSDPVecGetC(ytemp,&dsdp->cnorm);DSDPCHKERR(info); 00291 dsdp->cnorm=sqrt(dsdp->cnorm); 00292 info = DSDPVecSetR(ytemp,0);DSDPCHKERR(info); 00293 info = DSDPVecSetC(ytemp,0);DSDPCHKERR(info); 00294 info = DSDPVecNorm1(ytemp,&dsdp->anorm);DSDPCHKERR(info); 00295 dsdp->anorm=sqrt(dsdp->anorm); 00296 DSDPLogInfo(0,2,"Norm of data: %4.2e\n",dsdp->anorm); 00297 info=DSDPVecCopy(dsdp->b,ytemp);DSDPCHKERR(info); 00298 info = DSDPVecSetR(ytemp,0);DSDPCHKERR(info); 00299 info = DSDPVecSetC(ytemp,0);DSDPCHKERR(info); 00300 info = DSDPVecNorm2(ytemp,&dsdp->bnorm);DSDPCHKERR(info); 00301 DSDPFunctionReturn(0); 00302 } 00303 00304 #undef __FUNCT__ 00305 #define __FUNCT__ "DSDPScaleData" 00306 00311 int DSDPScaleData(DSDP dsdp){ 00312 int info; 00313 double scale; 00314 DSDPFunctionBegin; 00315 DSDPValid(dsdp); 00316 scale=1.0*dsdp->anorm; 00317 if (dsdp->bnorm){ scale/=dsdp->bnorm;} 00318 if (dsdp->cnorm){ scale/=dsdp->cnorm;} 00319 scale=DSDPMin(scale,1.0); 00320 scale=DSDPMax(scale,1.0e-6); 00321 if (dsdp->cnorm==0){ scale=1;} 00322 info=DSDPSetScale(dsdp,scale);DSDPCHKERR(info); 00323 DSDPFunctionReturn(0); 00324 } 00325 00341 #undef __FUNCT__ 00342 #define __FUNCT__ "DSDPSolve" 00343 int DSDPSolve(DSDP dsdp){ 00344 int info; 00345 DSDPFunctionBegin; 00346 info=DSDPEventLogBegin(dsdp->solvetime); 00347 dsdp->pdfeasible=DSDP_PDUNKNOWN; 00348 info=DSDPSetConvergenceFlag(dsdp,CONTINUE_ITERATING);DSDPCHKERR(info); 00349 info=DSDPInitializeVariables(dsdp);DSDPCHKERR(info); 00350 info=DSDPSolveDynamicRho(dsdp);DSDPCHKERR(info); 00351 if (dsdp->pstep==1){info=DSDPRefineStepDirection(dsdp,dsdp->xmakerrhs,dsdp->xmaker[0].dy);DSDPCHKERR(info);} 00352 if (dsdp->pdfeasible==DSDP_PDUNKNOWN) dsdp->pdfeasible=DSDP_PDFEASIBLE; 00353 info=DSDPEventLogEnd(dsdp->solvetime); 00354 DSDPFunctionReturn(0); 00355 } 00356 00357 00358 #undef __FUNCT__ 00359 #define __FUNCT__ "DSDPCallMonitors" 00360 00367 int DSDPCallMonitors(DSDP dsdp,DMonitor dmonitor[], int ndmonitors){ 00368 int i,info; 00369 DSDPFunctionBegin; 00370 for (i=0; i<ndmonitors;i++){ 00371 info=(dmonitor[i].monitor)(dsdp,dmonitor[i].monitorctx); DSDPCHKERR(info); 00372 } 00373 DSDPFunctionReturn(0); 00374 } 00375 /* ---------------------------------------------------------- */ 00376 #undef __FUNCT__ 00377 #define __FUNCT__ "DSDPCheckConvergence" 00378 00384 int DSDPCheckConvergence(DSDP dsdp,DSDPTerminationReason *reason){ 00385 int info; 00386 DSDPTruth unbounded; 00387 00388 DSDPFunctionBegin; 00389 info = DSDPGetConicDimension(dsdp,&dsdp->np);DSDPCHKERR(info); 00390 dsdp->rgap=(dsdp->ppobj-dsdp->ddobj)/(1.0+fabs(dsdp->ppobj)+fabs(dsdp->ddobj)); 00391 dsdp->pstepold=dsdp->pstep; 00392 if (dsdp->reason==CONTINUE_ITERATING){ 00393 if (dsdp->itnow>0){ 00394 info=DSDPCheckForUnboundedObjective(dsdp,&unbounded);DSDPCHKERR(info); 00395 if (unbounded==DSDP_TRUE){ 00396 dsdp->pdfeasible=DSDP_UNBOUNDED; 00397 info=DSDPSetConvergenceFlag(dsdp,DSDP_CONVERGED); DSDPCHKERR(info); 00398 } 00399 } 00400 if (dsdp->reason==CONTINUE_ITERATING){ 00401 if (dsdp->muold<dsdp->mutarget && dsdp->pstep==1 && dsdp->dstep==1 && dsdp->rgap<1e-5){ 00402 info=DSDPSetConvergenceFlag(dsdp,DSDP_NUMERICAL_ERROR); DSDPCHKERR(info); 00403 DSDPLogInfo(0,2,"DSDP Finished: Numerical issues: Increase in Barrier function. \n");} 00404 if (dsdp->itnow >= dsdp->maxiter){ 00405 info=DSDPSetConvergenceFlag(dsdp,DSDP_MAX_IT); DSDPCHKERR(info);} 00406 if (dsdp->Mshift>dsdp->maxschurshift){ 00407 info = DSDPSetConvergenceFlag(dsdp,DSDP_INDEFINITE_SCHUR_MATRIX); DSDPCHKERR(info); 00408 } 00409 } 00410 info=DSDPCallMonitors(dsdp,dsdp->dmonitor,dsdp->nmonitors);DSDPCHKERR(info); 00411 info=DSDPMonitorCones(dsdp,0); DSDPCHKERR(info); 00412 } 00413 dsdp->muold=dsdp->mutarget; 00414 info = DSDPStopReason(dsdp,reason); DSDPCHKERR(info); 00415 DSDPFunctionReturn(0); 00416 } 00417 00418 00419 00420 /* ---------------------------------------------------------- */ 00421 #undef __FUNCT__ 00422 #define __FUNCT__ "DSDPTakeDown" 00423 00428 int DSDPTakeDown(DSDP dsdp){ 00429 00430 int i,info; 00431 00432 DSDPFunctionBegin; 00433 DSDPValid(dsdp); 00434 info = DSDPVecDestroy(&dsdp->rhs);DSDPCHKERR(info); 00435 info = DSDPVecDestroy(&dsdp->rhs1);DSDPCHKERR(info); 00436 info = DSDPVecDestroy(&dsdp->rhs2);DSDPCHKERR(info); 00437 info = DSDPVecDestroy(&dsdp->rhstemp);DSDPCHKERR(info); 00438 info = DSDPVecDestroy(&dsdp->y);DSDPCHKERR(info); 00439 info = DSDPVecDestroy(&dsdp->ytemp);DSDPCHKERR(info); 00440 info = DSDPVecDestroy(&dsdp->dy1);DSDPCHKERR(info); 00441 info = DSDPVecDestroy(&dsdp->dy2);DSDPCHKERR(info); 00442 info = DSDPVecDestroy(&dsdp->dy);DSDPCHKERR(info); 00443 for (i=0;i<MAX_XMAKERS;i++){ 00444 info = DSDPVecDestroy(&dsdp->xmaker[i].y);DSDPCHKERR(info); 00445 info = DSDPVecDestroy(&dsdp->xmaker[i].dy);DSDPCHKERR(info); 00446 info = DSDPVecDestroy(&dsdp->xmaker[i].rhs);DSDPCHKERR(info); 00447 } 00448 info = DSDPVecDestroy(&dsdp->xmakerrhs);DSDPCHKERR(info); 00449 info = DSDPVecDestroy(&dsdp->y0);DSDPCHKERR(info); 00450 info = DSDPVecDestroy(&dsdp->b);DSDPCHKERR(info); 00451 00452 info = DSDPCGDestroy(&dsdp->sles);DSDPCHKERR(info); 00453 info = DSDPDestroyCones(dsdp);DSDPCHKERR(info); 00454 info = DSDPSchurMatDestroy(&dsdp->M);DSDPCHKERR(info); 00455 info = DSDPGetConicDimension(dsdp,&dsdp->np);DSDPCHKERR(info); 00456 dsdp->setupcalled=DSDP_FALSE; 00457 DSDPFunctionReturn(0); 00458 } 00459 00469 int DSDPSetDestroyRoutine(DSDP dsdp, int (*fd)(void*), void* ctx){ 00470 int nd=dsdp->ndroutines; 00471 if (nd<10){ 00472 dsdp->droutine[nd].f=fd; 00473 dsdp->droutine[nd].ptr=ctx; 00474 dsdp->ndroutines++; 00475 } else { 00476 printf("TOO MANY Destroy routines\n"); 00477 return 1; 00478 } 00479 return 0; 00480 } 00481 00482 00494 #undef __FUNCT__ 00495 #define __FUNCT__ "DSDPDestroy" 00496 int DSDPDestroy(DSDP dsdp){ 00497 int i,info; 00498 DSDPFunctionBegin; 00499 DSDPValid(dsdp); 00500 for (i=0;i<dsdp->ndroutines;i++){ 00501 info=(*dsdp->droutine[i].f)(dsdp->droutine[i].ptr);DSDPCHKERR(info); 00502 } 00503 info=DSDPTakeDown(dsdp);DSDPCHKERR(info); 00504 DSDPFREE(&dsdp,&info);DSDPCHKERR(info); 00505 DSDPFunctionReturn(0); 00506 }