source: branches/devel-1/PresolveSubst.cpp @ 33

Last change on this file since 33 was 31, checked in by forrest, 18 years ago

Presolve nearly ready

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 35.4 KB
Line 
1#include <stdio.h>
2#include <math.h>
3#include <strings.h>
4
5#include "PresolveMatrix.hpp"
6#include "PresolveEmpty.hpp"    // for DROP_COL/DROP_ROW
7#include "PresolveFixed.hpp"
8#include "PresolveZeros.hpp"
9#include "PresolveSubst.hpp"
10#include "ClpMessage.hpp"
11#include "CoinSort.hpp"
12
13inline void prepend_elem(int jcol, double coeff, int irow,
14                    int *mcstrt,
15                    double *colels,
16                    int *hrow,
17                    int *link, int *free_listp)
18{
19  int kk = *free_listp;
20  *free_listp = link[*free_listp];
21  check_free_list(*free_listp);
22
23  link[kk] = mcstrt[jcol];
24  mcstrt[jcol] = kk;
25  colels[kk] = coeff;
26  hrow[kk] = irow;
27}
28
29const char *subst_constraint_action::name() const
30{
31  return ("subst_constraint_action");
32}
33
34void compact_rep(double *elems, int *indices, int *starts, const int *lengths, int n,
35                 const presolvehlink *link);
36
37inline int min(int x, int y)
38{
39  return (x < y) ? x : y;
40}
41
42
43// copy of expand_col; have to rename params
44static void expand_row(int *mcstrt, 
45                    double *colels,
46                       int *hrow, // int *hcol,
47                       //int *hinrow,
48                       int *hincol,
49                    presolvehlink *clink, int ncols,
50
51                       int icolx
52                       ///, int icoly
53                       )
54{
55  /////int kcs = mcstrt[icoly];
56  /////int kce = kcs + hincol[icoly];
57  int kcsx = mcstrt[icolx];
58  int kcex = kcsx + hincol[icolx];
59
60  const int maxk = mcstrt[ncols];       // (22)
61
62  // update col rep - need to expand the column, though.
63  int nextcol = clink[icolx].suc;
64
65  // (22)
66  if (kcex + 1 < mcstrt[nextcol] || nextcol == ncols) {
67    if (! (kcex + 1 < mcstrt[nextcol])) {
68      // nextcol==ncols and no space - must compact
69      compact_rep(colels, hrow, mcstrt, hincol, ncols, clink);
70
71      // update vars
72      kcsx = mcstrt[icolx];
73      kcex = kcsx + hincol[icolx];
74
75      if (! (kcex + 1 < mcstrt[nextcol])) {
76        abort();
77      }
78    }
79  } else {
80    // this is not the last col
81    // fetch last non-empty col (presolve_make_memlists-1)
82    int lastcol = clink[ncols].pre;
83    // (clink[icolx].suc != ncols) ==> (icolx != lastcol)
84
85    // put it directly after the last column
86    int newkcsx = mcstrt[lastcol] + hincol[lastcol];
87
88    // well, pad it a bit
89    newkcsx += min(hincol[icolx], 5); // slack
90
91    //printf("EXPAND_ROW:  %d %d %d\n", newkcsx, maxk, icolx);
92
93    if (newkcsx + hincol[icolx] + 1 >= maxk) {
94      compact_rep(colels, hrow, mcstrt, hincol, ncols, clink);
95
96      // update vars
97      kcsx = mcstrt[icolx];
98      kcex = kcsx + hincol[icolx];
99
100      newkcsx = mcstrt[lastcol] + hincol[lastcol];
101
102      if (newkcsx + hincol[icolx] + 1 >= maxk) {
103        abort();
104      }
105      // have to adjust various induction variables
106      ////kcoly = mcstrt[icoly] + (kcoly - kcs);
107      /////kcs = mcstrt[icoly];                 // do this for ease of debugging
108      /////kce = mcstrt[icoly] + hincol[icoly];
109           
110      /////kcolx = mcstrt[icolx] + (kcolx - kcs);       // don't really need to do this
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#if 0
129    hincol[icolx]++;
130    kcex = newkcsx + hincol[icolx];
131
132    // move the column - 4:  add the new entry
133    hrow[kcex-1] = row;
134    colels[kcex-1] = colels[kcoly] * coeff_factor;
135#endif
136  }
137}
138
139// add coeff_factor * rowy to rowx
140void add_row(int *mrstrt, 
141             double *rlo, double *rup,
142             double *rowels,
143             int *hcol,
144             int *hinrow,
145             presolvehlink *rlink, int nrows,
146             double coeff_factor,
147             int irowx, int irowy,
148             int *x_to_y)
149{
150  int krs = mrstrt[irowy];
151  int kre = krs + hinrow[irowy];
152  int krsx = mrstrt[irowx];
153  int krex = krsx + hinrow[irowx];
154  const int maxk = mrstrt[nrows];       // (22)
155
156  // if irowx is very long, the searching gets very slow,
157  // so we always sort.
158  // whatever sorts rows should handle almost-sorted data efficiently
159  // (quicksort may not)
160  CoinSort_2(hcol+krsx,hcol+krsx+hinrow[irowx],rowels+krsx);
161  CoinSort_2(hcol+krs,hcol+krs+hinrow[irowy],rowels+krs);
162  //ekk_sort2(hcol+krsx, rowels+krsx, hinrow[irowx]);
163  //ekk_sort2(hcol+krs,  rowels+krs,  hinrow[irowy]);
164
165  //printf("%s x=%d y=%d cf=%g nx=%d ny=%d\n",
166  // "ADD_ROW:",
167  //  irowx, irowy, coeff_factor, hinrow[irowx], hinrow[irowy]);
168
169#if     DEBUG_PRESOLVE
170  printf("%s x=%d y=%d cf=%g nx=%d ycols=(",
171         "ADD_ROW:",
172          irowx, irowy, coeff_factor, hinrow[irowx]);
173#endif
174
175  // adjust row bounds of rowx;
176  // analogous to adjusting bounds info of colx in doubleton,
177  // or perhaps adjustment to rlo/rup in elim_doubleton
178  //
179  // I believe that since we choose a column that is implied free,
180  // no other column bounds need to be updated.
181  // This is what would happen in doubleton if y's bounds were implied free;
182  // in that case,
183  // lo1 would never improve clo[icolx] and
184  // up1 would never improve cup[icolx].
185  {
186    double rhsy = rlo[irowy];
187
188    // (1)
189    if (-PRESOLVE_INF < rlo[irowx]) {
190#if     DEBUG_PRESOLVE
191      if (rhsy * coeff_factor)
192        printf("ELIM_ROW RLO:  %g -> %g\n",
193               rlo[irowx],
194               rlo[irowx] + rhsy * coeff_factor);
195#endif
196      rlo[irowx] += rhsy * coeff_factor;
197    }
198    // (2)
199    if (rup[irowx] < PRESOLVE_INF) {
200#if     DEBUG_PRESOLVE
201      if (rhsy * coeff_factor)
202        printf("ELIM_ROW RUP:  %g -> %g\n",
203               rup[irowx],
204               rup[irowx] + rhsy * coeff_factor);
205#endif
206      rup[irowx] += rhsy * coeff_factor;
207    }
208  }
209
210  int kcolx = krsx;
211  int krex0 = krex;
212  int x_to_y_i = 0;
213
214  for (int krowy=krs; krowy<kre; krowy++) {
215    int jcol = hcol[krowy];
216
217    // even though these values are updated, they remain consistent
218    PRESOLVEASSERT(krex == krsx + hinrow[irowx]);
219
220    // see if row appears in colx
221    // do NOT look beyond the original elements of rowx
222    //int kcolx = presolve_find_row1(jcol, krsx, krex, hcol);
223    while (kcolx < krex0 && hcol[kcolx] < jcol)
224      kcolx++;
225
226#if     DEBUG_PRESOLVE
227    printf("%d%s ", jcol, (kcolx < krex0 && hcol[kcolx] == jcol) ? "+" : "");
228#endif
229
230    if (kcolx < krex0 && hcol[kcolx] == jcol) {
231      // before:  both x and y are in the jcol
232      // after:   only x is in the jcol
233      // so: number of elems in col x unchanged, and num elems in jcol is one less
234
235      // update row rep - just modify coefficent
236      // column y is deleted as a whole at the end of the loop
237#if     DEBUG_PRESOLVE
238      printf("CHANGING %g + %g -> %g\n",
239             rowels[kcolx],
240             rowels[krowy],
241             rowels[kcolx] + rowels[krowy] * coeff_factor);
242#endif
243      rowels[kcolx] += rowels[krowy] * coeff_factor;
244
245      // this is where this element in rowy ended up
246      x_to_y[x_to_y_i++] = kcolx - krsx;
247      kcolx++;
248    } else {
249      // before:  only y is in the jcol
250      // after:   only x is in the jcol
251      // so: number of elems in col x is one greater, but num elems in jcol remains same
252      {
253        expand_row(mrstrt, rowels, hcol, hinrow, rlink, nrows, irowx);
254        // this may force a compaction
255        // this will be called excessively if the rows are packed too tightly
256
257        // have to adjust various induction variables
258        krowy = mrstrt[irowy] + (krowy - krs);
259        krs = mrstrt[irowy];                    // do this for ease of debugging
260        kre = mrstrt[irowy] + hinrow[irowy];
261           
262        kcolx = mrstrt[irowx] + (kcolx - krsx); // don't really need to do this
263        krex0 = mrstrt[irowx] + (krex0 - krsx);
264        krsx = mrstrt[irowx];
265        krex = mrstrt[irowx] + hinrow[irowx];
266      }
267
268      // this is where this element in rowy ended up
269      x_to_y[x_to_y_i++] = krex - krsx;
270
271      // there is now an unused entry in the memory after the column - use it
272      // mrstrt[nrows] == penultimate index of arrays hcol/rowels
273      hcol[krex] = jcol;
274      rowels[krex] = rowels[krowy] * coeff_factor;
275      hinrow[irowx]++, krex++;  // expand the col
276
277      // do NOT increment kcolx
278    }
279  }
280
281#if     DEBUG_PRESOLVE
282  printf(")\n");
283#endif
284}
285
286
287// It is common in osl to copy from one representation to another
288// (say from a col rep to a row rep).
289// One such routine is ekkclcp.
290// This is similar, except that it does not assume that the
291// representation is packed, and it adds some slack space
292// in the target rep.
293// It assumes both hincol/hinrow are correct.
294// Note that such routines automatically sort the target rep by index,
295// because they sweep the rows in ascending order.
296void copyrep(const int * mrstrt, const int *hcol, const double *rowels, 
297             const int *hinrow, int nrows,
298             int *mcstrt, int *hrow, double *colels, 
299             int *hincol, int ncols)
300{
301  int pos = 0;
302  for (int j = 0; j < ncols; ++j) {
303    mcstrt[j] = pos;
304    pos += hincol[j];
305    pos += min(hincol[j], 10); // slack
306    hincol[j] = 0;
307  }
308
309  for (int i = 0; i < nrows; ++i) {
310    int krs = mrstrt[i];
311    int kre = krs + hinrow[i];
312    for (int kr = krs; kr < kre; ++kr) {
313      int icol = hcol[kr];
314      int iput = hincol[icol];
315      hincol[icol] = iput + 1;
316      iput += mcstrt[icol];
317     
318      hrow[iput] = i;
319      colels[iput] = rowels[kr];
320    }
321  }
322}
323
324#if 0
325// variant of function of same name in ekk_7.c
326// exactly the same, with some code #if'ed out
327int ekk_makeColumnOrdered(int * indexRow , int * indexColumn , double * element ,
328                          int * rowCount , int * columnCount , int * startColumn ,
329                          int numberRows , int numberColumns,
330                          int numberElements, double tolerance)
331{
332  int iColumn,i,k;
333#if 0
334  for (i=0;i<numberRows;i++) {
335    rowCount[i]=0;
336  }
337  for (i=0;i<numberColumns;i++) {
338    columnCount[i]=0;
339  }
340  for (i=0;i<numberElements;i++) {
341    int iRow=indexRow[i];
342    int iColumn=indexColumn[i];
343    rowCount[iRow]++;
344    columnCount[iColumn]++;
345  }
346#endif
347
348  i=0;
349  for (iColumn=0;iColumn<numberColumns;iColumn++) {
350    /* position after end of Column */
351    i+=columnCount[iColumn];
352    startColumn[iColumn]=i;
353  } /* endfor */
354  startColumn[iColumn]=i;
355
356  for (k=numberElements-1;k>=0;k--) {
357    iColumn=indexColumn[k];
358    if (iColumn>=0) {
359      /* pick up the entry with your right hand */
360      double value = element[k];
361      int iRow=indexRow[k];
362      int iColumnSave=0;
363      indexColumn[k]=-2;        /* the hole */
364
365      while (1) {
366        /* pick this up with your left */
367        int iLook=startColumn[iColumn]-1;
368        double valueSave=element[iLook];
369        int iColumnSave=indexColumn[iLook];
370        int iRowSave=indexRow[iLook];
371
372        /* put the right-hand entry where it wanted to go */
373        startColumn[iColumn]=iLook;
374        element[iLook]=value;
375        indexRow[iLook]=iRow;
376        indexColumn[iLook]=-1;  /* mark it as being where it wants to be */
377       
378        /* there was something there */
379        if (iColumnSave>=0) {
380          iColumn=iColumnSave;
381          value=valueSave;
382          iRow=iRowSave;
383        } else if (iColumnSave = -2)    /* that was the hole */
384          break;
385        else
386          ekkmesg_no(158);      /* should never happen */
387        /* endif */
388      } /* endwhile */
389    } /* endif */
390  } /* endfor */
391
392#if 0
393  /* now pack the elements and combine entries with the same row and column */
394  /* also, drop entries with "small" coefficients */
395  numberElements=0;
396  for (iColumn=0;iColumn<numberColumns;iColumn++) {
397    int start=startColumn[iColumn];
398    int end =startColumn[iColumn+1];
399    startColumn[iColumn]=numberElements;
400    if (end>start) {
401      int lastRow;
402      double lastValue;
403      CoinSort_2(indexRow+start,indexRow+end,element+start);
404      //ekk_sort2(indexRow+start,element+start,end-start);
405      lastRow=indexRow[start];
406      lastValue=element[start];
407      for (i=start+1;i<end;i++) {
408        int iRow=indexRow[i];
409        double value=element[i];
410        if (iRow>lastRow) {
411          if(fabs(lastValue)>tolerance) {
412            indexRow[numberElements]=lastRow;
413            element[numberElements]=lastValue;
414            numberElements++;
415          }
416          lastRow=iRow;
417          lastValue=value;
418        } else {
419          lastValue+=value;
420        } /* endif */
421      } /* endfor */
422      if(fabs(lastValue)>tolerance) {
423        indexRow[numberElements]=lastRow;
424        element[numberElements]=lastValue;
425        numberElements++;
426      }
427    }
428  } /* endfor */
429#endif
430
431  startColumn[numberColumns]=numberElements;
432  return numberElements;
433}
434#endif
435
436
437// add -x/y times row y to row x, thus cancelling out one column of rowx;
438// afterwards, that col will be singleton for rowy, so we drop the row.
439//
440// This no longer maintains the col rep as it goes along.
441// Instead, it reconstructs it from scratch afterward.
442//
443// This implements the functionality of ekkrdc3.
444const PresolveAction *subst_constraint_action::presolve(PresolveMatrix *prob,
445                                         char *implied_free,
446                                        const PresolveAction *next,
447                                                        int &try_fill_level)
448{
449  double *colels        = prob->colels_;
450  int *hrow     = prob->hrow_;
451  int *mcstrt   = prob->mcstrt_;
452  int *hincol   = prob->hincol_;
453  const int ncols       = prob->ncols_;
454
455  const double *clo     = prob->clo_;
456  const double *cup     = prob->cup_;
457
458  double *rowels        = prob->rowels_;
459  int *hcol     = prob->hcol_;
460  int *mrstrt   = prob->mrstrt_;
461  int *hinrow   = prob->hinrow_;
462  const int nrows       = prob->nrows_;
463
464  double *rlo   = prob->rlo_;
465  double *rup   = prob->rup_;
466
467  double *dcost         = prob->cost_;
468
469  presolvehlink *clink = prob->clink_;
470  presolvehlink *rlink = prob->rlink_;
471
472  const double tol = ZTOLDP;
473
474  action *actions       = new action [ncols];
475  int nactions = 0;
476
477  int *zerocols = new int[ncols];
478  int nzerocols = 0;
479
480  int *x_to_y   = new int[ncols];
481
482#if 0
483  // follmer.mps presents a challenge, since it has some very
484  // long rows.  I started experimenting with how to deal with it,
485  // but haven't yet finished.
486  // The idea was to space out the rows to add some padding between them.
487  // Ideally, we shouldn't have to do this just here, but could try to
488  // do it a little everywhere.
489
490  // sort the row rep by reconstructing from col rep
491  copyrep(mcstrt, hrow, colels, hincol, ncols,
492          mrstrt, hcol, rowels, hinrow, nrows);
493  presolve_make_memlists(mrstrt, hinrow, rlink, nrows);
494  // NEED SOME ASSERTION ABOUT NELEMS
495
496  copyrep(mrstrt, hcol, rowels, hinrow, nrows,
497          mcstrt, hrow, colels, hincol, ncols);
498  presolve_make_memlists(mcstrt, hincol, clink, ncols);
499#endif
500
501  // in the original presolve, I don't think the two representations were
502  // kept in sync.  It may be useful not to do that here, either,
503  // but rather just keep the columns with nfill_level rows accurate
504  // and resync at the end of the function.
505
506  // DEBUGGING
507  int nt = 0;
508  int ngood = 0;
509  int nsubst = 0;
510#ifdef  DEBUG_PRESOLVEx
511  int maxsubst = atoi(getenv("MAXSUBST"));
512#else
513  const int maxsubst = 1000000;
514#endif
515
516  // This loop does very nearly the same thing as
517  // the first loop in implied_free_action::presolve.
518  // We have to do it again in case constraints change while we process them (???).
519  int numberLook = prob->numberColsToDo_;
520  int iLook;
521  int * look = prob->colsToDo_;
522  int fill_level = try_fill_level;
523  int * look2 = NULL;
524  // if gone from 2 to 3 look at all
525  if (fill_level<0) {
526    fill_level=-fill_level;
527    try_fill_level=fill_level;
528    look2 = new int[ncols];
529    look=look2;
530    for (iLook=0;iLook<ncols;iLook++) 
531      look[iLook]=iLook;
532    numberLook=ncols;
533 }
534 
535
536  for (iLook=0;iLook<numberLook;iLook++) {
537    int jcoly=look[iLook];
538    if (hincol[jcoly] > 1 && hincol[jcoly] <= fill_level &&
539        implied_free[jcoly] == hincol[jcoly]) {
540      int kcs = mcstrt[jcoly];
541      int kce = kcs + hincol[jcoly];
542
543      int bestrowy_size = 0;
544      int bestrowy_row=-1;
545      int bestrowy_k=-1;
546      double bestrowy_coeff=0.0;
547
548      for (int k=kcs; k<kce; ++k) {
549        int row = hrow[k];
550        double coeffj = colels[k];
551
552        // we don't clean up zeros in the middle of the routine.
553        // if there is one, skip this candidate.
554        if (fabs(coeffj) <= ZTOLDP) {
555          bestrowy_size = 0;
556          break;
557        }
558
559        // if its row is an equality constraint...
560        if (hinrow[row] > 1 &&  // don't bother with singleton rows
561
562            fabs(rlo[row] - rup[row]) < tol) {
563          int krs = mrstrt[row];
564          int kre = krs + hinrow[row];
565
566          double maxup, maxdown, ilow, iup;
567
568          implied_bounds(rowels, clo, cup, hcol,
569                         krs, kre,
570                         &maxup, &maxdown,
571                         jcoly, rlo[row], rup[row], &ilow, &iup);
572
573          if (maxup < PRESOLVE_INF && maxup + tol < rlo[row]) {
574            prob->status_|= 1;
575            prob->originalModel_->messageHandler()->message(CLP_PRESOLVE_ROWINFEAS,
576                                             prob->originalModel_->messages())
577                                               <<row
578                                               <<rlo[row]
579                                               <<rup[row]
580                                               <<CoinMessageEol;
581            break;
582          } else if (-PRESOLVE_INF < maxdown && rup[row] < maxdown - tol) {
583            prob->status_|= 1;
584            prob->originalModel_->messageHandler()->message(CLP_PRESOLVE_ROWINFEAS,
585                                             prob->originalModel_->messages())
586                                               <<row
587                                               <<rlo[row]
588                                               <<rup[row]
589                                               <<CoinMessageEol;
590            break;
591          } else {
592            // the row has an implied upper or a lower bound
593
594            if (clo[jcoly] <= ilow && iup <= cup[jcoly]) {
595              // both column bounds implied by the constraint bounds
596
597              // we want coeffy to be smaller than x, BACKWARDS from in doubleton
598              if (bestrowy_size == 0 ||
599                  fabs(coeffj) > fabs(bestrowy_coeff) ||
600                  (fabs(coeffj) == fabs(bestrowy_coeff) &&
601                   hinrow[row] < bestrowy_size)) {
602                bestrowy_size = hinrow[row];
603                bestrowy_row = row;
604                bestrowy_coeff = coeffj;
605                bestrowy_k = k;
606              }
607            }
608          }
609        }
610      }
611
612      if (bestrowy_size == 0)
613        continue;
614
615      bool all_ok = true;
616      for (int k=kcs; k<kce; ++k) {
617        double coeff_factor = fabs(colels[k] / bestrowy_coeff);
618        if (fabs(coeff_factor) > 1.3)
619          all_ok = false;
620      }
621
622      // check fill-in
623      if (all_ok && hincol[jcoly] == 3) {
624        // compute fill-in
625        int row1 = -1;
626        int row2;
627        for (int kk=kcs; kk<kce; ++kk) 
628          if (kk != bestrowy_k) {
629            if (row1 == -1)
630              row1 = hrow[kk];
631            else
632              row2 = hrow[kk];
633          }
634
635
636        int krs = mrstrt[bestrowy_row];
637        int kre = krs + hinrow[bestrowy_row];
638        int krs1 = mrstrt[row1];
639        int kre1 = krs + hinrow[row1];
640        int krs2 = mrstrt[row2];
641        int kre2 = krs + hinrow[row2];
642
643        CoinSort_2(hcol+krs,hcol+krs+hinrow[bestrowy_row],rowels+krs);
644        CoinSort_2(hcol+krs1,hcol+krs1+hinrow[row1],rowels+krs1);
645        CoinSort_2(hcol+krs2,hcol+krs2+hinrow[row2],rowels+krs2);
646        //ekk_sort2(hcol+krs,  rowels+krs,  hinrow[bestrowy_row]);
647        //ekk_sort2(hcol+krs1, rowels+krs1, hinrow[row1]);
648        //ekk_sort2(hcol+krs2, rowels+krs2, hinrow[row2]);
649
650        int nfill = -hinrow[bestrowy_row];
651        int kcol1 = krs1;
652        for (int kk=krs; kk<kre; ++kk) {
653          int jcol = hcol[kk];
654
655          while (kcol1 < kre1 && hcol[kcol1] < jcol)
656            kcol1++;
657          if (! (kcol1 < kre1 && hcol[kcol1] == jcol))
658            nfill++;
659        }
660        int kcol2 = krs2;
661        for (int kk=krs; kk<kre; ++kk) {
662          int jcol = hcol[kk];
663
664          while (kcol2 < kre2 && hcol[kcol2] < jcol)
665            kcol2++;
666          if (! (kcol2 < kre2 && hcol[kcol2] == jcol))
667            nfill++;
668        }
669#if     DEBUG_PRESOLVE
670        printf("FILL:  %d\n", nfill);
671#endif
672
673#if 0
674        static int maxfill = atoi(getenv("MAXFILL"));
675
676        if (nfill > maxfill)
677          all_ok = false;
678#endif
679
680        // not too much
681        if (nfill <= 0)
682          ngood++;
683
684#if 0
685        static int nts = 0;
686        if (++nts > atoi(getenv("NTS")))
687          all_ok = false;
688        else
689          nt++;
690#endif
691      }
692
693      // probably never happens
694      if (all_ok && nzerocols + hinrow[bestrowy_row] >= ncols)
695        all_ok = false;
696
697      if (nsubst >= maxsubst) {
698        all_ok = false;
699      } 
700
701      if (all_ok) {
702        nsubst++;
703#if 0
704        // debug
705        if (numberLook<ncols&&iLook==numberLook-1) {
706          printf("found last one?? %d\n", jcoly);
707        }
708#endif
709
710        int kcs = mcstrt[jcoly];
711        int rowy = bestrowy_row;
712        double coeffy = bestrowy_coeff;
713
714        PRESOLVEASSERT(fabs(colels[kcs]) > ZTOLDP);
715        PRESOLVEASSERT(fabs(colels[kcs+1]) > ZTOLDP);
716
717        PRESOLVEASSERT(hinrow[rowy] > 1);
718
719        const bool nonzero_cost = (fabs(dcost[jcoly]) > tol);
720
721        double *costsx = nonzero_cost ? new double[hinrow[rowy]] : 0;
722
723        int ntotels = 0;
724        for (int k=kcs; k<kce; ++k) {
725          int irow = hrow[k];
726          ntotels += hinrow[irow];
727        }
728
729        {
730          action *ap = &actions[nactions++];
731          int nincol = hincol[jcoly];
732
733          ap->col = jcoly;
734          ap->rowy = rowy;
735
736          ap->nincol = nincol;
737          ap->rows = new int[nincol];
738          ap->rlos = new double[nincol];
739          ap->rups = new double[nincol];
740
741          // coefficients in deleted col
742          ap->coeffxs = new double[nincol];
743
744          ap->ninrowxs = new int[nincol];
745          ap->rowcolsxs = new int[ntotels];
746          ap->rowelsxs = new double[ntotels];
747
748          ap->costsx = costsx;
749
750          // copy all the rows for restoring later - wasteful
751          {
752            int nel = 0;
753            for (int k=kcs; k<kce; ++k) {
754              int irow = hrow[k];
755              int krs = mrstrt[irow];
756              int kre = krs + hinrow[irow];
757
758              prob->addRow(irow);
759              ap->rows[k-kcs] = irow;
760              ap->ninrowxs[k-kcs] = hinrow[irow];
761              ap->rlos[k-kcs] = rlo[irow];
762              ap->rups[k-kcs] = rup[irow];
763
764              ap->coeffxs[k-kcs] = colels[k];
765
766              bcopy(&hcol[krs], &ap->rowcolsxs[nel], hinrow[irow]*sizeof(int));
767              bcopy(&rowels[krs], &ap->rowelsxs[nel], hinrow[irow]*sizeof(double));
768              nel += hinrow[irow];
769            }
770          }
771        }
772
773        // rowy is supposed to be an equality row
774        PRESOLVEASSERT(fabs(rup[rowy] - rlo[rowy]) < ZTOLDP);
775
776        // now adjust for the implied free row - COPIED
777        if (nonzero_cost) {
778#ifdef  DEBUG_PRESOLVE
779          printf("NONZERO SUBST COST:  %d %g\n", jcoly, dcost[jcoly]);
780#endif
781          double *cost = dcost;
782          double *save_costs = costsx;
783          double coeffj = coeffy;
784          int krs = mrstrt[rowy];
785          int kre = krs + hinrow[rowy];
786
787          double rhs = rlo[rowy];
788          double costj = cost[jcoly];
789
790          for (int k=krs; k<kre; k++) {
791            int jcol = hcol[k];
792            prob->addCol(jcol);
793            save_costs[k-krs] = cost[jcol];
794
795            if (jcol != jcoly) {
796              double coeff = rowels[k];
797
798              /*
799               * Similar to eliminating doubleton:
800               *   cost1 x = cost1 (c - b y) / a = (c cost1)/a - (b cost1)/a
801               *   cost[icoly] += cost[icolx] * (-coeff2 / coeff1);
802               */
803              cost[jcol] += costj * (-coeff / coeffj);
804            }
805          }
806
807          // I'm not sure about this
808          prob->change_bias(costj * rhs / coeffj);
809
810          // ??
811          cost[jcoly] = 0.0;
812        }
813
814#if     DEBUG_PRESOLVE
815            if (hincol[jcoly] == 3) {
816              int krs = mrstrt[rowy];
817              int kre = krs + hinrow[rowy];
818              printf("HROW0 (%d):  ", rowy);
819              for (int k=krs; k<kre; ++k) {
820                int jcol = hcol[k];
821                double coeff = rowels[k];
822                printf("%d:%g (%d) ", jcol, coeff, hincol[jcol]);
823              }
824              printf("\n");
825            }
826#endif
827
828            if (hincol[jcoly] != 2) {
829              int krs = mrstrt[rowy];
830              int kre = krs + hinrow[rowy];
831              CoinSort_2(hcol+krs,hcol+krs+hinrow[rowy],rowels+krs);
832              //ekk_sort2(hcol+krs,  rowels+krs,  hinrow[rowy]);
833            }
834
835            // substitute away jcoly in the other rows
836            // Use ap as mcstrt etc may move if compacted
837            kce = hincol[jcoly];
838            int k;
839            action *ap = &actions[nactions-1];
840            for (k=0; k<kce; ++k) {
841              int rowx = ap->rows[k];
842              //assert(rowx==hrow[k+kcs]);
843              //assert(ap->coeffxs[k]==colels[k+kcs]);
844              if (rowx != rowy) {
845                double coeffx = ap->coeffxs[k];
846                double coeff_factor = -coeffx / coeffy; // backwards from doubleton
847
848#if     DEBUG_PRESOLVE
849                {
850                  int krs = mrstrt[rowx];
851                  int kre = krs + hinrow[rowx];
852                  printf("HROW (%d %d %d):  ", rowx, hinrow[rowx], jcoly);
853                  for (int k=krs; k<kre; ++k) {
854                    int jcol = hcol[k];
855                    double coeff = rowels[k];
856                    printf("%d ", jcol);
857                  }
858                  printf("\n");
859#if 0
860                  for (int k=krs; k<kre; ++k) {
861                    int jcol = hcol[k];
862                    prob->addCol(jcol);
863                    double coeff = rowels[k];
864                    printf("%g ", coeff);
865                  }
866                  printf("\n");
867#endif
868                }
869#endif
870                {
871                  int krsx = mrstrt[rowx];
872                  int krex = krsx + hinrow[rowx];
873                  int i;
874                  for (i=krsx;i<krex;i++) 
875                    prob->addCol(hcol[i]);
876                  if (hincol[jcoly] != 2) 
877                    CoinSort_2(hcol+krsx,hcol+krsx+hinrow[rowx],rowels+krsx);
878                  //ekk_sort2(hcol+krsx, rowels+krsx, hinrow[rowx]);
879                }
880               
881                // add (coeff_factor * <rowy>) to rowx
882                // does not affect rowy
883                // may introduce (or cancel) elements in rowx
884                add_row(mrstrt,
885                        rlo, rup,
886                        rowels, hcol,
887                        hinrow,
888                        rlink, nrows,
889                        coeff_factor,
890                        rowx, rowy,
891                        x_to_y);
892               
893                // update col rep of rowx from row rep:
894                // for every col in rowy, copy the elem for that col in rowx
895                // from the row rep to the col rep
896                {
897                  int krs = mrstrt[rowy];
898                  int kre = krs + hinrow[rowy];
899                  int niny = hinrow[rowy];
900                 
901                  int krsx = mrstrt[rowx];
902                  int krex = krsx + hinrow[rowx];
903                  for (int ki=0; ki<niny; ++ki) {
904                    int k = krs + ki;
905                    int jcol = hcol[k];
906                    prob->addCol(jcol);
907                    int kcs = mcstrt[jcol];
908                    int kce = kcs + hincol[jcol];
909                   
910                    //double coeff = rowels[presolve_find_row(jcol, krsx, krex, hcol)];
911                    if (hcol[krsx + x_to_y[ki]] != jcol)
912                      abort();
913                    double coeff = rowels[krsx + x_to_y[ki]];
914                   
915                    // see if rowx appeared in jcol in the col rep
916                    int k2 = presolve_find_row1(rowx, kcs, kce, hrow);
917                   
918                    //PRESOLVEASSERT(fabs(coeff) > ZTOLDP);
919                   
920                    if (k2 < kce) {
921                      // yes - just update the entry
922                      colels[k2] = coeff;
923                    } else {
924                      // no - make room, then append
925                      expand_row(mcstrt, colels, hrow, hincol, clink, ncols, jcol);
926                      kcs = mcstrt[jcol];
927                      kce = kcs + hincol[jcol];
928                     
929                      hrow[kce] = rowx;
930                      colels[kce] = coeff;
931                      hincol[jcol]++;
932                    }
933                  }
934                }
935                // now colels[k] == 0.0
936
937#if 1
938                // now remove jcoly from rowx in the row rep
939                // better if this were first
940                presolve_delete_from_row(rowx, jcoly, mrstrt, hinrow, hcol, rowels);
941#endif
942#if     DEBUG_PRESOLVE
943                {
944                  int krs = mrstrt[rowx];
945                  int kre = krs + hinrow[rowx];
946                  printf("HROW (%d %d %d):  ", rowx, hinrow[rowx], jcoly);
947                  for (int k=krs; k<kre; ++k) {
948                    int jcol = hcol[k];
949                    double coeff = rowels[k];
950                    printf("%d ", jcol);
951                  }
952                  printf("\n");
953#if 0
954                  for (int k=krs; k<kre; ++k) {
955                    int jcol = hcol[k];
956                    double coeff = rowels[k];
957                    printf("%g ", coeff);
958                  }
959                  printf("\n");
960#endif
961                }
962#endif
963               
964                // don't have to update col rep, since entire col deleted later
965              }
966            }
967
968#if     DEBUG_PRESOLVE
969            printf("\n");
970#endif
971
972            // the addition of rows may have created zero coefficients
973            bcopy(&hcol[mrstrt[rowy]], &zerocols[nzerocols], hinrow[rowy]*sizeof(int));
974            nzerocols += hinrow[rowy];
975           
976            // delete rowy in col rep
977            {
978              int krs = mrstrt[rowy];
979              int kre = krs + hinrow[rowy];
980              for (int k=krs; k<kre; ++k) {
981                int jcol = hcol[k];
982               
983                // delete rowy from the jcol
984                presolve_delete_from_row(jcol, rowy, mcstrt, hincol, hrow, colels);
985              }
986            }
987            // delete rowy in row rep
988            hinrow[rowy] = 0;
989           
990            // This last is entirely dual to doubleton, but for the cost adjustment
991           
992            // eliminate col entirely from the col rep
993            PRESOLVE_REMOVE_LINK(clink, jcoly);
994            hincol[jcoly] = 0;
995           
996            // eliminate rowy entirely from the row rep
997            PRESOLVE_REMOVE_LINK(rlink, rowy);
998            //cost[irowy] = 0.0;
999           
1000            rlo[rowy] = 0.0;
1001            rup[rowy] = 0.0;
1002           
1003#if     0 && DEBUG_PRESOLVE
1004            printf("ROWY COLS:  ");
1005            for (int k=0; k<save_ninrowy; ++k)
1006              if (rowycols[k] != col) {
1007                printf("%d ", rowycols[k]);
1008                (void)presolve_find_row(rowycols[k], mrstrt[rowx], mrstrt[rowx]+hinrow[rowx],
1009                                        hcol);
1010              }
1011            printf("\n");
1012#endif
1013#if 0
1014        presolve_links_ok(clink, mcstrt, hincol, ncols);
1015        presolve_links_ok(rlink, mrstrt, hinrow, nrows);
1016        prob->consistent();
1017#endif
1018      }
1019    }
1020  }
1021
1022  // general idea - only do doubletons until there are almost none left
1023  if (nactions < 30&&fill_level==2)
1024    try_fill_level = -3;
1025
1026  if (nactions) {
1027#if     PRESOLVE_SUMMARY
1028    printf("NSUBSTS:  %d\n", nactions);
1029    printf("NT: %d  NGOOD:  %d FILL_LEVEL:  %d\n", nt, ngood, fill_level);
1030#endif
1031    next = new subst_constraint_action(nactions, copyOfArray(actions,nactions), next);
1032
1033    next = drop_zero_coefficients_action::presolve(prob, zerocols, nzerocols, next);
1034  }
1035  delete [] look2;
1036  delete [] actions;
1037
1038  delete[]x_to_y;
1039  delete[]zerocols;
1040
1041  return (next);
1042}
1043
1044
1045void subst_constraint_action::postsolve(PostsolveMatrix *prob) const
1046{
1047  const action *const actions = actions_;
1048  const int nactions    = nactions_;
1049
1050  double *colels        = prob->colels_;
1051  int *hrow             = prob->hrow_;
1052  int *mcstrt           = prob->mcstrt_;
1053  int *hincol           = prob->hincol_;
1054  int *link             = prob->link_;
1055  int ncols             = prob->ncols_;
1056
1057  double *clo   = prob->clo_;
1058  double *cup   = prob->cup_;
1059
1060  double *rlo   = prob->rlo_;
1061  double *rup   = prob->rup_;
1062
1063  double *dcost = prob->cost_;
1064
1065  double *sol   = prob->sol_;
1066  double *rcosts        = prob->rcosts_;
1067
1068  double *acts  = prob->acts_;
1069  double *rowduals = prob->rowduals_;
1070
1071  char *cdone   = prob->cdone_;
1072  char *rdone   = prob->rdone_;
1073  int free_list = prob->free_list_;
1074
1075  const double ztolzb   = prob->ztolzb_;
1076  const double ztoldj   = prob->ztoldj_;
1077  const double maxmin = prob->maxmin_;
1078
1079  for (const action *f = &actions[nactions-1]; actions<=f; f--) {
1080    int icol = f->col;
1081
1082    int nincoly = f->nincol;
1083    double *rlos = f->rlos;
1084    double *rups = f->rups;
1085    int *rows = f->rows;
1086
1087    double *coeffxs = f->coeffxs;
1088
1089    int jrowy = f->rowy;
1090
1091    int *ninrowxs = f->ninrowxs;
1092    const int *rowcolsxs = f->rowcolsxs;
1093    const double *rowelsxs = f->rowelsxs;
1094   
1095    /* the row was in the reduced problem */
1096    for (int i=0; i<nincoly; ++i) {
1097      if (rows[i] != jrowy)
1098        PRESOLVEASSERT(rdone[rows[i]]);
1099    }
1100    PRESOLVEASSERT(cdone[icol]==DROP_COL);
1101    PRESOLVEASSERT(rdone[jrowy]==DROP_ROW);
1102
1103    // DEBUG CHECK
1104#if     0 && DEBUG_PRESOLVE
1105    {
1106      double actx = 0.0;
1107      for (int j=0; j<ncols; ++j)
1108        if (hincol[j] > 0 && cdone[j]) {
1109          int krow = presolve_find_row1(jrowx, mcstrt[j], mcstrt[j] + hincol[j], hrow);
1110          if (krow < mcstrt[j] + hincol[j])
1111            actx += colels[krow] * sol[j];
1112      }
1113      if (! (fabs(acts[jrowx] - actx) < 100*ztolzb))
1114        printf("BAD ACTSX:  acts[%d]==%g != %g\n",
1115               jrowx, acts[jrowx], actx);
1116      if (! (rlo[jrowx] - 100*ztolzb <= actx && actx <= rup[jrowx] + 100*ztolzb))
1117        printf("ACTSX NOT IN RANGE:  %d %g %g %g\n",
1118               jrowx, rlo[jrowx], actx, rup[jrowx]);
1119    }
1120#endif
1121
1122    int ninrowy=-1;
1123    const int *rowcolsy=NULL;
1124    const double *rowelsy=NULL;
1125    double coeffy=0.0;
1126
1127    double rloy;
1128    {
1129      int nel = 0;
1130      for (int i=0; i<nincoly; ++i) {
1131        int row = rows[i];
1132        rlo[row] = rlos[i];
1133        rup[row] = rups[i];
1134        if (row == jrowy) {
1135          ninrowy = ninrowxs[i];
1136          rowcolsy = &rowcolsxs[nel];
1137          rowelsy  = &rowelsxs[nel];
1138
1139          coeffy = coeffxs[i];
1140          rloy = rlo[row];
1141        }
1142        nel += ninrowxs[i];
1143      }
1144    }
1145    double rhsy = rloy;
1146
1147    // restore costs
1148    {
1149      const double *costs = f->costsx;
1150      if (costs)
1151        for (int i = 0; i<ninrowy; ++i) {
1152          dcost[rowcolsy[i]] = costs[i];
1153        }
1154    }
1155
1156    // solve for the equality to find the solution for the eliminated col
1157    // this is why we want coeffx < coeffy (55)
1158    {
1159      double sol0 = rloy;
1160      sol[icol] = 0.0;  // to avoid condition in loop
1161      for (int k = 0; k<ninrowy; ++k) {
1162        int jcolx = rowcolsy[k];
1163        double coeffx = rowelsy[k];
1164        sol0 -= coeffx * sol[jcolx];
1165      }
1166      sol[icol] = sol0 / coeffy;
1167
1168      if (! (sol[icol] > clo[icol] - ztolzb &&
1169             cup[icol] + ztolzb > sol[icol]))
1170        printf("NEW SOL OUT-OF-TOL:  %g %g %g\n", clo[icol], sol[icol], cup[icol]);
1171    }
1172
1173    // since this row is fixed
1174    acts[jrowy] = rloy;
1175
1176    // acts[irow] always ok, since slack is fixed
1177    prob->setRowStatus(jrowy,PrePostsolveMatrix::atLowerBound);
1178
1179    // remove old rowx from col rep
1180    // we don't explicitly store what the current rowx is;
1181    // however, after the presolve, rowx contains a col for every
1182    // col in either the original rowx or the original rowy.
1183    // If there were cancellations, those were handled in subsequent
1184    // presolves.
1185    {
1186      // erase those cols in the other rows that occur in rowy
1187      // (with the exception of icol, which was deleted);
1188      // the other rows *must* contain these cols
1189      for (int k = 0; k<ninrowy; ++k) {
1190        int col = rowcolsy[k];
1191
1192        // remove jrowx from col in the col rep
1193        // icol itself was deleted, so won't be there
1194        if (col != icol)
1195          for (int i = 0; i<nincoly; ++i) {
1196            if (rows[i] != jrowy)
1197              presolve_delete_from_row2(col, rows[i], mcstrt, hincol, hrow, colels, link, &free_list);
1198          }
1199      }
1200
1201      // initialize this for loops below
1202      hincol[icol] = 0;
1203
1204      // now restore the original rows (other than rowy).
1205      // those cols that were also in rowy were just removed;
1206      // otherwise, they *must* already be there.
1207      // This loop and the next automatically create the rep for the new col.
1208      {
1209        const int *rowcolsx = rowcolsxs;
1210        const double *rowelsx = rowelsxs;
1211
1212        for (int i = 0; i<nincoly; ++i) {
1213          int ninrowx = ninrowxs[i];
1214          int jrowx = rows[i];
1215
1216          if (jrowx != jrowy)
1217            for (int k = 0; k<ninrowx; ++k) {
1218              int col = rowcolsx[k];
1219              int kcolx = presolve_find_row3(jrowx, mcstrt[col], hincol[col], hrow, link);
1220
1221              if (kcolx != -1) {
1222                PRESOLVEASSERT(presolve_find_row1(col, 0, ninrowy, rowcolsy) == ninrowy);
1223                // overwrite the existing entry
1224                colels[kcolx] = rowelsx[k];
1225              } else {
1226                PRESOLVEASSERT(presolve_find_row1(col, 0, ninrowy, rowcolsy) < ninrowy);
1227
1228                {
1229                  int kk = free_list;
1230                  free_list = link[free_list];
1231                  check_free_list(free_list);
1232
1233                  link[kk] = mcstrt[col];
1234                  mcstrt[col] = kk;
1235                  colels[kk] = rowelsx[k];
1236                  hrow[kk] = jrowx;
1237                }
1238                ++hincol[col];
1239              }
1240            }
1241          rowcolsx += ninrowx;
1242          rowelsx += ninrowx;
1243        }
1244      }
1245
1246      // finally, add original rowy elements
1247      for (int k = 0; k<ninrowy; ++k) {
1248        int col = rowcolsy[k];
1249
1250        {
1251          prepend_elem(col, rowelsy[k], jrowy, mcstrt, colels, hrow, link, &free_list);
1252          ++hincol[col];
1253        }
1254      }
1255    }
1256
1257    // my guess is that the CLAIM in doubleton generalizes to
1258    // equations with more than one x-style variable.
1259    // Since I can't see how to distinguish among them,
1260    // I assume that any of them will do.
1261
1262    {
1263      int k;
1264      double dj = maxmin*dcost[icol];
1265      double bounds_factor = rhsy/coeffy;
1266      for (int i=0; i<nincoly; ++i)
1267        if (rows[i] != jrowy) {
1268          int row = rows[i];
1269          double coeff = coeffxs[i];
1270
1271          // PROBABLY DOESN'T MAKE SENSE
1272          acts[row] += coeff * bounds_factor;
1273
1274          dj -= rowduals[row] * coeff;
1275        }
1276
1277      // DEBUG CHECK
1278      double acty = 0.0;
1279      for (int k = 0; k<ninrowy; ++k) {
1280        int col = rowcolsy[k];
1281        acty += rowelsy[k] * sol[col];
1282      }
1283      PRESOLVEASSERT(fabs(acty - acts[jrowy]) < 100*ZTOLDP);
1284
1285      // RECOMPUTING
1286      {
1287        const int *rowcolsx = rowcolsxs;
1288        const double *rowelsx = rowelsxs;
1289
1290        for (int i=0; i<nincoly; ++i) {
1291          int ninrowx = ninrowxs[i];
1292
1293          if (rows[i] != jrowy) {
1294            int jrowx = rows[i];
1295
1296            double actx = 0.0;
1297            for (int k = 0; k<ninrowx; ++k) {
1298              int col = rowcolsx[k];
1299              actx += rowelsx[k] * sol[col];
1300            }
1301#if     DEBUG_PRESOLVE
1302            PRESOLVEASSERT(rlo[jrowx] - ztolzb <= actx && actx <= rup[jrowx] + ztolzb);
1303#endif
1304            acts[jrowx] = actx;
1305          }
1306          rowcolsx += ninrowx;
1307          rowelsx += ninrowx;
1308        }
1309      }
1310
1311      // this is the coefficient we need to force col y's reduced cost to 0.0;
1312      // for example, this is obviously true if y is a singleton column
1313      rowduals[jrowy] = dj / coeffy;
1314      rcosts[icol] = 0.0;
1315
1316      // furthermore, this should leave rcosts[colx] for all colx
1317      // in jrowx unchanged (????).
1318    }
1319
1320    // Unlike doubleton, there should never be a problem with keeping
1321    // the reduced costs the way they were, because the other
1322    // variable's bounds are never changed, since col was implied free.
1323    //rowstat[jrowy] = 0;
1324    prob->setColumnStatus(icol,PrePostsolveMatrix::basic);
1325
1326    cdone[icol] = SUBST_ROW;
1327    rdone[jrowy] = SUBST_ROW;
1328  }
1329
1330  prob->free_list_ = free_list;
1331}
1332
1333
1334
1335subst_constraint_action::~subst_constraint_action()
1336{
1337  const action *actions = actions_;
1338
1339  for (int i=0; i<nactions_; ++i) {
1340    delete[]actions[i].rows;
1341    delete[]actions[i].rlos;
1342    delete[]actions[i].rups;
1343    delete[]actions[i].coeffxs;
1344    delete[]actions[i].ninrowxs;
1345    delete[]actions[i].rowcolsxs;
1346    delete[]actions[i].rowelsxs;
1347    delete[]actions[i].costsx;
1348  }
1349
1350  delete[]actions;
1351
1352}
Note: See TracBrowser for help on using the repository browser.