source: trunk/PresolveDoubleton.cpp @ 63

Last change on this file since 63 was 63, checked in by jpfasano, 17 years ago

-modified to use CoinPragma?.h
-modified to compile w/MS Visual C++ and gcc (on cygwin) without errors or warnings.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 29.5 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include <stdio.h>
5#include <math.h>
6
7#include "CoinHelperFunctions.hpp"
8#include "PresolveMatrix.hpp"
9
10#include "PresolveEmpty.hpp"    // for DROP_COL/DROP_ROW
11#include "PresolveZeros.hpp"
12#include "PresolveFixed.hpp"
13#include "PresolveDoubleton.hpp"
14
15#include "PresolvePsdebug.hpp"
16#include "ClpMessage.hpp"
17
18void compact_rep(double *elems, int *indices, CoinBigIndex *starts, const int *lengths, int n,
19                        const presolvehlink *link)
20{
21#if     PRESOLVE_SUMMARY
22  printf("****COMPACTING****\n");
23#endif
24
25  // for now, just look for the first element of the list
26  int i = n;
27  while (link[i].pre != NO_LINK)
28    i = link[i].pre;
29
30  int j = 0;
31  for (; i != n; i = link[i].suc) {
32    CoinBigIndex s = starts[i];
33    CoinBigIndex e = starts[i] + lengths[i];
34
35    // because of the way link is organized, j <= s
36    starts[i] = j;
37    for (CoinBigIndex k = s; k < e; k++) {
38      elems[j] = elems[k];
39      indices[j] = indices[k];
40      j++;
41   }
42  }
43}
44
45// returns true if ran out of memory
46static bool expand_col(CoinBigIndex *mcstrt, 
47                       double *colels,
48                       int *hrow,
49                       int *hincol,
50                       presolvehlink *clink, int ncols,
51
52                       int icolx)
53{
54  CoinBigIndex kcsx = mcstrt[icolx];
55  CoinBigIndex kcex = kcsx + hincol[icolx];
56
57  const int maxk = mcstrt[ncols];       // (22)
58
59  // update col rep - need to expand the column, though.
60  int nextcol = clink[icolx].suc;
61
62  // (22)
63  if (kcex + 1 < mcstrt[nextcol] || nextcol == ncols) {
64    if (! (kcex + 1 < mcstrt[nextcol])) {
65      // nextcol==ncols and no space - must compact
66      compact_rep(colels, hrow, mcstrt, hincol, ncols, clink);
67
68      // update vars
69      kcsx = mcstrt[icolx];
70      kcex = kcsx + hincol[icolx];
71
72      if (! (kcex + 1 < mcstrt[nextcol])) {
73        return (true);
74      }
75    }
76  } else {
77    //printf("EXPAND_COL\n");
78
79    // this is not the last col
80    // fetch last non-empty col (presolve_make_memlists-1)
81    int lastcol = clink[ncols].pre;
82    // (clink[icolx].suc != ncols) ==> (icolx != lastcol)
83
84    // put it directly after the last column
85    int newkcsx = mcstrt[lastcol] + hincol[lastcol];
86
87    if (newkcsx + hincol[icolx] + 1 >= maxk) {
88      compact_rep(colels, hrow, mcstrt, hincol, ncols, clink);
89
90      // update vars
91      kcsx = mcstrt[icolx];
92      kcex = kcsx + hincol[icolx];
93
94      newkcsx = mcstrt[lastcol] + hincol[lastcol];
95
96      if (newkcsx + hincol[icolx] + 1 >= maxk) {
97        return (true);
98      }
99      // have to adjust various induction variables
100      kcsx = mcstrt[icolx];
101      kcex = mcstrt[icolx] + hincol[icolx];
102    }
103
104    // move the column - 1:  copy the entries
105    memcpy((void*)&hrow[newkcsx], (void*)&hrow[kcsx], hincol[icolx] * sizeof(int));
106    memcpy((void*)&colels[newkcsx], (void*)&colels[kcsx], hincol[icolx] * sizeof(double));
107
108    // move the column - 2:  update the memory-order linked list
109    PRESOLVE_REMOVE_LINK(clink, icolx);
110    PRESOLVE_INSERT_LINK(clink, icolx, lastcol);
111
112    // move the column - 3:  update loop variables to maintain invariant
113    mcstrt[icolx] = newkcsx;
114    kcsx = newkcsx;
115    kcex = newkcsx + hincol[icolx];
116  }
117
118  return (false);
119}
120
121#if 0
122static bool reject_doubleton(int *mcstrt,
123                             double *colels,
124                             int *hrow,
125                             int *hincol,
126                             double coeff_factor,
127                             double max_coeff_ratio,
128                             int row0, int icolx, int icoly)
129{
130  CoinBigIndex kcs = mcstrt[icoly];
131  CoinBigIndex kce = kcs + hincol[icoly];
132  CoinBigIndex kcsx = mcstrt[icolx];
133  CoinBigIndex kcex = kcsx + hincol[icolx];
134
135  for (CoinBigIndex kcoly=kcs; kcoly<kce; kcoly++) {
136    int row = hrow[kcoly];
137
138    if (row != row0) {
139      // see if row appears in colx
140      CoinBigIndex kcolx = presolve_find_row1(row, kcsx, kcex, hrow);
141
142      if (kcolx<kcex) {
143        // we will add these two terms
144        // they they are of different magnitudes,
145        // then their difference will be approximately the size of the largest.
146        double orig  = fabs(colels[kcoly] * coeff_factor);
147        double addin = fabs(colels[kcolx]);
148
149        if (max_coeff_ratio * min(orig,addin) < max(orig,addin)) {
150#if     DEBUG_PRESOLVE
151          printf("REJECTED %d %g %g\n", row0, orig, addin);
152#endif
153          return (true);
154        }
155
156        // cancellation is bad, too
157        double newval = fabs((colels[kcoly] * coeff_factor) + colels[kcolx]);
158
159        if (max_coeff_ratio * newval < orig) {
160#if     DEBUG_PRESOLVE
161          printf("REJECTED1 %d %g %g\n", row0,
162                 (colels[kcoly] * coeff_factor),
163                 colels[kcolx]);
164#endif
165          return (true);
166        }
167         
168      }
169    }
170  }
171  return (false);
172}
173#endif
174
175
176/*
177 * Substituting y away:
178 *
179 *       y = (c - a x) / b
180 *
181 * so adjust bounds by:   c/b
182 *           and x  by:  -a/b
183 *
184 * This affects both the row and col representations.
185 *
186 * mcstrt only modified if the column must be moved.
187 *
188 * for every row in icoly
189 *      if icolx is also has an entry for row
190 *              modify the icolx entry for row
191 *              drop the icoly entry from row and modify the icolx entry
192 *      else
193 *              add a new entry to icolx column
194 *              create a new icolx entry
195 *              (this may require moving the column in memory)
196 *              replace icoly entry from row and replace with icolx entry
197 *
198 * The row and column reps are inconsistent during the routine,
199 * because icolx in the column rep is updated, and the entries corresponding
200 * to icolx in the row rep are updated, but nothing concerning icoly
201 * in the col rep is changed.  icoly entries in the row rep are deleted,
202 * and icolx entries in both reps are consistent.
203 * At the end, we set the length of icoly to be zero, so the reps would
204 * be consistent if the row were deleted from the row rep.
205 * Both the row and icoly must be removed from both reps.
206 * In the col rep, icoly will be eliminated entirely, at the end of the routine;
207 * irow occurs in just two columns, one of which (icoly) is eliminated
208 * entirely, the other is icolx, which is not deleted here.
209 * In the row rep, irow will be eliminated entirely, but not here;
210 * icoly is removed from the rows it occurs in.
211 */
212static bool elim_doubleton(const char *msg,
213                           CoinBigIndex *mcstrt, 
214                           double *rlo, double * acts, double *rup,
215                           double *colels,
216                           int *hrow, int *hcol,
217                           int *hinrow, int *hincol,
218                           presolvehlink *clink, int ncols,
219                           CoinBigIndex *mrstrt, double *rowels,
220                           //double a, double b, double c,
221                           double coeff_factor,
222                           double bounds_factor,
223                           int row0, int icolx, int icoly)
224{
225  CoinBigIndex kcs = mcstrt[icoly];
226  CoinBigIndex kce = kcs + hincol[icoly];
227  CoinBigIndex kcsx = mcstrt[icolx];
228  CoinBigIndex kcex = kcsx + hincol[icolx];
229
230  //printf("%s %d x=%d y=%d cf=%g bf=%g nx=%d\n", msg,
231  // row0, icolx, icoly, coeff_factor, bounds_factor, hincol[icolx]);
232#if     DEBUG_PRESOLVE
233  printf("%s %d x=%d y=%d cf=%g bf=%g nx=%d yrows=(", msg,
234         row0, icolx, icoly, coeff_factor, bounds_factor, hincol[icolx]);
235#endif
236  for (CoinBigIndex kcoly=kcs; kcoly<kce; kcoly++) {
237    int row = hrow[kcoly];
238
239    // even though these values are updated, they remain consistent
240    PRESOLVEASSERT(kcex == kcsx + hincol[icolx]);
241
242    // we don't need to update the row being eliminated
243    if (row != row0/* && hinrow[row] > 0*/) {
244      // see if row appears in colx
245      CoinBigIndex kcolx = presolve_find_row1(row, kcsx, kcex, hrow);
246
247      if (bounds_factor != 0.0) {
248        // (1)
249        if (-PRESOLVE_INF < rlo[row])
250          rlo[row] -= colels[kcoly] * bounds_factor;
251
252        // (2)
253        if (rup[row] < PRESOLVE_INF)
254          rup[row] -= colels[kcoly] * bounds_factor;
255
256        // and solution
257        acts[row] -= colels[kcoly] * bounds_factor;
258      }
259
260#if     DEBUG_PRESOLVE
261      printf("%d%s ", row, (kcolx<kcex) ? "+" : "");
262#endif
263
264      if (kcolx<kcex) {
265        // before:  both x and y are in the row
266        // after:   only x is in the row
267        // so: number of elems in col x unchanged, and num elems in row is one less
268
269        // update col rep - just modify coefficent
270        // column y is deleted as a whole at the end of the loop
271        colels[kcolx] += colels[kcoly] * coeff_factor;
272
273        // update row rep
274        // first, copy new value for col x into proper place in rowels
275        CoinBigIndex k2 = presolve_find_row(icolx, mrstrt[row], mrstrt[row]+hinrow[row], hcol);
276        rowels[k2] = colels[kcolx];
277
278        // now delete col y from the row; this changes hinrow[row]
279        presolve_delete_from_row(row, icoly, mrstrt, hinrow, hcol, rowels);
280      } else {
281        // before:  only y is in the row
282        // after:   only x is in the row
283        // so: number of elems in col x is one greater, but num elems in row remains same
284        // update entry corresponding to icolx in row rep
285        // by just overwriting the icoly entry
286        {
287          CoinBigIndex k2 = presolve_find_row(icoly, mrstrt[row], mrstrt[row]+hinrow[row], hcol);
288          hcol[k2] = icolx;
289          rowels[k2] = colels[kcoly] * coeff_factor;
290        }
291
292        {
293          bool no_mem = expand_col(mcstrt, colels, hrow, hincol, clink, ncols,
294                                   icolx);
295          if (no_mem)
296            return (true);
297
298          // have to adjust various induction variables
299          kcoly = mcstrt[icoly] + (kcoly - kcs);
300          kcs = mcstrt[icoly];                  // do this for ease of debugging
301          kce = mcstrt[icoly] + hincol[icoly];
302           
303          kcolx = mcstrt[icolx] + (kcolx - kcs);        // don't really need to do this
304          kcsx = mcstrt[icolx];
305          kcex = mcstrt[icolx] + hincol[icolx];
306        }
307
308        // there is now an unused entry in the memory after the column - use it
309        // mcstrt[ncols] == penultimate index of arrays hrow/colels
310        hrow[kcex] = row;
311        colels[kcex] = colels[kcoly] * coeff_factor;    // y factor is 0.0
312        hincol[icolx]++, kcex++;        // expand the col
313      }
314    }
315  }
316
317#if     DEBUG_PRESOLVE
318  printf(")\n");
319#endif
320
321  // delete the whole column
322  hincol[icoly] = 0;
323
324  return (false);
325}
326
327
328void update_other_rep_quick(int col,
329                            const int *mcstrt, const int *hrow, const double *colels,
330                            const int *hincol,
331
332                            const int *mrstrt, int *hcol, double *rowels, int *hinrow)
333{
334  CoinBigIndex kcs = mcstrt[col];
335  CoinBigIndex kce = kcs + hincol[col];
336
337  for (CoinBigIndex k=kcs; k<kce; ++k) {
338    int row = hrow[k];
339    double coeff = colels[k];
340
341    CoinBigIndex krs = mrstrt[row];
342    CoinBigIndex kre = krs + hinrow[row];
343
344    // the "quick" refers to the assumption that there will be enough room,
345    // and that col does not already occur in the row.
346    hcol[kre] = col;
347    rowels[kre] = coeff;
348    ++hinrow[row];
349  }
350}   
351
352
353
354/*
355 * It is always the case that one of the variables of a doubleton
356 * will be (implied) free, but neither will necessarily be a singleton.
357 * Since in the case of a doubleton the number of non-zero entries
358 * will never increase, though, it makes sense to always eliminate them.
359 *
360 * The col rep and row rep must be consistent.
361 */
362const PresolveAction *doubleton_action::presolve(PresolveMatrix *prob,
363                                                  const PresolveAction *next)
364{
365  double *colels        = prob->colels_;
366  int *hrow             = prob->hrow_;
367  CoinBigIndex *mcstrt          = prob->mcstrt_;
368  int *hincol           = prob->hincol_;
369  int ncols             = prob->ncols_;
370
371  double *clo   = prob->clo_;
372  double *cup   = prob->cup_;
373
374  double *rowels        = prob->rowels_;
375  int *hcol             = prob->hcol_;
376  CoinBigIndex *mrstrt          = prob->mrstrt_;
377  int *hinrow           = prob->hinrow_;
378  int nrows             = prob->nrows_;
379
380  double *rlo   = prob->rlo_;
381  double *rup   = prob->rup_;
382
383  presolvehlink *clink = prob->clink_;
384  presolvehlink *rlink = prob->rlink_;
385
386  const char *integerType = prob->integerType_;
387
388  double *cost  = prob->cost_;
389
390  int numberLook = prob->numberRowsToDo_;
391  int iLook;
392  int * look = prob->rowsToDo_;
393  const double ztolzb   = prob->ztolzb_;
394
395  action * actions = new action [nrows];
396  int nactions = 0;
397
398  int *zeros    = new int[ncols];
399  int nzeros    = 0;
400
401  int *fixed    = new int[ncols];
402  int nfixed    = 0;
403
404  // If rowstat exists then all do
405  unsigned char *rowstat        = prob->rowstat_;
406  double *acts  = prob->acts_;
407  double * sol = prob->sol_;
408  //  unsigned char * colstat = prob->colstat_;
409
410
411#if     CHECK_CONSISTENCY
412  presolve_links_ok(clink, mcstrt, hincol, ncols);
413#endif
414
415  // wasfor (int irow=0; irow<nrows; irow++)
416  for (iLook=0;iLook<numberLook;iLook++) {
417    int irow = look[iLook];
418    if (hinrow[irow] == 2 &&
419        fabs(rup[irow] - rlo[irow]) <= ZTOLDP) {
420      double rhs = rlo[irow];
421      CoinBigIndex krs = mrstrt[irow];
422      CoinBigIndex kre = krs + hinrow[irow];
423      int icolx, icoly;
424      CoinBigIndex k;
425     
426      /* locate first column */
427      for (k=krs; k<kre; k++) {
428        if (hincol[hcol[k]] > 0) {
429          break;
430        }
431      }
432      PRESOLVEASSERT(k<kre);
433      if (fabs(rowels[k]) < ZTOLDP)
434        continue;
435      icolx = hcol[k];
436     
437      /* locate second column */
438      for (k++; k<kre; k++) {
439        if (hincol[hcol[k]] > 0) {
440          break;
441        }
442      }
443      PRESOLVEASSERT(k<kre);
444      if (fabs(rowels[k]) < ZTOLDP)
445        continue;
446      icoly = hcol[k];
447     
448      // don't bother with fixed variables
449      if (!(fabs(cup[icolx] - clo[icolx]) < ZTOLDP) &&
450          !(fabs(cup[icoly] - clo[icoly]) < ZTOLDP)) {
451        double coeffx, coeffy;
452        /* find this row in each of the columns */
453        CoinBigIndex krowx = presolve_find_row(irow, mcstrt[icolx], mcstrt[icolx] + hincol[icolx], hrow);
454        CoinBigIndex krowy = presolve_find_row(irow, mcstrt[icoly], mcstrt[icoly] + hincol[icoly], hrow);
455       
456        /* don't do if both integers for now - unless a variant
457           of x=y and 0-1 variables */
458        // 0 not, 1 x integer, 2 y integer, 3 both (but okay), -1 skip
459        int integerStatus=0;
460        if (integerType[icolx]) {
461          if (integerType[icoly]) {
462            // both integer
463            int good = 0;
464            double rhs2 = rhs;
465            double value;
466            value=colels[krowx];
467            if (value<0.0) {
468              value = - value;
469              rhs2 += 1;
470            }
471            if (cup[icolx]==1.0&&clo[icolx]==0.0&&fabs(value-1.0)<1.0e-7)
472              good =1;
473            value=colels[krowy];
474            if (value<0.0) {
475              value = - value;
476              rhs2 += 1;
477            }
478            if (cup[icoly]==1.0&&clo[icoly]==0.0&&fabs(value-1.0)<1.0e-7)
479              good  |= 2;
480            if (good==3&&fabs(rhs2-1.0)<1.0e-7)
481              integerStatus = 3;
482            else
483              integerStatus=-1;
484          } else {
485            integerStatus = 1;
486          }
487        } else if (integerType[icoly]) {
488          integerStatus = 2;
489        }
490        if (integerStatus<0)
491          continue;
492        if (integerStatus == 2) {
493          swap(icoly,icolx);
494          swap(krowy,krowx);
495        }
496       
497        // HAVE TO JIB WITH ABOVE swapS
498        // if x's coefficient is something like 1000, but y's only something like -1,
499        // then when we postsolve, if x's is close to being out of tolerance,
500        // then y is very likely to be (because y==1000x) . (55)
501        // It it interesting that the number of doubletons found may depend
502        // on which column is substituted away (this is true of baxter.mps).
503        if (!integerStatus) {
504          if (fabs(colels[krowy]) < fabs(colels[krowx])) {
505            swap(icoly,icolx);
506            swap(krowy,krowx);
507          }
508        }
509       
510#if 0
511        //?????
512        if (integerType[icolx] &&
513            clo[icoly] != -PRESOLVE_INF &&
514            cup[icoly] != PRESOLVE_INF) {
515          continue;
516        }
517#endif
518       
519        {
520          CoinBigIndex kcs = mcstrt[icoly];
521          CoinBigIndex kce = kcs + hincol[icoly];
522          for (k=kcs; k<kce; k++) {
523            if (hinrow[hrow[k]] == 1) {
524              break;
525            }
526          }
527          // let singleton rows be taken care of first
528          if (k<kce)
529            continue;
530        }
531       
532        coeffx = colels[krowx];
533        coeffy = colels[krowy];
534       
535        // it is possible that both x and y are singleton columns
536        // that can cause problems
537        if (hincol[icolx] == 1 && hincol[icoly] == 1)
538          continue;
539       
540        // BE CAUTIOUS and avoid very large relative differences
541        // if this is not done in baxter, then the computed solution isn't optimal,
542        // but gets it in 11995 iterations; the postsolve goes to iteration 16181.
543        // with this, the solution is optimal, but takes 18825 iters; postsolve 18871.
544#if 0
545        if (fabs(coeffx) * max_coeff_factor <= fabs(coeffy))
546          continue;
547#endif
548       
549#if 0
550        if (only_zero_rhs && rhs != 0)
551          continue;
552       
553        if (reject_doubleton(mcstrt, colels, hrow, hincol,
554                             -coeffx / coeffy,
555                             max_coeff_ratio,
556                             irow, icolx, icoly))
557          continue;
558#endif
559       
560        // common equations are of the form ax + by = 0, or x + y >= lo
561        {
562          action *s = &actions[nactions];         
563          nactions++;
564         
565          s->row = irow;
566          s->icolx = icolx;
567         
568          s->clox = clo[icolx];
569          s->cupx = cup[icolx];
570          s->costx = cost[icolx];
571         
572          s->icoly = icoly;
573          s->cloy = clo[icoly];
574          s->cupy = cup[icoly];
575          s->costy = cost[icoly];
576         
577          s->rlo = rlo[irow];
578          s->rup = rup[irow];
579         
580          s->coeffx = coeffx;
581         
582          s->coeffy = coeffy;
583         
584          s->ncolx      = hincol[icolx];
585          s->colx       = presolve_duparray(&colels[mcstrt[icolx]], hincol[icolx]);
586          s->indx       = presolve_duparray(&hrow[mcstrt[icolx]], hincol[icolx]);
587         
588          s->ncoly      = hincol[icoly];
589          s->coly       = presolve_duparray(&colels[mcstrt[icoly]], hincol[icoly]);
590          s->indy       = presolve_duparray(&hrow[mcstrt[icoly]], hincol[icoly]);
591        }
592       
593        /*
594         * This moves the bounds information for y onto x,
595         * making y free and allowing us to substitute it away.
596         *
597         * a x + b y = c
598         * l1 <= x <= u1
599         * l2 <= y <= u2        ==>
600         *
601         * l2 <= (c - a x) / b <= u2
602         * b/-a > 0 ==> (b l2 - c) / -a <= x <= (b u2 - c) / -a
603         * b/-a < 0 ==> (b u2 - c) / -a <= x <= (b l2 - c) / -a
604         */
605        {
606          double lo1 = -PRESOLVE_INF;
607          double up1 = PRESOLVE_INF;
608         
609          //PRESOLVEASSERT((coeffx < 0) == (coeffy/-coeffx < 0));
610          // (coeffy/-coeffx < 0) == (coeffy<0 == coeffx<0)
611          if (-PRESOLVE_INF < clo[icoly]) {
612            if ((coeffx < 0) != (coeffy < 0))
613              lo1 = (coeffy * clo[icoly] - rhs) / -coeffx;
614            else 
615              up1 = (coeffy * clo[icoly] - rhs) / -coeffx;
616          }
617         
618          if (cup[icoly] < PRESOLVE_INF) {
619            if ((coeffx < 0) != (coeffy < 0))
620              up1 = (coeffy * cup[icoly] - rhs) / -coeffx;
621            else 
622              lo1 = (coeffy * cup[icoly] - rhs) / -coeffx;
623          }
624         
625          // costy y = costy ((c - a x) / b) = (costy c)/b + x (costy -a)/b
626          // the effect of maxmin cancels out
627          cost[icolx] += cost[icoly] * (-coeffx / coeffy);
628         
629          prob->change_bias(cost[icoly] * rhs / coeffy);
630         
631          if (0    /*integerType[icolx]*/) {
632            abort();
633            /* no change possible for now */
634#if 0
635            lo1 = trunc(lo1);
636            up1 = trunc(up1);
637           
638            /* trunc(3.5) == 3.0 */
639            /* trunc(-3.5) == -3.0 */
640           
641            /* I think this is ok */
642            if (lo1 > clo[icolx]) {
643              (clo[icolx] <= 0.0)
644                clo[icolx] =  ? ilo
645               
646                clo[icolx] = ilo;
647              cup[icolx] = iup;
648            }
649#endif
650          } else {
651            double lo2 = max(clo[icolx], lo1);
652            double up2 = min(cup[icolx], up1);
653            if (lo2 > up2) {
654              if (lo2 <= up2 + prob->feasibilityTolerance_) {
655                // If close to integer then go there
656                double nearest = floor(lo2+0.5);
657                if (fabs(nearest-lo2)<2.0*prob->feasibilityTolerance_) {
658                  lo2 = nearest;
659                  up2 = nearest;
660                } else {
661                  lo2 = up2;
662                }
663              } else {
664                prob->status_ |= 1;
665                prob->messageHandler()->message(CLP_PRESOLVE_COLINFEAS,
666                                                         prob->messages())
667                                                           <<icolx
668                                                           <<lo2
669                                                           <<up2
670                                                           <<CoinMessageEol;
671                break;
672              }
673            }
674            clo[icolx] = lo2;
675            cup[icolx] = up2;
676
677            if (rowstat) {
678              // update solution and basis
679              int basisChoice=0;
680              int numberBasic=0;
681              if (prob->columnIsBasic(icolx))
682                numberBasic++;
683              if (prob->columnIsBasic(icoly))
684                numberBasic++;
685              if (prob->rowIsBasic(irow))
686                numberBasic++;
687              if (sol[icolx]<=lo2+ztolzb) {
688                sol[icolx] =lo2;
689                prob->setColumnStatus(icolx,PrePostsolveMatrix::atLowerBound);
690              } else if (sol[icolx]>=up2-ztolzb) {
691                sol[icolx] =up2;
692                prob->setColumnStatus(icolx,PrePostsolveMatrix::atUpperBound);
693              } else {
694                basisChoice=1;
695              }
696              if (numberBasic>1)
697                prob->setColumnStatus(icolx,PrePostsolveMatrix::basic);
698            }
699            if (lo2 == up2)
700              fixed[nfixed++] = icolx;
701          }
702        }
703       
704        // Update next set of actions
705        {
706          prob->addCol(icolx);
707          int i,kcs,kce;
708          kcs = mcstrt[icoly];
709          kce = kcs + hincol[icoly];
710          for (i=kcs;i<kce;i++) {
711            int row = hrow[i];
712            prob->addRow(row);
713          }
714          kcs = mcstrt[icolx];
715          kce = kcs + hincol[icolx];
716          for (i=kcs;i<kce;i++) {
717            int row = hrow[i];
718            prob->addRow(row);
719          }
720        }
721       
722        /* transfer the colx factors to coly */
723        bool no_mem = elim_doubleton("ELIMD",
724                                     mcstrt, rlo, acts, rup, colels,
725                                     hrow, hcol, hinrow, hincol,
726                                     clink, ncols, 
727                                     mrstrt, rowels,
728                                     -coeffx / coeffy,
729                                     rhs / coeffy,
730                                     irow, icolx, icoly);
731        if (no_mem) 
732          throwCoinError("out of memory",
733                         "doubleton_action::presolve");
734       
735        // now remove irow from icolx in the col rep
736        // better if this were first.
737        presolve_delete_from_row(icolx, irow, mcstrt, hincol, hrow, colels);
738       
739        // eliminate irow entirely from the row rep
740        hinrow[irow] = 0;
741       
742        // eliminate irow entirely from the row rep
743        PRESOLVE_REMOVE_LINK(rlink, irow);
744       
745        // eliminate coly entirely from the col rep
746        PRESOLVE_REMOVE_LINK(clink, icoly);
747        cost[icoly] = 0.0;
748       
749        rlo[irow] = 0.0;
750        rup[irow] = 0.0;
751       
752        zeros[nzeros++] = icolx;        // check for zeros
753       
754        // strictly speaking, since we didn't adjust {clo,cup}[icoly]
755        // or {rlo,rup}[irow], this col/row may be infeasible,
756        // because the solution/activity would be 0, whereas the
757        // bounds may be non-zero.
758      }
759     
760#if 0
761      presolve_links_ok(clink, mcstrt, ncols);
762      presolve_links_ok(rlink, mrstrt, nrows);
763      prob->consistent();
764#endif
765    }
766  }
767
768  if (nactions) {
769#if     PRESOLVE_SUMMARY
770    printf("NDOUBLETONS:  %d\n", nactions);
771#endif
772    action *actions1 = new action[nactions];
773    ClpDisjointCopyN(actions, nactions, actions1);
774
775    next = new doubleton_action(nactions, actions1, next);
776
777    if (nzeros)
778      next = drop_zero_coefficients_action::presolve(prob, zeros, nzeros, next);
779    if (nfixed)
780      next = remove_fixed_action::presolve(prob, fixed, nfixed, next);
781  }
782
783  delete[]zeros;
784  delete[]fixed;
785  deleteAction(actions,action*);
786
787  return (next);
788}
789
790
791void doubleton_action::postsolve(PostsolveMatrix *prob) const
792{
793  const action *const actions = actions_;
794  const int nactions = nactions_;
795
796  double *colels        = prob->colels_;
797  int *hrow             = prob->hrow_;
798  CoinBigIndex *mcstrt          = prob->mcstrt_;
799  int *hincol           = prob->hincol_;
800  int *link             = prob->link_;
801
802  double *clo   = prob->clo_;
803  double *cup   = prob->cup_;
804
805  double *rlo   = prob->rlo_;
806  double *rup   = prob->rup_;
807
808  double *dcost = prob->cost_;
809
810  double *sol   = prob->sol_;
811  double *rcosts        = prob->rcosts_;
812
813  double *acts  = prob->acts_;
814  double *rowduals = prob->rowduals_;
815
816  unsigned char *colstat        = prob->colstat_;
817  unsigned char *rowstat        = prob->rowstat_;
818
819  const double maxmin   = prob->maxmin_;
820
821  char *cdone   = prob->cdone_;
822  char *rdone   = prob->rdone_;
823
824  CoinBigIndex free_list = prob->free_list_;
825
826  const double ztolzb   = prob->ztolzb_;
827  const double ztoldj   = prob->ztoldj_;
828
829  for (const action *f = &actions[nactions-1]; actions<=f; f--) {
830    int irow = f->row;
831    double lo0 = f->clox;
832    double up0 = f->cupx;
833
834    // probably don't need this
835    double ylo0 = f->cloy;
836    double yup0 = f->cupy;
837
838    double coeffx = f->coeffx;
839    double coeffy = f->coeffy;
840    int jcolx = f->icolx;
841    int jcoly = f->icoly;
842
843    // needed?
844    double rhs = f->rlo;
845
846    /* the column was in the reduced problem */
847    PRESOLVEASSERT(cdone[jcolx] && rdone[irow]==DROP_ROW);
848    PRESOLVEASSERT(cdone[jcoly]==DROP_COL);
849
850    // probably don't need this
851    rlo[irow] = f->rlo;
852    rup[irow] = f->rup;
853
854    clo[jcolx] = lo0;
855    cup[jcolx] = up0;
856
857    // probably don't need this
858    clo[jcoly] = ylo0;
859    cup[jcoly] = yup0;
860
861    dcost[jcolx] = f->costx;
862    dcost[jcoly] = f->costy;
863
864#if     DEBUG_PRESOLVE
865    // I've forgotten what this is about
866    if ((rhs < 0) == ((coeffx * sol[jcolx]) < 0) &&
867        fabs(rhs - coeffx * sol[jcolx]) * 100 < rhs &&
868        fabs(rhs - coeffx * sol[jcolx]) * 100 < (coeffx * sol[jcolx]))
869      printf("DANGEROUS RHS??? %g %g %g\n",
870             rhs, coeffx * sol[jcolx],
871             (rhs - coeffx * sol[jcolx]));
872#endif
873    // this is why we want coeffx < coeffy (55)
874    sol[jcoly] = (rhs - coeffx * sol[jcolx]) / coeffy;
875         
876    // since this row is fixed
877    acts[irow] = rhs;
878
879    // acts[irow] always ok, since slack is fixed
880    if (rowstat)
881      prob->setRowStatus(irow,PrePostsolveMatrix::atLowerBound);
882
883
884
885    {
886      // find the tail
887      CoinBigIndex k=mcstrt[jcolx];
888      for (int i=0; i<hincol[jcolx]-1; ++i) {
889        k = link[k];
890      }
891      // append the old x column to the free list, freeing all links
892      link[k] = free_list;
893      free_list = mcstrt[jcolx];
894    }
895
896    // use create_row??
897    {
898      int xstart = NO_LINK;
899      for (int i=0; i<f->ncolx; ++i) {
900        CoinBigIndex k = free_list;
901        free_list = link[free_list];
902
903        check_free_list(free_list);
904
905        hrow[k] = f->indx[i];
906        PRESOLVEASSERT(rdone[hrow[k]] || hrow[k] == irow);
907        colels[k] = f->colx[i];
908        link[k] = xstart;
909        xstart = k;
910      }
911      mcstrt[jcolx] = xstart;
912    }
913    {
914      int ystart = NO_LINK;
915      for (int i=0; i<f->ncoly; ++i) {
916        CoinBigIndex k = free_list;
917        free_list = link[free_list];
918
919        check_free_list(free_list);
920
921        hrow[k] = f->indy[i];
922        PRESOLVEASSERT(rdone[hrow[k]] || hrow[k] == irow);
923        colels[k] = f->coly[i];
924        link[k] = ystart;
925        ystart = k;
926      }
927      mcstrt[jcoly] = ystart;
928    }
929
930    hincol[jcolx] = f->ncolx;
931    hincol[jcoly] = f->ncoly;
932
933
934
935    // CLAIM:
936    // if the new pi value is chosen to keep the reduced cost
937    // of col x at its prior value, then the reduced cost of
938    // col y will be 0.
939         
940    // also have to update row activities and bounds for rows affected by jcoly
941    //
942    // sol[jcolx] was found for coeffx that
943    // was += colels[kcoly] * coeff_factor;
944    // where coeff_factor == -coeffx / coeffy
945    //
946    // its contribution to activity was
947    // (colels[kcolx] + colels[kcoly] * coeff_factor) * sol[jcolx]      (1)
948    //
949    // After adjustment, the two columns contribute:
950    // colels[kcoly] * sol[jcoly] + colels[kcolx] * sol[jcolx]
951    // == colels[kcoly] * ((rhs - coeffx * sol[jcolx]) / coeffy) + colels[kcolx] * sol[jcolx]
952    // == colels[kcoly] * rhs/coeffy + colels[kcoly] * coeff_factor * sol[jcolx] + colels[kcolx] * sol[jcolx]
953    // colels[kcoly] * rhs/coeffy + the expression (1)
954    //
955    // therefore, we must increase the row bounds by colels[kcoly] * rhs/coeffy,
956    // which is similar to the bias
957    {
958      CoinBigIndex k = mcstrt[jcoly];
959      int ny = hincol[jcoly];
960      double dj = maxmin * dcost[jcoly];
961      double bounds_factor = rhs/coeffy;
962
963      // this probably doesn't work (???)
964           
965      for (int i=0; i<ny; ++i) {
966        int row = hrow[k];
967        double coeff = colels[k];
968        k = link[k];
969
970        if (row != irow) {
971          // are these tests always true???
972
973          // undo elim_doubleton(1)
974          if (-PRESOLVE_INF < rlo[row])
975            rlo[row] += coeff * bounds_factor;
976
977          // undo elim_doubleton(2)
978          if (rup[row] < PRESOLVE_INF)
979            rup[row] += coeff * bounds_factor;
980
981          acts[row] += coeff * bounds_factor;
982
983          dj -= rowduals[row] * coeff;
984        }
985      }
986
987      // this is the coefficient we need to force col y's reduced cost to 0.0;
988      // for example, this is obviously true if y is a singleton column
989      rowduals[irow] = dj / coeffy;
990      rcosts[jcoly] = 0.0;
991
992      // furthermore, this should leave rcosts[jcolx] unchanged.
993    }
994
995    // The only problem with keeping the reduced costs the way they were
996    // was that the variable's bound may have moved, requiring it
997    // to become basic.
998    if (colstat) {
999      if (prob->columnIsBasic(jcolx) ||
1000          (fabs(lo0 - sol[jcolx]) < ztolzb && rcosts[jcolx] >= -ztoldj) ||
1001          (fabs(up0 - sol[jcolx]) < ztolzb && rcosts[jcolx] <= ztoldj)) {
1002        // colx is fine as it is - make coly basic
1003       
1004        prob->setColumnStatus(jcoly,PrePostsolveMatrix::basic);
1005        // col y's reduced costs has already been set to 0.0
1006      } else {
1007        prob->setColumnStatus(jcolx,PrePostsolveMatrix::basic);
1008        prob->setColumnStatusUsingValue(jcoly);
1009
1010        if (rcosts[jcolx] != 0.0) {
1011          // change rowduals[jcolx] enough to cancel out rcosts[jcolx]
1012          rowduals[irow] += (rcosts[jcolx] / coeffx);
1013          rcosts[jcoly] -= ((rcosts[jcolx] / coeffx) * coeffy);
1014          rcosts[jcolx] = 0.0;
1015        }
1016      }
1017    }
1018
1019      // DEBUG CHECK
1020#if     DEBUG_PRESOLVE
1021      {
1022        CoinBigIndex k = mcstrt[jcolx];
1023        int nx = hincol[jcolx];
1024        double dj = maxmin * dcost[jcolx];
1025
1026        for (int i=0; i<nx; ++i) {
1027          int row = hrow[k];
1028          double coeff = colels[k];
1029          k = link[k];
1030
1031          dj -= rowduals[row] * coeff;
1032        }
1033        if (! (fabs(rcosts[jcolx] - dj) < 100*ZTOLDP))
1034          printf("BAD DOUBLE X DJ:  %d %d %g %g\n",
1035                 irow, jcolx, rcosts[jcolx], dj);
1036      }
1037      {
1038        CoinBigIndex k = mcstrt[jcoly];
1039        int ny = hincol[jcoly];
1040        double dj = maxmin * dcost[jcoly];
1041
1042        for (int i=0; i<ny; ++i) {
1043          int row = hrow[k];
1044          double coeff = colels[k];
1045          k = link[k];
1046
1047          dj -= rowduals[row] * coeff;
1048        }
1049        if (! (fabs(rcosts[jcoly] - dj) < 100*ZTOLDP))
1050          printf("BAD DOUBLE Y DJ:  %d %d %g %g\n",
1051                 irow, jcoly, rcosts[jcoly], dj);
1052      }
1053#endif
1054
1055    cdone[jcoly] = DOUBLETON;
1056    rdone[irow] = DOUBLETON;
1057  }
1058  prob->free_list_ = free_list;
1059}
1060
1061
1062doubleton_action::~doubleton_action()
1063{
1064  for (int i=nactions_-1; i>=0; i--) {
1065    delete[]actions_[i].colx;
1066    delete[]actions_[i].indx;
1067    delete[]actions_[i].coly;
1068    delete[]actions_[i].indy;
1069  }
1070  deleteAction(actions_,action*);
1071}
1072
1073
1074
1075static double *doubleton_mult;
1076static int *doubleton_id;
1077void check_doubletons(const PresolveAction * paction)
1078{
1079  const PresolveAction * paction0 = paction;
1080
1081  if (paction) {
1082    check_doubletons(paction->next);
1083   
1084    if (strcmp(paction0->name(), "doubleton_action") == 0) {
1085      const doubleton_action *daction = (const doubleton_action *)paction0;
1086      for (int i=daction->nactions_-1; i>=0; --i) {
1087        int icolx = daction->actions_[i].icolx;
1088        int icoly = daction->actions_[i].icoly;
1089        double coeffx = daction->actions_[i].coeffx;
1090        double coeffy = daction->actions_[i].coeffy;
1091
1092        doubleton_mult[icoly] = -coeffx/coeffy;
1093        doubleton_id[icoly] = icolx;
1094      }
1095    }
1096  }
1097}
1098
1099void check_doubletons1(const PresolveAction * paction,
1100                       int ncols)
1101{
1102#if     DEBUG_PRESOLVE
1103  doubleton_mult = new double[ncols];
1104  doubleton_id = new int[ncols];
1105  for (int i=0; i<ncols; ++i)
1106    doubleton_id[i] = i;
1107  check_doubletons(paction);
1108  double minmult = 1.0;
1109  int minid = -1;
1110  for (int i=0; i<ncols; ++i) {
1111    double mult = 1.0;
1112    int j = i;
1113    if (doubleton_id[j] != j) {
1114      printf("MULTS (%d):  ", j);
1115      while (doubleton_id[j] != j) {
1116        printf("%d %g, ", doubleton_id[j], doubleton_mult[j]);
1117        mult *= doubleton_mult[j];
1118        j = doubleton_id[j];
1119      }
1120      printf(" == %g\n", mult);
1121      if (minmult > fabs(mult)) {
1122        minmult = fabs(mult);
1123        minid = i;
1124      }
1125    }
1126  }
1127  if (minid != -1)
1128    printf("MIN MULT:  %d %g\n", minid, minmult);
1129#endif
1130}
Note: See TracBrowser for help on using the repository browser.