source: trunk/Clp/src/ClpDynamicMatrix.cpp

Last change on this file was 2385, checked in by unxusr, 10 months ago

formatting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 91.2 KB
Line 
1/* $Id: ClpDynamicMatrix.cpp 2385 2019-01-06 19:43:06Z stefan $ */
2// Copyright (C) 2004, 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 "ClpDynamicMatrix.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//-------------------------------------------------------------------
28ClpDynamicMatrix::ClpDynamicMatrix()
29  : ClpPackedMatrix()
30  , sumDualInfeasibilities_(0.0)
31  , sumPrimalInfeasibilities_(0.0)
32  , sumOfRelaxedDualInfeasibilities_(0.0)
33  , sumOfRelaxedPrimalInfeasibilities_(0.0)
34  , savedBestGubDual_(0.0)
35  , savedBestSet_(0)
36  , backToPivotRow_(NULL)
37  , keyVariable_(NULL)
38  , toIndex_(NULL)
39  , fromIndex_(NULL)
40  , numberSets_(0)
41  , numberActiveSets_(0)
42  , objectiveOffset_(0.0)
43  , lowerSet_(NULL)
44  , upperSet_(NULL)
45  , status_(NULL)
46  , model_(NULL)
47  , firstAvailable_(0)
48  , firstAvailableBefore_(0)
49  , firstDynamic_(0)
50  , lastDynamic_(0)
51  , numberStaticRows_(0)
52  , numberElements_(0)
53  , numberDualInfeasibilities_(0)
54  , numberPrimalInfeasibilities_(0)
55  , noCheck_(-1)
56  , infeasibilityWeight_(0.0)
57  , numberGubColumns_(0)
58  , maximumGubColumns_(0)
59  , maximumElements_(0)
60  , startSet_(NULL)
61  , next_(NULL)
62  , startColumn_(NULL)
63  , row_(NULL)
64  , element_(NULL)
65  , cost_(NULL)
66  , id_(NULL)
67  , dynamicStatus_(NULL)
68  , columnLower_(NULL)
69  , columnUpper_(NULL)
70{
71  setType(15);
72}
73
74//-------------------------------------------------------------------
75// Copy constructor
76//-------------------------------------------------------------------
77ClpDynamicMatrix::ClpDynamicMatrix(const ClpDynamicMatrix &rhs)
78  : ClpPackedMatrix(rhs)
79{
80  objectiveOffset_ = rhs.objectiveOffset_;
81  numberSets_ = rhs.numberSets_;
82  numberActiveSets_ = rhs.numberActiveSets_;
83  firstAvailable_ = rhs.firstAvailable_;
84  firstAvailableBefore_ = rhs.firstAvailableBefore_;
85  firstDynamic_ = rhs.firstDynamic_;
86  lastDynamic_ = rhs.lastDynamic_;
87  numberStaticRows_ = rhs.numberStaticRows_;
88  numberElements_ = rhs.numberElements_;
89  backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, lastDynamic_);
90  keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_);
91  toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_);
92  fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + 1 - numberStaticRows_);
93  lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
94  upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
95  status_ = ClpCopyOfArray(rhs.status_, static_cast< int >(2 * numberSets_ + 4 * sizeof(int)));
96  model_ = rhs.model_;
97  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
98  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
99  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
100  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
101  numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
102  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
103  savedBestGubDual_ = rhs.savedBestGubDual_;
104  savedBestSet_ = rhs.savedBestSet_;
105  noCheck_ = rhs.noCheck_;
106  infeasibilityWeight_ = rhs.infeasibilityWeight_;
107  // Now secondary data
108  numberGubColumns_ = rhs.numberGubColumns_;
109  maximumGubColumns_ = rhs.maximumGubColumns_;
110  maximumElements_ = rhs.maximumElements_;
111  startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_ + 1);
112  next_ = ClpCopyOfArray(rhs.next_, maximumGubColumns_);
113  startColumn_ = ClpCopyOfArray(rhs.startColumn_, maximumGubColumns_ + 1);
114  row_ = ClpCopyOfArray(rhs.row_, maximumElements_);
115  element_ = ClpCopyOfArray(rhs.element_, maximumElements_);
116  cost_ = ClpCopyOfArray(rhs.cost_, maximumGubColumns_);
117  id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_);
118  columnLower_ = ClpCopyOfArray(rhs.columnLower_, maximumGubColumns_);
119  columnUpper_ = ClpCopyOfArray(rhs.columnUpper_, maximumGubColumns_);
120  dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2 * maximumGubColumns_);
121}
122
123/* This is the real constructor*/
124ClpDynamicMatrix::ClpDynamicMatrix(ClpSimplex *model, int numberSets,
125  int numberGubColumns, const int *starts,
126  const double *lower, const double *upper,
127  const CoinBigIndex *startColumn, const int *row,
128  const double *element, const double *cost,
129  const double *columnLower, const double *columnUpper,
130  const unsigned char *status,
131  const unsigned char *dynamicStatus)
132  : ClpPackedMatrix()
133{
134  setType(15);
135  objectiveOffset_ = model->objectiveOffset();
136  model_ = model;
137  numberSets_ = numberSets;
138  numberGubColumns_ = numberGubColumns;
139  maximumGubColumns_ = numberGubColumns_;
140  if (numberGubColumns_)
141    maximumElements_ = startColumn[numberGubColumns_];
142  else
143    maximumElements_ = 0;
144  startSet_ = new int[numberSets_ + 1];
145  next_ = new int[maximumGubColumns_];
146  // fill in startSet and next
147  int iSet;
148  if (numberGubColumns_) {
149    for (iSet = 0; iSet < numberSets_; iSet++) {
150      int first = starts[iSet];
151      int last = starts[iSet + 1] - 1;
152      startSet_[iSet] = first;
153      for (int i = first; i < last; i++)
154        next_[i] = i + 1;
155      next_[last] = -iSet - 1;
156    }
157    startSet_[numberSets_] = starts[numberSets_];
158  }
159  int numberColumns = model->numberColumns();
160  int numberRows = model->numberRows();
161  numberStaticRows_ = numberRows;
162  savedBestGubDual_ = 0.0;
163  savedBestSet_ = 0;
164  // Number of columns needed
165  int frequency = model->factorizationFrequency();
166  int numberGubInSmall = numberRows + frequency + CoinMin(frequency, numberSets_) + 4;
167  // But we may have two per row + one for incoming (make it two)
168  numberGubInSmall = CoinMax(2 * numberRows + 2, numberGubInSmall);
169  // for small problems this could be too big
170  //numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_);
171  int numberNeeded = numberGubInSmall + numberColumns;
172  firstAvailable_ = numberColumns;
173  firstAvailableBefore_ = firstAvailable_;
174  firstDynamic_ = numberColumns;
175  lastDynamic_ = numberNeeded;
176  startColumn_ = ClpCopyOfArray(startColumn, numberGubColumns_ + 1);
177  if (!numberGubColumns_) {
178    if (!startColumn_)
179      startColumn_ = new CoinBigIndex[1];
180    startColumn_[0] = 0;
181  }
182  CoinBigIndex numberElements = startColumn_[numberGubColumns_];
183  row_ = ClpCopyOfArray(row, numberElements);
184  element_ = new double[numberElements];
185  CoinBigIndex i;
186  for (i = 0; i < numberElements; i++)
187    element_[i] = element[i];
188  cost_ = new double[numberGubColumns_];
189  for (i = 0; i < numberGubColumns_; i++) {
190    cost_[i] = cost[i];
191    // I don't think I need sorted but ...
192    CoinSort_2(row_ + startColumn_[i], row_ + startColumn_[i + 1], element_ + startColumn_[i]);
193  }
194  if (columnLower) {
195    columnLower_ = new double[numberGubColumns_];
196    for (i = 0; i < numberGubColumns_; i++)
197      columnLower_[i] = columnLower[i];
198  } else {
199    columnLower_ = NULL;
200  }
201  if (columnUpper) {
202    columnUpper_ = new double[numberGubColumns_];
203    for (i = 0; i < numberGubColumns_; i++)
204      columnUpper_[i] = columnUpper[i];
205  } else {
206    columnUpper_ = NULL;
207  }
208  lowerSet_ = new double[numberSets_];
209  for (i = 0; i < numberSets_; i++) {
210    if (lower[i] > -1.0e20)
211      lowerSet_[i] = lower[i];
212    else
213      lowerSet_[i] = -1.0e30;
214  }
215  upperSet_ = new double[numberSets_];
216  for (i = 0; i < numberSets_; i++) {
217    if (upper[i] < 1.0e20)
218      upperSet_[i] = upper[i];
219    else
220      upperSet_[i] = 1.0e30;
221  }
222  id_ = new int[numberGubInSmall];
223  for (i = 0; i < numberGubInSmall; i++)
224    id_[i] = -1;
225  ClpPackedMatrix *originalMatrixA = dynamic_cast< ClpPackedMatrix * >(model->clpMatrix());
226  assert(originalMatrixA);
227  CoinPackedMatrix *originalMatrix = originalMatrixA->getPackedMatrix();
228  originalMatrixA->setMatrixNull(); // so can be deleted safely
229  // guess how much space needed
230  double guess = numberElements;
231  guess /= static_cast< double >(numberColumns);
232  guess *= 2 * numberGubInSmall;
233  numberElements_ = static_cast< int >(guess);
234  numberElements_ = CoinMin(numberElements_, numberElements) + originalMatrix->getNumElements();
235  matrix_ = originalMatrix;
236  //delete originalMatrixA;
237  flags_ &= ~1;
238  // resize model (matrix stays same)
239  // modify frequency
240  if (frequency >= 50)
241    frequency = 50 + (frequency - 50) / 2;
242  int newRowSize = numberRows + CoinMin(numberSets_, frequency + numberRows) + 1;
243  model->resize(newRowSize, numberNeeded);
244  for (i = numberRows; i < newRowSize; i++)
245    model->setRowStatus(i, ClpSimplex::basic);
246  if (columnUpper_) {
247    // set all upper bounds so we have enough space
248    double *columnUpper = model->columnUpper();
249    for (i = firstDynamic_; i < lastDynamic_; i++)
250      columnUpper[i] = 1.0e10;
251  }
252  // resize matrix
253  // extra 1 is so can keep number of elements handy
254  originalMatrix->reserve(numberNeeded, numberElements_, true);
255  originalMatrix->reserve(numberNeeded + 1, numberElements_, false);
256  originalMatrix->getMutableVectorStarts()[numberColumns] = originalMatrix->getNumElements();
257  originalMatrix->setDimensions(newRowSize, -1);
258  numberActiveColumns_ = firstDynamic_;
259  // redo number of columns
260  numberColumns = matrix_->getNumCols();
261  backToPivotRow_ = new int[numberNeeded];
262  keyVariable_ = new int[numberSets_];
263  if (status) {
264    status_ = ClpCopyOfArray(status, static_cast< int >(2 * numberSets_ + 4 * sizeof(int)));
265    assert(dynamicStatus);
266    dynamicStatus_ = ClpCopyOfArray(dynamicStatus, 2 * numberGubColumns_);
267  } else {
268    status_ = new unsigned char[2 * numberSets_ + 4 * sizeof(int)];
269    memset(status_, 0, numberSets_);
270    int i;
271    for (i = 0; i < numberSets_; i++) {
272      // make slack key
273      setStatus(i, ClpSimplex::basic);
274    }
275    dynamicStatus_ = new unsigned char[2 * numberGubColumns_];
276    memset(dynamicStatus_, 0, numberGubColumns_); // for clarity
277    for (i = 0; i < numberGubColumns_; i++)
278      setDynamicStatus(i, atLowerBound);
279  }
280  toIndex_ = new int[numberSets_];
281  for (iSet = 0; iSet < numberSets_; iSet++)
282    toIndex_[iSet] = -1;
283  fromIndex_ = new int[newRowSize - numberStaticRows_ + 1];
284  numberActiveSets_ = 0;
285  rhsOffset_ = NULL;
286  if (numberGubColumns_) {
287    if (!status) {
288      gubCrash();
289    } else {
290      initialProblem();
291    }
292  }
293  noCheck_ = -1;
294  infeasibilityWeight_ = 0.0;
295}
296
297//-------------------------------------------------------------------
298// Destructor
299//-------------------------------------------------------------------
300ClpDynamicMatrix::~ClpDynamicMatrix()
301{
302  delete[] backToPivotRow_;
303  delete[] keyVariable_;
304  delete[] toIndex_;
305  delete[] fromIndex_;
306  delete[] lowerSet_;
307  delete[] upperSet_;
308  delete[] status_;
309  delete[] startSet_;
310  delete[] next_;
311  delete[] startColumn_;
312  delete[] row_;
313  delete[] element_;
314  delete[] cost_;
315  delete[] id_;
316  delete[] dynamicStatus_;
317  delete[] columnLower_;
318  delete[] columnUpper_;
319}
320
321//----------------------------------------------------------------
322// Assignment operator
323//-------------------------------------------------------------------
324ClpDynamicMatrix &
325ClpDynamicMatrix::operator=(const ClpDynamicMatrix &rhs)
326{
327  if (this != &rhs) {
328    ClpPackedMatrix::operator=(rhs);
329    delete[] backToPivotRow_;
330    delete[] keyVariable_;
331    delete[] toIndex_;
332    delete[] fromIndex_;
333    delete[] lowerSet_;
334    delete[] upperSet_;
335    delete[] status_;
336    delete[] startSet_;
337    delete[] next_;
338    delete[] startColumn_;
339    delete[] row_;
340    delete[] element_;
341    delete[] cost_;
342    delete[] id_;
343    delete[] dynamicStatus_;
344    delete[] columnLower_;
345    delete[] columnUpper_;
346    objectiveOffset_ = rhs.objectiveOffset_;
347    numberSets_ = rhs.numberSets_;
348    numberActiveSets_ = rhs.numberActiveSets_;
349    firstAvailable_ = rhs.firstAvailable_;
350    firstAvailableBefore_ = rhs.firstAvailableBefore_;
351    firstDynamic_ = rhs.firstDynamic_;
352    lastDynamic_ = rhs.lastDynamic_;
353    numberStaticRows_ = rhs.numberStaticRows_;
354    numberElements_ = rhs.numberElements_;
355    backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, lastDynamic_);
356    keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_);
357    toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_);
358    fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + 1 - numberStaticRows_);
359    lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
360    upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
361    status_ = ClpCopyOfArray(rhs.status_, static_cast< int >(2 * numberSets_ + 4 * sizeof(int)));
362    model_ = rhs.model_;
363    sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
364    sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
365    sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
366    sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
367    numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
368    numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
369    savedBestGubDual_ = rhs.savedBestGubDual_;
370    savedBestSet_ = rhs.savedBestSet_;
371    noCheck_ = rhs.noCheck_;
372    infeasibilityWeight_ = rhs.infeasibilityWeight_;
373    // Now secondary data
374    numberGubColumns_ = rhs.numberGubColumns_;
375    maximumGubColumns_ = rhs.maximumGubColumns_;
376    maximumElements_ = rhs.maximumElements_;
377    startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_ + 1);
378    next_ = ClpCopyOfArray(rhs.next_, maximumGubColumns_);
379    startColumn_ = ClpCopyOfArray(rhs.startColumn_, maximumGubColumns_ + 1);
380    row_ = ClpCopyOfArray(rhs.row_, maximumElements_);
381    element_ = ClpCopyOfArray(rhs.element_, maximumElements_);
382    cost_ = ClpCopyOfArray(rhs.cost_, maximumGubColumns_);
383    id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_);
384    columnLower_ = ClpCopyOfArray(rhs.columnLower_, maximumGubColumns_);
385    columnUpper_ = ClpCopyOfArray(rhs.columnUpper_, maximumGubColumns_);
386    dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2 * maximumGubColumns_);
387  }
388  return *this;
389}
390//-------------------------------------------------------------------
391// Clone
392//-------------------------------------------------------------------
393ClpMatrixBase *ClpDynamicMatrix::clone() const
394{
395  return new ClpDynamicMatrix(*this);
396}
397// Partial pricing
398void ClpDynamicMatrix::partialPricing(ClpSimplex *model, double startFraction, double endFraction,
399  int &bestSequence, int &numberWanted)
400{
401  numberWanted = currentWanted_;
402  assert(!model->rowScale());
403  if (numberSets_) {
404    // Do packed part before gub
405    // always???
406    //printf("normal packed price - start %d end %d (passed end %d, first %d)\n",
407    ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
408  } else {
409    // no gub
410    ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
411    return;
412  }
413  if (numberWanted > 0) {
414    // and do some proportion of full set
415    int startG2 = static_cast< int >(startFraction * numberSets_);
416    int endG2 = static_cast< int >(endFraction * numberSets_ + 0.1);
417    endG2 = CoinMin(endG2, numberSets_);
418    //printf("gub price - set start %d end %d\n",
419    //   startG2,endG2);
420    double tolerance = model->currentDualTolerance();
421    double *reducedCost = model->djRegion();
422    const double *duals = model->dualRowSolution();
423    double bestDj;
424    int numberRows = model->numberRows();
425    int slackOffset = lastDynamic_ + numberRows;
426    int structuralOffset = slackOffset + numberSets_;
427    // If nothing found yet can go all the way to end
428    int endAll = endG2;
429    if (bestSequence < 0 && !startG2)
430      endAll = numberSets_;
431    if (bestSequence >= 0) {
432      if (bestSequence != savedBestSequence_)
433        bestDj = fabs(reducedCost[bestSequence]); // dj from slacks or permanent
434      else
435        bestDj = savedBestDj_;
436    } else {
437      bestDj = tolerance;
438    }
439    int saveSequence = bestSequence;
440    double djMod = 0.0;
441    double bestDjMod = 0.0;
442    //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
443    //     startG2,endG2,numberWanted);
444    int bestSet = -1;
445#if 0
446          // make sure first available is clean (in case last iteration rejected)
447          cost[firstAvailable_] = 0.0;
448          length[firstAvailable_] = 0;
449          model->nonLinearCost()->setOne(firstAvailable_, 0.0, 0.0, COIN_DBL_MAX, 0.0);
450          model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
451          {
452               for (int i = firstAvailable_; i < lastDynamic_; i++)
453                    assert(!cost[i]);
454          }
455#endif
456    int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_;
457    int minNeg = minimumGoodReducedCosts_ < 0 ? 5 : minimumGoodReducedCosts_;
458    for (int iSet = startG2; iSet < endAll; iSet++) {
459      if (numberWanted + minNeg < originalWanted_ && iSet > startG2 + minSet) {
460        // give up
461        numberWanted = 0;
462        break;
463      } else if (iSet == endG2 && bestSequence >= 0) {
464        break;
465      }
466      int gubRow = toIndex_[iSet];
467      if (gubRow >= 0) {
468        djMod = duals[gubRow + numberStaticRows_]; // have I got sign right?
469      } else {
470        int iBasic = keyVariable_[iSet];
471        if (iBasic >= maximumGubColumns_) {
472          djMod = 0.0; // set not in
473        } else {
474          // get dj without
475          djMod = 0.0;
476          for (CoinBigIndex j = startColumn_[iBasic];
477               j < startColumn_[iBasic + 1]; j++) {
478            int jRow = row_[j];
479            djMod -= duals[jRow] * element_[j];
480          }
481          djMod += cost_[iBasic];
482          // See if gub slack possible - dj is djMod
483          if (getStatus(iSet) == ClpSimplex::atLowerBound) {
484            double value = -djMod;
485            if (value > tolerance) {
486              numberWanted--;
487              if (value > bestDj) {
488                // check flagged variable and correct dj
489                if (!flagged(iSet)) {
490                  bestDj = value;
491                  bestSequence = slackOffset + iSet;
492                  bestDjMod = djMod;
493                  bestSet = iSet;
494                } else {
495                  // just to make sure we don't exit before got something
496                  numberWanted++;
497                  abort();
498                }
499              }
500            }
501          } else if (getStatus(iSet) == ClpSimplex::atUpperBound) {
502            double value = djMod;
503            if (value > tolerance) {
504              numberWanted--;
505              if (value > bestDj) {
506                // check flagged variable and correct dj
507                if (!flagged(iSet)) {
508                  bestDj = value;
509                  bestSequence = slackOffset + iSet;
510                  bestDjMod = djMod;
511                  bestSet = iSet;
512                } else {
513                  // just to make sure we don't exit before got something
514                  numberWanted++;
515                  abort();
516                }
517              }
518            }
519          }
520        }
521      }
522      int iSequence = startSet_[iSet];
523      while (iSequence >= 0) {
524        DynamicStatus status = getDynamicStatus(iSequence);
525        if (status == atLowerBound || status == atUpperBound) {
526          double value = cost_[iSequence] - djMod;
527          for (CoinBigIndex j = startColumn_[iSequence];
528               j < startColumn_[iSequence + 1]; j++) {
529            int jRow = row_[j];
530            value -= duals[jRow] * element_[j];
531          }
532          // change sign if at lower bound
533          if (status == atLowerBound)
534            value = -value;
535          if (value > tolerance) {
536            numberWanted--;
537            if (value > bestDj) {
538              // check flagged variable and correct dj
539              if (!flagged(iSequence)) {
540                if (false /*status == atLowerBound
541                                                &&keyValue(iSet)<1.0e-7*/
542                ) {
543                  // can't come in because
544                  // of ones at ub
545                  numberWanted++;
546                } else {
547
548                  bestDj = value;
549                  bestSequence = structuralOffset + iSequence;
550                  bestDjMod = djMod;
551                  bestSet = iSet;
552                }
553              } else {
554                // just to make sure we don't exit before got something
555                numberWanted++;
556              }
557            }
558          }
559        }
560        iSequence = next_[iSequence]; //onto next in set
561      }
562      if (numberWanted <= 0) {
563        numberWanted = 0;
564        break;
565      }
566    }
567    if (bestSequence != saveSequence) {
568      savedBestGubDual_ = bestDjMod;
569      savedBestDj_ = bestDj;
570      savedBestSequence_ = bestSequence;
571      savedBestSet_ = bestSet;
572    }
573    // See if may be finished
574    if (!startG2 && bestSequence < 0)
575      infeasibilityWeight_ = model_->infeasibilityCost();
576    else if (bestSequence >= 0)
577      infeasibilityWeight_ = -1.0;
578  }
579  currentWanted_ = numberWanted;
580}
581/* Returns effective RHS if it is being used.  This is used for long problems
582   or big gub or anywhere where going through full columns is
583   expensive.  This may re-compute */
584double *
585ClpDynamicMatrix::rhsOffset(ClpSimplex *model, bool forceRefresh,
586  bool
587#ifdef CLP_DEBUG
588    check
589#endif
590)
591{
592  // forceRefresh=true;printf("take out forceRefresh\n");
593  if (!model_->numberIterations())
594    forceRefresh = true;
595    //check=false;
596#ifdef CLP_DEBUG
597  double *saveE = NULL;
598  if (rhsOffset_ && check) {
599    int numberRows = model->numberRows();
600    saveE = new double[numberRows];
601  }
602#endif
603  if (rhsOffset_) {
604#ifdef CLP_DEBUG
605    if (check) {
606      // no need - but check anyway
607      int numberRows = model->numberRows();
608      double *rhs = new double[numberRows];
609      int iRow;
610      int iSet;
611      CoinZeroN(rhs, numberRows);
612      // do ones at bounds before gub
613      const double *smallSolution = model->solutionRegion();
614      const double *element = matrix_->getElements();
615      const int *row = matrix_->getIndices();
616      const CoinBigIndex *startColumn = matrix_->getVectorStarts();
617      const int *length = matrix_->getVectorLengths();
618      int iColumn;
619      double objectiveOffset = 0.0;
620      for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
621        if (model->getStatus(iColumn) != ClpSimplex::basic) {
622          double value = smallSolution[iColumn];
623          for (CoinBigIndex j = startColumn[iColumn];
624               j < startColumn[iColumn] + length[iColumn]; j++) {
625            int jRow = row[j];
626            rhs[jRow] -= value * element[j];
627          }
628        }
629      }
630      if (columnLower_ || columnUpper_) {
631        double *solution = new double[numberGubColumns_];
632        for (iSet = 0; iSet < numberSets_; iSet++) {
633          int j = startSet_[iSet];
634          while (j >= 0) {
635            double value = 0.0;
636            if (getDynamicStatus(j) != inSmall) {
637              if (getDynamicStatus(j) == atLowerBound) {
638                if (columnLower_)
639                  value = columnLower_[j];
640              } else if (getDynamicStatus(j) == atUpperBound) {
641                value = columnUpper_[j];
642              } else if (getDynamicStatus(j) == soloKey) {
643                value = keyValue(iSet);
644              }
645              objectiveOffset += value * cost_[j];
646            }
647            solution[j] = value;
648            j = next_[j]; //onto next in set
649          }
650        }
651        // ones in gub and in small problem
652        for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
653          if (model_->getStatus(iColumn) != ClpSimplex::basic) {
654            int jFull = id_[iColumn - firstDynamic_];
655            solution[jFull] = smallSolution[iColumn];
656          }
657        }
658        for (iSet = 0; iSet < numberSets_; iSet++) {
659          int kRow = toIndex_[iSet];
660          if (kRow >= 0)
661            kRow += numberStaticRows_;
662          int j = startSet_[iSet];
663          while (j >= 0) {
664            double value = solution[j];
665            if (value) {
666              for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j + 1]; k++) {
667                int iRow = row_[k];
668                rhs[iRow] -= element_[k] * value;
669              }
670              if (kRow >= 0)
671                rhs[kRow] -= value;
672            }
673            j = next_[j]; //onto next in set
674          }
675        }
676        delete[] solution;
677      } else {
678        // bounds
679        ClpSimplex::Status iStatus;
680        for (int iSet = 0; iSet < numberSets_; iSet++) {
681          int kRow = toIndex_[iSet];
682          if (kRow < 0) {
683            int iColumn = keyVariable_[iSet];
684            if (iColumn < maximumGubColumns_) {
685              // key is not treated as basic
686              double b = 0.0;
687              // key is structural - where is slack
688              iStatus = getStatus(iSet);
689              assert(iStatus != ClpSimplex::basic);
690              if (iStatus == ClpSimplex::atLowerBound)
691                b = lowerSet_[iSet];
692              else
693                b = upperSet_[iSet];
694              if (b) {
695                objectiveOffset += b * cost_[iColumn];
696                for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn + 1]; j++) {
697                  int iRow = row_[j];
698                  rhs[iRow] -= element_[j] * b;
699                }
700              }
701            }
702          }
703        }
704      }
705      if (fabs(model->objectiveOffset() - objectiveOffset_ - objectiveOffset) > 1.0e-1)
706        printf("old offset %g, true %g\n", model->objectiveOffset(),
707          objectiveOffset_ - objectiveOffset);
708      for (iRow = 0; iRow < numberRows; iRow++) {
709        if (fabs(rhs[iRow] - rhsOffset_[iRow]) > 1.0e-3)
710          printf("** bad effective %d - true %g old %g\n", iRow, rhs[iRow], rhsOffset_[iRow]);
711      }
712      CoinMemcpyN(rhs, numberRows, saveE);
713      delete[] rhs;
714    }
715#endif
716    if (forceRefresh || (refreshFrequency_ && model->numberIterations() >= lastRefresh_ + refreshFrequency_)) {
717      int numberRows = model->numberRows();
718      int iSet;
719      CoinZeroN(rhsOffset_, numberRows);
720      // do ones at bounds before gub
721      const double *smallSolution = model->solutionRegion();
722      const double *element = matrix_->getElements();
723      const int *row = matrix_->getIndices();
724      const CoinBigIndex *startColumn = matrix_->getVectorStarts();
725      const int *length = matrix_->getVectorLengths();
726      int iColumn;
727      double objectiveOffset = 0.0;
728      for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
729        if (model->getStatus(iColumn) != ClpSimplex::basic) {
730          double value = smallSolution[iColumn];
731          for (CoinBigIndex j = startColumn[iColumn];
732               j < startColumn[iColumn] + length[iColumn]; j++) {
733            int jRow = row[j];
734            rhsOffset_[jRow] -= value * element[j];
735          }
736        }
737      }
738      if (columnLower_ || columnUpper_) {
739        double *solution = new double[numberGubColumns_];
740        for (iSet = 0; iSet < numberSets_; iSet++) {
741          int j = startSet_[iSet];
742          while (j >= 0) {
743            double value = 0.0;
744            if (getDynamicStatus(j) != inSmall) {
745              if (getDynamicStatus(j) == atLowerBound) {
746                if (columnLower_)
747                  value = columnLower_[j];
748              } else if (getDynamicStatus(j) == atUpperBound) {
749                value = columnUpper_[j];
750                assert(value < 1.0e30);
751              } else if (getDynamicStatus(j) == soloKey) {
752                value = keyValue(iSet);
753              }
754              objectiveOffset += value * cost_[j];
755            }
756            solution[j] = value;
757            j = next_[j]; //onto next in set
758          }
759        }
760        // ones in gub and in small problem
761        for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
762          if (model_->getStatus(iColumn) != ClpSimplex::basic) {
763            int jFull = id_[iColumn - firstDynamic_];
764            solution[jFull] = smallSolution[iColumn];
765          }
766        }
767        for (iSet = 0; iSet < numberSets_; iSet++) {
768          int kRow = toIndex_[iSet];
769          if (kRow >= 0)
770            kRow += numberStaticRows_;
771          int j = startSet_[iSet];
772          while (j >= 0) {
773            //? use DynamicStatus status = getDynamicStatus(j);
774            double value = solution[j];
775            if (value) {
776              for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j + 1]; k++) {
777                int iRow = row_[k];
778                rhsOffset_[iRow] -= element_[k] * value;
779              }
780              if (kRow >= 0)
781                rhsOffset_[kRow] -= value;
782            }
783            j = next_[j]; //onto next in set
784          }
785        }
786        delete[] solution;
787      } else {
788        // bounds
789        ClpSimplex::Status iStatus;
790        for (int iSet = 0; iSet < numberSets_; iSet++) {
791          int kRow = toIndex_[iSet];
792          if (kRow < 0) {
793            int iColumn = keyVariable_[iSet];
794            if (iColumn < maximumGubColumns_) {
795              // key is not treated as basic
796              double b = 0.0;
797              // key is structural - where is slack
798              iStatus = getStatus(iSet);
799              assert(iStatus != ClpSimplex::basic);
800              if (iStatus == ClpSimplex::atLowerBound)
801                b = lowerSet_[iSet];
802              else
803                b = upperSet_[iSet];
804              if (b) {
805                objectiveOffset += b * cost_[iColumn];
806                for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn + 1]; j++) {
807                  int iRow = row_[j];
808                  rhsOffset_[iRow] -= element_[j] * b;
809                }
810              }
811            }
812          }
813        }
814      }
815      model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
816#ifdef CLP_DEBUG
817      if (saveE) {
818        int iRow;
819        for (iRow = 0; iRow < numberRows; iRow++) {
820          if (fabs(saveE[iRow] - rhsOffset_[iRow]) > 1.0e-3)
821            printf("** %d - old eff %g new %g\n", iRow, saveE[iRow], rhsOffset_[iRow]);
822        }
823        delete[] saveE;
824      }
825#endif
826      lastRefresh_ = model->numberIterations();
827    }
828  }
829  return rhsOffset_;
830}
831/*
832  update information for a pivot (and effective rhs)
833*/
834int ClpDynamicMatrix::updatePivot(ClpSimplex *model, double oldInValue, double oldOutValue)
835{
836
837  // now update working model
838  int sequenceIn = model->sequenceIn();
839  int sequenceOut = model->sequenceOut();
840  int numberColumns = model->numberColumns();
841  if (sequenceIn != sequenceOut && sequenceIn < numberColumns)
842    backToPivotRow_[sequenceIn] = model->pivotRow();
843  if (sequenceIn >= firstDynamic_ && sequenceIn < numberColumns) {
844    int bigSequence = id_[sequenceIn - firstDynamic_];
845    if (getDynamicStatus(bigSequence) != inSmall) {
846      firstAvailable_++;
847      setDynamicStatus(bigSequence, inSmall);
848    }
849  }
850  // make sure slack is synchronized
851  if (sequenceIn >= numberColumns + numberStaticRows_) {
852    int iDynamic = sequenceIn - numberColumns - numberStaticRows_;
853    int iSet = fromIndex_[iDynamic];
854    setStatus(iSet, model->getStatus(sequenceIn));
855  }
856  if (sequenceOut >= numberColumns + numberStaticRows_) {
857    int iDynamic = sequenceOut - numberColumns - numberStaticRows_;
858    int iSet = fromIndex_[iDynamic];
859    // out may have gone through barrier - so check
860    double valueOut = model->lowerRegion()[sequenceOut];
861    if (fabs(valueOut - static_cast< double >(lowerSet_[iSet])) < fabs(valueOut - static_cast< double >(upperSet_[iSet])))
862      setStatus(iSet, ClpSimplex::atLowerBound);
863    else
864      setStatus(iSet, ClpSimplex::atUpperBound);
865    if (lowerSet_[iSet] == upperSet_[iSet])
866      setStatus(iSet, ClpSimplex::isFixed);
867#if 0
868          if (getStatus(iSet) != model->getStatus(sequenceOut))
869               printf("** set %d status %d, var status %d\n", iSet,
870                      getStatus(iSet), model->getStatus(sequenceOut));
871#endif
872  }
873  ClpMatrixBase::updatePivot(model, oldInValue, oldOutValue);
874#ifdef CLP_DEBUG
875  char *inSmall = new char[numberGubColumns_];
876  memset(inSmall, 0, numberGubColumns_);
877  const double *solution = model->solutionRegion();
878  for (int i = 0; i < numberGubColumns_; i++)
879    if (getDynamicStatus(i) == ClpDynamicMatrix::inSmall)
880      inSmall[i] = 1;
881  for (int i = firstDynamic_; i < firstAvailable_; i++) {
882    int k = id_[i - firstDynamic_];
883    inSmall[k] = 0;
884    //if (k>=23289&&k<23357&&solution[i])
885    //printf("var %d (in small %d) has value %g\n",k,i,solution[i]);
886  }
887  for (int i = 0; i < numberGubColumns_; i++)
888    assert(!inSmall[i]);
889  delete[] inSmall;
890#ifndef NDEBUG
891  for (int i = 0; i < numberActiveSets_; i++) {
892    int iSet = fromIndex_[i];
893    assert(toIndex_[iSet] == i);
894    //if (getStatus(iSet)!=model->getRowStatus(i+numberStaticRows_))
895    //printf("*** set %d status %d, var status %d\n",iSet,
896    //     getStatus(iSet),model->getRowStatus(i+numberStaticRows_));
897    //assert (model->getRowStatus(i+numberStaticRows_)==getStatus(iSet));
898    //if (iSet==1035) {
899    //printf("rhs for set %d (%d) is %g %g - cost %g\n",iSet,i,model->lowerRegion(0)[i+numberStaticRows_],
900    //     model->upperRegion(0)[i+numberStaticRows_],model->costRegion(0)[i+numberStaticRows_]);
901    //}
902  }
903#endif
904#endif
905  if (numberStaticRows_ + numberActiveSets_ < model->numberRows())
906    return 0;
907  else
908    return 1;
909}
910/*
911     utility dual function for dealing with dynamic constraints
912     mode=n see ClpGubMatrix.hpp for definition
913     Remember to update here when settled down
914*/
915void ClpDynamicMatrix::dualExpanded(ClpSimplex *model,
916  CoinIndexedVector * /*array*/,
917  double * /*other*/, int mode)
918{
919  switch (mode) {
920    // modify costs before transposeUpdate
921  case 0:
922    break;
923    // create duals for key variables (without check on dual infeasible)
924  case 1:
925    break;
926    // as 1 but check slacks and compute djs
927  case 2: {
928    // do pivots
929    int *pivotVariable = model->pivotVariable();
930    int numberRows = numberStaticRows_ + numberActiveSets_;
931    int numberColumns = model->numberColumns();
932    for (int iRow = 0; iRow < numberRows; iRow++) {
933      int iPivot = pivotVariable[iRow];
934      if (iPivot < numberColumns)
935        backToPivotRow_[iPivot] = iRow;
936    }
937    if (noCheck_ >= 0) {
938      if (infeasibilityWeight_ != model_->infeasibilityCost()) {
939        // don't bother checking
940        sumDualInfeasibilities_ = 100.0;
941        numberDualInfeasibilities_ = 1;
942        sumOfRelaxedDualInfeasibilities_ = 100.0;
943        return;
944      }
945    }
946    // In theory we should subtract out ones we have done but ....
947    // If key slack then dual 0.0
948    // If not then slack could be dual infeasible
949    // dj for key is zero so that defines dual on set
950    int i;
951    double *dual = model->dualRowSolution();
952    double dualTolerance = model->dualTolerance();
953    double relaxedTolerance = dualTolerance;
954    // we can't really trust infeasibilities if there is dual error
955    double error = CoinMin(1.0e-2, model->largestDualError());
956    // allow tolerance at least slightly bigger than standard
957    relaxedTolerance = relaxedTolerance + error;
958    // but we will be using difference
959    relaxedTolerance -= dualTolerance;
960    sumDualInfeasibilities_ = 0.0;
961    numberDualInfeasibilities_ = 0;
962    sumOfRelaxedDualInfeasibilities_ = 0.0;
963    for (i = 0; i < numberSets_; i++) {
964      double value = 0.0;
965      int gubRow = toIndex_[i];
966      if (gubRow < 0) {
967        int kColumn = keyVariable_[i];
968        if (kColumn < maximumGubColumns_) {
969          // dj without set
970          value = cost_[kColumn];
971          for (CoinBigIndex j = startColumn_[kColumn];
972               j < startColumn_[kColumn + 1]; j++) {
973            int iRow = row_[j];
974            value -= dual[iRow] * element_[j];
975          }
976          double infeasibility = 0.0;
977          if (getStatus(i) == ClpSimplex::atLowerBound) {
978            if (-value > dualTolerance)
979              infeasibility = -value - dualTolerance;
980          } else if (getStatus(i) == ClpSimplex::atUpperBound) {
981            if (value > dualTolerance)
982              infeasibility = value - dualTolerance;
983          }
984          if (infeasibility > 0.0) {
985            sumDualInfeasibilities_ += infeasibility;
986            if (infeasibility > relaxedTolerance)
987              sumOfRelaxedDualInfeasibilities_ += infeasibility;
988            numberDualInfeasibilities_++;
989          }
990        }
991      } else {
992        value = dual[gubRow + numberStaticRows_];
993      }
994      // Now subtract out from all
995      int k = startSet_[i];
996      while (k >= 0) {
997        if (getDynamicStatus(k) != inSmall) {
998          double djValue = cost_[k] - value;
999          for (CoinBigIndex j = startColumn_[k];
1000               j < startColumn_[k + 1]; j++) {
1001            int iRow = row_[j];
1002            djValue -= dual[iRow] * element_[j];
1003          }
1004          double infeasibility = 0.0;
1005          if (getDynamicStatus(k) == atLowerBound) {
1006            if (djValue < -dualTolerance)
1007              infeasibility = -djValue - dualTolerance;
1008          } else if (getDynamicStatus(k) == atUpperBound) {
1009            // at upper bound
1010            if (djValue > dualTolerance)
1011              infeasibility = djValue - dualTolerance;
1012          }
1013          if (infeasibility > 0.0) {
1014            sumDualInfeasibilities_ += infeasibility;
1015            if (infeasibility > relaxedTolerance)
1016              sumOfRelaxedDualInfeasibilities_ += infeasibility;
1017            numberDualInfeasibilities_++;
1018          }
1019        }
1020        k = next_[k]; //onto next in set
1021      }
1022    }
1023  }
1024    infeasibilityWeight_ = -1.0;
1025    break;
1026  // Report on infeasibilities of key variables
1027  case 3: {
1028    model->setSumDualInfeasibilities(model->sumDualInfeasibilities() + sumDualInfeasibilities_);
1029    model->setNumberDualInfeasibilities(model->numberDualInfeasibilities() + numberDualInfeasibilities_);
1030    model->setSumOfRelaxedDualInfeasibilities(model->sumOfRelaxedDualInfeasibilities() + sumOfRelaxedDualInfeasibilities_);
1031  } break;
1032  // modify costs before transposeUpdate for partial pricing
1033  case 4:
1034    break;
1035  }
1036}
1037/*
1038     general utility function for dealing with dynamic constraints
1039     mode=n see ClpGubMatrix.hpp for definition
1040     Remember to update here when settled down
1041*/
1042int ClpDynamicMatrix::generalExpanded(ClpSimplex *model, int mode, int &number)
1043{
1044  int returnCode = 0;
1045#if 0 //ndef NDEBUG
1046     {
1047       int numberColumns = model->numberColumns();
1048       int numberRows = model->numberRows();
1049       int * pivotVariable = model->pivotVariable();
1050       if (pivotVariable&&model->numberIterations()) {
1051         for (int i=numberStaticRows_+numberActiveSets_;i<numberRows;i++) {
1052           assert (pivotVariable[i]==i+numberColumns);
1053         }
1054       }
1055     }
1056#endif
1057  switch (mode) {
1058    // Fill in pivotVariable
1059  case 0: {
1060    // If no effective rhs - form it
1061    if (!rhsOffset_) {
1062      rhsOffset_ = new double[model->numberRows()];
1063      rhsOffset(model, true);
1064    }
1065    int i;
1066    int numberBasic = number;
1067    int numberColumns = model->numberColumns();
1068    // Use different array so can build from true pivotVariable_
1069    //int * pivotVariable = model->pivotVariable();
1070    int *pivotVariable = model->rowArray(0)->getIndices();
1071    for (i = 0; i < numberColumns; i++) {
1072      if (model->getColumnStatus(i) == ClpSimplex::basic)
1073        pivotVariable[numberBasic++] = i;
1074    }
1075    number = numberBasic;
1076  } break;
1077  // Do initial extra rows + maximum basic
1078  case 2: {
1079    number = model->numberRows();
1080  } break;
1081  // Before normal replaceColumn
1082  case 3: {
1083    if (numberActiveSets_ + numberStaticRows_ == model_->numberRows()) {
1084      // no space - re-factorize
1085      returnCode = 4;
1086      number = -1; // say no need for normal replaceColumn
1087    }
1088  } break;
1089  // To see if can dual or primal
1090  case 4: {
1091    returnCode = 1;
1092  } break;
1093  // save status
1094  case 5: {
1095    memcpy(status_ + numberSets_, status_, numberSets_);
1096    memcpy(status_ + 2 * numberSets_, &numberActiveSets_, sizeof(int));
1097    memcpy(dynamicStatus_ + maximumGubColumns_,
1098      dynamicStatus_, maximumGubColumns_);
1099  } break;
1100  // restore status
1101  case 6: {
1102    memcpy(status_, status_ + numberSets_, numberSets_);
1103    memcpy(&numberActiveSets_, status_ + 2 * numberSets_, sizeof(int));
1104    memcpy(dynamicStatus_, dynamicStatus_ + maximumGubColumns_,
1105      maximumGubColumns_);
1106    initialProblem();
1107  } break;
1108  // unflag all variables
1109  case 8: {
1110    for (int i = 0; i < numberGubColumns_; i++) {
1111      if (flagged(i)) {
1112        unsetFlagged(i);
1113        returnCode++;
1114      }
1115    }
1116  } break;
1117  // redo costs in primal
1118  case 9: {
1119    double *cost = model->costRegion();
1120    double *solution = model->solutionRegion();
1121    double *columnLower = model->lowerRegion();
1122    double *columnUpper = model->upperRegion();
1123    int i;
1124    bool doCosts = (number & 4) != 0;
1125    bool doBounds = (number & 1) != 0;
1126    for (i = firstDynamic_; i < firstAvailable_; i++) {
1127      int jColumn = id_[i - firstDynamic_];
1128      if (doBounds) {
1129        if (!columnLower_ && !columnUpper_) {
1130          columnLower[i] = 0.0;
1131          columnUpper[i] = COIN_DBL_MAX;
1132        } else {
1133          if (columnLower_)
1134            columnLower[i] = columnLower_[jColumn];
1135          else
1136            columnLower[i] = 0.0;
1137          if (columnUpper_)
1138            columnUpper[i] = columnUpper_[jColumn];
1139          else
1140            columnUpper[i] = COIN_DBL_MAX;
1141        }
1142      }
1143      if (doCosts) {
1144        cost[i] = cost_[jColumn];
1145        // Original bounds
1146        if (model->nonLinearCost())
1147          model->nonLinearCost()->setOne(i, solution[i],
1148            this->columnLower(jColumn),
1149            this->columnUpper(jColumn), cost_[jColumn]);
1150      }
1151    }
1152    // and active sets
1153    for (i = 0; i < numberActiveSets_; i++) {
1154      int iSet = fromIndex_[i];
1155      int iSequence = lastDynamic_ + numberStaticRows_ + i;
1156      if (doBounds) {
1157        if (lowerSet_[iSet] > -1.0e20)
1158          columnLower[iSequence] = lowerSet_[iSet];
1159        else
1160          columnLower[iSequence] = -COIN_DBL_MAX;
1161        if (upperSet_[iSet] < 1.0e20)
1162          columnUpper[iSequence] = upperSet_[iSet];
1163        else
1164          columnUpper[iSequence] = COIN_DBL_MAX;
1165      }
1166      if (doCosts) {
1167        if (model->nonLinearCost()) {
1168          double trueLower;
1169          if (lowerSet_[iSet] > -1.0e20)
1170            trueLower = lowerSet_[iSet];
1171          else
1172            trueLower = -COIN_DBL_MAX;
1173          double trueUpper;
1174          if (upperSet_[iSet] < 1.0e20)
1175            trueUpper = upperSet_[iSet];
1176          else
1177            trueUpper = COIN_DBL_MAX;
1178          model->nonLinearCost()->setOne(iSequence, solution[iSequence],
1179            trueLower, trueUpper, 0.0);
1180        }
1181      }
1182    }
1183  } break;
1184  // return 1 if there may be changing bounds on variable (column generation)
1185  case 10: {
1186    // return 1 as bounds on rhs will change
1187    returnCode = 1;
1188  } break;
1189  // make sure set is clean
1190  case 7: {
1191    // first flag
1192    if (number >= firstDynamic_ && number < lastDynamic_) {
1193      int sequence = id_[number - firstDynamic_];
1194      setFlagged(sequence);
1195      //model->clearFlagged(number);
1196    } else if (number >= model_->numberColumns() + numberStaticRows_) {
1197      // slack
1198      int iSet = fromIndex_[number - model_->numberColumns() - numberStaticRows_];
1199      setFlaggedSlack(iSet);
1200      //model->clearFlagged(number);
1201    }
1202  }
1203  case 11: {
1204    //int sequenceIn = model->sequenceIn();
1205    if (number >= firstDynamic_ && number < lastDynamic_) {
1206      //assert (number == model->sequenceIn());
1207      // take out variable (but leave key)
1208      double *cost = model->costRegion();
1209      double *columnLower = model->lowerRegion();
1210      double *columnUpper = model->upperRegion();
1211      double *solution = model->solutionRegion();
1212      int *length = matrix_->getMutableVectorLengths();
1213#ifndef NDEBUG
1214      if (length[number]) {
1215        int *row = matrix_->getMutableIndices();
1216        CoinBigIndex *startColumn = matrix_->getMutableVectorStarts();
1217        int iRow = row[startColumn[number] + length[number] - 1];
1218        int which = iRow - numberStaticRows_;
1219        assert(which >= 0);
1220        int iSet = fromIndex_[which];
1221        assert(toIndex_[iSet] == which);
1222      }
1223#endif
1224      // no need firstAvailable_--;
1225      solution[firstAvailable_] = 0.0;
1226      cost[firstAvailable_] = 0.0;
1227      length[firstAvailable_] = 0;
1228      model->nonLinearCost()->setOne(firstAvailable_, 0.0, 0.0, COIN_DBL_MAX, 0.0);
1229      model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
1230      columnLower[firstAvailable_] = 0.0;
1231      columnUpper[firstAvailable_] = COIN_DBL_MAX;
1232
1233      // not really in small problem
1234      int iBig = id_[number - firstDynamic_];
1235      if (model->getStatus(number) == ClpSimplex::atLowerBound) {
1236        setDynamicStatus(iBig, atLowerBound);
1237        if (columnLower_)
1238          modifyOffset(number, columnLower_[iBig]);
1239      } else {
1240        setDynamicStatus(iBig, atUpperBound);
1241        modifyOffset(number, columnUpper_[iBig]);
1242      }
1243    } else if (number >= model_->numberColumns() + numberStaticRows_) {
1244      // slack
1245      int iSet = fromIndex_[number - model_->numberColumns() - numberStaticRows_];
1246      printf("what now - set %d\n", iSet);
1247    }
1248  } break;
1249  default:
1250    break;
1251  }
1252  return returnCode;
1253}
1254/* Purely for column generation and similar ideas.  Allows
1255   matrix and any bounds or costs to be updated (sensibly).
1256   Returns non-zero if any changes.
1257*/
1258int ClpDynamicMatrix::refresh(ClpSimplex *model)
1259{
1260  // If at beginning then create initial problem
1261  if (firstDynamic_ == firstAvailable_) {
1262    initialProblem();
1263    return 1;
1264  } else if (!model->nonLinearCost()) {
1265    // will be same as last time
1266    return 1;
1267  }
1268#ifndef NDEBUG
1269  {
1270    int numberColumns = model->numberColumns();
1271    int numberRows = model->numberRows();
1272    int *pivotVariable = model->pivotVariable();
1273    for (int i = numberStaticRows_ + numberActiveSets_; i < numberRows; i++) {
1274      assert(pivotVariable[i] == i + numberColumns);
1275    }
1276  }
1277#endif
1278  // lookup array
1279  int *active = new int[numberActiveSets_];
1280  CoinZeroN(active, numberActiveSets_);
1281  int iColumn;
1282  int numberColumns = model->numberColumns();
1283  double *element = matrix_->getMutableElements();
1284  int *row = matrix_->getMutableIndices();
1285  CoinBigIndex *startColumn = matrix_->getMutableVectorStarts();
1286  int *length = matrix_->getMutableVectorLengths();
1287  double *cost = model->costRegion();
1288  double *columnLower = model->lowerRegion();
1289  double *columnUpper = model->upperRegion();
1290  CoinBigIndex numberElements = startColumn[firstDynamic_];
1291  // first just do lookup and basic stuff
1292  int currentNumber = firstAvailable_;
1293  firstAvailable_ = firstDynamic_;
1294  double objectiveChange = 0.0;
1295  double *solution = model->solutionRegion();
1296  int currentNumberActiveSets = numberActiveSets_;
1297  for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
1298    int iRow = row[startColumn[iColumn] + length[iColumn] - 1];
1299    int which = iRow - numberStaticRows_;
1300    assert(which >= 0);
1301    int iSet = fromIndex_[which];
1302    assert(toIndex_[iSet] == which);
1303    if (model->getStatus(iColumn) == ClpSimplex::basic) {
1304      active[which]++;
1305      // may as well make key
1306      keyVariable_[iSet] = id_[iColumn - firstDynamic_];
1307    }
1308  }
1309  int i;
1310  numberActiveSets_ = 0;
1311  int numberDeleted = 0;
1312  for (i = 0; i < currentNumberActiveSets; i++) {
1313    int iRow = i + numberStaticRows_;
1314    int numberActive = active[i];
1315    int iSet = fromIndex_[i];
1316    if (model->getRowStatus(iRow) == ClpSimplex::basic) {
1317      numberActive++;
1318      // may as well make key
1319      keyVariable_[iSet] = maximumGubColumns_ + iSet;
1320    }
1321    if (numberActive > 1) {
1322      // keep
1323      active[i] = numberActiveSets_;
1324      numberActiveSets_++;
1325    } else {
1326      active[i] = -1;
1327    }
1328  }
1329
1330  for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
1331    int iRow = row[startColumn[iColumn] + length[iColumn] - 1];
1332    int which = iRow - numberStaticRows_;
1333    int jColumn = id_[iColumn - firstDynamic_];
1334    if (model->getStatus(iColumn) == ClpSimplex::atLowerBound || model->getStatus(iColumn) == ClpSimplex::atUpperBound) {
1335      double value = solution[iColumn];
1336      bool toLowerBound = true;
1337      assert(jColumn >= 0);
1338      assert(iColumn >= 0);
1339      if (columnUpper_) {
1340        if (!columnLower_) {
1341          if (fabs(value - columnUpper_[jColumn]) < fabs(value))
1342            toLowerBound = false;
1343        } else if (fabs(value - columnUpper_[jColumn]) < fabs(value - columnLower_[jColumn])) {
1344          toLowerBound = false;
1345        }
1346      }
1347      if (toLowerBound) {
1348        // throw out to lower bound
1349        if (columnLower_) {
1350          setDynamicStatus(jColumn, atLowerBound);
1351          // treat solution as if exactly at a bound
1352          double value = columnLower[iColumn];
1353          objectiveChange += cost[iColumn] * value;
1354        } else {
1355          setDynamicStatus(jColumn, atLowerBound);
1356        }
1357      } else {
1358        // throw out to upper bound
1359        setDynamicStatus(jColumn, atUpperBound);
1360        double value = columnUpper[iColumn];
1361        objectiveChange += cost[iColumn] * value;
1362      }
1363      numberDeleted++;
1364    } else {
1365      assert(model->getStatus(iColumn) != ClpSimplex::superBasic); // deal with later
1366      int iPut = active[which];
1367      if (iPut >= 0) {
1368        // move
1369        id_[firstAvailable_ - firstDynamic_] = jColumn;
1370        int numberThis = startColumn_[jColumn + 1] - startColumn_[jColumn] + 1;
1371        length[firstAvailable_] = numberThis;
1372        cost[firstAvailable_] = cost[iColumn];
1373        columnLower[firstAvailable_] = columnLower[iColumn];
1374        columnUpper[firstAvailable_] = columnUpper[iColumn];
1375        model->nonLinearCost()->setOne(firstAvailable_, solution[iColumn], 0.0, COIN_DBL_MAX,
1376          cost_[jColumn]);
1377        CoinBigIndex base = startColumn_[jColumn];
1378        for (int j = 0; j < numberThis - 1; j++) {
1379          row[numberElements] = row_[base + j];
1380          element[numberElements++] = element_[base + j];
1381        }
1382        row[numberElements] = iPut + numberStaticRows_;
1383        element[numberElements++] = 1.0;
1384        model->setStatus(firstAvailable_, model->getStatus(iColumn));
1385        solution[firstAvailable_] = solution[iColumn];
1386        firstAvailable_++;
1387        startColumn[firstAvailable_] = numberElements;
1388      }
1389    }
1390  }
1391  // zero solution
1392  CoinZeroN(solution + firstAvailable_, currentNumber - firstAvailable_);
1393  // zero cost
1394  CoinZeroN(cost + firstAvailable_, currentNumber - firstAvailable_);
1395  // zero lengths
1396  CoinZeroN(length + firstAvailable_, currentNumber - firstAvailable_);
1397  for (iColumn = firstAvailable_; iColumn < currentNumber; iColumn++) {
1398    model->nonLinearCost()->setOne(iColumn, 0.0, 0.0, COIN_DBL_MAX, 0.0);
1399    model->setStatus(iColumn, ClpSimplex::atLowerBound);
1400    columnLower[iColumn] = 0.0;
1401    columnUpper[iColumn] = COIN_DBL_MAX;
1402  }
1403  // move rhs and set rest to infinite
1404  numberActiveSets_ = 0;
1405  for (i = 0; i < currentNumberActiveSets; i++) {
1406    int iSet = fromIndex_[i];
1407    assert(toIndex_[iSet] == i);
1408    int iRow = i + numberStaticRows_;
1409    int iPut = active[i];
1410    if (iPut >= 0) {
1411      int in = iRow + numberColumns;
1412      int out = iPut + numberColumns + numberStaticRows_;
1413      solution[out] = solution[in];
1414      columnLower[out] = columnLower[in];
1415      columnUpper[out] = columnUpper[in];
1416      cost[out] = cost[in];
1417      double trueLower;
1418      if (lowerSet_[iSet] > -1.0e20)
1419        trueLower = lowerSet_[iSet];
1420      else
1421        trueLower = -COIN_DBL_MAX;
1422      double trueUpper;
1423      if (upperSet_[iSet] < 1.0e20)
1424        trueUpper = upperSet_[iSet];
1425      else
1426        trueUpper = COIN_DBL_MAX;
1427      model->nonLinearCost()->setOne(out, solution[out],
1428        trueLower, trueUpper, 0.0);
1429      model->setStatus(out, model->getStatus(in));
1430      toIndex_[iSet] = numberActiveSets_;
1431      rhsOffset_[numberActiveSets_ + numberStaticRows_] = rhsOffset_[i + numberStaticRows_];
1432      fromIndex_[numberActiveSets_++] = fromIndex_[i];
1433    } else {
1434      // adjust offset
1435      // put out as key
1436      int jColumn = keyVariable_[iSet];
1437      toIndex_[iSet] = -1;
1438      if (jColumn < maximumGubColumns_) {
1439        setDynamicStatus(jColumn, soloKey);
1440        double value = keyValue(iSet);
1441        objectiveChange += cost_[jColumn] * value;
1442        modifyOffset(jColumn, -value);
1443      }
1444    }
1445  }
1446  model->setObjectiveOffset(model->objectiveOffset() - objectiveChange);
1447  delete[] active;
1448  for (i = numberActiveSets_; i < currentNumberActiveSets; i++) {
1449    int iSequence = i + numberStaticRows_ + numberColumns;
1450    solution[iSequence] = 0.0;
1451    columnLower[iSequence] = -COIN_DBL_MAX;
1452    columnUpper[iSequence] = COIN_DBL_MAX;
1453    cost[iSequence] = 0.0;
1454    model->nonLinearCost()->setOne(iSequence, solution[iSequence],
1455      columnLower[iSequence],
1456      columnUpper[iSequence], 0.0);
1457    model->setStatus(iSequence, ClpSimplex::basic);
1458    rhsOffset_[i + numberStaticRows_] = 0.0;
1459  }
1460  if (currentNumberActiveSets != numberActiveSets_ || numberDeleted) {
1461    // deal with pivotVariable
1462    int *pivotVariable = model->pivotVariable();
1463    int sequence = firstDynamic_;
1464    int iRow = 0;
1465    int base1 = firstDynamic_;
1466    int base2 = lastDynamic_;
1467    int base3 = numberColumns + numberStaticRows_;
1468    int numberRows = numberStaticRows_ + currentNumberActiveSets;
1469    while (sequence < firstAvailable_) {
1470      int iPivot = pivotVariable[iRow];
1471      while (iPivot < base1 || (iPivot >= base2 && iPivot < base3)) {
1472        iPivot = pivotVariable[++iRow];
1473      }
1474      pivotVariable[iRow++] = sequence;
1475      sequence++;
1476    }
1477    // move normal basic ones down
1478    int iPut = iRow;
1479    for (; iRow < numberRows; iRow++) {
1480      int iPivot = pivotVariable[iRow];
1481      if (iPivot < base1 || (iPivot >= base2 && iPivot < base3)) {
1482        pivotVariable[iPut++] = iPivot;
1483      }
1484    }
1485    // and basic others
1486    for (i = 0; i < numberActiveSets_; i++) {
1487      if (model->getRowStatus(i + numberStaticRows_) == ClpSimplex::basic) {
1488        pivotVariable[iPut++] = i + base3;
1489      }
1490    }
1491    if (iPut < numberStaticRows_ + numberActiveSets_) {
1492      printf("lost %d sets\n",
1493        numberStaticRows_ + numberActiveSets_ - iPut);
1494      iPut = numberStaticRows_ + numberActiveSets_;
1495    }
1496    for (i = numberActiveSets_; i < currentNumberActiveSets; i++) {
1497      pivotVariable[iPut++] = i + base3;
1498    }
1499    //assert (iPut == numberRows);
1500  }
1501#ifdef CLP_DEBUG
1502#if 0
1503     printf("row for set 244 is %d, row status %d value %g ", toIndex_[244], status_[244],
1504            keyValue(244));
1505     int jj = startSet_[244];
1506     while (jj >= 0) {
1507          printf("var %d status %d ", jj, dynamicStatus_[jj]);
1508          jj = next_[jj];
1509     }
1510     printf("\n");
1511#endif
1512#ifdef CLP_DEBUG
1513  {
1514    // debug coding to analyze set
1515    int i;
1516    int kSet = -6;
1517    if (kSet >= 0) {
1518      int *back = new int[numberGubColumns_];
1519      for (i = 0; i < numberGubColumns_; i++)
1520        back[i] = -1;
1521      for (i = firstDynamic_; i < firstAvailable_; i++)
1522        back[id_[i - firstDynamic_]] = i;
1523      int sequence = startSet_[kSet];
1524      if (toIndex_[kSet] < 0) {
1525        printf("Not in - Set %d - slack status %d, key %d\n", kSet, status_[kSet], keyVariable_[kSet]);
1526        while (sequence >= 0) {
1527          printf("( %d status %d ) ", sequence, getDynamicStatus(sequence));
1528          sequence = next_[sequence];
1529        }
1530      } else {
1531        int iRow = numberStaticRows_ + toIndex_[kSet];
1532        printf("In - Set %d - slack status %d, key %d offset %g slack %g\n", kSet, status_[kSet], keyVariable_[kSet],
1533          rhsOffset_[iRow], model->solutionRegion(0)[iRow]);
1534        while (sequence >= 0) {
1535          int iBack = back[sequence];
1536          printf("( %d status %d value %g) ", sequence, getDynamicStatus(sequence), model->solutionRegion()[iBack]);
1537          sequence = next_[sequence];
1538        }
1539      }
1540      printf("\n");
1541      delete[] back;
1542    }
1543  }
1544#endif
1545  int n = numberActiveSets_;
1546  for (i = 0; i < numberSets_; i++) {
1547    if (toIndex_[i] < 0) {
1548      //assert(keyValue(i)>=lowerSet_[i]&&keyValue(i)<=upperSet_[i]);
1549      n++;
1550    }
1551    int k = 0;
1552    for (int j = startSet_[i]; j < startSet_[i + 1]; j++) {
1553      if (getDynamicStatus(j) == soloKey)
1554        k++;
1555    }
1556    assert(k < 2);
1557  }
1558  assert(n == numberSets_);
1559#endif
1560  return 1;
1561}
1562void ClpDynamicMatrix::times(double scalar,
1563  const double *x, double *y) const
1564{
1565  if (model_->specialOptions() != 16) {
1566    ClpPackedMatrix::times(scalar, x, y);
1567  } else {
1568    int iRow;
1569    const double *element = matrix_->getElements();
1570    const int *row = matrix_->getIndices();
1571    const CoinBigIndex *startColumn = matrix_->getVectorStarts();
1572    const int *length = matrix_->getVectorLengths();
1573    int *pivotVariable = model_->pivotVariable();
1574    for (iRow = 0; iRow < numberStaticRows_ + numberActiveSets_; iRow++) {
1575      y[iRow] -= scalar * rhsOffset_[iRow];
1576      int iColumn = pivotVariable[iRow];
1577      if (iColumn < lastDynamic_) {
1578        CoinBigIndex j;
1579        double value = scalar * x[iColumn];
1580        if (value) {
1581          for (j = startColumn[iColumn];
1582               j < startColumn[iColumn] + length[iColumn]; j++) {
1583            int jRow = row[j];
1584            y[jRow] += value * element[j];
1585          }
1586        }
1587      }
1588    }
1589  }
1590}
1591// Modifies rhs offset
1592void ClpDynamicMatrix::modifyOffset(int sequence, double amount)
1593{
1594  if (amount) {
1595    assert(rhsOffset_);
1596    CoinBigIndex j;
1597    for (j = startColumn_[sequence]; j < startColumn_[sequence + 1]; j++) {
1598      int iRow = row_[j];
1599      rhsOffset_[iRow] += amount * element_[j];
1600    }
1601  }
1602}
1603// Gets key value when none in small
1604double
1605ClpDynamicMatrix::keyValue(int iSet) const
1606{
1607  double value = 0.0;
1608  if (toIndex_[iSet] < 0) {
1609    int key = keyVariable_[iSet];
1610    if (key < maximumGubColumns_) {
1611      if (getStatus(iSet) == ClpSimplex::atLowerBound)
1612        value = lowerSet_[iSet];
1613      else
1614        value = upperSet_[iSet];
1615      int numberKey = 0;
1616      int j = startSet_[iSet];
1617      while (j >= 0) {
1618        DynamicStatus status = getDynamicStatus(j);
1619        assert(status != inSmall);
1620        if (status == soloKey) {
1621          numberKey++;
1622        } else if (status == atUpperBound) {
1623          value -= columnUpper_[j];
1624        } else if (columnLower_) {
1625          value -= columnLower_[j];
1626        }
1627        j = next_[j]; //onto next in set
1628      }
1629      assert(numberKey == 1);
1630    } else {
1631      int j = startSet_[iSet];
1632      while (j >= 0) {
1633        DynamicStatus status = getDynamicStatus(j);
1634        assert(status != inSmall);
1635        assert(status != soloKey);
1636        if (status == atUpperBound) {
1637          value += columnUpper_[j];
1638        } else if (columnLower_) {
1639          value += columnLower_[j];
1640        }
1641        j = next_[j]; //onto next in set
1642      }
1643#if 0
1644               // slack must be feasible
1645               double oldValue=value;
1646               value = CoinMax(value,lowerSet_[iSet]);
1647               value = CoinMin(value,upperSet_[iSet]);
1648               if (value!=oldValue)
1649                 printf("using %g (not %g) for slack on set %d (%g,%g)\n",
1650                        value,oldValue,iSet,lowerSet_[iSet],upperSet_[iSet]);
1651#endif
1652    }
1653  }
1654  return value;
1655}
1656// Switches off dj checking each factorization (for BIG models)
1657void ClpDynamicMatrix::switchOffCheck()
1658{
1659  noCheck_ = 0;
1660  infeasibilityWeight_ = 0.0;
1661}
1662/* Creates a variable.  This is called after partial pricing and may modify matrix.
1663   May update bestSequence.
1664*/
1665void ClpDynamicMatrix::createVariable(ClpSimplex *model, int &bestSequence)
1666{
1667  int numberRows = model->numberRows();
1668  int slackOffset = lastDynamic_ + numberRows;
1669  int structuralOffset = slackOffset + numberSets_;
1670  int bestSequence2 = savedBestSequence_ - structuralOffset;
1671  if (bestSequence >= slackOffset) {
1672    double *columnLower = model->lowerRegion();
1673    double *columnUpper = model->upperRegion();
1674    double *solution = model->solutionRegion();
1675    double *reducedCost = model->djRegion();
1676    const double *duals = model->dualRowSolution();
1677    if (toIndex_[savedBestSet_] < 0) {
1678      // need to put key into basis
1679      int newRow = numberActiveSets_ + numberStaticRows_;
1680      model->dualRowSolution()[newRow] = savedBestGubDual_;
1681      double valueOfKey = keyValue(savedBestSet_); // done before toIndex_ set
1682      toIndex_[savedBestSet_] = numberActiveSets_;
1683      fromIndex_[numberActiveSets_++] = savedBestSet_;
1684      int iSequence = lastDynamic_ + newRow;
1685      // we need to get lower and upper correct
1686      double shift = 0.0;
1687      int j = startSet_[savedBestSet_];
1688      while (j >= 0) {
1689        if (getDynamicStatus(j) == atUpperBound)
1690          shift += columnUpper_[j];
1691        else if (getDynamicStatus(j) == atLowerBound && columnLower_)
1692          shift += columnLower_[j];
1693        j = next_[j]; //onto next in set
1694      }
1695      if (lowerSet_[savedBestSet_] > -1.0e20)
1696        columnLower[iSequence] = lowerSet_[savedBestSet_];
1697      else
1698        columnLower[iSequence] = -COIN_DBL_MAX;
1699      if (upperSet_[savedBestSet_] < 1.0e20)
1700        columnUpper[iSequence] = upperSet_[savedBestSet_];
1701      else
1702        columnUpper[iSequence] = COIN_DBL_MAX;
1703#ifdef CLP_DEBUG
1704      if (model_->logLevel() == 63) {
1705        printf("%d in in set %d, key is %d rhs %g %g - keyvalue %g\n",
1706          bestSequence2, savedBestSet_, keyVariable_[savedBestSet_],
1707          columnLower[iSequence], columnUpper[iSequence], valueOfKey);
1708        int j = startSet_[savedBestSet_];
1709        while (j >= 0) {
1710          if (getDynamicStatus(j) == atUpperBound)
1711            printf("%d atup ", j);
1712          else if (getDynamicStatus(j) == atLowerBound)
1713            printf("%d atlo ", j);
1714          else if (getDynamicStatus(j) == soloKey)
1715            printf("%d solo ", j);
1716          else
1717            abort();
1718          j = next_[j]; //onto next in set
1719        }
1720        printf("\n");
1721      }
1722#endif
1723      if (keyVariable_[savedBestSet_] < maximumGubColumns_) {
1724        // slack not key
1725        model_->pivotVariable()[newRow] = firstAvailable_;
1726        backToPivotRow_[firstAvailable_] = newRow;
1727        model->setStatus(iSequence, getStatus(savedBestSet_));
1728        model->djRegion()[iSequence] = savedBestGubDual_;
1729        solution[iSequence] = valueOfKey;
1730        // create variable and pivot in
1731        int key = keyVariable_[savedBestSet_];
1732        setDynamicStatus(key, inSmall);
1733        double *element = matrix_->getMutableElements();
1734        int *row = matrix_->getMutableIndices();
1735        CoinBigIndex *startColumn = matrix_->getMutableVectorStarts();
1736        int *length = matrix_->getMutableVectorLengths();
1737        CoinBigIndex numberElements = startColumn[firstAvailable_];
1738        int numberThis = startColumn_[key + 1] - startColumn_[key] + 1;
1739        if (numberElements + numberThis > numberElements_) {
1740          // need to redo
1741          numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
1742          matrix_->reserve(lastDynamic_, numberElements_);
1743          element = matrix_->getMutableElements();
1744          row = matrix_->getMutableIndices();
1745          // these probably okay but be safe
1746          startColumn = matrix_->getMutableVectorStarts();
1747          length = matrix_->getMutableVectorLengths();
1748        }
1749        // already set startColumn[firstAvailable_]=numberElements;
1750        length[firstAvailable_] = numberThis;
1751        model->costRegion()[firstAvailable_] = cost_[key];
1752        CoinBigIndex base = startColumn_[key];
1753        for (int j = 0; j < numberThis - 1; j++) {
1754          row[numberElements] = row_[base + j];
1755          element[numberElements++] = element_[base + j];
1756        }
1757        row[numberElements] = newRow;
1758        element[numberElements++] = 1.0;
1759        id_[firstAvailable_ - firstDynamic_] = key;
1760        model->setObjectiveOffset(model->objectiveOffset() + cost_[key] * valueOfKey);
1761        model->solutionRegion()[firstAvailable_] = valueOfKey;
1762        model->setStatus(firstAvailable_, ClpSimplex::basic);
1763        // ***** need to adjust effective rhs
1764        if (!columnLower_ && !columnUpper_) {
1765          columnLower[firstAvailable_] = 0.0;
1766          columnUpper[firstAvailable_] = COIN_DBL_MAX;
1767        } else {
1768          if (columnLower_)
1769            columnLower[firstAvailable_] = columnLower_[key];
1770          else
1771            columnLower[firstAvailable_] = 0.0;
1772          if (columnUpper_)
1773            columnUpper[firstAvailable_] = columnUpper_[key];
1774          else
1775            columnUpper[firstAvailable_] = COIN_DBL_MAX;
1776        }
1777        model->nonLinearCost()->setOne(firstAvailable_, solution[firstAvailable_],
1778          columnLower[firstAvailable_],
1779          columnUpper[firstAvailable_], cost_[key]);
1780        startColumn[firstAvailable_ + 1] = numberElements;
1781        reducedCost[firstAvailable_] = 0.0;
1782        modifyOffset(key, valueOfKey);
1783        rhsOffset_[newRow] = -shift; // sign?
1784#ifdef CLP_DEBUG
1785        model->rowArray(1)->checkClear();
1786#endif
1787        // now pivot in
1788        unpack(model, model->rowArray(1), firstAvailable_);
1789        model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(1));
1790        double alpha = model->rowArray(1)->denseVector()[newRow];
1791        int updateStatus = model->factorization()->replaceColumn(model,
1792          model->rowArray(2),
1793          model->rowArray(1),
1794          newRow, alpha);
1795        model->rowArray(1)->clear();
1796        if (updateStatus) {
1797          if (updateStatus == 3) {
1798            // out of memory
1799            // increase space if not many iterations
1800            if (model->factorization()->pivots() < 0.5 * model->factorization()->maximumPivots() && model->factorization()->pivots() < 400)
1801              model->factorization()->areaFactor(
1802                model->factorization()->areaFactor() * 1.1);
1803          } else {
1804            printf("Bad returncode %d from replaceColumn\n", updateStatus);
1805          }
1806          bestSequence = -1;
1807          return;
1808        }
1809        // firstAvailable_ only finally updated if good pivot (in updatePivot)
1810        // otherwise it reverts to firstAvailableBefore_
1811        firstAvailable_++;
1812      } else {
1813        // slack key
1814        model->setStatus(iSequence, ClpSimplex::basic);
1815        model->djRegion()[iSequence] = 0.0;
1816        solution[iSequence] = valueOfKey + shift;
1817        rhsOffset_[newRow] = -shift; // sign?
1818      }
1819      // correct slack
1820      model->costRegion()[iSequence] = 0.0;
1821      model->nonLinearCost()->setOne(iSequence, solution[iSequence], columnLower[iSequence],
1822        columnUpper[iSequence], 0.0);
1823    }
1824    if (savedBestSequence_ >= structuralOffset) {
1825      // recompute dj and create
1826      double value = cost_[bestSequence2] - savedBestGubDual_;
1827      for (CoinBigIndex jBigIndex = startColumn_[bestSequence2];
1828           jBigIndex < startColumn_[bestSequence2 + 1]; jBigIndex++) {
1829        int jRow = row_[jBigIndex];
1830        value -= duals[jRow] * element_[jBigIndex];
1831      }
1832      int gubRow = toIndex_[savedBestSet_] + numberStaticRows_;
1833      double *element = matrix_->getMutableElements();
1834      int *row = matrix_->getMutableIndices();
1835      CoinBigIndex *startColumn = matrix_->getMutableVectorStarts();
1836      int *length = matrix_->getMutableVectorLengths();
1837      CoinBigIndex numberElements = startColumn[firstAvailable_];
1838      int numberThis = startColumn_[bestSequence2 + 1] - startColumn_[bestSequence2] + 1;
1839      if (numberElements + numberThis > numberElements_) {
1840        // need to redo
1841        numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
1842        matrix_->reserve(lastDynamic_, numberElements_);
1843        element = matrix_->getMutableElements();
1844        row = matrix_->getMutableIndices();
1845        // these probably okay but be safe
1846        startColumn = matrix_->getMutableVectorStarts();
1847        length = matrix_->getMutableVectorLengths();
1848      }
1849      // already set startColumn[firstAvailable_]=numberElements;
1850      length[firstAvailable_] = numberThis;
1851      model->costRegion()[firstAvailable_] = cost_[bestSequence2];
1852      CoinBigIndex base = startColumn_[bestSequence2];
1853      for (int j = 0; j < numberThis - 1; j++) {
1854        row[numberElements] = row_[base + j];
1855        element[numberElements++] = element_[base + j];
1856      }
1857      row[numberElements] = gubRow;
1858      element[numberElements++] = 1.0;
1859      id_[firstAvailable_ - firstDynamic_] = bestSequence2;
1860      //printf("best %d\n",bestSequence2);
1861      model->solutionRegion()[firstAvailable_] = 0.0;
1862      model->clearFlagged(firstAvailable_);
1863      if (!columnLower_ && !columnUpper_) {
1864        model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
1865        columnLower[firstAvailable_] = 0.0;
1866        columnUpper[firstAvailable_] = COIN_DBL_MAX;
1867      } else {
1868        DynamicStatus status = getDynamicStatus(bestSequence2);
1869        if (columnLower_)
1870          columnLower[firstAvailable_] = columnLower_[bestSequence2];
1871        else
1872          columnLower[firstAvailable_] = 0.0;
1873        if (columnUpper_)
1874          columnUpper[firstAvailable_] = columnUpper_[bestSequence2];
1875        else
1876          columnUpper[firstAvailable_] = COIN_DBL_MAX;
1877        if (status == atLowerBound) {
1878          solution[firstAvailable_] = columnLower[firstAvailable_];
1879          model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
1880        } else {
1881          solution[firstAvailable_] = columnUpper[firstAvailable_];
1882          model->setStatus(firstAvailable_, ClpSimplex::atUpperBound);
1883        }
1884      }
1885      model->setObjectiveOffset(model->objectiveOffset() + cost_[bestSequence2] * solution[firstAvailable_]);
1886      model->nonLinearCost()->setOne(firstAvailable_, solution[firstAvailable_],
1887        columnLower[firstAvailable_],
1888        columnUpper[firstAvailable_], cost_[bestSequence2]);
1889      bestSequence = firstAvailable_;
1890      // firstAvailable_ only updated if good pivot (in updatePivot)
1891      startColumn[firstAvailable_ + 1] = numberElements;
1892      //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,savedBestGubDual_);
1893      reducedCost[bestSequence] = value;
1894    } else {
1895      // slack - row must just have been created
1896      assert(toIndex_[savedBestSet_] == numberActiveSets_ - 1);
1897      int newRow = numberStaticRows_ + numberActiveSets_ - 1;
1898      bestSequence = lastDynamic_ + newRow;
1899      reducedCost[bestSequence] = savedBestGubDual_;
1900    }
1901  }
1902  // clear for next iteration
1903  savedBestSequence_ = -1;
1904}
1905// Returns reduced cost of a variable
1906double
1907ClpDynamicMatrix::reducedCost(ClpSimplex *model, int sequence) const
1908{
1909  int numberRows = model->numberRows();
1910  int slackOffset = lastDynamic_ + numberRows;
1911  if (sequence < slackOffset)
1912    return model->djRegion()[sequence];
1913  else
1914    return savedBestDj_;
1915}
1916// Does gub crash
1917void ClpDynamicMatrix::gubCrash()
1918{
1919  // Do basis - cheapest or slack if feasible
1920  int longestSet = 0;
1921  int iSet;
1922  for (iSet = 0; iSet < numberSets_; iSet++) {
1923    int n = 0;
1924    int j = startSet_[iSet];
1925    while (j >= 0) {
1926      n++;
1927      j = next_[j];
1928    }
1929    longestSet = CoinMax(longestSet, n);
1930  }
1931  double *upper = new double[longestSet + 1];
1932  double *cost = new double[longestSet + 1];
1933  double *lower = new double[longestSet + 1];
1934  double *solution = new double[longestSet + 1];
1935  int *back = new int[longestSet + 1];
1936  double tolerance = model_->primalTolerance();
1937  double objectiveOffset = 0.0;
1938  for (iSet = 0; iSet < numberSets_; iSet++) {
1939    int iBasic = -1;
1940    double value = 0.0;
1941    // find cheapest
1942    int numberInSet = 0;
1943    int j = startSet_[iSet];
1944    while (j >= 0) {
1945      if (!columnLower_)
1946        lower[numberInSet] = 0.0;
1947      else
1948        lower[numberInSet] = columnLower_[j];
1949      if (!columnUpper_)
1950        upper[numberInSet] = COIN_DBL_MAX;
1951      else
1952        upper[numberInSet] = columnUpper_[j];
1953      back[numberInSet++] = j;
1954      j = next_[j];
1955    }
1956    CoinFillN(solution, numberInSet, 0.0);
1957    // and slack
1958    iBasic = numberInSet;
1959    solution[iBasic] = -value;
1960    lower[iBasic] = -upperSet_[iSet];
1961    upper[iBasic] = -lowerSet_[iSet];
1962    int kphase;
1963    if (value < lowerSet_[iSet] - tolerance || value > upperSet_[iSet] + tolerance) {
1964      // infeasible
1965      kphase = 0;
1966      // remember bounds are flipped so opposite to natural
1967      if (value < lowerSet_[iSet] - tolerance)
1968        cost[iBasic] = 1.0;
1969      else
1970        cost[iBasic] = -1.0;
1971      CoinZeroN(cost, numberInSet);
1972      double dualTolerance = model_->dualTolerance();
1973      for (int iphase = kphase; iphase < 2; iphase++) {
1974        if (iphase) {
1975          cost[numberInSet] = 0.0;
1976          for (int j = 0; j < numberInSet; j++)
1977            cost[j] = cost_[back[j]];
1978        }
1979        // now do one row lp
1980        bool improve = true;
1981        while (improve) {
1982          improve = false;
1983          double dual = cost[iBasic];
1984          int chosen = -1;
1985          double best = dualTolerance;
1986          int way = 0;
1987          for (int i = 0; i <= numberInSet; i++) {
1988            double dj = cost[i] - dual;
1989            double improvement = 0.0;
1990            if (iphase || i < numberInSet)
1991              assert(solution[i] >= lower[i] && solution[i] <= upper[i]);
1992            if (dj > dualTolerance)
1993              improvement = dj * (solution[i] - lower[i]);
1994            else if (dj < -dualTolerance)
1995              improvement = dj * (solution[i] - upper[i]);
1996            if (improvement > best) {
1997              best = improvement;
1998              chosen = i;
1999              if (dj < 0.0) {
2000                way = 1;
2001              } else {
2002                way = -1;
2003              }
2004            }
2005          }
2006          if (chosen >= 0) {
2007            improve = true;
2008            // now see how far
2009            if (way > 0) {
2010              // incoming increasing so basic decreasing
2011              // if phase 0 then go to nearest bound
2012              double distance = upper[chosen] - solution[chosen];
2013              double basicDistance;
2014              if (!iphase) {
2015                assert(iBasic == numberInSet);
2016                assert(solution[iBasic] > upper[iBasic]);
2017                basicDistance = solution[iBasic] - upper[iBasic];
2018              } else {
2019                basicDistance = solution[iBasic] - lower[iBasic];
2020              }
2021              // need extra coding for unbounded
2022              assert(CoinMin(distance, basicDistance) < 1.0e20);
2023              if (distance > basicDistance) {
2024                // incoming becomes basic
2025                solution[chosen] += basicDistance;
2026                if (!iphase)
2027                  solution[iBasic] = upper[iBasic];
2028                else
2029                  solution[iBasic] = lower[iBasic];
2030                iBasic = chosen;
2031              } else {
2032                // flip
2033                solution[chosen] = upper[chosen];
2034                solution[iBasic] -= distance;
2035              }
2036            } else {
2037              // incoming decreasing so basic increasing
2038              // if phase 0 then go to nearest bound
2039              double distance = solution[chosen] - lower[chosen];
2040              double basicDistance;
2041              if (!iphase) {
2042                assert(iBasic == numberInSet);
2043                assert(solution[iBasic] < lower[iBasic]);
2044                basicDistance = lower[iBasic] - solution[iBasic];
2045              } else {
2046                basicDistance = upper[iBasic] - solution[iBasic];
2047              }
2048              // need extra coding for unbounded - for now just exit
2049              if (CoinMin(distance, basicDistance) > 1.0e20) {
2050                printf("unbounded on set %d\n", iSet);
2051                iphase = 1;
2052                iBasic = numberInSet;
2053                break;
2054              }
2055              if (distance > basicDistance) {
2056                // incoming becomes basic
2057                solution[chosen] -= basicDistance;
2058                if (!iphase)
2059                  solution[iBasic] = lower[iBasic];
2060                else
2061                  solution[iBasic] = upper[iBasic];
2062                iBasic = chosen;
2063              } else {
2064                // flip
2065                solution[chosen] = lower[chosen];
2066                solution[iBasic] += distance;
2067              }
2068            }
2069            if (!iphase) {
2070              if (iBasic < numberInSet)
2071                break; // feasible
2072              else if (solution[iBasic] >= lower[iBasic] && solution[iBasic] <= upper[iBasic])
2073                break; // feasible (on flip)
2074            }
2075          }
2076        }
2077      }
2078    }
2079    // do solution i.e. bounds
2080    if (columnLower_ || columnUpper_) {
2081      for (int j = 0; j < numberInSet; j++) {
2082        if (j != iBasic) {
2083          objectiveOffset += solution[j] * cost[j];
2084          if (columnLower_ && columnUpper_) {
2085            if (fabs(solution[j] - columnLower_[back[j]]) > fabs(solution[j] - columnUpper_[back[j]]))
2086              setDynamicStatus(back[j], atUpperBound);
2087          } else if (columnUpper_ && solution[j] > 0.0) {
2088            setDynamicStatus(back[j], atUpperBound);
2089          } else {
2090            setDynamicStatus(back[j], atLowerBound);
2091            assert(!solution[j]);
2092          }
2093        }
2094      }
2095    }
2096    // convert iBasic back and do bounds
2097    if (iBasic == numberInSet) {
2098      // slack basic
2099      setStatus(iSet, ClpSimplex::basic);
2100      iBasic = iSet + maximumGubColumns_;
2101    } else {
2102      iBasic = back[iBasic];
2103      setDynamicStatus(iBasic, soloKey);
2104      // remember bounds flipped
2105      if (upper[numberInSet] == lower[numberInSet])
2106        setStatus(iSet, ClpSimplex::isFixed);
2107      else if (solution[numberInSet] == upper[numberInSet])
2108        setStatus(iSet, ClpSimplex::atLowerBound);
2109      else if (solution[numberInSet] == lower[numberInSet])
2110        setStatus(iSet, ClpSimplex::atUpperBound);
2111      else
2112        abort();
2113    }
2114    keyVariable_[iSet] = iBasic;
2115  }
2116  model_->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
2117  delete[] lower;
2118  delete[] solution;
2119  delete[] upper;
2120  delete[] cost;
2121  delete[] back;
2122  // make sure matrix is in good shape
2123  matrix_->orderMatrix();
2124}
2125// Populates initial matrix from dynamic status
2126void ClpDynamicMatrix::initialProblem()
2127{
2128  int iSet;
2129  double *element = matrix_->getMutableElements();
2130  int *row = matrix_->getMutableIndices();
2131  CoinBigIndex *startColumn = matrix_->getMutableVectorStarts();
2132  int *length = matrix_->getMutableVectorLengths();
2133  double *cost = model_->objective();
2134  double *solution = model_->primalColumnSolution();
2135  double *columnLower = model_->columnLower();
2136  double *columnUpper = model_->columnUpper();
2137  double *rowSolution = model_->primalRowSolution();
2138  double *rowLower = model_->rowLower();
2139  double *rowUpper = model_->rowUpper();
2140  CoinBigIndex numberElements = startColumn[firstDynamic_];
2141
2142  firstAvailable_ = firstDynamic_;
2143  numberActiveSets_ = 0;
2144  for (iSet = 0; iSet < numberSets_; iSet++) {
2145    toIndex_[iSet] = -1;
2146    int numberActive = 0;
2147    int whichKey = -1;
2148    if (getStatus(iSet) == ClpSimplex::basic) {
2149      whichKey = maximumGubColumns_ + iSet;
2150      numberActive = 1;
2151    } else {
2152      whichKey = -1;
2153    }
2154    int j = startSet_[iSet];
2155    while (j >= 0) {
2156      assert(getDynamicStatus(j) != soloKey || whichKey == -1);
2157      if (getDynamicStatus(j) == inSmall) {
2158        numberActive++;
2159      } else if (getDynamicStatus(j) == soloKey) {
2160        whichKey = j;
2161        numberActive++;
2162      }
2163      j = next_[j]; //onto next in set
2164    }
2165    if (numberActive > 1) {
2166      int iRow = numberActiveSets_ + numberStaticRows_;
2167      rowSolution[iRow] = 0.0;
2168      double lowerValue;
2169      if (lowerSet_[iSet] > -1.0e20)
2170        lowerValue = lowerSet_[iSet];
2171      else
2172        lowerValue = -COIN_DBL_MAX;
2173      double upperValue;
2174      if (upperSet_[iSet] < 1.0e20)
2175        upperValue = upperSet_[iSet];
2176      else
2177        upperValue = COIN_DBL_MAX;
2178      rowLower[iRow] = lowerValue;
2179      rowUpper[iRow] = upperValue;
2180      if (getStatus(iSet) == ClpSimplex::basic) {
2181        model_->setRowStatus(iRow, ClpSimplex::basic);
2182        rowSolution[iRow] = 0.0;
2183      } else if (getStatus(iSet) == ClpSimplex::atLowerBound) {
2184        model_->setRowStatus(iRow, ClpSimplex::atLowerBound);
2185        rowSolution[iRow] = lowerValue;
2186      } else {
2187        model_->setRowStatus(iRow, ClpSimplex::atUpperBound);
2188        rowSolution[iRow] = upperValue;
2189      }
2190      j = startSet_[iSet];
2191      while (j >= 0) {
2192        DynamicStatus status = getDynamicStatus(j);
2193        if (status == inSmall) {
2194          int numberThis = startColumn_[j + 1] - startColumn_[j] + 1;
2195          if (numberElements + numberThis > numberElements_) {
2196            // need to redo
2197            numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
2198            matrix_->reserve(lastDynamic_, numberElements_);
2199            element = matrix_->getMutableElements();
2200            row = matrix_->getMutableIndices();
2201            // these probably okay but be safe
2202            startColumn = matrix_->getMutableVectorStarts();
2203            length = matrix_->getMutableVectorLengths();
2204          }
2205          length[firstAvailable_] = numberThis;
2206          cost[firstAvailable_] = cost_[j];
2207          CoinBigIndex base = startColumn_[j];
2208          for (int k = 0; k < numberThis - 1; k++) {
2209            row[numberElements] = row_[base + k];
2210            element[numberElements++] = element_[base + k];
2211          }
2212          row[numberElements] = iRow;
2213          element[numberElements++] = 1.0;
2214          id_[firstAvailable_ - firstDynamic_] = j;
2215          solution[firstAvailable_] = 0.0;
2216          model_->setStatus(firstAvailable_, ClpSimplex::basic);
2217          if (!columnLower_ && !columnUpper_) {
2218            columnLower[firstAvailable_] = 0.0;
2219            columnUpper[firstAvailable_] = COIN_DBL_MAX;
2220          } else {
2221            if (columnLower_)
2222              columnLower[firstAvailable_] = columnLower_[j];
2223            else
2224              columnLower[firstAvailable_] = 0.0;
2225            if (columnUpper_)
2226              columnUpper[firstAvailable_] = columnUpper_[j];
2227            else
2228              columnUpper[firstAvailable_] = COIN_DBL_MAX;
2229            if (status != atUpperBound) {
2230              solution[firstAvailable_] = columnLower[firstAvailable_];
2231            } else {
2232              solution[firstAvailable_] = columnUpper[firstAvailable_];
2233            }
2234          }
2235          firstAvailable_++;
2236          startColumn[firstAvailable_] = numberElements;
2237        }
2238        j = next_[j]; //onto next in set
2239      }
2240      model_->setRowStatus(numberActiveSets_ + numberStaticRows_, getStatus(iSet));
2241      toIndex_[iSet] = numberActiveSets_;
2242      fromIndex_[numberActiveSets_++] = iSet;
2243    } else {
2244      // solo key
2245      bool needKey = false;
2246      if (numberActive) {
2247        if (whichKey < maximumGubColumns_) {
2248          // structural - assume ok
2249          needKey = false;
2250        } else {
2251          // slack
2252          keyVariable_[iSet] = maximumGubColumns_ + iSet;
2253          double value = keyValue(iSet);
2254          if (value < lowerSet_[iSet] - 1.0e-8 || value > upperSet_[iSet] + 1.0e-8)
2255            needKey = true;
2256        }
2257      } else {
2258        needKey = true;
2259      }
2260      if (needKey) {
2261        // all to lb then up some (slack/null if possible)
2262        int length = 99999999;
2263        int which = -1;
2264        double sum = 0.0;
2265        for (int iColumn = startSet_[iSet]; iColumn < startSet_[iSet + 1]; iColumn++) {
2266          setDynamicStatus(iColumn, atLowerBound);
2267          sum += columnLower_[iColumn];
2268          if (length > startColumn_[iColumn + 1] - startColumn_[iColumn]) {
2269            which = iColumn;
2270            length = startColumn_[iColumn + 1] - startColumn_[iColumn];
2271          }
2272        }
2273        if (sum > lowerSet_[iSet] - 1.0e-8) {
2274          // slack can be basic
2275          setStatus(iSet, ClpSimplex::basic);
2276          keyVariable_[iSet] = maximumGubColumns_ + iSet;
2277        } else {
2278          // use shortest
2279          setDynamicStatus(which, soloKey);
2280          keyVariable_[iSet] = which;
2281          setStatus(iSet, ClpSimplex::atLowerBound);
2282        }
2283      }
2284    }
2285    assert(toIndex_[iSet] >= 0 || whichKey >= 0);
2286    keyVariable_[iSet] = whichKey;
2287  }
2288  // clean up pivotVariable
2289  int numberColumns = model_->numberColumns();
2290  int numberRows = model_->numberRows();
2291  int *pivotVariable = model_->pivotVariable();
2292  if (pivotVariable) {
2293    for (int i = 0; i < numberStaticRows_ + numberActiveSets_; i++) {
2294      if (model_->getRowStatus(i) != ClpSimplex::basic)
2295        pivotVariable[i] = -1;
2296      else
2297        pivotVariable[i] = numberColumns + i;
2298    }
2299    for (int i = numberStaticRows_ + numberActiveSets_; i < numberRows; i++) {
2300      pivotVariable[i] = i + numberColumns;
2301    }
2302    int put = -1;
2303    for (int i = 0; i < numberColumns; i++) {
2304      if (model_->getColumnStatus(i) == ClpSimplex::basic) {
2305        while (put < numberRows) {
2306          put++;
2307          if (pivotVariable[put] == -1) {
2308            pivotVariable[put] = i;
2309            break;
2310          }
2311        }
2312      }
2313    }
2314    for (int i = CoinMax(put, 0); i < numberRows; i++) {
2315      if (pivotVariable[i] == -1)
2316        pivotVariable[i] = i + numberColumns;
2317    }
2318  }
2319  if (rhsOffset_ && model_->costRegion()) {
2320    double *cost = model_->costRegion();
2321    double *columnLower = model_->lowerRegion();
2322    double *columnUpper = model_->upperRegion();
2323    double *solution = model_->solutionRegion();
2324    int numberRows = model_->numberRows();
2325    for (int i = numberActiveSets_; i < numberRows - numberStaticRows_; i++) {
2326      int iSequence = i + numberStaticRows_ + numberColumns;
2327      solution[iSequence] = 0.0;
2328      columnLower[iSequence] = -COIN_DBL_MAX;
2329      columnUpper[iSequence] = COIN_DBL_MAX;
2330      cost[iSequence] = 0.0;
2331      model_->nonLinearCost()->setOne(iSequence, solution[iSequence],
2332        columnLower[iSequence],
2333        columnUpper[iSequence], 0.0);
2334      model_->setStatus(iSequence, ClpSimplex::basic);
2335      rhsOffset_[i + numberStaticRows_] = 0.0;
2336    }
2337#if 0
2338       for (int i=0;i<numberStaticRows_;i++)
2339         printf("%d offset %g\n",
2340                i,rhsOffset_[i]);
2341#endif
2342  }
2343  numberActiveColumns_ = firstAvailable_;
2344#if 0
2345     for (iSet = 0; iSet < numberSets_; iSet++) {
2346       for (int j=startSet_[iSet];j<startSet_[iSet+1];j++) {
2347         if (getDynamicStatus(j)==soloKey)
2348           printf("%d in set %d is solo key\n",j,iSet);
2349         else if (getDynamicStatus(j)==inSmall)
2350           printf("%d in set %d is in small\n",j,iSet);
2351       }
2352     }
2353#endif
2354  return;
2355}
2356// Writes out model (without names)
2357void ClpDynamicMatrix::writeMps(const char *name)
2358{
2359  int numberTotalRows = numberStaticRows_ + numberSets_;
2360  int numberTotalColumns = firstDynamic_ + numberGubColumns_;
2361  // over estimate
2362  int numberElements = getNumElements() + startColumn_[numberGubColumns_]
2363    + numberGubColumns_;
2364  double *columnLower = new double[numberTotalColumns];
2365  double *columnUpper = new double[numberTotalColumns];
2366  double *cost = new double[numberTotalColumns];
2367  double *rowLower = new double[numberTotalRows];
2368  double *rowUpper = new double[numberTotalRows];
2369  CoinBigIndex *start = new CoinBigIndex[numberTotalColumns + 1];
2370  int *row = new int[numberElements];
2371  double *element = new double[numberElements];
2372  // Fill in
2373  const CoinBigIndex *startA = getVectorStarts();
2374  const int *lengthA = getVectorLengths();
2375  const int *rowA = getIndices();
2376  const double *elementA = getElements();
2377  const double *columnLowerA = model_->columnLower();
2378  const double *columnUpperA = model_->columnUpper();
2379  const double *costA = model_->objective();
2380  const double *rowLowerA = model_->rowLower();
2381  const double *rowUpperA = model_->rowUpper();
2382  start[0] = 0;
2383  numberElements = 0;
2384  for (int i = 0; i < firstDynamic_; i++) {
2385    columnLower[i] = columnLowerA[i];
2386    columnUpper[i] = columnUpperA[i];
2387    cost[i] = costA[i];
2388    for (CoinBigIndex j = startA[i]; j < startA[i] + lengthA[i]; j++) {
2389      row[numberElements] = rowA[j];
2390      element[numberElements++] = elementA[j];
2391    }
2392    start[i + 1] = numberElements;
2393  }
2394  for (int i = 0; i < numberStaticRows_; i++) {
2395    rowLower[i] = rowLowerA[i];
2396    rowUpper[i] = rowUpperA[i];
2397  }
2398  int putC = firstDynamic_;
2399  int putR = numberStaticRows_;
2400  for (int i = 0; i < numberSets_; i++) {
2401    rowLower[putR] = lowerSet_[i];
2402    rowUpper[putR] = upperSet_[i];
2403    for (CoinBigIndex k = startSet_[i]; k < startSet_[i + 1]; k++) {
2404      columnLower[putC] = columnLower_[k];
2405      columnUpper[putC] = columnUpper_[k];
2406      cost[putC] = cost_[k];
2407      putC++;
2408      for (CoinBigIndex j = startColumn_[k]; j < startColumn_[k + 1]; j++) {
2409        row[numberElements] = row_[j];
2410        element[numberElements++] = element_[j];
2411      }
2412      row[numberElements] = putR;
2413      element[numberElements++] = 1.0;
2414      start[putC] = numberElements;
2415    }
2416    putR++;
2417  }
2418
2419  assert(putR == numberTotalRows);
2420  assert(putC == numberTotalColumns);
2421  ClpSimplex modelOut;
2422  modelOut.loadProblem(numberTotalColumns, numberTotalRows,
2423    start, row, element,
2424    columnLower, columnUpper, cost,
2425    rowLower, rowUpper);
2426  modelOut.writeMps(name);
2427  delete[] columnLower;
2428  delete[] columnUpper;
2429  delete[] cost;
2430  delete[] rowLower;
2431  delete[] rowUpper;
2432  delete[] start;
2433  delete[] row;
2434  delete[] element;
2435}
2436// Adds in a column to gub structure (called from descendant)
2437int ClpDynamicMatrix::addColumn(CoinBigIndex numberEntries, const int *row, const double *element,
2438  double cost, double lower, double upper, int iSet,
2439  DynamicStatus status)
2440{
2441  // check if already in
2442  int j = startSet_[iSet];
2443  while (j >= 0) {
2444    if (startColumn_[j + 1] - startColumn_[j] == numberEntries) {
2445      const int *row2 = row_ + startColumn_[j];
2446      const double *element2 = element_ + startColumn_[j];
2447      bool same = true;
2448      for (int k = 0; k < numberEntries; k++) {
2449        if (row[k] != row2[k] || element[k] != element2[k]) {
2450          same = false;
2451          break;
2452        }
2453      }
2454      if (same) {
2455        bool odd = false;
2456        if (cost != cost_[j])
2457          odd = true;
2458        if (columnLower_ && lower != columnLower_[j])
2459          odd = true;
2460        if (columnUpper_ && upper != columnUpper_[j])
2461          odd = true;
2462        if (odd) {
2463          printf("seems odd - same els but cost,lo,up are %g,%g,%g and %g,%g,%g\n",
2464            cost, lower, upper, cost_[j],
2465            columnLower_ ? columnLower_[j] : 0.0,
2466            columnUpper_ ? columnUpper_[j] : 1.0e100);
2467        } else {
2468          setDynamicStatus(j, status);
2469          return j;
2470        }
2471      }
2472    }
2473    j = next_[j];
2474  }
2475
2476  if (numberGubColumns_ == maximumGubColumns_ || startColumn_[numberGubColumns_] + numberEntries > maximumElements_) {
2477    CoinBigIndex j;
2478    int i;
2479    int put = 0;
2480    int numberElements = 0;
2481    CoinBigIndex start = 0;
2482    // compress - leave ones at ub and basic
2483    int *which = new int[numberGubColumns_];
2484    for (i = 0; i < numberGubColumns_; i++) {
2485      CoinBigIndex end = startColumn_[i + 1];
2486      // what about ubs if column generation?
2487      if (getDynamicStatus(i) != atLowerBound) {
2488        // keep in
2489        for (j = start; j < end; j++) {
2490          row_[numberElements] = row_[j];
2491          element_[numberElements++] = element_[j];
2492        }
2493        startColumn_[put + 1] = numberElements;
2494        cost_[put] = cost_[i];
2495        if (columnLower_)
2496          columnLower_[put] = columnLower_[i];
2497        if (columnUpper_)
2498          columnUpper_[put] = columnUpper_[i];
2499        dynamicStatus_[put] = dynamicStatus_[i];
2500        id_[put] = id_[i];
2501        which[i] = put;
2502        put++;
2503      } else {
2504        // out
2505        which[i] = -1;
2506      }
2507      start = end;
2508    }
2509    // now redo startSet_ and next_
2510    int *newNext = new int[maximumGubColumns_];
2511    for (int jSet = 0; jSet < numberSets_; jSet++) {
2512      int sequence = startSet_[jSet];
2513      while (which[sequence] < 0) {
2514        // out
2515        assert(next_[sequence] >= 0);
2516        sequence = next_[sequence];
2517      }
2518      startSet_[jSet] = which[sequence];
2519      int last = which[sequence];
2520      while (next_[sequence] >= 0) {
2521        sequence = next_[sequence];
2522        if (which[sequence] >= 0) {
2523          // keep
2524          newNext[last] = which[sequence];
2525          last = which[sequence];
2526        }
2527      }
2528      newNext[last] = -jSet - 1;
2529    }
2530    delete[] next_;
2531    next_ = newNext;
2532    delete[] which;
2533    abort();
2534  }
2535  CoinBigIndex start = startColumn_[numberGubColumns_];
2536  CoinMemcpyN(row, numberEntries, row_ + start);
2537  CoinMemcpyN(element, numberEntries, element_ + start);
2538  startColumn_[numberGubColumns_ + 1] = start + numberEntries;
2539  cost_[numberGubColumns_] = cost;
2540  if (columnLower_)
2541    columnLower_[numberGubColumns_] = lower;
2542  else
2543    assert(!lower);
2544  if (columnUpper_)
2545    columnUpper_[numberGubColumns_] = upper;
2546  else
2547    assert(upper > 1.0e20);
2548  setDynamicStatus(numberGubColumns_, status);
2549  // Do next_
2550  j = startSet_[iSet];
2551  startSet_[iSet] = numberGubColumns_;
2552  next_[numberGubColumns_] = j;
2553  numberGubColumns_++;
2554  return numberGubColumns_ - 1;
2555}
2556// Returns which set a variable is in
2557int ClpDynamicMatrix::whichSet(int sequence) const
2558{
2559  while (next_[sequence] >= 0)
2560    sequence = next_[sequence];
2561  int iSet = -next_[sequence] - 1;
2562  return iSet;
2563}
2564
2565/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
2566*/
Note: See TracBrowser for help on using the repository browser.