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

Last change on this file since 29 was 29, checked in by forrest, 18 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: 35.0 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        for (int k=kcs; k<kce; ++k)
837          if (hrow[k] != rowy) {
838            int rowx = hrow[k];
839            double coeffx = colels[k];
840            double coeff_factor = -coeffx / coeffy;     // backwards from doubleton
841
842#if     DEBUG_PRESOLVE
843            {
844              int krs = mrstrt[rowx];
845              int kre = krs + hinrow[rowx];
846              printf("HROW (%d %d %d):  ", rowx, hinrow[rowx], jcoly);
847              for (int k=krs; k<kre; ++k) {
848                int jcol = hcol[k];
849                double coeff = rowels[k];
850                printf("%d ", jcol);
851              }
852              printf("\n");
853#if 0
854              for (int k=krs; k<kre; ++k) {
855                int jcol = hcol[k];
856                prob->addCol(jcol);
857                double coeff = rowels[k];
858                printf("%g ", coeff);
859              }
860              printf("\n");
861#endif
862            }
863#endif
864            {
865                int krsx = mrstrt[rowx];
866                int krex = krsx + hinrow[rowx];
867                int i;
868                for (i=krsx;i<krex;i++) 
869                  prob->addCol(hcol[i]);
870                if (hincol[jcoly] != 2) 
871                  CoinSort_2(hcol+krsx,hcol+krsx+hinrow[rowx],rowels+krsx);
872                //ekk_sort2(hcol+krsx, rowels+krsx, hinrow[rowx]);
873            }
874
875            // add (coeff_factor * <rowy>) to rowx
876            // does not affect rowy
877            // may introduce (or cancel) elements in rowx
878            add_row(mrstrt,
879                    rlo, rup,
880                    rowels, hcol,
881                    hinrow,
882                    rlink, nrows,
883                    coeff_factor,
884                    rowx, rowy,
885                    x_to_y);
886
887            // update col rep of rowx from row rep:
888            // for every col in rowy, copy the elem for that col in rowx
889            // from the row rep to the col rep
890            {
891              int krs = mrstrt[rowy];
892              int kre = krs + hinrow[rowy];
893              int niny = hinrow[rowy];
894
895              int krsx = mrstrt[rowx];
896              int krex = krsx + hinrow[rowx];
897              for (int ki=0; ki<niny; ++ki) {
898                int k = krs + ki;
899                int jcol = hcol[k];
900                prob->addCol(jcol);
901                int kcs = mcstrt[jcol];
902                int kce = kcs + hincol[jcol];
903
904                //double coeff = rowels[presolve_find_row(jcol, krsx, krex, hcol)];
905                if (hcol[krsx + x_to_y[ki]] != jcol)
906                  abort();
907                double coeff = rowels[krsx + x_to_y[ki]];
908
909                // see if rowx appeared in jcol in the col rep
910                int k2 = presolve_find_row1(rowx, kcs, kce, hrow);
911
912                //PRESOLVEASSERT(fabs(coeff) > ZTOLDP);
913
914                if (k2 < kce) {
915                  // yes - just update the entry
916                  colels[k2] = coeff;
917                } else {
918                  // no - make room, then append
919                  expand_row(mcstrt, colels, hrow, hincol, clink, ncols, jcol);
920                  kcs = mcstrt[jcol];
921                  kce = kcs + hincol[jcol];
922
923                  hrow[kce] = rowx;
924                  colels[kce] = coeff;
925                  hincol[jcol]++;
926                }
927              }
928            }
929            // now colels[k] == 0.0
930
931#if 1
932            // now remove jcoly from rowx in the row rep
933            // better if this were first
934            presolve_delete_from_row(rowx, jcoly, mrstrt, hinrow, hcol, rowels);
935#endif
936#if     DEBUG_PRESOLVE
937{
938              int krs = mrstrt[rowx];
939              int kre = krs + hinrow[rowx];
940              printf("HROW (%d %d %d):  ", rowx, hinrow[rowx], jcoly);
941              for (int k=krs; k<kre; ++k) {
942                int jcol = hcol[k];
943                double coeff = rowels[k];
944                printf("%d ", jcol);
945              }
946              printf("\n");
947#if 0
948              for (int k=krs; k<kre; ++k) {
949                int jcol = hcol[k];
950                double coeff = rowels[k];
951                printf("%g ", coeff);
952              }
953              printf("\n");
954#endif
955            }
956#endif
957
958            // don't have to update col rep, since entire col deleted later
959          }
960
961#if     DEBUG_PRESOLVE
962        printf("\n");
963#endif
964
965        // the addition of rows may have created zero coefficients
966        bcopy(&hcol[mrstrt[rowy]], &zerocols[nzerocols], hinrow[rowy]*sizeof(int));
967        nzerocols += hinrow[rowy];
968
969        // delete rowy in col rep
970        {
971          int krs = mrstrt[rowy];
972          int kre = krs + hinrow[rowy];
973          for (int k=krs; k<kre; ++k) {
974            int jcol = hcol[k];
975
976            // delete rowy from the jcol
977            presolve_delete_from_row(jcol, rowy, mcstrt, hincol, hrow, colels);
978          }
979        }
980        // delete rowy in row rep
981        hinrow[rowy] = 0;
982   
983        // This last is entirely dual to doubleton, but for the cost adjustment
984
985        // eliminate col entirely from the col rep
986        PRESOLVE_REMOVE_LINK(clink, jcoly);
987        hincol[jcoly] = 0;
988
989        // eliminate rowy entirely from the row rep
990        PRESOLVE_REMOVE_LINK(rlink, rowy);
991        //cost[irowy] = 0.0;
992
993        rlo[rowy] = 0.0;
994        rup[rowy] = 0.0;
995
996#if     0 && DEBUG_PRESOLVE
997        printf("ROWY COLS:  ");
998        for (int k=0; k<save_ninrowy; ++k)
999          if (rowycols[k] != col) {
1000            printf("%d ", rowycols[k]);
1001            (void)presolve_find_row(rowycols[k], mrstrt[rowx], mrstrt[rowx]+hinrow[rowx],
1002                           hcol);
1003          }
1004        printf("\n");
1005#endif
1006
1007#if 0
1008        presolve_links_ok(clink, mcstrt, hincol, ncols);
1009        presolve_links_ok(rlink, mrstrt, hinrow, nrows);
1010        prob->consistent();
1011#endif
1012      }
1013    }
1014  }
1015
1016  // general idea - only do doubletons until there are almost none left
1017  if (nactions < 30&&fill_level==2)
1018    try_fill_level = -3;
1019
1020  if (nactions) {
1021#if     PRESOLVE_SUMMARY
1022    printf("NSUBSTS:  %d\n", nactions);
1023    printf("NT: %d  NGOOD:  %d FILL_LEVEL:  %d\n", nt, ngood, fill_level);
1024#endif
1025    next = new subst_constraint_action(nactions, copyOfArray(actions,nactions), next);
1026
1027    next = drop_zero_coefficients_action::presolve(prob, zerocols, nzerocols, next);
1028  }
1029  delete [] look2;
1030  delete [] actions;
1031
1032  delete[]x_to_y;
1033  delete[]zerocols;
1034
1035  return (next);
1036}
1037
1038
1039void subst_constraint_action::postsolve(PostsolveMatrix *prob) const
1040{
1041  const action *const actions = actions_;
1042  const int nactions    = nactions_;
1043
1044  double *colels        = prob->colels_;
1045  int *hrow             = prob->hrow_;
1046  int *mcstrt           = prob->mcstrt_;
1047  int *hincol           = prob->hincol_;
1048  int *link             = prob->link_;
1049  int ncols             = prob->ncols_;
1050
1051  double *clo   = prob->clo_;
1052  double *cup   = prob->cup_;
1053
1054  double *rlo   = prob->rlo_;
1055  double *rup   = prob->rup_;
1056
1057  double *dcost = prob->cost_;
1058
1059  double *sol   = prob->sol_;
1060  double *rcosts        = prob->rcosts_;
1061
1062  double *acts  = prob->acts_;
1063  double *rowduals = prob->rowduals_;
1064
1065  char *cdone   = prob->cdone_;
1066  char *rdone   = prob->rdone_;
1067  int free_list = prob->free_list_;
1068
1069  const double ztolzb   = prob->ztolzb_;
1070  const double ztoldj   = prob->ztoldj_;
1071
1072  for (const action *f = &actions[nactions-1]; actions<=f; f--) {
1073    int icol = f->col;
1074
1075    int nincoly = f->nincol;
1076    double *rlos = f->rlos;
1077    double *rups = f->rups;
1078    int *rows = f->rows;
1079
1080    double *coeffxs = f->coeffxs;
1081
1082    int jrowy = f->rowy;
1083
1084    int *ninrowxs = f->ninrowxs;
1085    const int *rowcolsxs = f->rowcolsxs;
1086    const double *rowelsxs = f->rowelsxs;
1087   
1088    /* the row was in the reduced problem */
1089    for (int i=0; i<nincoly; ++i) {
1090      if (rows[i] != jrowy)
1091        PRESOLVEASSERT(rdone[rows[i]]);
1092    }
1093    PRESOLVEASSERT(cdone[icol]==DROP_COL);
1094    PRESOLVEASSERT(rdone[jrowy]==DROP_ROW);
1095
1096    // DEBUG CHECK
1097#if     0 && DEBUG_PRESOLVE
1098    {
1099      double actx = 0.0;
1100      for (int j=0; j<ncols; ++j)
1101        if (hincol[j] > 0 && cdone[j]) {
1102          int krow = presolve_find_row1(jrowx, mcstrt[j], mcstrt[j] + hincol[j], hrow);
1103          if (krow < mcstrt[j] + hincol[j])
1104            actx += colels[krow] * sol[j];
1105      }
1106      if (! (fabs(acts[jrowx] - actx) < 100*ztolzb))
1107        printf("BAD ACTSX:  acts[%d]==%g != %g\n",
1108               jrowx, acts[jrowx], actx);
1109      if (! (rlo[jrowx] - 100*ztolzb <= actx && actx <= rup[jrowx] + 100*ztolzb))
1110        printf("ACTSX NOT IN RANGE:  %d %g %g %g\n",
1111               jrowx, rlo[jrowx], actx, rup[jrowx]);
1112    }
1113#endif
1114
1115    int ninrowy=-1;
1116    const int *rowcolsy=NULL;
1117    const double *rowelsy=NULL;
1118    double coeffy=0.0;
1119
1120    double rloy;
1121    {
1122      int nel = 0;
1123      for (int i=0; i<nincoly; ++i) {
1124        int row = rows[i];
1125        rlo[row] = rlos[i];
1126        rup[row] = rups[i];
1127        if (row == jrowy) {
1128          ninrowy = ninrowxs[i];
1129          rowcolsy = &rowcolsxs[nel];
1130          rowelsy  = &rowelsxs[nel];
1131
1132          coeffy = coeffxs[i];
1133          rloy = rlo[row];
1134        }
1135        nel += ninrowxs[i];
1136      }
1137    }
1138    double rhsy = rloy;
1139
1140    // restore costs
1141    {
1142      const double *costs = f->costsx;
1143      if (costs)
1144        for (int i = 0; i<ninrowy; ++i) {
1145          dcost[rowcolsy[i]] = costs[i];
1146        }
1147    }
1148
1149    // solve for the equality to find the solution for the eliminated col
1150    // this is why we want coeffx < coeffy (55)
1151    {
1152      double sol0 = rloy;
1153      sol[icol] = 0.0;  // to avoid condition in loop
1154      for (int k = 0; k<ninrowy; ++k) {
1155        int jcolx = rowcolsy[k];
1156        double coeffx = rowelsy[k];
1157        sol0 -= coeffx * sol[jcolx];
1158      }
1159      sol[icol] = sol0 / coeffy;
1160
1161      if (! (sol[icol] > clo[icol] - ztolzb &&
1162             cup[icol] + ztolzb > sol[icol]))
1163        printf("NEW SOL OUT-OF-TOL:  %g %g %g\n", clo[icol], sol[icol], cup[icol]);
1164    }
1165
1166    // since this row is fixed
1167    acts[jrowy] = rloy;
1168
1169    // acts[irow] always ok, since slack is fixed
1170    prob->setRowStatus(jrowy,PrePostsolveMatrix::atLowerBound);
1171
1172    // remove old rowx from col rep
1173    // we don't explicitly store what the current rowx is;
1174    // however, after the presolve, rowx contains a col for every
1175    // col in either the original rowx or the original rowy.
1176    // If there were cancellations, those were handled in subsequent
1177    // presolves.
1178    {
1179      // erase those cols in the other rows that occur in rowy
1180      // (with the exception of icol, which was deleted);
1181      // the other rows *must* contain these cols
1182      for (int k = 0; k<ninrowy; ++k) {
1183        int col = rowcolsy[k];
1184
1185        // remove jrowx from col in the col rep
1186        // icol itself was deleted, so won't be there
1187        if (col != icol)
1188          for (int i = 0; i<nincoly; ++i) {
1189            if (rows[i] != jrowy)
1190              presolve_delete_from_row2(col, rows[i], mcstrt, hincol, hrow, colels, link, &free_list);
1191          }
1192      }
1193
1194      // initialize this for loops below
1195      hincol[icol] = 0;
1196
1197      // now restore the original rows (other than rowy).
1198      // those cols that were also in rowy were just removed;
1199      // otherwise, they *must* already be there.
1200      // This loop and the next automatically create the rep for the new col.
1201      {
1202        const int *rowcolsx = rowcolsxs;
1203        const double *rowelsx = rowelsxs;
1204
1205        for (int i = 0; i<nincoly; ++i) {
1206          int ninrowx = ninrowxs[i];
1207          int jrowx = rows[i];
1208
1209          if (jrowx != jrowy)
1210            for (int k = 0; k<ninrowx; ++k) {
1211              int col = rowcolsx[k];
1212              int kcolx = presolve_find_row3(jrowx, mcstrt[col], hincol[col], hrow, link);
1213
1214              if (kcolx != -1) {
1215                PRESOLVEASSERT(presolve_find_row1(col, 0, ninrowy, rowcolsy) == ninrowy);
1216                // overwrite the existing entry
1217                colels[kcolx] = rowelsx[k];
1218              } else {
1219                PRESOLVEASSERT(presolve_find_row1(col, 0, ninrowy, rowcolsy) < ninrowy);
1220
1221                {
1222                  int kk = free_list;
1223                  free_list = link[free_list];
1224                  check_free_list(free_list);
1225
1226                  link[kk] = mcstrt[col];
1227                  mcstrt[col] = kk;
1228                  colels[kk] = rowelsx[k];
1229                  hrow[kk] = jrowx;
1230                }
1231                ++hincol[col];
1232              }
1233            }
1234          rowcolsx += ninrowx;
1235          rowelsx += ninrowx;
1236        }
1237      }
1238
1239      // finally, add original rowy elements
1240      for (int k = 0; k<ninrowy; ++k) {
1241        int col = rowcolsy[k];
1242
1243        {
1244          prepend_elem(col, rowelsy[k], jrowy, mcstrt, colels, hrow, link, &free_list);
1245          ++hincol[col];
1246        }
1247      }
1248    }
1249
1250    // my guess is that the CLAIM in doubleton generalizes to
1251    // equations with more than one x-style variable.
1252    // Since I can't see how to distinguish among them,
1253    // I assume that any of them will do.
1254
1255    {
1256      int k;
1257      double dj = dcost[icol];
1258      double bounds_factor = rhsy/coeffy;
1259      for (int i=0; i<nincoly; ++i)
1260        if (rows[i] != jrowy) {
1261          int row = rows[i];
1262          double coeff = coeffxs[i];
1263
1264          // PROBABLY DOESN'T MAKE SENSE
1265          acts[row] += coeff * bounds_factor;
1266
1267          dj -= rowduals[row] * coeff;
1268        }
1269
1270      // DEBUG CHECK
1271      double acty = 0.0;
1272      for (int k = 0; k<ninrowy; ++k) {
1273        int col = rowcolsy[k];
1274        acty += rowelsy[k] * sol[col];
1275      }
1276      PRESOLVEASSERT(fabs(acty - acts[jrowy]) < 100*ZTOLDP);
1277
1278      // RECOMPUTING
1279      {
1280        const int *rowcolsx = rowcolsxs;
1281        const double *rowelsx = rowelsxs;
1282
1283        for (int i=0; i<nincoly; ++i) {
1284          int ninrowx = ninrowxs[i];
1285
1286          if (rows[i] != jrowy) {
1287            int jrowx = rows[i];
1288
1289            double actx = 0.0;
1290            for (int k = 0; k<ninrowx; ++k) {
1291              int col = rowcolsx[k];
1292              actx += rowelsx[k] * sol[col];
1293            }
1294#if     DEBUG_PRESOLVE
1295            PRESOLVEASSERT(rlo[jrowx] - ztolzb <= actx && actx <= rup[jrowx] + ztolzb);
1296#endif
1297            acts[jrowx] = actx;
1298          }
1299          rowcolsx += ninrowx;
1300          rowelsx += ninrowx;
1301        }
1302      }
1303
1304      // this is the coefficient we need to force col y's reduced cost to 0.0;
1305      // for example, this is obviously true if y is a singleton column
1306      rowduals[jrowy] = dj / coeffy;
1307      rcosts[icol] = 0.0;
1308
1309      // furthermore, this should leave rcosts[colx] for all colx
1310      // in jrowx unchanged (????).
1311    }
1312
1313    // Unlike doubleton, there should never be a problem with keeping
1314    // the reduced costs the way they were, because the other
1315    // variable's bounds are never changed, since col was implied free.
1316    //rowstat[jrowy] = 0;
1317    prob->setColumnStatus(icol,PrePostsolveMatrix::basic);
1318
1319    cdone[icol] = SUBST_ROW;
1320    rdone[jrowy] = SUBST_ROW;
1321  }
1322
1323  prob->free_list_ = free_list;
1324}
1325
1326
1327
1328subst_constraint_action::~subst_constraint_action()
1329{
1330  const action *actions = actions_;
1331
1332  for (int i=0; i<nactions_; ++i) {
1333    delete[]actions[i].rows;
1334    delete[]actions[i].rlos;
1335    delete[]actions[i].rups;
1336    delete[]actions[i].coeffxs;
1337    delete[]actions[i].ninrowxs;
1338    delete[]actions[i].rowcolsxs;
1339    delete[]actions[i].rowelsxs;
1340    delete[]actions[i].costsx;
1341  }
1342
1343  delete[]actions;
1344
1345}
Note: See TracBrowser for help on using the repository browser.