1/* $Id: CoinPresolveEmpty.cpp 1448 2011-06-19 15:34:41Z stefan $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others. All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include <stdio.h>
7#include <math.h>
8
9#include "CoinFinite.hpp"
10#include "CoinHelperFunctions.hpp"
11#include "CoinPresolveMatrix.hpp"
12
13#include "CoinPresolveEmpty.hpp"
14#include "CoinMessage.hpp"
15
16#if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
17#include "CoinPresolvePsdebug.hpp"
18#endif
19
20
21/* \file
22
23 Routines to remove/reinsert empty columns and rows.
24*/
25
26
27/*
28 Physically remove empty columns, compressing mcstrt and hincol. The major
29 side effect is that columns are renumbered, thus clink_ is no longer valid
30 and must be rebuilt.
31
32 It's necessary to rebuild clink_ in order to do direct conversion of a
33 CoinPresolveMatrix to a CoinPostsolveMatrix by transferring the data arrays.
34 Without clink_, it's impractical to build link_ to match the transferred bulk
35 storage.
36*/
37const CoinPresolveAction
38 *drop_empty_cols_action::presolve (CoinPresolveMatrix *prob,
39 int *ecols,
40 int necols,
41 const CoinPresolveAction *next)
42{
43 int ncols = prob->ncols_;
44 CoinBigIndex *mcstrt = prob->mcstrt_;
45 int *hincol = prob->hincol_;
46
47 presolvehlink *clink = prob->clink_;
48
49 double *clo = prob->clo_;
50 double *cup = prob->cup_;
51 double *dcost = prob->cost_;
52
53 const double ztoldj = prob->ztoldj_;
54
55 unsigned char * integerType = prob->integerType_;
56 int * originalColumn = prob->originalColumn_;
57
58 const double maxmin = prob->maxmin_;
59
60 double * sol = prob->sol_;
61 unsigned char * colstat = prob->colstat_;
62
63 action *actions = new action[necols];
64 int * colmapping = new int [ncols+1];
65 bool fixInfeasibility = (prob->presolveOptions_&16384)!=0;
66
67 CoinZeroN(colmapping,ncols);
68 int i;
69 for (i=necols-1; i>=0; i--) {
70 int jcol = ecols[i];
71 colmapping[jcol]=-1;
72 action &e = actions[i];
73
74 e.jcol = jcol;
75 e.clo = clo[jcol];
76 e.cup = cup[jcol];
77 // adjust if integer
78 if (integerType[jcol]) {
79 e.clo = ceil(e.clo-1.0e-9);
80 e.cup = floor(e.cup+1.0e-9);
81 if (e.clo>e.cup&&!fixInfeasibility) {
82 // infeasible
83 prob->status_|= 1;
84 prob->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
85 prob->messages())
86 <<jcol
87 <<e.clo
88 <<e.cup
89 <<CoinMessageEol;
90 }
91 }
92 e.cost = dcost[jcol];
93
94 // there are no more constraints on this variable,
95 // so we had better be able to compute the answer now
96 if (fabs(dcost[jcol])<ztoldj)
97 dcost[jcol]=0.0;
98 if (dcost[jcol] * maxmin == 0.0) {
99 // hopefully, we can make this non-basic
100 // what does OSL currently do in this case???
101 e.sol = (-PRESOLVE_INF < e.clo
102 ? e.clo
103 : e.cup < PRESOLVE_INF
104 ? e.cup
105 : 0.0);
106 } else if (dcost[jcol] * maxmin > 0.0) {
107 if (-PRESOLVE_INF < e.clo)
108 e.sol = e.clo;
109 else {
110 prob->messageHandler()->message(COIN_PRESOLVE_COLUMNBOUNDB,
111 prob->messages())
112 <<jcol
113 <<CoinMessageEol;
114 prob->status_ |= 2;
115 break;
116 }
117 } else {
118 if (e.cup < PRESOLVE_INF)
119 e.sol = e.cup;
120 else {
121 prob->messageHandler()->message(COIN_PRESOLVE_COLUMNBOUNDA,
122 prob->messages())
123 <<jcol
124 <<CoinMessageEol;
125 prob->status_ |= 2;
126 break;
127 }
128 }
129
130#if PRESOLVE_DEBUG
131 if (e.sol * dcost[jcol]) {
132 //printf("\a>>>NON-ZERO COST FOR EMPTY COL %d: %g\n", jcol, dcost[jcol]);
133 }
134#endif
135 prob->change_bias(e.sol * dcost[jcol]);
136
137
138 }
139 int ncols2=0;
140
141 // now move remaining ones down
142 for (i=0;i<ncols;i++) {
143 if (!colmapping[i]) {
144 mcstrt[ncols2] = mcstrt[i];
145 hincol[ncols2] = hincol[i];
146
147 clo[ncols2] = clo[i];
148 cup[ncols2] = cup[i];
149
150 dcost[ncols2] = dcost[i];
151 if (sol) {
152 sol[ncols2] = sol[i];
153 colstat[ncols2] = colstat[i];
154 }
155
156 integerType[ncols2] = integerType[i];
157 originalColumn[ncols2] = originalColumn[i];
158 colmapping[i] = ncols2++;
159 }
160 }
161 mcstrt[ncols2] = mcstrt[ncols] ;
162 colmapping[ncols] = ncols2 ;
163/*
164 Rebuild clink_. At this point, all empty columns are linked out, so the
165 only columns left are columns that are to be saved, hence available in
166 colmapping. All we need to do is walk clink_ and write the new entries
167 into a new array.
168*/
169
170 { presolvehlink *newclink = new presolvehlink [ncols2+1] ;
171 for (int oldj = ncols ; oldj >= 0 ; oldj = clink[oldj].pre)
172 { presolvehlink &oldlnk = clink[oldj] ;
173 int newj = colmapping[oldj] ;
174 assert(newj >= 0 && newj <= ncols2) ;
175 presolvehlink &newlnk = newclink[newj] ;
176 if (oldlnk.suc >= 0)
177 { newlnk.suc = colmapping[oldlnk.suc] ; }
178 else
179 { newlnk.suc = NO_LINK ; }
180 if (oldlnk.pre >= 0)
181 { newlnk.pre = colmapping[oldlnk.pre] ; }
182 else
183 { newlnk.pre = NO_LINK ; } }
184 delete [] clink ;
185 prob->clink_ = newclink ; }
186
187 delete [] colmapping;
188 prob->ncols_ = ncols2;
189
190 return (new drop_empty_cols_action(necols, actions, next));
191}
192
193
194const CoinPresolveAction *drop_empty_cols_action::presolve(CoinPresolveMatrix *prob,
195 const CoinPresolveAction *next)
196{
197 const int *hincol = prob->hincol_;
198 // const double *clo = prob->clo_;
199 // const double *cup = prob->cup_;
200 int ncols = prob->ncols_;
201 int i;
202 int nempty = 0;
203 int * empty = new int [ncols];
204 CoinBigIndex nelems2 = 0 ;
205
206 // count empty cols
207 for (i=0; i<ncols; i++)
208 { nelems2 += hincol[i] ;
209 if (hincol[i] == 0) {
210#if PRESOLVE_DEBUG
211 if (nempty==0)
212 printf("UNUSED COLS: ");
213 else
214 if (i < 100 && nempty%25 == 0)
215 printf("\n") ;
216 else
217 if (i >= 100 && i < 1000 && nempty%19 == 0)
218 printf("\n") ;
219 else
220 if (i >= 1000 && nempty%15 == 0)
221 printf("\n") ;
222 printf("%d ", i);
223#endif
224 empty[nempty++] = i;
225 }
226 }
227 prob->nelems_ = nelems2 ;
228
229 if (nempty) {
230#if PRESOLVE_DEBUG
231 printf("\ndropped %d cols\n", nempty);
232#endif
233 next = drop_empty_cols_action::presolve(prob, empty, nempty, next);
234 }
235 delete [] empty;
236 return (next);
237}
238
239
240void drop_empty_cols_action::postsolve(CoinPostsolveMatrix *prob) const
241{
242 const int nactions = nactions_;
243 const action *const actions = actions_;
244
245 int ncols = prob->ncols_;
246
247 CoinBigIndex *mcstrt = prob->mcstrt_;
248 int *hincol = prob->hincol_;
249 // int *hrow = prob->hrow_;
250
251 double *clo = prob->clo_;
252 double *cup = prob->cup_;
253
254 double *sol = prob->sol_;
255 double *cost = prob->cost_;
256 double *rcosts = prob->rcosts_;
257 unsigned char *colstat = prob->colstat_;
258 const double maxmin = prob->maxmin_;
259
260 int ncols2 = ncols+nactions;
261 int * colmapping = new int [ncols2];
262
263 CoinZeroN(colmapping,ncols2);
264# if PRESOLVE_DEBUG
265 char *cdone = prob->cdone_;
266# endif
267 int action_i;
268 for (action_i = 0; action_i < nactions; action_i++) {
269 const action *e = &actions[action_i];
270 int jcol = e->jcol;
271 colmapping[jcol]=-1;
272 }
273
274 int i;
275
276 // now move remaining ones up
277 for (i=ncols2-1;i>=0;i--) {
278 if (!colmapping[i]) {
279 ncols--;
280 mcstrt[i] = mcstrt[ncols];
281 hincol[i] = hincol[ncols];
282
283 clo[i] = clo[ncols];
284 cup[i] = cup[ncols];
285
286 cost[i] = cost[ncols];
287
288 if (sol)
289 sol[i] = sol[ncols];
290
291 if (rcosts)
292 rcosts[i] = rcosts[ncols];
293
294 if (colstat)
295 colstat[i] = colstat[ncols];
296# if PRESOLVE_DEBUG
297 cdone[i] = cdone[ncols];
298# endif
299 }
300 }
301 assert (!ncols);
302
303 delete [] colmapping;
304
305 for (action_i = 0; action_i < nactions; action_i++) {
306 const action *e = &actions[action_i];
307 int jcol = e->jcol;
308
309 // now recreate jcol
310 clo[jcol] = e->clo;
311 cup[jcol] = e->cup;
312 if (sol)
313 sol[jcol] = e->sol;
314 cost[jcol] = e->cost;
315
316 if (rcosts)
317 rcosts[jcol] = maxmin*cost[jcol];
318
319 hincol[jcol] = 0;
320 mcstrt[jcol] = NO_LINK ;
321# if PRESOLVE_DEBUG
322 cdone[jcol] = DROP_COL;
323# endif
324 if (colstat)
325 prob->setColumnStatusUsingValue(jcol);
326 }
327
328 prob->ncols_ += nactions;
329
330# if PRESOLVE_CONSISTENCY
331 presolve_check_threads(prob) ;
332# endif
333
334}
335
336
337
338const CoinPresolveAction *drop_empty_rows_action::presolve(CoinPresolveMatrix *prob,
339 const CoinPresolveAction *next)
340{
341 int ncols = prob->ncols_;
342 CoinBigIndex *mcstrt = prob->mcstrt_;
343 int *hincol = prob->hincol_;
344 int *hrow = prob->hrow_;
345
346 int nrows = prob->nrows_;
347 // This is done after row copy needed
348 //int *mrstrt = prob->mrstrt_;
349 int *hinrow = prob->hinrow_;
350 //int *hcol = prob->hcol_;
351
352 double *rlo = prob->rlo_;
353 double *rup = prob->rup_;
354
355 unsigned char *rowstat = prob->rowstat_;
356 double *acts = prob->acts_;
357 int * originalRow = prob->originalRow_;
358
359 //presolvehlink *rlink = prob->rlink_;
360 bool fixInfeasibility = (prob->presolveOptions_&16384)!=0;
361 // Relax tolerance
362 double tolerance = 10.0*prob->feasibilityTolerance_;
363
364
365 int i;
366 int nactions = 0;
367 for (i=0; i<nrows; i++)
368 if (hinrow[i] == 0)
369 nactions++;
370
371 if (nactions == 0)
372 return next;
373 else {
374 action *actions = new action[nactions];
375 int * rowmapping = new int [nrows];
376
377 nactions = 0;
378 int nrows2=0;
379 for (i=0; i<nrows; i++) {
380 if (hinrow[i] == 0) {
381 action &e = actions[nactions];
382
383# if PRESOLVE_DEBUG
384 if (nactions == 0)
385 printf("UNUSED ROWS: ");
386 else
387 if (i < 100 && nactions%25 == 0)
388 printf("\n") ;
389 else
390 if (i >= 100 && i < 1000 && nactions%19 == 0)
391 printf("\n") ;
392 else
393 if (i >= 1000 && nactions%15 == 0)
394 printf("\n") ;
395 printf("%d ", i);
396# endif
397
398 nactions++;
399 if (rlo[i] > 0.0 || rup[i] < 0.0) {
400 if ((rlo[i]<=tolerance &&
401 rup[i]>=-tolerance)||fixInfeasibility) {
402 rlo[i]=0.0;
403 rup[i]=0.0;
404 } else {
405 prob->status_|= 1;
406 prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
407 prob->messages())
408 <<i
409 <<rlo[i]
410 <<rup[i]
411 <<CoinMessageEol;
412 break;
413 }
414 }
415 e.row = i;
416 e.rlo = rlo[i];
417 e.rup = rup[i];
418 rowmapping[i]=-1;
419
420 } else {
421 // move down - we want to preserve order
422 rlo[nrows2]=rlo[i];
423 rup[nrows2]=rup[i];
424 originalRow[nrows2]=i;
425 if (acts) {
426 acts[nrows2]=acts[i];
427 rowstat[nrows2]=rowstat[i];
428 }
429 rowmapping[i]=nrows2++;
430 }
431 }
432
433 // remap matrix
434 for (i=0;i<ncols;i++) {
435 int j;
436 for (j=mcstrt[i];j<mcstrt[i]+hincol[i];j++)
437 hrow[j] = rowmapping[hrow[j]];
438 }
439 delete [] rowmapping;
440
441 prob->nrows_ = nrows2;
442
443#if PRESOLVE_DEBUG
444 presolve_check_nbasic(prob) ;
445 if (nactions)
446 printf("\ndropped %d rows\n", nactions);
447#endif
448 return (new drop_empty_rows_action(nactions, actions, next));
449 }
450}
451
452void drop_empty_rows_action::postsolve(CoinPostsolveMatrix *prob) const
453{
454 const int nactions = nactions_;
455 const action *const actions = actions_;
456
457 int ncols = prob->ncols_;
458 CoinBigIndex *mcstrt = prob->mcstrt_;
459 int *hincol = prob->hincol_;
460
461 int *hrow = prob->hrow_;
462
463 double *rlo = prob->rlo_;
464 double *rup = prob->rup_;
465 unsigned char *rowstat = prob->rowstat_;
466 double *rowduals = prob->rowduals_;
467 double *acts = prob->acts_;
468# if PRESOLVE_DEBUG
469 char *rdone = prob->rdone_;
470# endif
471
472 int nrows0 = prob->nrows0_;
473 int nrows = prob->nrows_;
474
475 int * rowmapping = new int [nrows0];
476 CoinZeroN(rowmapping,nrows0) ;
477
478 int i, action_i;
479 for (action_i = 0; action_i<nactions; action_i++) {
480 const action *e = &actions[action_i];
481 int hole = e->row;
482 rowmapping[hole]=-1;
483 }
484
485 // move data
486 for (i=nrows0-1; i>=0; i--) {
487 if (!rowmapping[i]) {
488 // not a hole
489 nrows--;
490 rlo[i]=rlo[nrows];
491 rup[i]=rup[nrows];
492 acts[i]=acts[nrows];
493 rowduals[i]=rowduals[nrows];
494 if (rowstat)
495 rowstat[i] = rowstat[nrows];
496# if PRESOLVE_DEBUG
497 rdone[i] = rdone[nrows] ;
498# endif
499 }
500 }
501 assert (!nrows);
502 // set up mapping for matrix
503 for (i=0;i<nrows0;i++) {
504 if (!rowmapping[i])
505 rowmapping[nrows++]=i;
506 }
507
508 for (int j=0; j<ncols; j++) {
509 CoinBigIndex start = mcstrt[j];
510 CoinBigIndex end = start + hincol[j];
511
512 for (CoinBigIndex k=start; k<end; ++k) {
513 hrow[k] = rowmapping[hrow[k]];
514 }
515 }
516
517
518 delete [] rowmapping;
519
520 for (action_i = 0; action_i < nactions; action_i++) {
521 const action *e = &actions[action_i];
522 int irow = e->row;
523
524 // Now recreate irow
525 rlo[irow] = e->rlo;
526 rup[irow] = e->rup;
527
528 if (rowstat)
529 prob->setRowStatus(irow,CoinPrePostsolveMatrix::basic);
530 rowduals[irow] = 0.0; // ???
531 acts[irow] = 0.0;
532
533# if PRESOLVE_DEBUG
534 rdone[irow] = DROP_ROW;
535# endif
536 }
537
538 prob->nrows_ = prob->nrows_+nactions;
539
540# if PRESOLVE_DEBUG
541 presolve_check_threads(prob) ;
542# endif
543
544}
545