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

Last change on this file since 29 was 29, checked in by forrest, 17 years ago

Presolve (no changes to Makefile)

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