source: trunk/PresolveSubst.cpp @ 173

Last change on this file since 173 was 173, checked in by forrest, 16 years ago

Improvements to presolve

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