1 /* $Id$
2 *
3 * Name: CouenneLPtightenBoundsCLP.cpp
4 * Authors: Pietro Belotti, Carnegie Mellon University
5 * Purpose: adaptation of OsiClpSolverInterface::tightenBounds()
6 *
7 * (C) Carnegie-Mellon University, 2009.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11 #include "CouennePrecisions.hpp"
12 #include "CouenneProblem.hpp"
13 #include "CouenneCutGenerator.hpp"
14 #include "CouenneExprVar.hpp"
15
16 namespace Couenne {
17
18 //#define COIN_DEVELOP 4
19
20 // Tighten bounds. Returns -1 if infeasible, otherwise number of
21 // variables tightened.
22 template <class T>
tightenBoundsCLP(int lightweight)23 int CouenneSolverInterface<T>::tightenBoundsCLP (int lightweight) {
24
25 // Copied from OsiClpSolverInterface::tightenBounds
26
27 int
28 numberRows = T::getNumRows(),
29 numberColumns = T::getNumCols(),
30 iRow, iColumn;
31
32 const double * columnUpper = T::getColUpper();
33 const double * columnLower = T::getColLower();
34 const double * rowUpper = T::getRowUpper();
35 const double * rowLower = T::getRowLower();
36
37 // Column copy of matrix
38 const double * element = T::getMatrixByCol()->getElements();
39 const int * row = T::getMatrixByCol()->getIndices();
40 const CoinBigIndex * columnStart = T::getMatrixByCol()->getVectorStarts();
41 const int * columnLength = T::getMatrixByCol()->getVectorLengths();
42 const double *objective = T::getObjCoefficients() ;
43
44 double direction = T::getObjSense();
45 double * down = new double [numberRows];
46
47 if (lightweight)
48 return tightenBoundsCLP_Light (lightweight);
49
50 // NOT LIGHTWEIGHT /////////////////////////////////////////////////////
51
52 double * up = new double [numberRows];
53 double * sum = new double [numberRows];
54 int * type = new int [numberRows];
55 CoinZeroN(down,numberRows);
56 CoinZeroN(up,numberRows);
57 CoinZeroN(sum,numberRows);
58 CoinZeroN(type,numberRows);
59 double infinity = T::getInfinity();
60
61 for (iColumn=0;iColumn<numberColumns;iColumn++) {
62 CoinBigIndex start = columnStart[iColumn];
63 CoinBigIndex end = start + columnLength[iColumn];
64 double lower = columnLower[iColumn];
65 double upper = columnUpper[iColumn];
66 if (lower==upper) {
67 for (CoinBigIndex j=start;j<end;j++) {
68 int iRow = row[j];
69 double value = element[j];
70 sum[iRow]+=2.0*fabs(value*lower);
71 if ((type[iRow]&1)==0)
72 down[iRow] += value*lower;
73 if ((type[iRow]&2)==0)
74 up[iRow] += value*lower;
75 }
76 } else {
77 for (CoinBigIndex j=start;j<end;j++) {
78 int iRow = row[j];
79 double value = element[j];
80 if (value>0.0) {
81 if ((type[iRow]&1)==0) {
82 if (lower!=-infinity) {
83 down[iRow] += value*lower;
84 sum[iRow]+=fabs(value*lower);
85 } else {
86 type[iRow] |= 1;
87 }
88 }
89 if ((type[iRow]&2)==0) {
90 if (upper!=infinity) {
91 up[iRow] += value*upper;
92 sum[iRow]+=fabs(value*upper);
93 } else {
94 type[iRow] |= 2;
95 }
96 }
97 } else {
98 if ((type[iRow]&1)==0) {
99 if (upper!=infinity) {
100 down[iRow] += value*upper;
101 sum[iRow]+=fabs(value*upper);
102 } else {
103 type[iRow] |= 1;
104 }
105 }
106 if ((type[iRow]&2)==0) {
107 if (lower!=-infinity) {
108 up[iRow] += value*lower;
109 sum[iRow]+=fabs(value*lower);
110 } else {
111 type[iRow] |= 2;
112 }
113 }
114 }
115 }
116 }
117 }
118
119 int nTightened = 0;
120 double tolerance = 1.0e-6;
121
122 for (iRow=0;iRow<numberRows;iRow++) {
123 if ((type[iRow]&1)!=0)
124 down[iRow]=-infinity;
125 if (down[iRow]>rowUpper[iRow]) {
126 if (down[iRow]>rowUpper[iRow]+tolerance+1.0e-8*sum[iRow]) {
127 // infeasible
128 #ifdef COIN_DEVELOP
129 printf("infeasible on row %d\n",iRow);
130 #endif
131 nTightened=-1;
132 break;
133 } else {
134 down[iRow]=rowUpper[iRow];
135 }
136 }
137 if ((type[iRow]&2)!=0)
138 up[iRow]=infinity;
139 if (up[iRow]<rowLower[iRow]) {
140 if (up[iRow]<rowLower[iRow]-tolerance-1.0e-8*sum[iRow]) {
141 // infeasible
142 #ifdef COIN_DEVELOP
143 printf("infeasible on row %d\n",iRow);
144 #endif
145 nTightened=-1;
146 break;
147 } else {
148 up[iRow]=rowLower[iRow];
149 }
150 }
151 }
152
153 if (nTightened)
154 numberColumns = 0; // so will skip
155
156 for (iColumn=0;iColumn<numberColumns;iColumn++) {
157 double lower = columnLower[iColumn];
158 double upper = columnUpper[iColumn];
159 double gap = upper-lower;
160
161 if (!gap)
162 continue;
163
164 int canGo=0;
165
166 CoinBigIndex
167 start = columnStart [iColumn],
168 end = start + columnLength [iColumn];
169
170 if (lower < -1.0e8 && upper > 1.0e8)
171 continue; // Could do severe damage to accuracy
172
173
174 // there was an ifInteger condition here. We do like tightened
175 // bounds for continuous variables too, so we don't test for
176 // integrality.
177
178 std::vector <exprVar *> &vars = cutgen_ -> Problem () -> Variables ();
179
180 {
181 if (vars [iColumn] -> isInteger ()) {
182
183 if (lower < ceil (lower - COUENNE_EPS) - COUENNE_EPS) {
184 #ifdef COIN_DEVELOP
185 printf("increasing lower bound on %d from %e to %e\n",iColumn,
186 lower,ceil(lower - COUENNE_EPS));
187 #endif
188 lower=ceil(lower - COUENNE_EPS);
189 gap=upper-lower;
190 T::setColLower(iColumn,lower);
191 }
192
193 if (upper > floor(upper + COUENNE_EPS) + COUENNE_EPS) {
194 #ifdef COIN_DEVELOP
195 printf("decreasing upper bound on %d from %e to %e\n",iColumn,
196 upper,floor(upper + COUENNE_EPS));
197 #endif
198 upper=floor(upper + COUENNE_EPS);
199 gap=upper-lower;
200 T::setColUpper(iColumn,upper);
201 }
202 }
203
204 double newLower=lower;
205 double newUpper=upper;
206
207 for (CoinBigIndex j=start;j<end;j++) {
208 int iRow = row[j];
209 double value = element[j];
210 if (value>0.0) {
211 if ((type[iRow]&1)==0) {
212 // has to be at most something
213 if (down[iRow] + value*gap > rowUpper[iRow]+tolerance) {
214 double newGap = (rowUpper[iRow]-down[iRow])/value;
215 // adjust
216 newGap += 1.0e-10*sum[iRow];
217 if (vars [iColumn] -> isInteger ())
218 newGap = floor(newGap);
219 if (lower+newGap<newUpper)
220 newUpper=lower+newGap;
221 }
222 }
223 if (down[iRow]<rowLower[iRow])
224 canGo |=1; // can't go down without affecting result
225 if ((type[iRow]&2)==0) {
226 // has to be at least something
227 if (up[iRow] - value*gap < rowLower[iRow]-tolerance) {
228 double newGap = (up[iRow]-rowLower[iRow])/value;
229 // adjust
230 newGap += 1.0e-10*sum[iRow];
231 if (vars [iColumn] -> isInteger ())
232 newGap = floor(newGap);
233 if (upper-newGap>newLower)
234 newLower=upper-newGap;
235 }
236 }
237 if (up[iRow]>rowUpper[iRow])
238 canGo |=2; // can't go up without affecting result
239 } else {
240 if ((type[iRow]&1)==0) {
241 // has to be at least something
242 if (down[iRow] - value*gap > rowUpper[iRow]+tolerance) {
243 double newGap = -(rowUpper[iRow]-down[iRow])/value;
244 // adjust
245 newGap += 1.0e-10*sum[iRow];
246 if (vars [iColumn] -> isInteger ())
247 newGap = floor(newGap);
248 if (upper-newGap>newLower)
249 newLower=upper-newGap;
250 }
251 }
252 if (up[iRow]>rowUpper[iRow])
253 canGo |=1; // can't go down without affecting result
254 if ((type[iRow]&2)==0) {
255 // has to be at most something
256 if (up[iRow] + value*gap < rowLower[iRow]-tolerance) {
257 double newGap = -(up[iRow]-rowLower[iRow])/value;
258 // adjust
259 newGap += 1.0e-10*sum[iRow];
260 if (vars [iColumn] -> isInteger ())
261 newGap = floor(newGap);
262 if (lower+newGap<newUpper)
263 newUpper=lower+newGap;
264 }
265 }
266 if (down[iRow]<rowLower[iRow])
267 canGo |=2; // can't go up without affecting result
268 }
269 }
270
271 if (newUpper<upper || newLower>lower) {
272 nTightened++;
273 if (newLower>newUpper) {
274 // infeasible
275 #if COIN_DEVELOP>1
276 printf("infeasible on column %d\n",iColumn);
277 #endif
278 nTightened=-1;
279 break;
280 } else {
281 T::setColLower(iColumn,newLower);
282 T::setColUpper(iColumn,newUpper);
283 }
284 for (CoinBigIndex j=start;j<end;j++) {
285 int iRow = row[j];
286 double value = element[j];
287 if (value>0.0) {
288 if ((type[iRow]&1)==0) down [iRow] += value*(newLower-lower);
289 if ((type[iRow]&2)==0) up [iRow] += value*(newUpper-upper);
290 } else {
291 if ((type[iRow]&1)==0) down [iRow] += value*(newUpper-upper);
292 if ((type[iRow]&2)==0) up [iRow] += value*(newLower-lower);
293 }
294 }
295 } else {
296
297 if (canGo!=3) {
298
299 double objValue = direction*objective[iColumn];
300
301 if (objValue>=0.0&&(canGo&1)==0) {
302 #if COIN_DEVELOP>2
303 printf("dual fix down on column %d\n",iColumn);
304 #endif
305 nTightened++;
306 T::setColUpper(iColumn,lower);
307 } else if (objValue<=0.0 && (canGo&2)==0) {
308 #if COIN_DEVELOP>2
309 printf("dual fix up on column %d\n",iColumn);
310 #endif
311 nTightened++;
312 T::setColLower(iColumn,upper);
313 }
314 }
315 }
316 }
317
318 // else {
319
320 // // CONTINUOUS //////////////////////////////////////////
321
322 // // just do dual tests
323 // for (CoinBigIndex j=start;j<end;j++) {
324 // int iRow = row[j];
325 // double value = element[j];
326 // if (value>0.0) {
327 // if (down [iRow] < rowLower [iRow]) canGo |=1; // can't go down without affecting result
328 // if (up [iRow] > rowUpper [iRow]) canGo |=2; // can't go up without affecting result
329 // } else {
330 // if (up [iRow] > rowUpper [iRow]) canGo |=1; // can't go down without affecting result
331 // if (down [iRow] < rowLower [iRow]) canGo |=2; // can't go up without affecting result
332 // }
333 // }
334
335 // if (canGo!=3) {
336 // double objValue = direction*objective[iColumn];
337 // if (objValue>=0.0&&(canGo&1)==0) {
338 // #if COIN_DEVELOP>2
339 // printf("dual fix down on continuous column %d lower %g\n",
340 // iColumn,lower);
341 // #endif
342 // // Only if won't cause numerical problems
343 // if (lower>-1.0e10) {
344 // nTightened++;;
345 // setColUpper(iColumn,lower);
346 // }
347 // } else if (objValue<=0.0&&(canGo&2)==0) {
348 // #if COIN_DEVELOP>2
349 // printf("dual fix up on continuous column %d upper %g\n",
350 // iColumn,upper);
351 // #endif
352 // // Only if won't cause numerical problems
353 // if (upper<1.0e10) {
354 // nTightened++;;
355 // setColLower(iColumn,upper);
356 // }
357 // }
358 // }
359 // }
360
361 }
362
363 delete [] type;
364 delete [] down;
365 delete [] up;
366 delete [] sum;
367
368 return nTightened;
369 }
370
371 }
372