source: trunk/Clp/src/ClpDynamicMatrix.cpp @ 1732

Last change on this file since 1732 was 1732, checked in by forrest, 8 years ago

various fixes plus slightly weighted pricing plus lagomory switches

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