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
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
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
00066 info = DSDPSchurMatInitialize(&dsdp->M);DSDPCHKERR(info);
00067 info = DSDPSetDefaultSchurMatrixStructure(dsdp); DSDPCHKERR(info);
00068 info = DSDPCGInitialize(&dsdp->sles); DSDPCHKERR(info);
00069
00070
00071
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
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
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
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
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 }