source: trunk/Clp/src/ClpGubDynamicMatrix.cpp @ 2470

Last change on this file since 2470 was 2385, checked in by unxusr, 9 months ago

formatting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 76.8 KB
Line 
1/* $Id: ClpGubDynamicMatrix.cpp 2385 2019-01-06 19:43:06Z stefan $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include <cstdio>
7
8#include "CoinPragma.hpp"
9#include "CoinIndexedVector.hpp"
10#include "CoinHelperFunctions.hpp"
11
12#include "ClpSimplex.hpp"
13#include "ClpFactorization.hpp"
14#include "ClpQuadraticObjective.hpp"
15#include "ClpNonLinearCost.hpp"
16// at end to get min/max!
17#include "ClpGubDynamicMatrix.hpp"
18#include "ClpMessage.hpp"
19//#define CLP_DEBUG
20//#define CLP_DEBUG_PRINT
21//#############################################################################
22// Constructors / Destructor / Assignment
23//#############################################################################
24
25//-------------------------------------------------------------------
26// Default Constructor
27//-------------------------------------------------------------------
28ClpGubDynamicMatrix::ClpGubDynamicMatrix()
29  : ClpGubMatrix()
30  , objectiveOffset_(0.0)
31  , startColumn_(NULL)
32  , row_(NULL)
33  , element_(NULL)
34  , cost_(NULL)
35  , fullStart_(NULL)
36  , id_(NULL)
37  , dynamicStatus_(NULL)
38  , lowerColumn_(NULL)
39  , upperColumn_(NULL)
40  , lowerSet_(NULL)
41  , upperSet_(NULL)
42  , numberGubColumns_(0)
43  , firstAvailable_(0)
44  , savedFirstAvailable_(0)
45  , firstDynamic_(0)
46  , lastDynamic_(0)
47  , numberElements_(0)
48{
49  setType(13);
50}
51
52//-------------------------------------------------------------------
53// Copy constructor
54//-------------------------------------------------------------------
55ClpGubDynamicMatrix::ClpGubDynamicMatrix(const ClpGubDynamicMatrix &rhs)
56  : ClpGubMatrix(rhs)
57{
58  objectiveOffset_ = rhs.objectiveOffset_;
59  numberGubColumns_ = rhs.numberGubColumns_;
60  firstAvailable_ = rhs.firstAvailable_;
61  savedFirstAvailable_ = rhs.savedFirstAvailable_;
62  firstDynamic_ = rhs.firstDynamic_;
63  lastDynamic_ = rhs.lastDynamic_;
64  numberElements_ = rhs.numberElements_;
65  startColumn_ = ClpCopyOfArray(rhs.startColumn_, numberGubColumns_ + 1);
66  CoinBigIndex numberElements = startColumn_[numberGubColumns_];
67  row_ = ClpCopyOfArray(rhs.row_, numberElements);
68  ;
69  element_ = ClpCopyOfArray(rhs.element_, numberElements);
70  ;
71  cost_ = ClpCopyOfArray(rhs.cost_, numberGubColumns_);
72  fullStart_ = ClpCopyOfArray(rhs.fullStart_, numberSets_ + 1);
73  id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_);
74  lowerColumn_ = ClpCopyOfArray(rhs.lowerColumn_, numberGubColumns_);
75  upperColumn_ = ClpCopyOfArray(rhs.upperColumn_, numberGubColumns_);
76  dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, numberGubColumns_);
77  lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
78  upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
79}
80
81/* This is the real constructor*/
82ClpGubDynamicMatrix::ClpGubDynamicMatrix(ClpSimplex *model, int numberSets,
83  int numberGubColumns, const int *starts,
84  const double *lower, const double *upper,
85  const CoinBigIndex *startColumn, const int *row,
86  const double *element, const double *cost,
87  const double *lowerColumn, const double *upperColumn,
88  const unsigned char *status)
89  : ClpGubMatrix()
90{
91  objectiveOffset_ = model->objectiveOffset();
92  model_ = model;
93  numberSets_ = numberSets;
94  numberGubColumns_ = numberGubColumns;
95  fullStart_ = ClpCopyOfArray(starts, numberSets_ + 1);
96  lower_ = ClpCopyOfArray(lower, numberSets_);
97  upper_ = ClpCopyOfArray(upper, numberSets_);
98  int numberColumns = model->numberColumns();
99  int numberRows = model->numberRows();
100  // Number of columns needed
101  int numberGubInSmall = numberSets_ + numberRows + 2 * model->factorizationFrequency() + 2;
102  // for small problems this could be too big
103  //numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_);
104  int numberNeeded = numberGubInSmall + numberColumns;
105  firstAvailable_ = numberColumns;
106  savedFirstAvailable_ = numberColumns;
107  firstDynamic_ = numberColumns;
108  lastDynamic_ = numberNeeded;
109  startColumn_ = ClpCopyOfArray(startColumn, numberGubColumns_ + 1);
110  CoinBigIndex numberElements = startColumn_[numberGubColumns_];
111  row_ = ClpCopyOfArray(row, numberElements);
112  element_ = new double[numberElements];
113  CoinBigIndex i;
114  for (i = 0; i < numberElements; i++)
115    element_[i] = element[i];
116  cost_ = new double[numberGubColumns_];
117  for (i = 0; i < numberGubColumns_; i++) {
118    cost_[i] = cost[i];
119    // need sorted
120    CoinSort_2(row_ + startColumn_[i], row_ + startColumn_[i + 1], element_ + startColumn_[i]);
121  }
122  if (lowerColumn) {
123    lowerColumn_ = new double[numberGubColumns_];
124    for (i = 0; i < numberGubColumns_; i++)
125      lowerColumn_[i] = lowerColumn[i];
126  } else {
127    lowerColumn_ = NULL;
128  }
129  if (upperColumn) {
130    upperColumn_ = new double[numberGubColumns_];
131    for (i = 0; i < numberGubColumns_; i++)
132      upperColumn_[i] = upperColumn[i];
133  } else {
134    upperColumn_ = NULL;
135  }
136  if (upperColumn || lowerColumn) {
137    lowerSet_ = new double[numberSets_];
138    for (i = 0; i < numberSets_; i++) {
139      if (lower[i] > -1.0e20)
140        lowerSet_[i] = lower[i];
141      else
142        lowerSet_[i] = -1.0e30;
143    }
144    upperSet_ = new double[numberSets_];
145    for (i = 0; i < numberSets_; i++) {
146      if (upper[i] < 1.0e20)
147        upperSet_[i] = upper[i];
148      else
149        upperSet_[i] = 1.0e30;
150    }
151  } else {
152    lowerSet_ = NULL;
153    upperSet_ = NULL;
154  }
155  start_ = NULL;
156  end_ = NULL;
157  dynamicStatus_ = NULL;
158  id_ = new int[numberGubInSmall];
159  for (i = 0; i < numberGubInSmall; i++)
160    id_[i] = -1;
161  ClpPackedMatrix *originalMatrixA = dynamic_cast< ClpPackedMatrix * >(model->clpMatrix());
162  assert(originalMatrixA);
163  CoinPackedMatrix *originalMatrix = originalMatrixA->getPackedMatrix();
164  originalMatrixA->setMatrixNull(); // so can be deleted safely
165  // guess how much space needed
166  double guess = originalMatrix->getNumElements() + 10;
167  guess /= static_cast< double >(numberColumns);
168  guess *= 2 * numberGubColumns_;
169  numberElements_ = static_cast< int >(CoinMin(guess, 10000000.0));
170  numberElements_ = CoinMin(numberElements_, numberElements) + originalMatrix->getNumElements();
171  matrix_ = originalMatrix;
172  flags_ &= ~1;
173  // resize model (matrix stays same)
174  model->resize(numberRows, numberNeeded);
175  if (upperColumn_) {
176    // set all upper bounds so we have enough space
177    double *columnUpper = model->columnUpper();
178    for (i = firstDynamic_; i < lastDynamic_; i++)
179      columnUpper[i] = 1.0e10;
180  }
181  // resize matrix
182  // extra 1 is so can keep number of elements handy
183  originalMatrix->reserve(numberNeeded, numberElements_, true);
184  originalMatrix->reserve(numberNeeded + 1, numberElements_, false);
185  originalMatrix->getMutableVectorStarts()[numberColumns] = originalMatrix->getNumElements();
186  // redo number of columns
187  numberColumns = matrix_->getNumCols();
188  backward_ = new int[numberNeeded];
189  backToPivotRow_ = new int[numberNeeded];
190  // We know a bit better
191  delete[] changeCost_;
192  changeCost_ = new double[numberRows + numberSets_];
193  keyVariable_ = new int[numberSets_];
194  // signal to need new ordering
195  next_ = NULL;
196  for (int iColumn = 0; iColumn < numberNeeded; iColumn++)
197    backward_[iColumn] = -1;
198
199  firstGub_ = firstDynamic_;
200  lastGub_ = lastDynamic_;
201  if (!lowerColumn_ && !upperColumn_)
202    gubType_ = 8;
203  if (status) {
204    status_ = ClpCopyOfArray(status, numberSets_);
205  } else {
206    status_ = new unsigned char[numberSets_];
207    memset(status_, 0, numberSets_);
208    int i;
209    for (i = 0; i < numberSets_; i++) {
210      // make slack key
211      setStatus(i, ClpSimplex::basic);
212    }
213  }
214  saveStatus_ = new unsigned char[numberSets_];
215  memset(saveStatus_, 0, numberSets_);
216  savedKeyVariable_ = new int[numberSets_];
217  memset(savedKeyVariable_, 0, numberSets_ * sizeof(int));
218}
219
220//-------------------------------------------------------------------
221// Destructor
222//-------------------------------------------------------------------
223ClpGubDynamicMatrix::~ClpGubDynamicMatrix()
224{
225  delete[] startColumn_;
226  delete[] row_;
227  delete[] element_;
228  delete[] cost_;
229  delete[] fullStart_;
230  delete[] id_;
231  delete[] dynamicStatus_;
232  delete[] lowerColumn_;
233  delete[] upperColumn_;
234  delete[] lowerSet_;
235  delete[] upperSet_;
236}
237
238//----------------------------------------------------------------
239// Assignment operator
240//-------------------------------------------------------------------
241ClpGubDynamicMatrix &
242ClpGubDynamicMatrix::operator=(const ClpGubDynamicMatrix &rhs)
243{
244  if (this != &rhs) {
245    ClpGubMatrix::operator=(rhs);
246    delete[] startColumn_;
247    delete[] row_;
248    delete[] element_;
249    delete[] cost_;
250    delete[] fullStart_;
251    delete[] id_;
252    delete[] dynamicStatus_;
253    delete[] lowerColumn_;
254    delete[] upperColumn_;
255    delete[] lowerSet_;
256    delete[] upperSet_;
257    objectiveOffset_ = rhs.objectiveOffset_;
258    numberGubColumns_ = rhs.numberGubColumns_;
259    firstAvailable_ = rhs.firstAvailable_;
260    savedFirstAvailable_ = rhs.savedFirstAvailable_;
261    firstDynamic_ = rhs.firstDynamic_;
262    lastDynamic_ = rhs.lastDynamic_;
263    numberElements_ = rhs.numberElements_;
264    startColumn_ = ClpCopyOfArray(rhs.startColumn_, numberGubColumns_ + 1);
265    int numberElements = startColumn_[numberGubColumns_];
266    row_ = ClpCopyOfArray(rhs.row_, numberElements);
267    ;
268    element_ = ClpCopyOfArray(rhs.element_, numberElements);
269    ;
270    cost_ = ClpCopyOfArray(rhs.cost_, numberGubColumns_);
271    fullStart_ = ClpCopyOfArray(rhs.fullStart_, numberSets_ + 1);
272    id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_);
273    lowerColumn_ = ClpCopyOfArray(rhs.lowerColumn_, numberGubColumns_);
274    upperColumn_ = ClpCopyOfArray(rhs.upperColumn_, numberGubColumns_);
275    dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, numberGubColumns_);
276    lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
277    upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
278  }
279  return *this;
280}
281//-------------------------------------------------------------------
282// Clone
283//-------------------------------------------------------------------
284ClpMatrixBase *ClpGubDynamicMatrix::clone() const
285{
286  return new ClpGubDynamicMatrix(*this);
287}
288// Partial pricing
289void ClpGubDynamicMatrix::partialPricing(ClpSimplex *model, double startFraction, double endFraction,
290  int &bestSequence, int &numberWanted)
291{
292  assert(!model->rowScale());
293  numberWanted = currentWanted_;
294  if (!numberSets_) {
295    // no gub
296    ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
297    return;
298  } else {
299    // and do some proportion of full set
300    int startG2 = static_cast< int >(startFraction * numberSets_);
301    int endG2 = static_cast< int >(endFraction * numberSets_ + 0.1);
302    endG2 = CoinMin(endG2, numberSets_);
303    //printf("gub price - set start %d end %d\n",
304    //   startG2,endG2);
305    double tolerance = model->currentDualTolerance();
306    double *reducedCost = model->djRegion();
307    const double *duals = model->dualRowSolution();
308    double *cost = model->costRegion();
309    double bestDj;
310    int numberRows = model->numberRows();
311    int numberColumns = lastDynamic_;
312    // If nothing found yet can go all the way to end
313    int endAll = endG2;
314    if (bestSequence < 0 && !startG2)
315      endAll = numberSets_;
316    if (bestSequence >= 0)
317      bestDj = fabs(reducedCost[bestSequence]);
318    else
319      bestDj = tolerance;
320    int saveSequence = bestSequence;
321    double djMod = 0.0;
322    double infeasibilityCost = model->infeasibilityCost();
323    double bestDjMod = 0.0;
324    //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
325    //     startG2,endG2,numberWanted);
326    int bestType = -1;
327    int bestSet = -1;
328    const double *element = matrix_->getElements();
329    const int *row = matrix_->getIndices();
330    const CoinBigIndex *startColumn = matrix_->getVectorStarts();
331    int *length = matrix_->getMutableVectorLengths();
332#if 0
333          // make sure first available is clean (in case last iteration rejected)
334          cost[firstAvailable_] = 0.0;
335          length[firstAvailable_] = 0;
336          model->nonLinearCost()->setOne(firstAvailable_, 0.0, 0.0, COIN_DBL_MAX, 0.0);
337          model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
338          {
339               for (int i = firstAvailable_; i < lastDynamic_; i++)
340                    assert(!cost[i]);
341          }
342#endif
343#ifdef CLP_DEBUG
344    {
345      for (int i = firstDynamic_; i < firstAvailable_; i++) {
346        assert(getDynamicStatus(id_[i - firstDynamic_]) == inSmall);
347      }
348    }
349#endif
350    int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_;
351    int minNeg = minimumGoodReducedCosts_ < 0 ? 5 : minimumGoodReducedCosts_;
352    for (int iSet = startG2; iSet < endAll; iSet++) {
353      if (numberWanted + minNeg < originalWanted_ && iSet > startG2 + minSet) {
354        // give up
355        numberWanted = 0;
356        break;
357      } else if (iSet == endG2 && bestSequence >= 0) {
358        break;
359      }
360      CoinBigIndex j;
361      int iBasic = keyVariable_[iSet];
362      if (iBasic >= numberColumns) {
363        djMod = -weight(iSet) * infeasibilityCost;
364      } else {
365        // get dj without
366        assert(model->getStatus(iBasic) == ClpSimplex::basic);
367        djMod = 0.0;
368
369        for (j = startColumn[iBasic];
370             j < startColumn[iBasic] + length[iBasic]; j++) {
371          int jRow = row[j];
372          djMod -= duals[jRow] * element[j];
373        }
374        djMod += cost[iBasic];
375        // See if gub slack possible - dj is djMod
376        if (getStatus(iSet) == ClpSimplex::atLowerBound) {
377          double value = -djMod;
378          if (value > tolerance) {
379            numberWanted--;
380            if (value > bestDj) {
381              // check flagged variable and correct dj
382              if (!flagged(iSet)) {
383                bestDj = value;
384                bestSequence = numberRows + numberColumns + iSet;
385                bestDjMod = djMod;
386                bestType = 0;
387                bestSet = iSet;
388              } else {
389                // just to make sure we don't exit before got something
390                numberWanted++;
391                abort();
392              }
393            }
394          }
395        } else if (getStatus(iSet) == ClpSimplex::atUpperBound) {
396          double value = djMod;
397          if (value > tolerance) {
398            numberWanted--;
399            if (value > bestDj) {
400              // check flagged variable and correct dj
401              if (!flagged(iSet)) {
402                bestDj = value;
403                bestSequence = numberRows + numberColumns + iSet;
404                bestDjMod = djMod;
405                bestType = 0;
406                bestSet = iSet;
407              } else {
408                // just to make sure we don't exit before got something
409                numberWanted++;
410                abort();
411              }
412            }
413          }
414        }
415      }
416      for (int iSequence = fullStart_[iSet]; iSequence < fullStart_[iSet + 1]; iSequence++) {
417        DynamicStatus status = getDynamicStatus(iSequence);
418        if (status != inSmall) {
419          double value = cost_[iSequence] - djMod;
420          for (j = startColumn_[iSequence];
421               j < startColumn_[iSequence + 1]; j++) {
422            int jRow = row_[j];
423            value -= duals[jRow] * element_[j];
424          }
425          // change sign if at lower bound
426          if (status == atLowerBound)
427            value = -value;
428          if (value > tolerance) {
429            numberWanted--;
430            if (value > bestDj) {
431              // check flagged variable and correct dj
432              if (!flagged(iSequence)) {
433                bestDj = value;
434                bestSequence = iSequence;
435                bestDjMod = djMod;
436                bestType = 1;
437                bestSet = iSet;
438              } else {
439                // just to make sure we don't exit before got something
440                numberWanted++;
441              }
442            }
443          }
444        }
445      }
446      if (numberWanted <= 0) {
447        numberWanted = 0;
448        break;
449      }
450    }
451    // Do packed part before gub and small gub - but lightly
452    int saveMinNeg = minimumGoodReducedCosts_;
453    int saveSequence2 = bestSequence;
454    if (bestSequence >= 0)
455      minimumGoodReducedCosts_ = -2;
456    int saveLast = lastGub_;
457    lastGub_ = firstAvailable_;
458    currentWanted_ = numberWanted;
459    ClpGubMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
460    minimumGoodReducedCosts_ = saveMinNeg;
461    lastGub_ = saveLast;
462    if (bestSequence != saveSequence2) {
463      bestType = -1; // in normal or small gub part
464      saveSequence = bestSequence;
465    }
466    if (bestSequence != saveSequence || bestType >= 0) {
467      double *lowerColumn = model->lowerRegion();
468      double *upperColumn = model->upperRegion();
469      double *solution = model->solutionRegion();
470      if (bestType > 0) {
471        // recompute dj and create
472        double value = cost_[bestSequence] - bestDjMod;
473        for (CoinBigIndex jBigIndex = startColumn_[bestSequence];
474             jBigIndex < startColumn_[bestSequence + 1]; jBigIndex++) {
475          int jRow = row_[jBigIndex];
476          value -= duals[jRow] * element_[jBigIndex];
477        }
478        double *element = matrix_->getMutableElements();
479        int *row = matrix_->getMutableIndices();
480        CoinBigIndex *startColumn = matrix_->getMutableVectorStarts();
481        int *length = matrix_->getMutableVectorLengths();
482        CoinBigIndex numberElements = startColumn[firstAvailable_];
483        int numberThis = startColumn_[bestSequence + 1] - startColumn_[bestSequence];
484        if (numberElements + numberThis > numberElements_) {
485          // need to redo
486          numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
487          matrix_->reserve(numberColumns, numberElements_);
488          element = matrix_->getMutableElements();
489          row = matrix_->getMutableIndices();
490          // these probably okay but be safe
491          startColumn = matrix_->getMutableVectorStarts();
492          length = matrix_->getMutableVectorLengths();
493        }
494        // already set startColumn[firstAvailable_]=numberElements;
495        length[firstAvailable_] = numberThis;
496        model->costRegion()[firstAvailable_] = cost_[bestSequence];
497        CoinBigIndex base = startColumn_[bestSequence];
498        for (int j = 0; j < numberThis; j++) {
499          row[numberElements] = row_[base + j];
500          element[numberElements++] = element_[base + j];
501        }
502        id_[firstAvailable_ - firstDynamic_] = bestSequence;
503        //printf("best %d\n",bestSequence);
504        backward_[firstAvailable_] = bestSet;
505        model->solutionRegion()[firstAvailable_] = 0.0;
506        if (!lowerColumn_ && !upperColumn_) {
507          model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
508          lowerColumn[firstAvailable_] = 0.0;
509          upperColumn[firstAvailable_] = COIN_DBL_MAX;
510        } else {
511          DynamicStatus status = getDynamicStatus(bestSequence);
512          if (lowerColumn_)
513            lowerColumn[firstAvailable_] = lowerColumn_[bestSequence];
514          else
515            lowerColumn[firstAvailable_] = 0.0;
516          if (upperColumn_)
517            upperColumn[firstAvailable_] = upperColumn_[bestSequence];
518          else
519            upperColumn[firstAvailable_] = COIN_DBL_MAX;
520          if (status == atLowerBound) {
521            solution[firstAvailable_] = lowerColumn[firstAvailable_];
522            model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
523          } else {
524            solution[firstAvailable_] = upperColumn[firstAvailable_];
525            model->setStatus(firstAvailable_, ClpSimplex::atUpperBound);
526          }
527        }
528        model->nonLinearCost()->setOne(firstAvailable_, solution[firstAvailable_],
529          lowerColumn[firstAvailable_],
530          upperColumn[firstAvailable_], cost_[bestSequence]);
531        bestSequence = firstAvailable_;
532        // firstAvailable_ only updated if good pivot (in updatePivot)
533        startColumn[firstAvailable_ + 1] = numberElements;
534        //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,bestDjMod);
535        reducedCost[bestSequence] = value;
536        gubSlackIn_ = -1;
537      } else {
538        // slack - make last column
539        gubSlackIn_ = bestSequence - numberRows - numberColumns;
540        bestSequence = numberColumns + 2 * numberRows;
541        reducedCost[bestSequence] = bestDjMod;
542        //printf("price slack %d - gubpi %g\n",gubSlackIn_,bestDjMod);
543        model->setStatus(bestSequence, getStatus(gubSlackIn_));
544        if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound)
545          solution[bestSequence] = upper_[gubSlackIn_];
546        else
547          solution[bestSequence] = lower_[gubSlackIn_];
548        lowerColumn[bestSequence] = lower_[gubSlackIn_];
549        upperColumn[bestSequence] = upper_[gubSlackIn_];
550        model->costRegion()[bestSequence] = 0.0;
551        model->nonLinearCost()->setOne(bestSequence, solution[bestSequence], lowerColumn[bestSequence],
552          upperColumn[bestSequence], 0.0);
553      }
554      savedBestSequence_ = bestSequence;
555      savedBestDj_ = reducedCost[savedBestSequence_];
556    }
557    // See if may be finished
558    if (!startG2 && bestSequence < 0)
559      infeasibilityWeight_ = model_->infeasibilityCost();
560    else if (bestSequence >= 0)
561      infeasibilityWeight_ = -1.0;
562  }
563  currentWanted_ = numberWanted;
564}
565// This is local to Gub to allow synchronization when status is good
566int ClpGubDynamicMatrix::synchronize(ClpSimplex *model, int mode)
567{
568  int returnNumber = 0;
569  switch (mode) {
570  case 0: {
571#ifdef CLP_DEBUG
572    {
573      for (int i = 0; i < numberSets_; i++)
574        assert(toIndex_[i] == -1);
575    }
576#endif
577    // lookup array
578    int *lookup = new int[lastDynamic_];
579    int iColumn;
580    int numberColumns = model->numberColumns();
581    double *element = matrix_->getMutableElements();
582    int *row = matrix_->getMutableIndices();
583    CoinBigIndex *startColumn = matrix_->getMutableVectorStarts();
584    int *length = matrix_->getMutableVectorLengths();
585    double *cost = model->costRegion();
586    double *lowerColumn = model->lowerRegion();
587    double *upperColumn = model->upperRegion();
588    int *pivotVariable = model->pivotVariable();
589    CoinBigIndex numberElements = startColumn[firstDynamic_];
590    // first just do lookup and basic stuff
591    int currentNumber = firstAvailable_;
592    firstAvailable_ = firstDynamic_;
593    int numberToDo = 0;
594    double objectiveChange = 0.0;
595    double *solution = model->solutionRegion();
596    for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
597      int iSet = backward_[iColumn];
598      if (toIndex_[iSet] < 0) {
599        toIndex_[iSet] = 0;
600        fromIndex_[numberToDo++] = iSet;
601      }
602      if (model->getStatus(iColumn) == ClpSimplex::basic || iColumn == keyVariable_[iSet]) {
603        lookup[iColumn] = firstAvailable_;
604        if (iColumn != keyVariable_[iSet]) {
605          int iPivot = backToPivotRow_[iColumn];
606          backToPivotRow_[firstAvailable_] = iPivot;
607          pivotVariable[iPivot] = firstAvailable_;
608        }
609        firstAvailable_++;
610      } else {
611        int jColumn = id_[iColumn - firstDynamic_];
612        setDynamicStatus(jColumn, atLowerBound);
613        if (lowerColumn_ || upperColumn_) {
614          if (model->getStatus(iColumn) == ClpSimplex::atUpperBound)
615            setDynamicStatus(jColumn, atUpperBound);
616          // treat solution as if exactly at a bound
617          double value = solution[iColumn];
618          if (fabs(value - lowerColumn[iColumn]) < fabs(value - upperColumn[iColumn]))
619            value = lowerColumn[iColumn];
620          else
621            value = upperColumn[iColumn];
622          objectiveChange += cost[iColumn] * value;
623          // redo lower and upper on sets
624          double shift = value;
625          if (lowerSet_[iSet] > -1.0e20)
626            lower_[iSet] = lowerSet_[iSet] - shift;
627          if (upperSet_[iSet] < 1.0e20)
628            upper_[iSet] = upperSet_[iSet] - shift;
629        }
630        lookup[iColumn] = -1;
631      }
632    }
633    model->setObjectiveOffset(model->objectiveOffset() + objectiveChange);
634    firstAvailable_ = firstDynamic_;
635    for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
636      if (lookup[iColumn] >= 0) {
637        // move
638        int jColumn = id_[iColumn - firstDynamic_];
639        id_[firstAvailable_ - firstDynamic_] = jColumn;
640        int numberThis = startColumn_[jColumn + 1] - startColumn_[jColumn];
641        length[firstAvailable_] = numberThis;
642        cost[firstAvailable_] = cost[iColumn];
643        lowerColumn[firstAvailable_] = lowerColumn[iColumn];
644        upperColumn[firstAvailable_] = upperColumn[iColumn];
645        double originalLower = lowerColumn_ ? lowerColumn_[jColumn] : 0.0;
646        double originalUpper = upperColumn_ ? upperColumn_[jColumn] : COIN_DBL_MAX;
647        if (originalUpper > 1.0e30)
648          originalUpper = COIN_DBL_MAX;
649        model->nonLinearCost()->setOne(firstAvailable_, solution[iColumn],
650          originalLower, originalUpper,
651          cost_[jColumn]);
652        CoinBigIndex base = startColumn_[jColumn];
653        for (int j = 0; j < numberThis; j++) {
654          row[numberElements] = row_[base + j];
655          element[numberElements++] = element_[base + j];
656        }
657        model->setStatus(firstAvailable_, model->getStatus(iColumn));
658        backward_[firstAvailable_] = backward_[iColumn];
659        solution[firstAvailable_] = solution[iColumn];
660        firstAvailable_++;
661        startColumn[firstAvailable_] = numberElements;
662      }
663    }
664    // clean up next_
665    int *temp = new int[firstAvailable_];
666    for (int jSet = 0; jSet < numberToDo; jSet++) {
667      int iSet = fromIndex_[jSet];
668      toIndex_[iSet] = -1;
669      int last = keyVariable_[iSet];
670      int j = next_[last];
671      bool setTemp = true;
672      if (last < lastDynamic_) {
673        last = lookup[last];
674        assert(last >= 0);
675        keyVariable_[iSet] = last;
676      } else if (j >= 0) {
677        int newJ = lookup[j];
678        assert(newJ >= 0);
679        j = next_[j];
680        next_[last] = newJ;
681        last = newJ;
682      } else {
683        next_[last] = -(iSet + numberColumns + 1);
684        setTemp = false;
685      }
686      while (j >= 0) {
687        int newJ = lookup[j];
688        assert(newJ >= 0);
689        temp[last] = newJ;
690        last = newJ;
691        j = next_[j];
692      }
693      if (setTemp)
694        temp[last] = -(keyVariable_[iSet] + 1);
695      if (lowerSet_) {
696        // we only need to get lower_ and upper_ correct
697        double shift = 0.0;
698        for (int j = fullStart_[iSet]; j < fullStart_[iSet + 1]; j++)
699          if (getDynamicStatus(j) == atUpperBound)
700            shift += upperColumn_[j];
701          else if (getDynamicStatus(j) == atLowerBound && lowerColumn_)
702            shift += lowerColumn_[j];
703        if (lowerSet_[iSet] > -1.0e20)
704          lower_[iSet] = lowerSet_[iSet] - shift;
705        if (upperSet_[iSet] < 1.0e20)
706          upper_[iSet] = upperSet_[iSet] - shift;
707      }
708    }
709    // move to next_
710    CoinMemcpyN(temp + firstDynamic_, (firstAvailable_ - firstDynamic_), next_ + firstDynamic_);
711    // if odd iterations may be one out so adjust currentNumber
712    currentNumber = CoinMin(currentNumber + 1, lastDynamic_);
713    // zero solution
714    CoinZeroN(solution + firstAvailable_, currentNumber - firstAvailable_);
715    // zero cost
716    CoinZeroN(cost + firstAvailable_, currentNumber - firstAvailable_);
717    // zero lengths
718    CoinZeroN(length + firstAvailable_, currentNumber - firstAvailable_);
719    for (iColumn = firstAvailable_; iColumn < currentNumber; iColumn++) {
720      model->nonLinearCost()->setOne(iColumn, 0.0, 0.0, COIN_DBL_MAX, 0.0);
721      model->setStatus(iColumn, ClpSimplex::atLowerBound);
722      backward_[iColumn] = -1;
723    }
724    delete[] lookup;
725    delete[] temp;
726    // make sure fromIndex clean
727    fromIndex_[0] = -1;
728    //#define CLP_DEBUG
729#ifdef CLP_DEBUG
730    // debug
731    {
732      int i;
733      int numberRows = model->numberRows();
734      char *xxxx = new char[numberColumns];
735      memset(xxxx, 0, numberColumns);
736      for (i = 0; i < numberRows; i++) {
737        int iPivot = pivotVariable[i];
738        assert(model->getStatus(iPivot) == ClpSimplex::basic);
739        if (iPivot < numberColumns && backward_[iPivot] >= 0)
740          xxxx[iPivot] = 1;
741      }
742      for (i = 0; i < numberSets_; i++) {
743        int key = keyVariable_[i];
744        int iColumn = next_[key];
745        int k = 0;
746        while (iColumn >= 0) {
747          k++;
748          assert(k < 100);
749          assert(backward_[iColumn] == i);
750          iColumn = next_[iColumn];
751        }
752        int stop = -(key + 1);
753        while (iColumn != stop) {
754          assert(iColumn < 0);
755          iColumn = -iColumn - 1;
756          k++;
757          assert(k < 100);
758          assert(backward_[iColumn] == i);
759          iColumn = next_[iColumn];
760        }
761        iColumn = next_[key];
762        while (iColumn >= 0) {
763          assert(xxxx[iColumn]);
764          xxxx[iColumn] = 0;
765          iColumn = next_[iColumn];
766        }
767      }
768      for (i = 0; i < numberColumns; i++) {
769        if (i < numberColumns && backward_[i] >= 0) {
770          assert(!xxxx[i] || i == keyVariable_[backward_[i]]);
771        }
772      }
773      delete[] xxxx;
774    }
775    {
776      for (int i = 0; i < numberSets_; i++)
777        assert(toIndex_[i] == -1);
778    }
779#endif
780    savedFirstAvailable_ = firstAvailable_;
781  } break;
782  // flag a variable
783  case 1: {
784    // id will be sitting at firstAvailable
785    int sequence = id_[firstAvailable_ - firstDynamic_];
786    assert(!flagged(sequence));
787    setFlagged(sequence);
788    model->clearFlagged(firstAvailable_);
789  } break;
790  // unflag all variables
791  case 2: {
792    for (int i = 0; i < numberGubColumns_; i++) {
793      if (flagged(i)) {
794        unsetFlagged(i);
795        returnNumber++;
796      }
797    }
798  } break;
799  //  just reset costs and bounds (primal)
800  case 3: {
801    double *cost = model->costRegion();
802    double *solution = model->solutionRegion();
803    double *lowerColumn = model->columnLower();
804    double *upperColumn = model->columnUpper();
805    for (int i = firstDynamic_; i < firstAvailable_; i++) {
806      int jColumn = id_[i - firstDynamic_];
807      cost[i] = cost_[jColumn];
808      if (!lowerColumn_ && !upperColumn_) {
809        lowerColumn[i] = 0.0;
810        upperColumn[i] = COIN_DBL_MAX;
811      } else {
812        if (lowerColumn_)
813          lowerColumn[i] = lowerColumn_[jColumn];
814        else
815          lowerColumn[i] = 0.0;
816        if (upperColumn_)
817          upperColumn[i] = upperColumn_[jColumn];
818        else
819          upperColumn[i] = COIN_DBL_MAX;
820      }
821      if (model->nonLinearCost())
822        model->nonLinearCost()->setOne(i, solution[i],
823          lowerColumn[i],
824          upperColumn[i], cost_[jColumn]);
825    }
826    if (!model->numberIterations() && rhsOffset_) {
827      lastRefresh_ = -refreshFrequency_; // force refresh
828    }
829  } break;
830  // and get statistics for column generation
831  case 4: {
832    // In theory we should subtract out ones we have done but ....
833    // If key slack then dual 0.0
834    // If not then slack could be dual infeasible
835    // dj for key is zero so that defines dual on set
836    int i;
837    int numberColumns = model->numberColumns();
838    double *dual = model->dualRowSolution();
839    double infeasibilityCost = model->infeasibilityCost();
840    double dualTolerance = model->dualTolerance();
841    double relaxedTolerance = dualTolerance;
842    // we can't really trust infeasibilities if there is dual error
843    double error = CoinMin(1.0e-2, model->largestDualError());
844    // allow tolerance at least slightly bigger than standard
845    relaxedTolerance = relaxedTolerance + error;
846    // but we will be using difference
847    relaxedTolerance -= dualTolerance;
848    double objectiveOffset = 0.0;
849    for (i = 0; i < numberSets_; i++) {
850      int kColumn = keyVariable_[i];
851      double value = 0.0;
852      if (kColumn < numberColumns) {
853        kColumn = id_[kColumn - firstDynamic_];
854        // dj without set
855        value = cost_[kColumn];
856        for (CoinBigIndex j = startColumn_[kColumn];
857             j < startColumn_[kColumn + 1]; j++) {
858          int iRow = row_[j];
859          value -= dual[iRow] * element_[j];
860        }
861        double infeasibility = 0.0;
862        if (getStatus(i) == ClpSimplex::atLowerBound) {
863          if (-value > dualTolerance)
864            infeasibility = -value - dualTolerance;
865        } else if (getStatus(i) == ClpSimplex::atUpperBound) {
866          if (value > dualTolerance)
867            infeasibility = value - dualTolerance;
868        }
869        if (infeasibility > 0.0) {
870          sumDualInfeasibilities_ += infeasibility;
871          if (infeasibility > relaxedTolerance)
872            sumOfRelaxedDualInfeasibilities_ += infeasibility;
873          numberDualInfeasibilities_++;
874        }
875      } else {
876        // slack key - may not be feasible
877        assert(getStatus(i) == ClpSimplex::basic);
878        // negative as -1.0 for slack
879        value = -weight(i) * infeasibilityCost;
880      }
881      // Now subtract out from all
882      for (CoinBigIndex k = fullStart_[i]; k < fullStart_[i + 1]; k++) {
883        if (getDynamicStatus(k) != inSmall) {
884          double djValue = cost_[k] - value;
885          for (CoinBigIndex j = startColumn_[k];
886               j < startColumn_[k + 1]; j++) {
887            int iRow = row_[j];
888            djValue -= dual[iRow] * element_[j];
889          }
890          double infeasibility = 0.0;
891          double shift = 0.0;
892          if (getDynamicStatus(k) == atLowerBound) {
893            if (lowerColumn_)
894              shift = lowerColumn_[k];
895            if (djValue < -dualTolerance)
896              infeasibility = -djValue - dualTolerance;
897          } else {
898            // at upper bound
899            shift = upperColumn_[k];
900            if (djValue > dualTolerance)
901              infeasibility = djValue - dualTolerance;
902          }
903          objectiveOffset += shift * cost_[k];
904          if (infeasibility > 0.0) {
905            sumDualInfeasibilities_ += infeasibility;
906            if (infeasibility > relaxedTolerance)
907              sumOfRelaxedDualInfeasibilities_ += infeasibility;
908            numberDualInfeasibilities_++;
909          }
910        }
911      }
912    }
913    model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
914  } break;
915  // see if time to re-factorize
916  case 5: {
917    if (firstAvailable_ > numberSets_ + model->numberRows() + model->factorizationFrequency())
918      returnNumber = 4;
919  } break;
920  // return 1 if there may be changing bounds on variable (column generation)
921  case 6: {
922    returnNumber = (lowerColumn_ != NULL || upperColumn_ != NULL) ? 1 : 0;
923#if 0
924          if (!returnNumber) {
925               // may be gub slacks
926               for (int i = 0; i < numberSets_; i++) {
927                    if (upper_[i] > lower_[i]) {
928                         returnNumber = 1;
929                         break;
930                    }
931               }
932          }
933#endif
934  } break;
935  // restore firstAvailable_
936  case 7: {
937    int iColumn;
938    int *length = matrix_->getMutableVectorLengths();
939    double *cost = model->costRegion();
940    double *solution = model->solutionRegion();
941    int currentNumber = firstAvailable_;
942    firstAvailable_ = savedFirstAvailable_;
943    // zero solution
944    CoinZeroN(solution + firstAvailable_, currentNumber - firstAvailable_);
945    // zero cost
946    CoinZeroN(cost + firstAvailable_, currentNumber - firstAvailable_);
947    // zero lengths
948    CoinZeroN(length + firstAvailable_, currentNumber - firstAvailable_);
949    for (iColumn = firstAvailable_; iColumn < currentNumber; iColumn++) {
950      model->nonLinearCost()->setOne(iColumn, 0.0, 0.0, COIN_DBL_MAX, 0.0);
951      model->setStatus(iColumn, ClpSimplex::atLowerBound);
952      backward_[iColumn] = -1;
953    }
954  } break;
955  // make sure set is clean
956  case 8: {
957    int sequenceIn = model->sequenceIn();
958    if (sequenceIn < model->numberColumns()) {
959      int iSet = backward_[sequenceIn];
960      if (iSet >= 0 && lowerSet_) {
961        // we only need to get lower_ and upper_ correct
962        double shift = 0.0;
963        for (int j = fullStart_[iSet]; j < fullStart_[iSet + 1]; j++)
964          if (getDynamicStatus(j) == atUpperBound)
965            shift += upperColumn_[j];
966          else if (getDynamicStatus(j) == atLowerBound && lowerColumn_)
967            shift += lowerColumn_[j];
968        if (lowerSet_[iSet] > -1.0e20)
969          lower_[iSet] = lowerSet_[iSet] - shift;
970        if (upperSet_[iSet] < 1.0e20)
971          upper_[iSet] = upperSet_[iSet] - shift;
972      }
973      if (sequenceIn == firstAvailable_) {
974        // not really in small problem
975        int iBig = id_[sequenceIn - firstDynamic_];
976        if (model->getStatus(sequenceIn) == ClpSimplex::atLowerBound)
977          setDynamicStatus(iBig, atLowerBound);
978        else
979          setDynamicStatus(iBig, atUpperBound);
980      }
981    }
982  } break;
983  // adjust lower,upper
984  case 9: {
985    int sequenceIn = model->sequenceIn();
986    if (sequenceIn >= firstDynamic_ && sequenceIn < lastDynamic_ && lowerSet_) {
987      int iSet = backward_[sequenceIn];
988      assert(iSet >= 0);
989      int inBig = id_[sequenceIn - firstDynamic_];
990      const double *solution = model->solutionRegion();
991      setDynamicStatus(inBig, inSmall);
992      if (lowerSet_[iSet] > -1.0e20)
993        lower_[iSet] += solution[sequenceIn];
994      if (upperSet_[iSet] < 1.0e20)
995        upper_[iSet] += solution[sequenceIn];
996      model->setObjectiveOffset(model->objectiveOffset() - solution[sequenceIn] * cost_[inBig]);
997    }
998  }
999  }
1000  return returnNumber;
1001}
1002// Add a new variable to a set
1003void ClpGubDynamicMatrix::insertNonBasic(int sequence, int iSet)
1004{
1005  int last = keyVariable_[iSet];
1006  int j = next_[last];
1007  while (j >= 0) {
1008    last = j;
1009    j = next_[j];
1010  }
1011  next_[last] = -(sequence + 1);
1012  next_[sequence] = j;
1013}
1014// Sets up an effective RHS and does gub crash if needed
1015void ClpGubDynamicMatrix::useEffectiveRhs(ClpSimplex *model, bool cheapest)
1016{
1017  // Do basis - cheapest or slack if feasible (unless cheapest set)
1018  int longestSet = 0;
1019  int iSet;
1020  for (iSet = 0; iSet < numberSets_; iSet++)
1021    longestSet = CoinMax(longestSet, fullStart_[iSet + 1] - fullStart_[iSet]);
1022
1023  double *upper = new double[longestSet + 1];
1024  double *cost = new double[longestSet + 1];
1025  double *lower = new double[longestSet + 1];
1026  double *solution = new double[longestSet + 1];
1027  assert(!next_);
1028  delete[] next_;
1029  int numberColumns = model->numberColumns();
1030  next_ = new int[numberColumns + numberSets_ + CoinMax(2 * longestSet, lastDynamic_ - firstDynamic_)];
1031  char *mark = new char[numberColumns];
1032  memset(mark, 0, numberColumns);
1033  for (int iColumn = 0; iColumn < numberColumns; iColumn++)
1034    next_[iColumn] = COIN_INT_MAX;
1035  int i;
1036  int *keys = new int[numberSets_];
1037  int *back = new int[numberGubColumns_];
1038  CoinFillN(back, numberGubColumns_, -1);
1039  for (i = 0; i < numberSets_; i++)
1040    keys[i] = COIN_INT_MAX;
1041  delete[] dynamicStatus_;
1042  dynamicStatus_ = new unsigned char[numberGubColumns_];
1043  memset(dynamicStatus_, 0, numberGubColumns_); // for clarity
1044  for (i = 0; i < numberGubColumns_; i++)
1045    setDynamicStatus(i, atLowerBound);
1046  // set up chains
1047  for (i = firstDynamic_; i < lastDynamic_; i++) {
1048    if (id_[i - firstDynamic_] >= 0) {
1049      if (model->getStatus(i) == ClpSimplex::basic)
1050        mark[i] = 1;
1051      int iSet = backward_[i];
1052      assert(iSet >= 0);
1053      int iNext = keys[iSet];
1054      next_[i] = iNext;
1055      keys[iSet] = i;
1056      back[id_[i - firstDynamic_]] = i;
1057    } else {
1058      model->setStatus(i, ClpSimplex::atLowerBound);
1059      backward_[i] = -1;
1060    }
1061  }
1062  double *columnSolution = model->solutionRegion();
1063  int numberRows = getNumRows();
1064  toIndex_ = new int[numberSets_];
1065  for (iSet = 0; iSet < numberSets_; iSet++)
1066    toIndex_[iSet] = -1;
1067  fromIndex_ = new int[numberRows + numberSets_];
1068  double tolerance = model->primalTolerance();
1069  double *element = matrix_->getMutableElements();
1070  int *row = matrix_->getMutableIndices();
1071  CoinBigIndex *startColumn = matrix_->getMutableVectorStarts();
1072  int *length = matrix_->getMutableVectorLengths();
1073  double objectiveOffset = 0.0;
1074  for (iSet = 0; iSet < numberSets_; iSet++) {
1075    int j;
1076    int numberBasic = 0;
1077    int iBasic = -1;
1078    int iStart = fullStart_[iSet];
1079    int iEnd = fullStart_[iSet + 1];
1080    // find one with smallest length
1081    int smallest = numberRows + 1;
1082    double value = 0.0;
1083    j = keys[iSet];
1084    while (j != COIN_INT_MAX) {
1085      if (model->getStatus(j) == ClpSimplex::basic) {
1086        if (length[j] < smallest) {
1087          smallest = length[j];
1088          iBasic = j;
1089        }
1090        numberBasic++;
1091      }
1092      value += columnSolution[j];
1093      j = next_[j];
1094    }
1095    bool done = false;
1096    if (numberBasic > 1 || (numberBasic == 1 && getStatus(iSet) == ClpSimplex::basic)) {
1097      if (getStatus(iSet) == ClpSimplex::basic)
1098        iBasic = iSet + numberColumns; // slack key - use
1099      done = true;
1100    } else if (numberBasic == 1) {
1101      // see if can be key
1102      double thisSolution = columnSolution[iBasic];
1103      if (thisSolution < 0.0) {
1104        value -= thisSolution;
1105        thisSolution = 0.0;
1106        columnSolution[iBasic] = thisSolution;
1107      }
1108      // try setting slack to a bound
1109      assert(upper_[iSet] < 1.0e20 || lower_[iSet] > -1.0e20);
1110      double cost1 = COIN_DBL_MAX;
1111      int whichBound = -1;
1112      if (upper_[iSet] < 1.0e20) {
1113        // try slack at ub
1114        double newBasic = thisSolution + upper_[iSet] - value;
1115        if (newBasic >= -tolerance) {
1116          // can go
1117          whichBound = 1;
1118          cost1 = newBasic * cost_[iBasic];
1119          // But if exact then may be good solution
1120          if (fabs(upper_[iSet] - value) < tolerance)
1121            cost1 = -COIN_DBL_MAX;
1122        }
1123      }
1124      if (lower_[iSet] > -1.0e20) {
1125        // try slack at lb
1126        double newBasic = thisSolution + lower_[iSet] - value;
1127        if (newBasic >= -tolerance) {
1128          // can go but is it cheaper
1129          double cost2 = newBasic * cost_[iBasic];
1130          // But if exact then may be good solution
1131          if (fabs(lower_[iSet] - value) < tolerance)
1132            cost2 = -COIN_DBL_MAX;
1133          if (cost2 < cost1)
1134            whichBound = 0;
1135        }
1136      }
1137      if (whichBound != -1) {
1138        // key
1139        done = true;
1140        if (whichBound) {
1141          // slack to upper
1142          columnSolution[iBasic] = thisSolution + upper_[iSet] - value;
1143          setStatus(iSet, ClpSimplex::atUpperBound);
1144        } else {
1145          // slack to lower
1146          columnSolution[iBasic] = thisSolution + lower_[iSet] - value;
1147          setStatus(iSet, ClpSimplex::atLowerBound);
1148        }
1149      }
1150    }
1151    if (!done) {
1152      if (!cheapest) {
1153        // see if slack can be key
1154        if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
1155          done = true;
1156          setStatus(iSet, ClpSimplex::basic);
1157          iBasic = iSet + numberColumns;
1158        }
1159      }
1160      if (!done) {
1161        // set non basic if there was one
1162        if (iBasic >= 0)
1163          model->setStatus(iBasic, ClpSimplex::atLowerBound);
1164        // find cheapest
1165        int numberInSet = iEnd - iStart;
1166        if (!lowerColumn_) {
1167          CoinZeroN(lower, numberInSet);
1168        } else {
1169          for (int j = 0; j < numberInSet; j++)
1170            lower[j] = lowerColumn_[j + iStart];
1171        }
1172        if (!upperColumn_) {
1173          CoinFillN(upper, numberInSet, COIN_DBL_MAX);
1174        } else {
1175          for (int j = 0; j < numberInSet; j++)
1176            upper[j] = upperColumn_[j + iStart];
1177        }
1178        CoinFillN(solution, numberInSet, 0.0);
1179        // and slack
1180        iBasic = numberInSet;
1181        solution[iBasic] = -value;
1182        lower[iBasic] = -upper_[iSet];
1183        upper[iBasic] = -lower_[iSet];
1184        int kphase;
1185        if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
1186          // feasible
1187          kphase = 1;
1188          cost[iBasic] = 0.0;
1189          for (int j = 0; j < numberInSet; j++)
1190            cost[j] = cost_[j + iStart];
1191        } else {
1192          // infeasible
1193          kphase = 0;
1194          // remember bounds are flipped so opposite to natural
1195          if (value < lower_[iSet] - tolerance)
1196            cost[iBasic] = 1.0;
1197          else
1198            cost[iBasic] = -1.0;
1199          CoinZeroN(cost, numberInSet);
1200        }
1201        double dualTolerance = model->dualTolerance();
1202        for (int iphase = kphase; iphase < 2; iphase++) {
1203          if (iphase) {
1204            cost[numberInSet] = 0.0;
1205            for (int j = 0; j < numberInSet; j++)
1206              cost[j] = cost_[j + iStart];
1207          }
1208          // now do one row lp
1209          bool improve = true;
1210          while (improve) {
1211            improve = false;
1212            double dual = cost[iBasic];
1213            int chosen = -1;
1214            double best = dualTolerance;
1215            int way = 0;
1216            for (int i = 0; i <= numberInSet; i++) {
1217              double dj = cost[i] - dual;
1218              double improvement = 0.0;
1219              if (iphase || i < numberInSet)
1220                assert(solution[i] >= lower[i] && solution[i] <= upper[i]);
1221              if (dj > dualTolerance)
1222                improvement = dj * (solution[i] - lower[i]);
1223              else if (dj < -dualTolerance)
1224                improvement = dj * (solution[i] - upper[i]);
1225              if (improvement > best) {
1226                best = improvement;
1227                chosen = i;
1228                if (dj < 0.0) {
1229                  way = 1;
1230                } else {
1231                  way = -1;
1232                }
1233              }
1234            }
1235            if (chosen >= 0) {
1236              improve = true;
1237              // now see how far
1238              if (way > 0) {
1239                // incoming increasing so basic decreasing
1240                // if phase 0 then go to nearest bound
1241                double distance = upper[chosen] - solution[chosen];
1242                double basicDistance;
1243                if (!iphase) {
1244                  assert(iBasic == numberInSet);
1245                  assert(solution[iBasic] > upper[iBasic]);
1246                  basicDistance = solution[iBasic] - upper[iBasic];
1247                } else {
1248                  basicDistance = solution[iBasic] - lower[iBasic];
1249                }
1250                // need extra coding for unbounded
1251                assert(CoinMin(distance, basicDistance) < 1.0e20);
1252                if (distance > basicDistance) {
1253                  // incoming becomes basic
1254                  solution[chosen] += basicDistance;
1255                  if (!iphase)
1256                    solution[iBasic] = upper[iBasic];
1257                  else
1258                    solution[iBasic] = lower[iBasic];
1259                  iBasic = chosen;
1260                } else {
1261                  // flip
1262                  solution[chosen] = upper[chosen];
1263                  solution[iBasic] -= distance;
1264                }
1265              } else {
1266                // incoming decreasing so basic increasing
1267                // if phase 0 then go to nearest bound
1268                double distance = solution[chosen] - lower[chosen];
1269                double basicDistance;
1270                if (!iphase) {
1271                  assert(iBasic == numberInSet);
1272                  assert(solution[iBasic] < lower[iBasic]);
1273                  basicDistance = lower[iBasic] - solution[iBasic];
1274                } else {
1275                  basicDistance = upper[iBasic] - solution[iBasic];
1276                }
1277                // need extra coding for unbounded - for now just exit
1278                if (CoinMin(distance, basicDistance) > 1.0e20) {
1279                  printf("unbounded on set %d\n", iSet);
1280                  iphase = 1;
1281                  iBasic = numberInSet;
1282                  break;
1283                }
1284                if (distance > basicDistance) {
1285                  // incoming becomes basic
1286                  solution[chosen] -= basicDistance;
1287                  if (!iphase)
1288                    solution[iBasic] = lower[iBasic];
1289                  else
1290                    solution[iBasic] = upper[iBasic];
1291                  iBasic = chosen;
1292                } else {
1293                  // flip
1294                  solution[chosen] = lower[chosen];
1295                  solution[iBasic] += distance;
1296                }
1297              }
1298              if (!iphase) {
1299                if (iBasic < numberInSet)
1300                  break; // feasible
1301                else if (solution[iBasic] >= lower[iBasic] && solution[iBasic] <= upper[iBasic])
1302                  break; // feasible (on flip)
1303              }
1304            }
1305          }
1306        }
1307        // do solution i.e. bounds
1308        if (lowerColumn_ || upperColumn_) {
1309          for (int j = 0; j < numberInSet; j++) {
1310            if (j != iBasic) {
1311              objectiveOffset += solution[j] * cost[j];
1312              if (lowerColumn_ && upperColumn_) {
1313                if (fabs(solution[j] - lowerColumn_[j + iStart]) > fabs(solution[j] - upperColumn_[j + iStart]))
1314                  setDynamicStatus(j + iStart, atUpperBound);
1315              } else if (upperColumn_ && solution[j] > 0.0) {
1316                setDynamicStatus(j + iStart, atUpperBound);
1317              } else {
1318                setDynamicStatus(j + iStart, atLowerBound);
1319              }
1320            }
1321          }
1322        }
1323        // convert iBasic back and do bounds
1324        if (iBasic == numberInSet) {
1325          // slack basic
1326          setStatus(iSet, ClpSimplex::basic);
1327          iBasic = iSet + numberColumns;
1328        } else {
1329          iBasic += fullStart_[iSet];
1330          if (back[iBasic] >= 0) {
1331            // exists
1332            iBasic = back[iBasic];
1333          } else {
1334            // create
1335            CoinBigIndex numberElements = startColumn[firstAvailable_];
1336            int numberThis = startColumn_[iBasic + 1] - startColumn_[iBasic];
1337            if (numberElements + numberThis > numberElements_) {
1338              // need to redo
1339              numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
1340              matrix_->reserve(numberColumns, numberElements_);
1341              element = matrix_->getMutableElements();
1342              row = matrix_->getMutableIndices();
1343              // these probably okay but be safe
1344              startColumn = matrix_->getMutableVectorStarts();
1345              length = matrix_->getMutableVectorLengths();
1346            }
1347            length[firstAvailable_] = numberThis;
1348            model->costRegion()[firstAvailable_] = cost_[iBasic];
1349            if (lowerColumn_)
1350              model->lowerRegion()[firstAvailable_] = lowerColumn_[iBasic];
1351            else
1352              model->lowerRegion()[firstAvailable_] = 0.0;
1353            if (upperColumn_)
1354              model->upperRegion()[firstAvailable_] = upperColumn_[iBasic];
1355            else
1356              model->upperRegion()[firstAvailable_] = COIN_DBL_MAX;
1357            columnSolution[firstAvailable_] = solution[iBasic - fullStart_[iSet]];
1358            CoinBigIndex base = startColumn_[iBasic];
1359            for (int j = 0; j < numberThis; j++) {
1360              row[numberElements] = row_[base + j];
1361              element[numberElements++] = element_[base + j];
1362            }
1363            // already set startColumn[firstAvailable_]=numberElements;
1364            id_[firstAvailable_ - firstDynamic_] = iBasic;
1365            setDynamicStatus(iBasic, inSmall);
1366            backward_[firstAvailable_] = iSet;
1367            iBasic = firstAvailable_;
1368            firstAvailable_++;
1369            startColumn[firstAvailable_] = numberElements;
1370          }
1371          model->setStatus(iBasic, ClpSimplex::basic);
1372          // remember bounds flipped
1373          if (upper[numberInSet] == lower[numberInSet])
1374            setStatus(iSet, ClpSimplex::isFixed);
1375          else if (solution[numberInSet] == upper[numberInSet])
1376            setStatus(iSet, ClpSimplex::atLowerBound);
1377          else if (solution[numberInSet] == lower[numberInSet])
1378            setStatus(iSet, ClpSimplex::atUpperBound);
1379          else
1380            abort();
1381        }
1382        for (j = iStart; j < iEnd; j++) {
1383          int iBack = back[j];
1384          if (iBack >= 0) {
1385            if (model->getStatus(iBack) != ClpSimplex::basic) {
1386              int inSet = j - iStart;
1387              columnSolution[iBack] = solution[inSet];
1388              if (upper[inSet] == lower[inSet])
1389                model->setStatus(iBack, ClpSimplex::isFixed);
1390              else if (solution[inSet] == upper[inSet])
1391                model->setStatus(iBack, ClpSimplex::atUpperBound);
1392              else if (solution[inSet] == lower[inSet])
1393                model->setStatus(iBack, ClpSimplex::atLowerBound);
1394            }
1395          }
1396        }
1397      }
1398    }
1399    keyVariable_[iSet] = iBasic;
1400  }
1401  model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
1402  delete[] lower;
1403  delete[] solution;
1404  delete[] upper;
1405  delete[] cost;
1406  // make sure matrix is in good shape
1407  matrix_->orderMatrix();
1408  // create effective rhs
1409  delete[] rhsOffset_;
1410  rhsOffset_ = new double[numberRows];
1411  // and redo chains
1412  memset(mark, 0, numberColumns);
1413  for (int iColumnX = 0; iColumnX < firstAvailable_; iColumnX++)
1414    next_[iColumnX] = COIN_INT_MAX;
1415  for (i = 0; i < numberSets_; i++) {
1416    keys[i] = COIN_INT_MAX;
1417    int iKey = keyVariable_[i];
1418    if (iKey < numberColumns)
1419      model->setStatus(iKey, ClpSimplex::basic);
1420  }
1421  // set up chains
1422  for (i = 0; i < firstAvailable_; i++) {
1423    if (model->getStatus(i) == ClpSimplex::basic)
1424      mark[i] = 1;
1425    int iSet = backward_[i];
1426    if (iSet >= 0) {
1427      int iNext = keys[iSet];
1428      next_[i] = iNext;
1429      keys[iSet] = i;
1430    }
1431  }
1432  for (i = 0; i < numberSets_; i++) {
1433    if (keys[i] != COIN_INT_MAX) {
1434      // something in set
1435      int j;
1436      if (getStatus(i) != ClpSimplex::basic) {
1437        // make sure fixed if it is
1438        if (upper_[i] == lower_[i])
1439          setStatus(i, ClpSimplex::isFixed);
1440        // slack not key - choose one with smallest length
1441        int smallest = numberRows + 1;
1442        int key = -1;
1443        j = keys[i];
1444        while (1) {
1445          if (mark[j] && length[j] < smallest) {
1446            key = j;
1447            smallest = length[j];
1448          }
1449          if (next_[j] != COIN_INT_MAX) {
1450            j = next_[j];
1451          } else {
1452            // correct end
1453            next_[j] = -(keys[i] + 1);
1454            break;
1455          }
1456        }
1457        if (key >= 0) {
1458          keyVariable_[i] = key;
1459        } else {
1460          // nothing basic - make slack key
1461          //((ClpGubMatrix *)this)->setStatus(i,ClpSimplex::basic);
1462          // fudge to avoid const problem
1463          status_[i] = 1;
1464        }
1465      } else {
1466        // slack key
1467        keyVariable_[i] = numberColumns + i;
1468        int j;
1469        double sol = 0.0;
1470        j = keys[i];
1471        while (1) {
1472          sol += columnSolution[j];
1473          if (next_[j] != COIN_INT_MAX) {
1474            j = next_[j];
1475          } else {
1476            // correct end
1477            next_[j] = -(keys[i] + 1);
1478            break;
1479          }
1480        }
1481        if (sol > upper_[i] + tolerance) {
1482          setAbove(i);
1483        } else if (sol < lower_[i] - tolerance) {
1484          setBelow(i);
1485        } else {
1486          setFeasible(i);
1487        }
1488      }
1489      // Create next_
1490      int key = keyVariable_[i];
1491      redoSet(model, key, keys[i], i);
1492    } else {
1493      // nothing in set!
1494      next_[i + numberColumns] = -(i + numberColumns + 1);
1495      keyVariable_[i] = numberColumns + i;
1496      double sol = 0.0;
1497      if (sol > upper_[i] + tolerance) {
1498        setAbove(i);
1499      } else if (sol < lower_[i] - tolerance) {
1500        setBelow(i);
1501      } else {
1502        setFeasible(i);
1503      }
1504    }
1505  }
1506  delete[] keys;
1507  delete[] mark;
1508  delete[] back;
1509  rhsOffset(model, true);
1510}
1511/* Returns effective RHS if it is being used.  This is used for long problems
1512   or big gub or anywhere where going through full columns is
1513   expensive.  This may re-compute */
1514double *
1515ClpGubDynamicMatrix::rhsOffset(ClpSimplex *model, bool forceRefresh,
1516  bool
1517#ifdef CLP_DEBUG
1518    check
1519#endif
1520)
1521{
1522  //forceRefresh=true;
1523  //check=false;
1524#ifdef CLP_DEBUG
1525  double *saveE = NULL;
1526  if (rhsOffset_ && check) {
1527    int numberRows = model->numberRows();
1528    saveE = new double[numberRows];
1529  }
1530#endif
1531  if (rhsOffset_) {
1532#ifdef CLP_DEBUG
1533    if (check) {
1534      // no need - but check anyway
1535      int numberRows = model->numberRows();
1536      double *rhs = new double[numberRows];
1537      int numberColumns = model->numberColumns();
1538      int iRow;
1539      CoinZeroN(rhs, numberRows);
1540      // do ones at bounds before gub
1541      const double *smallSolution = model->solutionRegion();
1542      const double *element = matrix_->getElements();
1543      const int *row = matrix_->getIndices();
1544      const CoinBigIndex *startColumn = matrix_->getVectorStarts();
1545      const int *length = matrix_->getVectorLengths();
1546      int iColumn;
1547      for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
1548        if (model->getStatus(iColumn) != ClpSimplex::basic) {
1549          double value = smallSolution[iColumn];
1550          for (CoinBigIndex j = startColumn[iColumn];
1551               j < startColumn[iColumn] + length[iColumn]; j++) {
1552            int jRow = row[j];
1553            rhs[jRow] -= value * element[j];
1554          }
1555        }
1556      }
1557      if (lowerColumn_ || upperColumn_) {
1558        double *solution = new double[numberGubColumns_];
1559        for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
1560          double value = 0.0;
1561          if (getDynamicStatus(iColumn) == atUpperBound)
1562            value = upperColumn_[iColumn];
1563          else if (lowerColumn_)
1564            value = lowerColumn_[iColumn];
1565          solution[iColumn] = value;
1566        }
1567        // ones at bounds in small and gub
1568        for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
1569          int jFull = id_[iColumn - firstDynamic_];
1570          solution[jFull] = smallSolution[iColumn];
1571        }
1572        // zero all basic in small model
1573        int *pivotVariable = model->pivotVariable();
1574        for (iRow = 0; iRow < numberRows; iRow++) {
1575          int iColumn = pivotVariable[iRow];
1576          if (iColumn >= firstDynamic_ && iColumn < lastDynamic_) {
1577            int iSequence = id_[iColumn - firstDynamic_];
1578            solution[iSequence] = 0.0;
1579          }
1580        }
1581        // and now compute value to use for key
1582        ClpSimplex::Status iStatus;
1583        for (int iSet = 0; iSet < numberSets_; iSet++) {
1584          iColumn = keyVariable_[iSet];
1585          if (iColumn < numberColumns) {
1586            int iSequence = id_[iColumn - firstDynamic_];
1587            solution[iSequence] = 0.0;
1588            double b = 0.0;
1589            // key is structural - where is slack
1590            iStatus = getStatus(iSet);
1591            assert(iStatus != ClpSimplex::basic);
1592            if (iStatus == ClpSimplex::atLowerBound)
1593              b = lowerSet_[iSet];
1594            else
1595              b = upperSet_[iSet];
1596            // subtract out others at bounds
1597            for (int j = fullStart_[iSet]; j < fullStart_[iSet + 1]; j++)
1598              b -= solution[j];
1599            solution[iSequence] = b;
1600          }
1601        }
1602        for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
1603          double value = solution[iColumn];
1604          if (value) {
1605            for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn + 1]; j++) {
1606              int iRow = row_[j];
1607              rhs[iRow] -= element_[j] * value;
1608            }
1609          }
1610        }
1611        // now do lower and upper bounds on sets
1612        for (int iSet = 0; iSet < numberSets_; iSet++) {
1613          iColumn = keyVariable_[iSet];
1614          double shift = 0.0;
1615          for (int j = fullStart_[iSet]; j < fullStart_[iSet + 1]; j++) {
1616            if (getDynamicStatus(j) != inSmall && j != iColumn) {
1617              if (getDynamicStatus(j) == atLowerBound) {
1618                if (lowerColumn_)
1619                  shift += lowerColumn_[j];
1620              } else {
1621                shift += upperColumn_[j];
1622              }
1623            }
1624          }
1625          if (lowerSet_[iSet] > -1.0e20)
1626            assert(fabs(lower_[iSet] - (lowerSet_[iSet] - shift)) < 1.0e-3);
1627          if (upperSet_[iSet] < 1.0e20)
1628            assert(fabs(upper_[iSet] - (upperSet_[iSet] - shift)) < 1.0e-3);
1629        }
1630        delete[] solution;
1631      } else {
1632        // no bounds
1633        ClpSimplex::Status iStatus;
1634        for (int iSet = 0; iSet < numberSets_; iSet++) {
1635          int iColumn = keyVariable_[iSet];
1636          if (iColumn < numberColumns) {
1637            int iSequence = id_[iColumn - firstDynamic_];
1638            double b = 0.0;
1639            // key is structural - where is slack
1640            iStatus = getStatus(iSet);
1641            assert(iStatus != ClpSimplex::basic);
1642            if (iStatus == ClpSimplex::atLowerBound)
1643              b = lower_[iSet];
1644            else
1645              b = upper_[iSet];
1646            if (b) {
1647              for (CoinBigIndex j = startColumn_[iSequence]; j < startColumn_[iSequence + 1]; j++) {
1648                int iRow = row_[j];
1649                rhs[iRow] -= element_[j] * b;
1650              }
1651            }
1652          }
1653        }
1654      }
1655      for (iRow = 0; iRow < numberRows; iRow++) {
1656        if (fabs(rhs[iRow] - rhsOffset_[iRow]) > 1.0e-3)
1657          printf("** bad effective %d - true %g old %g\n", iRow, rhs[iRow], rhsOffset_[iRow]);
1658      }
1659      CoinMemcpyN(rhs, numberRows, saveE);
1660      delete[] rhs;
1661    }
1662#endif
1663    if (forceRefresh || (refreshFrequency_ && model->numberIterations() >= lastRefresh_ + refreshFrequency_)) {
1664      int numberRows = model->numberRows();
1665      int numberColumns = model->numberColumns();
1666      int iRow;
1667      CoinZeroN(rhsOffset_, numberRows);
1668      // do ones at bounds before gub
1669      const double *smallSolution = model->solutionRegion();
1670      const double *element = matrix_->getElements();
1671      const int *row = matrix_->getIndices();
1672      const CoinBigIndex *startColumn = matrix_->getVectorStarts();
1673      const int *length = matrix_->getVectorLengths();
1674      int iColumn;
1675      for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
1676        if (model->getStatus(iColumn) != ClpSimplex::basic) {
1677          double value = smallSolution[iColumn];
1678          for (CoinBigIndex j = startColumn[iColumn];
1679               j < startColumn[iColumn] + length[iColumn]; j++) {
1680            int jRow = row[j];
1681            rhsOffset_[jRow] -= value * element[j];
1682          }
1683        }
1684      }
1685      if (lowerColumn_ || upperColumn_) {
1686        double *solution = new double[numberGubColumns_];
1687        for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
1688          double value = 0.0;
1689          if (getDynamicStatus(iColumn) == atUpperBound)
1690            value = upperColumn_[iColumn];
1691          else if (lowerColumn_)
1692            value = lowerColumn_[iColumn];
1693          solution[iColumn] = value;
1694        }
1695        // ones in gub and in small problem
1696        for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
1697          int jFull = id_[iColumn - firstDynamic_];
1698          solution[jFull] = smallSolution[iColumn];
1699        }
1700        // zero all basic in small model
1701        int *pivotVariable = model->pivotVariable();
1702        for (iRow = 0; iRow < numberRows; iRow++) {
1703          int iColumn = pivotVariable[iRow];
1704          if (iColumn >= firstDynamic_ && iColumn < lastDynamic_) {
1705            int iSequence = id_[iColumn - firstDynamic_];
1706            solution[iSequence] = 0.0;
1707          }
1708        }
1709        // and now compute value to use for key
1710        ClpSimplex::Status iStatus;
1711        int iSet;
1712        for (iSet = 0; iSet < numberSets_; iSet++) {
1713          iColumn = keyVariable_[iSet];
1714          if (iColumn < numberColumns) {
1715            int iSequence = id_[iColumn - firstDynamic_];
1716            solution[iSequence] = 0.0;
1717            double b = 0.0;
1718            // key is structural - where is slack
1719            iStatus = getStatus(iSet);
1720            assert(iStatus != ClpSimplex::basic);
1721            if (iStatus == ClpSimplex::atLowerBound)
1722              b = lowerSet_[iSet];
1723            else
1724              b = upperSet_[iSet];
1725            // subtract out others at bounds
1726            for (int j = fullStart_[iSet]; j < fullStart_[iSet + 1]; j++)
1727              b -= solution[j];
1728            solution[iSequence] = b;
1729          }
1730        }
1731        for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
1732          double value = solution[iColumn];
1733          if (value) {
1734            for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn + 1]; j++) {
1735              int iRow = row_[j];
1736              rhsOffset_[iRow] -= element_[j] * value;
1737            }
1738          }
1739        }
1740        // now do lower and upper bounds on sets
1741        // and offset
1742        double objectiveOffset = 0.0;
1743        for (iSet = 0; iSet < numberSets_; iSet++) {
1744          iColumn = keyVariable_[iSet];
1745          double shift = 0.0;
1746          for (CoinBigIndex j = fullStart_[iSet]; j < fullStart_[iSet + 1]; j++) {
1747            if (getDynamicStatus(j) != inSmall) {
1748              double value = 0.0;
1749              if (getDynamicStatus(j) == atLowerBound) {
1750                if (lowerColumn_)
1751                  value = lowerColumn_[j];
1752              } else {
1753                value = upperColumn_[j];
1754              }
1755              if (j != iColumn)
1756                shift += value;
1757              objectiveOffset += value * cost_[j];
1758            }
1759          }
1760          if (lowerSet_[iSet] > -1.0e20)
1761            lower_[iSet] = lowerSet_[iSet] - shift;
1762          if (upperSet_[iSet] < 1.0e20)
1763            upper_[iSet] = upperSet_[iSet] - shift;
1764        }
1765        delete[] solution;
1766        model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
1767      } else {
1768        // no bounds
1769        ClpSimplex::Status iStatus;
1770        for (int iSet = 0; iSet < numberSets_; iSet++) {
1771          int iColumn = keyVariable_[iSet];
1772          if (iColumn < numberColumns) {
1773            int iSequence = id_[iColumn - firstDynamic_];
1774            double b = 0.0;
1775            // key is structural - where is slack
1776            iStatus = getStatus(iSet);
1777            assert(iStatus != ClpSimplex::basic);
1778            if (iStatus == ClpSimplex::atLowerBound)
1779              b = lower_[iSet];
1780            else
1781              b = upper_[iSet];
1782            if (b) {
1783              for (CoinBigIndex j = startColumn_[iSequence]; j < startColumn_[iSequence + 1]; j++) {
1784                int iRow = row_[j];
1785                rhsOffset_[iRow] -= element_[j] * b;
1786              }
1787            }
1788          }
1789        }
1790      }
1791#ifdef CLP_DEBUG
1792      if (saveE) {
1793        for (iRow = 0; iRow < numberRows; iRow++) {
1794          if (fabs(saveE[iRow] - rhsOffset_[iRow]) > 1.0e-3)
1795            printf("** %d - old eff %g new %g\n", iRow, saveE[iRow], rhsOffset_[iRow]);
1796        }
1797        delete[] saveE;
1798      }
1799#endif
1800      lastRefresh_ = model->numberIterations();
1801    }
1802  }
1803  return rhsOffset_;
1804}
1805/*
1806  update information for a pivot (and effective rhs)
1807*/
1808int ClpGubDynamicMatrix::updatePivot(ClpSimplex *model, double oldInValue, double oldOutValue)
1809{
1810
1811  // now update working model
1812  int sequenceIn = model->sequenceIn();
1813  int sequenceOut = model->sequenceOut();
1814  bool doPrinting = (model->messageHandler()->logLevel() == 63);
1815  bool print = false;
1816  int iSet;
1817  int trueIn = -1;
1818  int trueOut = -1;
1819  int numberRows = model->numberRows();
1820  int numberColumns = model->numberColumns();
1821  if (sequenceIn == firstAvailable_) {
1822    if (doPrinting)
1823      printf("New variable ");
1824    if (sequenceIn != sequenceOut) {
1825      insertNonBasic(firstAvailable_, backward_[firstAvailable_]);
1826      setDynamicStatus(id_[sequenceIn - firstDynamic_], inSmall);
1827      firstAvailable_++;
1828    } else {
1829      int bigSequence = id_[sequenceIn - firstDynamic_];
1830      if (model->getStatus(sequenceIn) == ClpSimplex::atUpperBound)
1831        setDynamicStatus(bigSequence, atUpperBound);
1832      else
1833        setDynamicStatus(bigSequence, atLowerBound);
1834    }
1835    synchronize(model, 8);
1836  }
1837  if (sequenceIn < lastDynamic_) {
1838    iSet = backward_[sequenceIn];
1839    if (iSet >= 0) {
1840      int bigSequence = id_[sequenceIn - firstDynamic_];
1841      trueIn = bigSequence + numberRows + numberColumns + numberSets_;
1842      if (doPrinting)
1843        printf(" incoming set %d big seq %d", iSet, bigSequence);
1844      print = true;
1845    }
1846  } else if (sequenceIn >= numberRows + numberColumns) {
1847    trueIn = numberRows + numberColumns + gubSlackIn_;
1848  }
1849  if (sequenceOut < lastDynamic_) {
1850    iSet = backward_[sequenceOut];
1851    if (iSet >= 0) {
1852      int bigSequence = id_[sequenceOut - firstDynamic_];
1853      trueOut = bigSequence + firstDynamic_;
1854      if (getDynamicStatus(bigSequence) != inSmall) {
1855        if (model->getStatus(sequenceOut) == ClpSimplex::atUpperBound)
1856          setDynamicStatus(bigSequence, atUpperBound);
1857        else
1858          setDynamicStatus(bigSequence, atLowerBound);
1859      }
1860      if (doPrinting)
1861        printf(" ,outgoing set %d big seq %d,", iSet, bigSequence);
1862      print = true;
1863      model->setSequenceIn(sequenceOut);
1864      synchronize(model, 8);
1865      model->setSequenceIn(sequenceIn);
1866    }
1867  }
1868  if (print && doPrinting)
1869    printf("\n");
1870  ClpGubMatrix::updatePivot(model, oldInValue, oldOutValue);
1871  // Redo true in and out
1872  if (trueIn >= 0)
1873    trueSequenceIn_ = trueIn;
1874  if (trueOut >= 0)
1875    trueSequenceOut_ = trueOut;
1876  if (doPrinting && 0) {
1877    for (int i = 0; i < numberSets_; i++) {
1878      printf("set %d key %d lower %g upper %g\n", i, keyVariable_[i], lower_[i], upper_[i]);
1879      for (int j = fullStart_[i]; j < fullStart_[i + 1]; j++)
1880        if (getDynamicStatus(j) == atUpperBound) {
1881          bool print = true;
1882          for (int k = firstDynamic_; k < firstAvailable_; k++) {
1883            if (id_[k - firstDynamic_] == j)
1884              print = false;
1885            if (id_[k - firstDynamic_] == j)
1886              assert(getDynamicStatus(j) == inSmall);
1887          }
1888          if (print)
1889            printf("variable %d at ub\n", j);
1890        }
1891    }
1892  }
1893#ifdef CLP_DEBUG
1894  char *inSmall = new char[numberGubColumns_];
1895  memset(inSmall, 0, numberGubColumns_);
1896  for (int i = 0; i < numberGubColumns_; i++)
1897    if (getDynamicStatus(i) == ClpGubDynamicMatrix::inSmall)
1898      inSmall[i] = 1;
1899  for (int i = firstDynamic_; i < firstAvailable_; i++) {
1900    int k = id_[i - firstDynamic_];
1901    inSmall[k] = 0;
1902  }
1903  for (int i = 0; i < numberGubColumns_; i++)
1904    assert(!inSmall[i]);
1905  delete[] inSmall;
1906#endif
1907  return 0;
1908}
1909void ClpGubDynamicMatrix::times(double scalar,
1910  const double *x, double *y) const
1911{
1912  if (model_->specialOptions() != 16) {
1913    ClpPackedMatrix::times(scalar, x, y);
1914  } else {
1915    int iRow;
1916    int numberColumns = model_->numberColumns();
1917    int numberRows = model_->numberRows();
1918    const double *element = matrix_->getElements();
1919    const int *row = matrix_->getIndices();
1920    const CoinBigIndex *startColumn = matrix_->getVectorStarts();
1921    const int *length = matrix_->getVectorLengths();
1922    int *pivotVariable = model_->pivotVariable();
1923    int numberToDo = 0;
1924    for (iRow = 0; iRow < numberRows; iRow++) {
1925      y[iRow] -= scalar * rhsOffset_[iRow];
1926      int iColumn = pivotVariable[iRow];
1927      if (iColumn < numberColumns) {
1928        int iSet = backward_[iColumn];
1929        if (iSet >= 0 && toIndex_[iSet] < 0) {
1930          toIndex_[iSet] = 0;
1931          fromIndex_[numberToDo++] = iSet;
1932        }
1933        CoinBigIndex j;
1934        double value = scalar * x[iColumn];
1935        if (value) {
1936          for (j = startColumn[iColumn];
1937               j < startColumn[iColumn] + length[iColumn]; j++) {
1938            int jRow = row[j];
1939            y[jRow] += value * element[j];
1940          }
1941        }
1942      }
1943    }
1944    // and gubs which are interacting
1945    for (int jSet = 0; jSet < numberToDo; jSet++) {
1946      int iSet = fromIndex_[jSet];
1947      toIndex_[iSet] = -1;
1948      int iKey = keyVariable_[iSet];
1949      if (iKey < numberColumns) {
1950        double valueKey;
1951        if (getStatus(iSet) == ClpSimplex::atLowerBound)
1952          valueKey = lower_[iSet];
1953        else
1954          valueKey = upper_[iSet];
1955        double value = scalar * (x[iKey] - valueKey);
1956        if (value) {
1957          for (CoinBigIndex j = startColumn[iKey];
1958               j < startColumn[iKey] + length[iKey]; j++) {
1959            int jRow = row[j];
1960            y[jRow] += value * element[j];
1961          }
1962        }
1963      }
1964    }
1965  }
1966}
1967/* Just for debug - may be extended to other matrix types later.
1968   Returns number and sum of primal infeasibilities.
1969*/
1970int ClpGubDynamicMatrix::checkFeasible(ClpSimplex * /*model*/, double &sum) const
1971{
1972  int numberRows = model_->numberRows();
1973  double *rhs = new double[numberRows];
1974  int numberColumns = model_->numberColumns();
1975  int iRow;
1976  CoinZeroN(rhs, numberRows);
1977  // do ones at bounds before gub
1978  const double *smallSolution = model_->solutionRegion();
1979  const double *element = matrix_->getElements();
1980  const int *row = matrix_->getIndices();
1981  const CoinBigIndex *startColumn = matrix_->getVectorStarts();
1982  const int *length = matrix_->getVectorLengths();
1983  int iColumn;
1984  int numberInfeasible = 0;
1985  const double *rowLower = model_->rowLower();
1986  const double *rowUpper = model_->rowUpper();
1987  sum = 0.0;
1988  for (iRow = 0; iRow < numberRows; iRow++) {
1989    double value = smallSolution[numberColumns + iRow];
1990    if (value < rowLower[iRow] - 1.0e-5 || value > rowUpper[iRow] + 1.0e-5) {
1991      //printf("row %d %g %g %g\n",
1992      //     iRow,rowLower[iRow],value,rowUpper[iRow]);
1993      numberInfeasible++;
1994      sum += CoinMax(rowLower[iRow] - value, value - rowUpper[iRow]);
1995    }
1996    rhs[iRow] = value;
1997  }
1998  const double *columnLower = model_->columnLower();
1999  const double *columnUpper = model_->columnUpper();
2000  for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
2001    double value = smallSolution[iColumn];
2002    if (value < columnLower[iColumn] - 1.0e-5 || value > columnUpper[iColumn] + 1.0e-5) {
2003      //printf("column %d %g %g %g\n",
2004      //     iColumn,columnLower[iColumn],value,columnUpper[iColumn]);
2005      numberInfeasible++;
2006      sum += CoinMax(columnLower[iColumn] - value, value - columnUpper[iColumn]);
2007    }
2008    for (CoinBigIndex j = startColumn[iColumn];
2009         j < startColumn[iColumn] + length[iColumn]; j++) {
2010      int jRow = row[j];
2011      rhs[jRow] -= value * element[j];
2012    }
2013  }
2014  double *solution = new double[numberGubColumns_];
2015  for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
2016    double value = 0.0;
2017    if (getDynamicStatus(iColumn) == atUpperBound)
2018      value = upperColumn_[iColumn];
2019    else if (lowerColumn_)
2020      value = lowerColumn_[iColumn];
2021    solution[iColumn] = value;
2022  }
2023  // ones in small and gub
2024  for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
2025    int jFull = id_[iColumn - firstDynamic_];
2026    solution[jFull] = smallSolution[iColumn];
2027  }
2028  // fill in all basic in small model
2029  int *pivotVariable = model_->pivotVariable();
2030  for (iRow = 0; iRow < numberRows; iRow++) {
2031    int iColumn = pivotVariable[iRow];
2032    if (iColumn >= firstDynamic_ && iColumn < lastDynamic_) {
2033      int iSequence = id_[iColumn - firstDynamic_];
2034      solution[iSequence] = smallSolution[iColumn];
2035    }
2036  }
2037  // and now compute value to use for key
2038  ClpSimplex::Status iStatus;
2039  for (int iSet = 0; iSet < numberSets_; iSet++) {
2040    iColumn = keyVariable_[iSet];
2041    if (iColumn < numberColumns) {
2042      int iSequence = id_[iColumn - firstDynamic_];
2043      solution[iSequence] = 0.0;
2044      double b = 0.0;
2045      // key is structural - where is slack
2046      iStatus = getStatus(iSet);
2047      assert(iStatus != ClpSimplex::basic);
2048      if (iStatus == ClpSimplex::atLowerBound)
2049        b = lower_[iSet];
2050      else
2051        b = upper_[iSet];
2052      // subtract out others at bounds
2053      for (int j = fullStart_[iSet]; j < fullStart_[iSet + 1]; j++)
2054        b -= solution[j];
2055      solution[iSequence] = b;
2056    }
2057  }
2058  for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
2059    double value = solution[iColumn];
2060    if ((lowerColumn_ && value < lowerColumn_[iColumn] - 1.0e-5) || (!lowerColumn_ && value < -1.0e-5) || (upperColumn_ && value > upperColumn_[iColumn] + 1.0e-5)) {
2061      //printf("column %d %g %g %g\n",
2062      //     iColumn,lowerColumn_[iColumn],value,upperColumn_[iColumn]);
2063      numberInfeasible++;
2064    }
2065    if (value) {
2066      for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn + 1]; j++) {
2067        int iRow = row_[j];
2068        rhs[iRow] -= element_[j] * value;
2069      }
2070    }
2071  }
2072  for (iRow = 0; iRow < numberRows; iRow++) {
2073    if (fabs(rhs[iRow]) > 1.0e-5)
2074      printf("rhs mismatch %d %g\n", iRow, rhs[iRow]);
2075  }
2076  delete[] solution;
2077  delete[] rhs;
2078  return numberInfeasible;
2079}
2080// Cleans data after setWarmStart
2081void ClpGubDynamicMatrix::cleanData(ClpSimplex *model)
2082{
2083  // and redo chains
2084  int numberColumns = model->numberColumns();
2085  int iColumn;
2086  // do backward
2087  int *mark = new int[numberGubColumns_];
2088  for (iColumn = 0; iColumn < numberGubColumns_; iColumn++)
2089    mark[iColumn] = -1;
2090  int i;
2091  for (i = 0; i < firstDynamic_; i++) {
2092    assert(backward_[i] == -1);
2093    next_[i] = -1;
2094  }
2095  for (i = firstDynamic_; i < firstAvailable_; i++) {
2096    iColumn = id_[i - firstDynamic_];
2097    mark[iColumn] = i;
2098  }
2099  for (i = 0; i < numberSets_; i++) {
2100    int iKey = keyVariable_[i];
2101    int lastNext = -1;
2102    int firstNext = -1;
2103    for (CoinBigIndex k = fullStart_[i]; k < fullStart_[i + 1]; k++) {
2104      iColumn = mark[k];
2105      if (iColumn >= 0) {
2106        if (iColumn != iKey) {
2107          if (lastNext >= 0)
2108            next_[lastNext] = iColumn;
2109          else
2110            firstNext = iColumn;
2111          lastNext = iColumn;
2112        }
2113        backward_[iColumn] = i;
2114      }
2115    }
2116    setFeasible(i);
2117    if (firstNext >= 0) {
2118      // others
2119      next_[iKey] = firstNext;
2120      next_[lastNext] = -(iKey + 1);
2121    } else if (iKey < numberColumns) {
2122      next_[iKey] = -(iKey + 1);
2123    }
2124  }
2125  delete[] mark;
2126  // fill matrix
2127  double *element = matrix_->getMutableElements();
2128  int *row = matrix_->getMutableIndices();
2129  CoinBigIndex *startColumn = matrix_->getMutableVectorStarts();
2130  int *length = matrix_->getMutableVectorLengths();
2131  CoinBigIndex numberElements = startColumn[firstDynamic_];
2132  for (i = firstDynamic_; i < firstAvailable_; i++) {
2133    int iColumn = id_[i - firstDynamic_];
2134    int numberThis = startColumn_[iColumn + 1] - startColumn_[iColumn];
2135    length[i] = numberThis;
2136    for (CoinBigIndex jBigIndex = startColumn_[iColumn];
2137         jBigIndex < startColumn_[iColumn + 1]; jBigIndex++) {
2138      row[numberElements] = row_[jBigIndex];
2139      element[numberElements++] = element_[jBigIndex];
2140    }
2141    startColumn[i + 1] = numberElements;
2142  }
2143}
2144
2145/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
2146*/
Note: See TracBrowser for help on using the repository browser.