source: trunk/PresolveMatrix.cpp @ 109

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

To help OsiSimplexInterface?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.7 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4//#define       CHECK_CONSISTENCY       1
5
6#include <stdio.h>
7
8#include <assert.h>
9#include <iostream>
10
11#include "CoinHelperFunctions.hpp"
12
13#include "CoinPackedMatrix.hpp"
14#include "ClpSimplex.hpp"
15
16#include "PresolveMatrix.hpp"
17#include "Presolve.hpp"
18
19void PresolveAction::throwCoinError(const char *error,
20                                       const char *ps_routine)
21{
22  throw CoinError(error, ps_routine, "Presolve");
23}
24
25
26  class presolvehlink;
27
28void presolve_make_memlists(CoinBigIndex *starts, int *lengths,
29                       presolvehlink *link,
30                       int n);
31
32/*static*/ void presolve_prefix(const int *starts, int *diffs, int n, int limit);
33
34/*static*/ void presolve_prefix(const CoinBigIndex *starts, int *diffs, int n, int limit)
35{
36  int i;
37
38  for (i=0; i<n; i++) {
39    diffs[i] = starts[i+1] - starts[i];
40  }
41}
42
43// afterward, link[i].pre is the largest index less than i st lengths[i]!=0
44//  (or -1 if all such lengths are 0).
45// link[i].suc is the opposite.
46// That is, it is the same thing as setting link[i].pre = i-1 and link[i].suc = i+1
47// and then deleting all the empty elements.
48// This list is maintained together with hrow/hcol so that as we relocate
49// columns or rows, we can quickly determine what column/row precedes a given
50// column/row in the memory region hrow/hcol.
51// Note that n itself has a pre and a suc.
52void presolve_make_memlists(CoinBigIndex *starts, int *lengths, presolvehlink *link, int n)
53{
54  int i;
55  int pre = NO_LINK;
56 
57  for (i=0; i<n; i++) {
58    if (lengths[i]) {
59      link[i].pre = pre;
60      if (pre != NO_LINK)
61        link[pre].suc = i;
62      pre = i;
63    }
64    else {
65      link[i].pre = NO_LINK;
66      link[i].suc = NO_LINK;
67    }
68  }
69  if (pre != NO_LINK)
70    link[pre].suc = n;
71
72  // (1) Arbitrarily place the last non-empty entry in link[n].pre
73  link[n].pre = pre;
74
75  link[n].suc = NO_LINK;
76}
77
78
79
80double *presolve_duparray(const double *d, int n, int n2)
81{
82  double *d1 = new double[n2];
83  memcpy(d1, d, n*sizeof(double));
84  return (d1);
85}
86double *presolve_duparray(const double *d, int n)
87{
88  return presolve_duparray(d, n, n);
89}
90
91int *presolve_duparray(const int *d, int n, int n2)
92{
93  int *d1 = new int[n2];
94  memcpy(d1, d, n*sizeof(int));
95  return (d1);
96}
97int *presolve_duparray(const int *d, int n)
98{
99  return presolve_duparray(d, n, n);
100}
101
102char *presolve_duparray(const char *d, int n, int n2)
103{
104  char *d1 = new char[n2];
105  memcpy(d1, d, n*sizeof(char));
106  return (d1);
107}
108char *presolve_duparray(const char *d, int n)
109{
110  return presolve_duparray(d, n, n);
111}
112
113#if 0
114double *presolve_duparray(const double *d, int n, char **end_mmapp)
115{
116  double *d1 = (double*)*end_mmapp;
117  memcpy(d1, d, n*sizeof(double));
118  *end_mmapp += ALIGN_DOUBLE(n*sizeof(double));
119  return (d1);
120}
121int *presolve_duparray(const int *d, int n, char **end_mmapp)
122{
123  int *d1 = (int*)*end_mmapp;
124  memcpy(d1, d, n*sizeof(int));
125  *end_mmapp += ALIGN_DOUBLE(n*sizeof(int));
126  return (d1);
127}
128
129double *presolve_duparray(const double *d, int n, int n2, char **end_mmapp)
130{
131  double *d1 = (double*)*end_mmapp;
132  memcpy(d1, d, n*sizeof(double));
133  *end_mmapp += ALIGN_DOUBLE(n2*sizeof(double));
134  return (d1);
135}
136int *presolve_duparray(const int *d, int n, int n2, char **end_mmapp)
137{
138  int *d1 = (int*)*end_mmapp;
139  memcpy(d1, d, n*sizeof(int));
140  *end_mmapp += ALIGN_DOUBLE(n2*sizeof(int));
141  return (d1);
142}
143#endif
144
145
146int presolve_find_row(int row, CoinBigIndex kcs, CoinBigIndex kce, const int *hrow)
147{
148  CoinBigIndex k;
149  for (k=kcs; k<kce; k++)
150    if (hrow[k] == row)
151      return (k);
152  DIE("FIND_ROW");
153  return (0);
154}
155
156int presolve_find_row1(int row, CoinBigIndex kcs, CoinBigIndex kce, const int *hrow)
157{
158  CoinBigIndex k;
159  for (k=kcs; k<kce; k++)
160    if (hrow[k] == row)
161      return (k);
162  return (kce);
163}
164
165int presolve_find_row2(int irow, CoinBigIndex ks, int nc, const int *hrow, const int *link)
166{
167  for (int i=0; i<nc; ++i) {
168    if (hrow[ks] == irow)
169      return (ks);
170    ks = link[ks];
171  }
172  abort();
173  return -1;
174}
175
176int presolve_find_row3(int irow, CoinBigIndex ks, int nc, const int *hrow, const int *link)
177{
178  for (int i=0; i<nc; ++i) {
179    if (hrow[ks] == irow)
180      return (ks);
181    ks = link[ks];
182  }
183  return (-1);
184}
185
186
187// delete the entry for col from row
188// this keeps the row loosely packed
189// if we didn't want to maintain that property, we could just decrement hinrow[row].
190void presolve_delete_from_row(int row, int col /* thing to delete */,
191                     const CoinBigIndex *mrstrt,
192                     int *hinrow, int *hcol, double *dels)
193{
194  CoinBigIndex krs = mrstrt[row];
195  CoinBigIndex kre = krs + hinrow[row];
196
197  CoinBigIndex kcol = presolve_find_row(col, krs, kre, hcol);
198
199  swap(hcol[kcol], hcol[kre-1]);
200  swap(dels[kcol], dels[kre-1]);
201  hinrow[row]--;
202}
203
204
205void presolve_delete_from_row2(int row, int col /* thing to delete */,
206                      CoinBigIndex *mrstrt,
207                     int *hinrow, int *hcol, double *dels, int *link, 
208                               CoinBigIndex *free_listp)
209{
210  CoinBigIndex k = mrstrt[row];
211
212  if (hcol[k] == col) {
213    mrstrt[row] = link[k];
214    link[k] = *free_listp;
215    *free_listp = k;
216    check_free_list(k);
217    hinrow[row]--;
218  } else {
219    int n = hinrow[row] - 1;
220    CoinBigIndex k0 = k;
221    k = link[k];
222    for (int i=0; i<n; ++i) {
223      if (hcol[k] == col) {
224        link[k0] = link[k];
225        link[k] = *free_listp;
226        *free_listp = k;
227        check_free_list(k);
228        hinrow[row]--;
229        return;
230      }
231      k0 = k;
232      k = link[k];
233    }
234    abort();
235  }
236}
237
238static inline double getTolerance(const ClpSimplex &si, ClpDblParam key)
239{
240  double tol;
241  if (! si.getDblParam(key, tol)) {
242    PresolveAction::throwCoinError("getDblParam failed",
243                                      "PrePostsolveMatrix::PrePostsolveMatrix");
244  }
245  return (tol);
246}
247
248#if 0
249PresolveAction::~PresolveAction()
250{
251  cout << "OOPS\n";
252}
253#endif
254
255// Assumptions:
256// 1. nrows>=m.getNumRows()
257// 2. ncols>=m.getNumCols()
258//
259// In presolve, these values are equal.
260// In postsolve, they may be inequal, since the reduced problem
261// may be smaller, but we need room for the large problem.
262// ncols may be larger than si.getNumCols() in postsolve,
263// this at that point si will be the reduced problem,
264// but we need to reserve enough space for the original problem.
265PrePostsolveMatrix::PrePostsolveMatrix(const ClpSimplex& si,
266                                             int ncols_in,
267                                             int nrows_in,
268                                             CoinBigIndex nelems_in) :
269  ncols_(si.getNumCols()),
270  ncols0_(ncols_in),
271  nelems_(si.getNumElements()),
272
273  mcstrt_(new CoinBigIndex[ncols_in+1]),
274  hincol_(new int[ncols_in+1]),
275  hrow_  (new int   [2*nelems_in]),
276  colels_(new double[2*nelems_in]),
277
278  cost_(new double[ncols_in]),
279  clo_(new double[ncols_in]),
280  cup_(new double[ncols_in]),
281  rlo_(new double[nrows_in]),
282  rup_(new double[nrows_in]),
283  originalColumn_(new int[ncols_in]),
284
285  ztolzb_(getTolerance(si, ClpPrimalTolerance)),
286  ztoldj_(getTolerance(si, ClpDualTolerance)),
287
288  maxmin_(si.getObjSense())
289
290{
291  si.getDblParam(ClpObjOffset,originalOffset_);
292  int ncols = si.getNumCols();
293  int nrows = si.getNumRows();
294
295  ClpDisjointCopyN(si.getColLower(), ncols, clo_);
296  ClpDisjointCopyN(si.getColUpper(), ncols, cup_);
297  ClpDisjointCopyN(si.getObjCoefficients(), ncols, cost_);
298  ClpDisjointCopyN(si.getRowLower(), nrows,  rlo_);
299  ClpDisjointCopyN(si.getRowUpper(), nrows,  rup_);
300  int i;
301  for (i=0;i<ncols_in;i++) 
302    originalColumn_[i]=i;
303  sol_=NULL;
304  rowduals_=NULL;
305  acts_=NULL;
306
307  rcosts_=NULL;
308  colstat_=NULL;
309  rowstat_=NULL;
310}
311
312
313PrePostsolveMatrix::~PrePostsolveMatrix()
314{
315  delete[]mcstrt_;
316  delete[]hrow_;
317  delete[]colels_;
318  delete[]hincol_;
319
320  delete[]cost_;
321  delete[]clo_;
322  delete[]cup_;
323  delete[]rlo_;
324  delete[]rup_;
325  delete[]originalColumn_;
326  delete[]rowduals_;
327
328  delete[]rcosts_;
329}
330
331// Sets status (non -basic ) using value
332void 
333PrePostsolveMatrix::setRowStatusUsingValue(int iRow)
334{
335  double value = acts_[iRow];
336  double lower = rlo_[iRow];
337  double upper = rup_[iRow];
338  if (lower<-1.0e20&&upper>1.0e20) {
339    setRowStatus(iRow,isFree);
340  } else if (fabs(lower-value)<=ztolzb_) {
341    setRowStatus(iRow,atLowerBound);
342  } else if (fabs(upper-value)<=ztolzb_) {
343    setRowStatus(iRow,atUpperBound);
344  } else {
345    setRowStatus(iRow,superBasic);
346  }
347}
348void 
349PrePostsolveMatrix::setColumnStatusUsingValue(int iColumn)
350{
351  double value = sol_[iColumn];
352  double lower = clo_[iColumn];
353  double upper = cup_[iColumn];
354  if (lower<-1.0e20&&upper>1.0e20) {
355    setColumnStatus(iColumn,isFree);
356  } else if (fabs(lower-value)<=ztolzb_) {
357    setColumnStatus(iColumn,atLowerBound);
358  } else if (fabs(upper-value)<=ztolzb_) {
359    setColumnStatus(iColumn,atUpperBound);
360  } else {
361    setColumnStatus(iColumn,superBasic);
362  }
363}
364
365// I am not familiar enough with CoinPackedMatrix to be confident
366// that I will implement a row-ordered version of toColumnOrderedGapFree
367// properly.
368static bool isGapFree(const CoinPackedMatrix& matrix)
369{
370  const CoinBigIndex * start = matrix.getVectorStarts();
371  const int * length = matrix.getVectorLengths();
372  int i;
373  for (i = matrix.getSizeVectorLengths() - 1; i >= 0; --i) {
374    if (start[i+1] - start[i] != length[i])
375      break;
376  }
377  return (! (i >= 0));
378}
379
380
381#if     DEBUG_PRESOLVE
382static void matrix_bounds_ok(const double *lo, const double *up, int n)
383{
384  int i;
385  for (i=0; i<n; i++) {
386    PRESOLVEASSERT(lo[i] <= up[i]);
387    PRESOLVEASSERT(lo[i] < PRESOLVE_INF);
388    PRESOLVEASSERT(-PRESOLVE_INF < up[i]);
389  }
390}
391#endif
392
393PresolveMatrix::PresolveMatrix(int ncols0_in,
394                                     double maxmin_,
395                                     // end prepost members
396
397                                     ClpSimplex &si,
398
399                                     // rowrep
400                                     int nrows_in,
401                                     CoinBigIndex nelems_in,
402                               bool doStatus) :
403
404  PrePostsolveMatrix(si,
405                        ncols0_in, nrows_in, nelems_in),
406  clink_(new presolvehlink[ncols0_in+1]),
407  rlink_(new presolvehlink[nrows_in+1]),
408
409  dobias_(0.0),
410
411  nrows_(si.getNumRows()),
412
413  // temporary init
414  mrstrt_(new CoinBigIndex[nrows_in+1]),
415  hinrow_(new int[nrows_in+1]),
416  rowels_(new double[2*nelems_in]),
417  hcol_(new int[2*nelems_in]),
418  integerType_(new char[ncols0_in]),
419  feasibilityTolerance_(0.0),
420  status_(-1),
421  rowsToDo_(new int [nrows_in]),
422  numberRowsToDo_(0),
423  nextRowsToDo_(new int[nrows_in]),
424  numberNextRowsToDo_(0),
425  colsToDo_(new int [ncols0_in]),
426  numberColsToDo_(0),
427  nextColsToDo_(new int[ncols0_in]),
428  numberNextColsToDo_(0)
429
430{
431  const int bufsize = 2*nelems_in;
432
433  // Set up change bits
434  rowChanged_ = new unsigned int[(nrows_+31)>>5];
435  memset(rowChanged_,0,((nrows_+31)>>5)*sizeof(unsigned int));
436  colChanged_ = new unsigned int[(ncols_+31)>>5];
437  memset(colChanged_,0,((ncols_+31)>>5)*sizeof(unsigned int));
438  CoinPackedMatrix * m = si.matrix();
439
440  // The coefficient matrix is a big hunk of stuff.
441  // Do the copy here to try to avoid running out of memory.
442
443  const CoinBigIndex * start = m->getVectorStarts();
444  const int * length = m->getVectorLengths();
445  const int * row = m->getIndices();
446  const double * element = m->getElements();
447  int icol,nel=0;
448  mcstrt_[0]=0;
449  for (icol=0;icol<ncols_;icol++) {
450    int j;
451    for (j=start[icol];j<start[icol]+length[icol];j++) {
452      hrow_[nel]=row[j];
453      colels_[nel++]=element[j];
454    }
455    mcstrt_[icol+1]=nel;
456  }
457  assert(mcstrt_[ncols_] == nelems_);
458  ClpDisjointCopyN(m->getVectorLengths(),ncols_,  hincol_);
459
460  // same thing for row rep
461  m = new CoinPackedMatrix();
462  m->reverseOrderedCopyOf(*si.matrix());
463  m->removeGaps();
464
465
466  ClpDisjointCopyN(m->getVectorStarts(),  nrows_,  mrstrt_);
467  mrstrt_[nrows_] = nelems_;
468  ClpDisjointCopyN(m->getVectorLengths(), nrows_,  hinrow_);
469  ClpDisjointCopyN(m->getIndices(),       nelems_, hcol_);
470  ClpDisjointCopyN(m->getElements(),      nelems_, rowels_);
471
472  delete m;
473  if (si.integerInformation()) {
474    memcpy(integerType_,si.integerInformation(),ncols_*sizeof(char));
475  } else {
476    ClpFillN<char>(integerType_, ncols_, 0);
477  }
478
479  if (doStatus) {
480    // allow for status and solution
481    sol_ = new double[ncols_];
482    memcpy(sol_,si.primalColumnSolution(),ncols_*sizeof(double));;
483    acts_ = new double [nrows_];
484    memcpy(acts_,si.primalRowSolution(),nrows_*sizeof(double));
485    if (!si.statusArray())
486      si.createStatus();
487    colstat_ = new unsigned char [nrows_+ncols_];
488    memcpy(colstat_,si.statusArray(),
489           (nrows_+ncols_)*sizeof(unsigned char));
490    rowstat_ = colstat_+ncols_;
491  }
492
493  // the original model's fields are now unneeded - free them
494 
495  si.resize(0,0);
496
497#if     DEBUG_PRESOLVE
498  matrix_bounds_ok(rlo_, rup_, nrows_);
499  matrix_bounds_ok(clo_, cup_, ncols_);
500#endif
501
502#if 0
503  for (i=0; i<nrows; ++i)
504    printf("NR: %6d\n", hinrow[i]);
505  for (int i=0; i<ncols; ++i)
506    printf("NC: %6d\n", hincol[i]);
507#endif
508
509  presolve_make_memlists(mcstrt_, hincol_, clink_, ncols_);
510  presolve_make_memlists(mrstrt_, hinrow_, rlink_, nrows_);
511
512  // this allows last col/row to expand up to bufsize-1 (22);
513  // this must come after the calls to presolve_prefix
514  mcstrt_[ncols_] = bufsize-1;
515  mrstrt_[nrows_] = bufsize-1;
516
517#if     CHECK_CONSISTENCY
518  consistent(false);
519#endif
520}
521
522PresolveMatrix::~PresolveMatrix()
523{
524  delete[]clink_;
525  delete[]rlink_;
526 
527  delete[]mrstrt_;
528  delete[]hinrow_;
529  delete[]rowels_;
530  delete[]hcol_;
531
532  delete[]integerType_;
533  delete[] rowChanged_;
534  delete[] rowsToDo_;
535  delete[] nextRowsToDo_;
536  delete[] colChanged_;
537  delete[] colsToDo_;
538  delete[] nextColsToDo_;
539}
540
541
542void PresolveMatrix::update_model(ClpSimplex& si,
543                                     int nrows0,
544                                     int ncols0,
545                                     CoinBigIndex nelems0)
546{
547#if 0
548  // out as it was hardware error !!
549  {
550    // temporary sanity check (for funny nw04 bug)
551    int icol;
552    for (icol=0;icol<ncols_;icol++) {
553      int j;
554      for (j=mcstrt_[icol];j<mcstrt_[icol]+hincol_[icol];j++) {
555        assert(hrow_[j]>=0&&hrow_[j]<nrows_);
556      }
557    }
558  }
559#endif
560  si.loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
561                 clo_, cup_, cost_, rlo_, rup_);
562
563  delete [] si.integerInformation();
564  int numberIntegers=0;
565  for (int i=0; i<ncols_; i++) {
566    if (integerType_[i])
567      numberIntegers++;
568  }
569  if (numberIntegers) 
570    si.copyInIntegerInformation(integerType_);
571  else
572    si.copyInIntegerInformation(NULL);
573
574#if     PRESOLVE_SUMMARY
575  printf("NEW NCOL/NROW/NELS:  %d(-%d) %d(-%d) %d(-%d)\n",
576         ncols_, ncols0-ncols_,
577         nrows_, nrows0-nrows_,
578         si.getNumElements(), nelems0-si.getNumElements());
579#endif
580  si.setDblParam(ClpObjOffset,originalOffset_-dobias_);
581
582}
583
584#if     CHECK_CONSISTENCY
585
586// The matrix is represented redundantly in both row and column format,
587// in what I call "loosely packed" format.
588// "packed" format would be as in normal OSL:  a vector of column starts,
589// together with two vectors that give the row indices and coefficients.
590//
591// This format is "loosely packed" because some of the elements may correspond
592// to empty rows.  This is so that we can quickly delete rows without having
593// to update the column rep and vice versa.
594//
595// Checks whether an entry is in the col rep iff it is also in the row rep,
596// and also that their values are the same (if testvals is non-zero).
597//
598// Note that there may be entries in a row that correspond to empty columns
599// and vice-versa.  --- HUH???
600static void matrix_consistent(const CoinBigIndex *mrstrt, const int *hinrow, const int *hcol,
601                              const CoinBigIndex *mcstrt, const int *hincol, const int *hrow,
602                              const double *rowels,
603                              const double *colels,
604                              int nrows, int testvals,
605                              const char *ROW, const char *COL)
606{
607  for (int irow=0; irow<nrows; irow++) {
608    if (hinrow[irow] > 0) {
609      CoinBigIndex krs = mrstrt[irow];
610      CoinBigIndex kre = krs + hinrow[irow];
611
612      for (CoinBigIndex k=krs; k<kre; k++) {
613        int jcol = hcol[k];
614        CoinBigIndex kcs = mcstrt[jcol];
615        CoinBigIndex kce = kcs + hincol[jcol];
616
617        CoinBigIndex kk = presolve_find_row1(irow, kcs, kce, hrow);
618        if (kk == kce) {
619          printf("MATRIX INCONSISTENT:  can't find %s %d in %s %d\n",
620                 ROW, irow, COL, jcol);
621          fflush(stdout);
622          abort();
623        }
624        if (testvals && colels[kk] != rowels[k]) {
625          printf("MATRIX INCONSISTENT:  value for %s %d and %s %d\n",
626                 ROW, irow, COL, jcol);
627          fflush(stdout);
628          abort();
629        }
630      }
631    }
632  }
633}
634#endif
635
636
637void PresolveMatrix::consistent(bool testvals)
638{
639#if     CHECK_CONSISTENCY
640  matrix_consistent(mrstrt_, hinrow_, hcol_,
641                    mcstrt_, hincol_, hrow_,
642                    rowels_, colels_,
643                    nrows_, testvals,
644                    "row", "col");
645  matrix_consistent(mcstrt_, hincol_, hrow_,
646                    mrstrt_, hinrow_, hcol_,
647                    colels_, rowels_, 
648                    ncols_, testvals,
649                    "col", "row");
650#endif
651}
652
653
654
655
656
657
658
659
660
661
662
663////////////////  POSTSOLVE
664
665PostsolveMatrix::PostsolveMatrix(ClpSimplex& si,
666                                       int ncols0_in,
667                                       int nrows0_in,
668                                       CoinBigIndex nelems0,
669                                   
670                                       double maxmin_,
671                                       // end prepost members
672
673                                       double *sol_in,
674                                       double *acts_in,
675
676                                       unsigned char *colstat_in,
677                                       unsigned char *rowstat_in) :
678  PrePostsolveMatrix(si,
679                        ncols0_in, nrows0_in, nelems0),
680
681  free_list_(0),
682  // link, free_list, maxlink
683  maxlink_(2*nelems0),
684  link_(new int[/*maxlink*/ 2*nelems0]),
685     
686  cdone_(new char[ncols0_]),
687  rdone_(new char[nrows0_in]),
688
689  nrows_(si.getNumRows()),
690  nrows0_(nrows0_in)
691{
692
693  sol_=sol_in;
694  rowduals_=NULL;
695  acts_=acts_in;
696
697  rcosts_=NULL;
698  colstat_=colstat_in;
699  rowstat_=rowstat_in;
700
701  // this is the *reduced* model, which is probably smaller
702  int ncols1 = si.getNumCols();
703  int nrows1 = si.getNumRows();
704
705  const CoinPackedMatrix * m = si.matrix();
706
707  if (! isGapFree(*m)) {
708    PresolveAction::throwCoinError("Matrix not gap free",
709                                      "PostsolveMatrix");
710  }
711
712  const CoinBigIndex nelemsr = m->getNumElements();
713
714  ClpDisjointCopyN(m->getVectorStarts(), ncols1, mcstrt_);
715  mcstrt_[ncols_] = nelems0;    // ??
716  ClpDisjointCopyN(m->getVectorLengths(),ncols1,  hincol_);
717  ClpDisjointCopyN(m->getIndices(),      nelemsr, hrow_);
718  ClpDisjointCopyN(m->getElements(),     nelemsr, colels_);
719
720
721#if     0 && DEBUG_PRESOLVE
722  presolve_check_costs(model, &colcopy);
723#endif
724
725  // This determines the size of the data structure that contains
726  // the matrix being postsolved.  Links are taken from the free_list
727  // to recreate matrix entries that were presolved away,
728  // and links are added to the free_list when entries created during
729  // presolve are discarded.  There is never a need to gc this list.
730  // Naturally, it should contain
731  // exactly nelems0 entries "in use" when postsolving is done,
732  // but I don't know whether the matrix could temporarily get
733  // larger during postsolving.  Substitution into more than two
734  // rows could do that, in principle.  I am being very conservative
735  // here by reserving much more than the amount of space I probably need.
736  // If this guess is wrong, check_free_list may be called.
737  //  int bufsize = 2*nelems0;
738
739  memset(cdone_, -1, ncols0_);
740  memset(rdone_, -1, nrows0_);
741
742  rowduals_ = new double[nrows0_];
743  ClpDisjointCopyN(si.getRowPrice(), nrows1, rowduals_);
744
745  rcosts_ = new double[ncols0_];
746  ClpDisjointCopyN(si.getReducedCost(), ncols1, rcosts_);
747  if (maxmin_<0.0) {
748    // change so will look as if minimize
749    int i;
750    for (i=0;i<nrows1;i++)
751      rowduals_[i] = - rowduals_[i];
752    for (i=0;i<ncols1;i++) {
753      rcosts_[i] = - rcosts_[i];
754    }
755  }
756
757  //ClpDisjointCopyN(si.getRowUpper(), nrows1, rup_);
758  //ClpDisjointCopyN(si.getRowLower(), nrows1, rlo_);
759
760  ClpDisjointCopyN(si.getColSolution(), ncols1, sol_);
761  si.setDblParam(ClpObjOffset,originalOffset_);
762
763  for (int j=0; j<ncols1; j++) {
764    CoinBigIndex kcs = mcstrt_[j];
765    CoinBigIndex kce = kcs + hincol_[j];
766    for (CoinBigIndex k=kcs; k<kce; ++k) {
767      link_[k] = k+1;
768    }
769  }
770  {
771    int ml = maxlink_;
772    for (CoinBigIndex k=nelemsr; k<ml; ++k)
773      link_[k] = k+1;
774    link_[ml-1] = NO_LINK;
775  }
776  free_list_ = nelemsr;
777}
778
779PostsolveMatrix::~PostsolveMatrix()
780{
781  delete[]link_;
782
783  delete[]cdone_;
784  delete[]rdone_;
785}
786
787
788void PostsolveMatrix::check_nbasic()
789{
790  int nbasic = 0;
791
792  int i;
793  for (i=0; i<ncols_; i++)
794    if (columnIsBasic(i))
795      nbasic++;
796
797  for (i=0; i<nrows_; i++)
798    if (rowIsBasic(i))
799      nbasic++;
800
801  if (nbasic != nrows_) {
802    printf("WRONG NUMBER NBASIC:  is:  %d  should be:  %d\n",
803           nbasic, nrows_);
804    fflush(stdout);
805  }
806}
807
808
809
810
811
812
Note: See TracBrowser for help on using the repository browser.