00001 #include "simplex.h"
00002 #include "error.h"
00003 #include "geom.h"
00004 #include "defines.h"
00005
00006 #include <float.h>
00007 #include <math.h>
00008 #include <string.h>
00009
00022 void SimplexCreate(double epsilon,unsigned int ncols,TSimplex *s)
00023 {
00024 double ct;
00025 unsigned int i;
00026
00027
00028 s->lp=Clp_newModel();
00029
00030
00031 Clp_resize(s->lp,0,ncols);
00032
00033 {
00034
00035 double *v;
00036 unsigned int *e;
00037
00038 double rl[1]={-INF};
00039 double ru[1]={ INF};
00040 int rs[2]={0,ncols};
00041
00042 NEW(v,ncols,double);
00043 NEW(e,ncols,unsigned int);
00044 for(i=0;i<ncols;i++)
00045 {
00046 e[i]=i;
00047 v[i]=0;
00048 }
00049
00050 Clp_addRows(s->lp,1,rl,ru,rs,(signed int *)e,v);
00051
00052 free(v);
00053 free(e);
00054
00055 s->fakeRows=1;
00056 }
00057
00058 NEW(s->lower,ncols,double);
00059 NEW(s->upper,ncols,double);
00060 NEW(s->obj,ncols,double);
00061 for(i=0;i<ncols;i++)
00062 {
00063 s->lower[i]=-INF;
00064 s->upper[i]=+INF;
00065 s->obj[i]=0;
00066 }
00067
00068
00069 #if (_DEBUG<5)
00070 Clp_setLogLevel(s->lp,0);
00071 #endif
00072
00073
00074 Clp_scaling(s->lp,0);
00075
00076 ct=Clp_primalTolerance(s->lp);
00077 if (ct>(epsilon*1e-3))
00078 Clp_setPrimalTolerance(s->lp,epsilon*1e-3);
00079
00080 ct=Clp_dualTolerance(s->lp);
00081 if (ct>(epsilon*1e-3))
00082 Clp_setDualTolerance(s->lp,epsilon*1e-3);
00083
00084 ct=Clp_getSmallElementValue(s->lp);
00085 if (ct>epsilon)
00086 Clp_setSmallElementValue(s->lp,epsilon);
00087
00088
00089 Clp_setOptimizationDirection(s->lp,1);
00090
00091
00092 Clp_setMaximumSeconds(s->lp,ncols*SIMPLEX_TIMEOUT);
00093 }
00094
00095
00096 void ResetSimplex(TSimplex *s)
00097 {
00098
00099
00100 }
00101
00102 unsigned int SimplexNColumns(TSimplex *s)
00103 {
00104 return((unsigned int)Clp_numberColumns(s->lp));
00105 }
00106
00107 unsigned int SimplexNRows(TSimplex *s)
00108 {
00109 return((unsigned int)Clp_numberRows(s->lp)-s->fakeRows);
00110 }
00111
00112
00113 void SimplexSetColBounds(unsigned int ncol,Tinterval *i,TSimplex *s)
00114 {
00115 s->lower[ncol]=LowerLimit(i);
00116 s->upper[ncol]=UpperLimit(i);
00117
00118 Clp_chgColumnLower(s->lp,s->lower);
00119 Clp_chgColumnUpper(s->lp,s->upper);
00120 }
00121
00122
00123 void SimplexGetColBounds(unsigned int ncol,Tinterval *i,TSimplex *s)
00124 {
00125 double *colUpper,*colLower;
00126
00127 colLower=Clp_columnLower(s->lp);
00128 colUpper=Clp_columnUpper(s->lp);
00129
00130 NewInterval(colLower[ncol],colUpper[ncol],i);
00131 }
00132
00133 void SimplexGetColConstraint(unsigned int ncol,TLinearConstraint *lc,TSimplex *s)
00134 {
00135 Tinterval error;
00136 const double *elements;
00137 const int *rowIndices;
00138 const int *columnLengths;
00139 const CoinBigIndex *columnStarts;
00140 CoinBigIndex a,b,i;
00141
00142 rowIndices=Clp_getIndices(s->lp);
00143 elements=Clp_getElements(s->lp);
00144
00145 columnLengths=Clp_getVectorLengths(s->lp);
00146 columnStarts=Clp_getVectorStarts(s->lp);
00147
00148 a=columnStarts[ncol];
00149 b=a+columnLengths[ncol];
00150
00151 InitLinearConstraint(lc);
00152 for(i=a;i<b;++i)
00153 AddTerm2LinearConstraint((unsigned int)rowIndices[i],elements[i],lc);
00154
00155 SimplexGetColBounds(ncol,&error,s);
00156 SetLinearConstraintError(&error,lc);
00157 }
00158
00159 boolean SimplexColEmpty(unsigned int ncol,TSimplex *s)
00160 {
00161 const int *columnLengths;
00162
00163 columnLengths=Clp_getVectorLengths(s->lp);
00164
00165 return(columnLengths[ncol]==0);
00166 }
00167
00168 double SimplexGetColPrimal(unsigned int ncol,TSimplex *s)
00169 {
00170 double *obj;
00171
00172 obj=Clp_primalColumnSolution(s->lp);
00173
00174 return(obj[ncol]);
00175 }
00176
00177 double SimplexGetColDual(unsigned int ncol,TSimplex *s)
00178 {
00179 double *obj;
00180
00181 obj=Clp_dualColumnSolution(s->lp);
00182
00183 return(obj[ncol]);
00184 }
00185
00186
00187 void SimplexSetRowBounds(unsigned int nrow,Tinterval *i,TSimplex *s)
00188 {
00189 Error("SimplexSetRowBounds is not implemented with CLP");
00190
00191 #if (0)
00192 double *rowUpper,*rowLower;
00193
00194 rowLower=Clp_rowLower(s->lp);
00195 rowUpper=Clp_rowUpper(s->lp);
00196
00197 rowLower[nrow+s->fakeRows]=LowerLimit(i);
00198 rowUpper[nrow+s->fakeRows]=UpperLimit(i);
00199 #endif
00200 }
00201
00202
00203 void SimplexGetRowBounds(unsigned int nrow,Tinterval *i,TSimplex *s)
00204 {
00205 double *rowUpper,*rowLower;
00206
00207 rowLower=Clp_rowLower(s->lp);
00208 rowUpper=Clp_rowUpper(s->lp);
00209
00210 NewInterval(rowLower[nrow+s->fakeRows],rowUpper[nrow+s->fakeRows],i);
00211 }
00212
00213 void SimplexGetRowConstraint(unsigned int nrow,TLinearConstraint *lc,TSimplex *s)
00214 {
00215
00216 Tinterval error;
00217 const double *elements;
00218 const int *rowIndices;
00219 const int *columnLengths;
00220 const CoinBigIndex *columnStarts;
00221 CoinBigIndex a,b,i;
00222 unsigned int j,m,n;
00223
00224 m=SimplexNColumns(s);
00225
00226 rowIndices=Clp_getIndices(s->lp);
00227 elements=Clp_getElements(s->lp);
00228
00229 columnLengths=Clp_getVectorLengths(s->lp);
00230 columnStarts=Clp_getVectorStarts(s->lp);
00231
00232 n=nrow+s->fakeRows;
00233
00234 InitLinearConstraint(lc);
00235 for(j=0;j<m;j++)
00236 {
00237 a=columnStarts[j];
00238 b=a+columnLengths[j];
00239 for(i=a;i<b;++i)
00240 {
00241 if (rowIndices[i]==n)
00242 AddTerm2LinearConstraint(j,elements[i],lc);
00243 }
00244 }
00245
00246 SimplexGetRowBounds(nrow,&error,s);
00247 SetLinearConstraintError(&error,lc);
00248 }
00249
00250
00251 double SimplexGetRowPrimal(unsigned int nrow,TSimplex *s)
00252 {
00253 double *obj;
00254
00255 obj=Clp_primalRowSolution(s->lp);
00256
00257 return(obj[nrow+s->fakeRows]);
00258 }
00259
00260 double SimplexGetRowDual(unsigned int nrow,TSimplex *s)
00261 {
00262 double *obj;
00263
00264 obj=Clp_dualRowSolution(s->lp);
00265
00266 return(obj[nrow+s->fakeRows]);
00267 }
00268
00269 void SimplexAddNewConstraintRaw(TLinearConstraint *lc,TSimplex *s)
00270 {
00271 Tinterval bounds;
00272 double rl[1];
00273 double ru[1];
00274 int rs[2];
00275
00276 GetLinearConstraintError(&bounds,lc);
00277
00278 rl[0]=LowerLimit(&bounds);
00279 ru[0]=UpperLimit(&bounds);
00280
00281 rs[0]=0;
00282 rs[1]=GetNumTermsInLinearConstraint(lc);
00283
00284 Clp_addRows(s->lp,1,rl,ru,rs,
00285 (signed int *)GetLinearConstraintVariables(lc),
00286 GetLinearConstraintCoefficients(lc));
00287
00288 #if (_DEBUG>4)
00289 printf(" Setting Row %u with range ",SimplexNRows(s));
00290 PrintInterval(stdout,&bounds);
00291 printf(" [%g %g]\n ",rl[0],ru[0]);
00292 PrintLinearConstraint(stdout,TRUE,NULL,lc);
00293 #endif
00294 }
00295
00296 void SimplexSetOptimizationFunction(TLinearConstraint *obj,TSimplex *s)
00297 {
00298 unsigned int i,n;
00299
00300 n=SimplexNColumns(s);
00301
00302 for(i=0;i<n;i++)
00303 s->obj[i]=0.0;
00304
00305 n=GetNumTermsInLinearConstraint(obj);
00306 for(i=0;i<n;i++)
00307 {
00308 #if (_DEBUG>4)
00309 printf("+%f*v[%u]",
00310 GetLinearConstraintCoefficient(i,obj),
00311 GetLinearConstraintVariable(i,obj));
00312 #endif
00313
00314 s->obj[GetLinearConstraintVariable(i,obj)]=GetLinearConstraintCoefficient(i,obj);
00315 }
00316 #if (_DEBUG>4)
00317 printf("\n");
00318 #endif
00319 Clp_chgObjCoefficients(s->lp,s->obj);
00320 }
00321
00322 void SimplexGetOptimizationFunction(TLinearConstraint *obj,TSimplex *s)
00323 {
00324 unsigned int i,n;
00325 double *obj_c;
00326
00327 obj_c=Clp_objective(s->lp);
00328
00329 InitLinearConstraint(obj);
00330
00331 n=SimplexNColumns(s);
00332
00333 for(i=0;i<n;i++)
00334 AddTerm2LinearConstraint(i,obj_c[i],obj);
00335 }
00336
00337 double SimplexGetOptimalValueRaw(TSimplex *s)
00338 {
00339 return(Clp_objectiveValue(s->lp));
00340 }
00341
00342 unsigned int SimplexSolve(TSimplex *s)
00343 {
00344 int status;
00345 unsigned int e;
00346
00347 if (s->fakeRows>0)
00348 {
00349 const int w[1]={0};
00350
00351 Clp_deleteRows(s->lp,1,w);
00352
00353 s->fakeRows=0;
00354 }
00355
00356 Clp_primal(s->lp,1);
00357
00358 status=Clp_status(s->lp);
00359
00360 switch(status)
00361 {
00362 case 0:
00363 e=REDUCED_BOX;
00364 break;
00365 case 1:
00366 e=EMPTY_BOX;
00367 break;
00368 case 2:
00369 e=UNBOUNDED_BOX;
00370 break;
00371 default:
00372 e=ERROR_IN_PROCESS;
00373 }
00374 return(e);
00375 }
00376
00377 void SimplexDelete(TSimplex *s)
00378 {
00379 if (SimplexNRows(s)>0)
00380 Clp_deleteModel(s->lp);
00381
00382 free(s->lower);
00383 free(s->upper);
00384 free(s->obj);
00385 }
00386
00387