source: branches/devel-1/PresolveDoubleton.cpp @ 2351

Last change on this file since 2351 was 49, checked in by ladanyi, 17 years ago

everything compiles

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