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

Last change on this file since 1726 was 1726, checked in by stefan, 8 years ago

fix compiler warnings, including one that pointed to a bug

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 112.2 KB
Line 
1/* $Id: ClpDynamicMatrix.cpp 1726 2011-05-02 08:58:39Z stefan $ */
2// Copyright (C) 2004, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6
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                              double distance = 0.0;
2024                              if (iphase || i < numberInSet)
2025                                   assert (solution[i] >= lower[i] && solution[i] <= upper[i]);
2026                              if (dj > dualTolerance)
2027                                   improvement = dj * (solution[i] - lower[i]);
2028                              else if (dj < -dualTolerance)
2029                                   improvement = dj * (solution[i] - upper[i]);
2030                              if (improvement > best) {
2031                                   best = improvement;
2032                                   chosen = i;
2033                                   if (dj < 0.0) {
2034                                        way = 1;
2035                                        distance = upper[i] - solution[i];
2036                                   } else {
2037                                        way = -1;
2038                                        distance = solution[i] - lower[i];
2039                                   }
2040                              }
2041                         }
2042                         if (chosen >= 0) {
2043                              improve = true;
2044                              // now see how far
2045                              if (way > 0) {
2046                                   // incoming increasing so basic decreasing
2047                                   // if phase 0 then go to nearest bound
2048                                   double distance = upper[chosen] - solution[chosen];
2049                                   double basicDistance;
2050                                   if (!iphase) {
2051                                        assert (iBasic == numberInSet);
2052                                        assert (solution[iBasic] > upper[iBasic]);
2053                                        basicDistance = solution[iBasic] - upper[iBasic];
2054                                   } else {
2055                                        basicDistance = solution[iBasic] - lower[iBasic];
2056                                   }
2057                                   // need extra coding for unbounded
2058                                   assert (CoinMin(distance, basicDistance) < 1.0e20);
2059                                   if (distance > basicDistance) {
2060                                        // incoming becomes basic
2061                                        solution[chosen] += basicDistance;
2062                                        if (!iphase)
2063                                             solution[iBasic] = upper[iBasic];
2064                                        else
2065                                             solution[iBasic] = lower[iBasic];
2066                                        iBasic = chosen;
2067                                   } else {
2068                                        // flip
2069                                        solution[chosen] = upper[chosen];
2070                                        solution[iBasic] -= distance;
2071                                   }
2072                              } else {
2073                                   // incoming decreasing so basic increasing
2074                                   // if phase 0 then go to nearest bound
2075                                   double distance = solution[chosen] - lower[chosen];
2076                                   double basicDistance;
2077                                   if (!iphase) {
2078                                        assert (iBasic == numberInSet);
2079                                        assert (solution[iBasic] < lower[iBasic]);
2080                                        basicDistance = lower[iBasic] - solution[iBasic];
2081                                   } else {
2082                                        basicDistance = upper[iBasic] - solution[iBasic];
2083                                   }
2084                                   // need extra coding for unbounded - for now just exit
2085                                   if (CoinMin(distance, basicDistance) > 1.0e20) {
2086                                        printf("unbounded on set %d\n", iSet);
2087                                        iphase = 1;
2088                                        iBasic = numberInSet;
2089                                        break;
2090                                   }
2091                                   if (distance > basicDistance) {
2092                                        // incoming becomes basic
2093                                        solution[chosen] -= basicDistance;
2094                                        if (!iphase)
2095                                             solution[iBasic] = lower[iBasic];
2096                                        else
2097                                             solution[iBasic] = upper[iBasic];
2098                                        iBasic = chosen;
2099                                   } else {
2100                                        // flip
2101                                        solution[chosen] = lower[chosen];
2102                                        solution[iBasic] += distance;
2103                                   }
2104                              }
2105                              if (!iphase) {
2106                                   if(iBasic < numberInSet)
2107                                        break; // feasible
2108                                   else if (solution[iBasic] >= lower[iBasic] &&
2109                                             solution[iBasic] <= upper[iBasic])
2110                                        break; // feasible (on flip)
2111                              }
2112                         }
2113                    }
2114               }
2115          }
2116          // do solution i.e. bounds
2117          if (columnLower_ || columnUpper_) {
2118               for (int j = 0; j < numberInSet; j++) {
2119                    if (j != iBasic) {
2120                         objectiveOffset += solution[j] * cost[j];
2121                         if (columnLower_ && columnUpper_) {
2122                              if (fabs(solution[j] - columnLower_[back[j]]) >
2123                                        fabs(solution[j] - columnUpper_[back[j]]))
2124                                   setDynamicStatus(back[j], atUpperBound);
2125                         } else if (columnUpper_ && solution[j] > 0.0) {
2126                              setDynamicStatus(back[j], atUpperBound);
2127                         } else {
2128                              setDynamicStatus(back[j], atLowerBound);
2129                              assert(!solution[j]);
2130                         }
2131                    }
2132               }
2133          }
2134          // convert iBasic back and do bounds
2135          if (iBasic == numberInSet) {
2136               // slack basic
2137               setStatus(iSet, ClpSimplex::basic);
2138               iBasic = iSet + maximumGubColumns_;
2139          } else {
2140               iBasic = back[iBasic];
2141               setDynamicStatus(iBasic, soloKey);
2142               // remember bounds flipped
2143               if (upper[numberInSet] == lower[numberInSet])
2144                    setStatus(iSet, ClpSimplex::isFixed);
2145               else if (solution[numberInSet] == upper[numberInSet])
2146                    setStatus(iSet, ClpSimplex::atLowerBound);
2147               else if (solution[numberInSet] == lower[numberInSet])
2148                    setStatus(iSet, ClpSimplex::atUpperBound);
2149               else
2150                    abort();
2151          }
2152          keyVariable_[iSet] = iBasic;
2153     }
2154     model_->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
2155     delete [] lower;
2156     delete [] solution;
2157     delete [] upper;
2158     delete [] cost;
2159     delete [] back;
2160     // make sure matrix is in good shape
2161     matrix_->orderMatrix();
2162}
2163// Populates initial matrix from dynamic status
2164void
2165ClpDynamicMatrix::initialProblem()
2166{
2167     int iSet;
2168     double * element =  matrix_->getMutableElements();
2169     int * row = matrix_->getMutableIndices();
2170     CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
2171     int * length = matrix_->getMutableVectorLengths();
2172     double * cost = model_->objective();
2173     double * solution = model_->primalColumnSolution();
2174     double * columnLower = model_->columnLower();
2175     double * columnUpper = model_->columnUpper();
2176     double * rowSolution = model_->primalRowSolution();
2177     double * rowLower = model_->rowLower();
2178     double * rowUpper = model_->rowUpper();
2179     CoinBigIndex numberElements = startColumn[firstDynamic_];
2180
2181     firstAvailable_ = firstDynamic_;
2182     numberActiveSets_ = 0;
2183     for (iSet = 0; iSet < numberSets_; iSet++) {
2184          toIndex_[iSet] = -1;
2185          int numberActive = 0;
2186          int whichKey = -1;
2187          if (getStatus(iSet) == ClpSimplex::basic) {
2188               whichKey = maximumGubColumns_ + iSet;
2189               numberActive = 1;
2190          } else {
2191               whichKey = -1;
2192          }
2193          int j = startSet_[iSet];
2194          while (j >= 0) {
2195               assert (getDynamicStatus(j) != soloKey || whichKey == -1);
2196               if (getDynamicStatus(j) == inSmall) {
2197                    numberActive++;
2198               } else if (getDynamicStatus(j) == soloKey) {
2199                    whichKey = j;
2200                    numberActive++;
2201               }
2202               j = next_[j]; //onto next in set
2203          }
2204          if (numberActive > 1) {
2205               int iRow = numberActiveSets_ + numberStaticRows_;
2206               rowSolution[iRow] = 0.0;
2207               double lowerValue;
2208               if (lowerSet_[iSet] > -1.0e20)
2209                    lowerValue = lowerSet_[iSet];
2210               else
2211                    lowerValue = -COIN_DBL_MAX;
2212               double upperValue;
2213               if (upperSet_[iSet] < 1.0e20)
2214                    upperValue = upperSet_[iSet];
2215               else
2216                    upperValue = COIN_DBL_MAX;
2217               rowLower[iRow] = lowerValue;
2218               rowUpper[iRow] = upperValue;
2219               if (getStatus(iSet) == ClpSimplex::basic) {
2220                    model_->setRowStatus(iRow, ClpSimplex::basic);
2221                    rowSolution[iRow] = 0.0;
2222               } else if (getStatus(iSet) == ClpSimplex::atLowerBound) {
2223                    model_->setRowStatus(iRow, ClpSimplex::atLowerBound);
2224                    rowSolution[iRow] = lowerValue;
2225               } else {
2226                    model_->setRowStatus(iRow, ClpSimplex::atUpperBound);
2227                    rowSolution[iRow] = upperValue;
2228               }
2229               j = startSet_[iSet];
2230               while (j >= 0) {
2231                    DynamicStatus status = getDynamicStatus(j);
2232                    if (status == inSmall) {
2233                         int numberThis = startColumn_[j+1] - startColumn_[j] + 1;
2234                         if (numberElements + numberThis > numberElements_) {
2235                              // need to redo
2236                              numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
2237                              matrix_->reserve(lastDynamic_, numberElements_);
2238                              element =  matrix_->getMutableElements();
2239                              row = matrix_->getMutableIndices();
2240                              // these probably okay but be safe
2241                              startColumn = matrix_->getMutableVectorStarts();
2242                              length = matrix_->getMutableVectorLengths();
2243                         }
2244                         length[firstAvailable_] = numberThis;
2245                         cost[firstAvailable_] = cost_[j];
2246                         CoinBigIndex base = startColumn_[j];
2247                         for (int k = 0; k < numberThis - 1; k++) {
2248                              row[numberElements] = row_[base+k];
2249                              element[numberElements++] = element_[base+k];
2250                         }
2251                         row[numberElements] = iRow;
2252                         element[numberElements++] = 1.0;
2253                         id_[firstAvailable_-firstDynamic_] = j;
2254                         solution[firstAvailable_] = 0.0;
2255                         model_->setStatus(firstAvailable_, ClpSimplex::basic);
2256                         if (!columnLower_ && !columnUpper_) {
2257                              columnLower[firstAvailable_] = 0.0;
2258                              columnUpper[firstAvailable_] = COIN_DBL_MAX;
2259                         }  else {
2260                              if (columnLower_)
2261                                   columnLower[firstAvailable_] = columnLower_[j];
2262                              else
2263                                   columnLower[firstAvailable_] = 0.0;
2264                              if (columnUpper_)
2265                                   columnUpper[firstAvailable_] = columnUpper_[j];
2266                              else
2267                                   columnUpper[firstAvailable_] = COIN_DBL_MAX;
2268                              if (status != atUpperBound) {
2269                                   solution[firstAvailable_] = columnLower[firstAvailable_];
2270                              } else {
2271                                   solution[firstAvailable_] = columnUpper[firstAvailable_];
2272                              }
2273                         }
2274                         firstAvailable_++;
2275                         startColumn[firstAvailable_] = numberElements;
2276                    }
2277                    j = next_[j]; //onto next in set
2278               }
2279               model_->setRowStatus(numberActiveSets_ + numberStaticRows_, getStatus(iSet));
2280               toIndex_[iSet] = numberActiveSets_;
2281               fromIndex_[numberActiveSets_++] = iSet;
2282          } else {
2283            // solo key
2284            bool needKey=false;
2285            if (numberActive) {
2286              if (whichKey<maximumGubColumns_) {
2287                // structural - assume ok
2288                needKey = false;
2289              } else {
2290                // slack
2291                keyVariable_[iSet] = maximumGubColumns_ + iSet;
2292                double value = keyValue(iSet);
2293                if (value<lowerSet_[iSet]-1.0e-8||
2294                    value>upperSet_[iSet]+1.0e-8)
2295                  needKey=true;
2296              }
2297            } else {
2298              needKey = true;
2299            }
2300            if (needKey) {
2301              // all to lb then up some (slack/null if possible)
2302              int length=99999999;
2303              int which=-1;
2304              double sum=0.0;
2305              for (int iColumn=startSet_[iSet];iColumn<startSet_[iSet+1];iColumn++) {
2306                setDynamicStatus(iColumn,atLowerBound);
2307                sum += columnLower_[iColumn];
2308                if (length>startColumn_[iColumn+1]-startColumn_[iColumn]) {
2309                  which=iColumn;
2310                  length=startColumn_[iColumn+1]-startColumn_[iColumn];
2311                }
2312              }
2313              if (sum>lowerSet_[iSet]-1.0e-8) {
2314                // slack can be basic
2315                setStatus(iSet,ClpSimplex::basic);
2316                keyVariable_[iSet] = maximumGubColumns_ + iSet;
2317              } else {
2318                // use shortest
2319                setDynamicStatus(which,soloKey);
2320                keyVariable_[iSet] = which;
2321                setStatus(iSet,ClpSimplex::atLowerBound);
2322              }
2323            }
2324          }
2325          assert (toIndex_[iSet] >= 0 || whichKey >= 0);
2326          keyVariable_[iSet] = whichKey;
2327     }
2328     // clean up pivotVariable
2329     int numberColumns = model_->numberColumns();
2330     int numberRows = model_->numberRows();
2331     int * pivotVariable = model_->pivotVariable();
2332     if (pivotVariable) {
2333       for (int i=0; i<numberStaticRows_+numberActiveSets_;i++) {
2334         if (model_->getRowStatus(i)!=ClpSimplex::basic)
2335           pivotVariable[i]=-1;
2336         else
2337           pivotVariable[i]=numberColumns+i;
2338       }
2339       for (int i=numberStaticRows_+numberActiveSets_;i<numberRows;i++) {
2340         pivotVariable[i]=i+numberColumns;
2341       }
2342       int put=-1;
2343       for (int i=0;i<numberColumns;i++) {
2344         if (model_->getColumnStatus(i)==ClpSimplex::basic) {
2345           while(put<numberRows) {
2346             put++;
2347             if (pivotVariable[put]==-1) {
2348               pivotVariable[put]=i;
2349               break;
2350             }
2351           }
2352         }
2353       }
2354       for (int i=CoinMax(put,0);i<numberRows;i++) {
2355         if (pivotVariable[i]==-1) 
2356           pivotVariable[i]=i+numberColumns;
2357       }
2358     }
2359     if (rhsOffset_) {
2360       double * cost = model_->costRegion();
2361       double * columnLower = model_->lowerRegion();
2362       double * columnUpper = model_->upperRegion();
2363       double * solution = model_->solutionRegion();
2364       int numberRows = model_->numberRows();
2365       for (int i = numberActiveSets_; i < numberRows-numberStaticRows_; i++) {
2366         int iSequence = i + numberStaticRows_ + numberColumns;
2367         solution[iSequence] = 0.0;
2368         columnLower[iSequence] = -COIN_DBL_MAX;
2369         columnUpper[iSequence] = COIN_DBL_MAX;
2370         cost[iSequence] = 0.0;
2371         model_->nonLinearCost()->setOne(iSequence, solution[iSequence],
2372                                        columnLower[iSequence],
2373                                        columnUpper[iSequence], 0.0);
2374         model_->setStatus(iSequence, ClpSimplex::basic);
2375         rhsOffset_[i+numberStaticRows_] = 0.0;
2376       }
2377#if 0
2378       for (int i=0;i<numberStaticRows_;i++)
2379         printf("%d offset %g\n",
2380                i,rhsOffset_[i]);
2381#endif
2382     }
2383     numberActiveColumns_ = firstAvailable_;
2384#if 0
2385     for (iSet = 0; iSet < numberSets_; iSet++) {
2386       for (int j=startSet_[iSet];j<startSet_[iSet+1];j++) {
2387         if (getDynamicStatus(j)==soloKey)
2388           printf("%d in set %d is solo key\n",j,iSet);
2389         else if (getDynamicStatus(j)==inSmall)
2390           printf("%d in set %d is in small\n",j,iSet);
2391       }
2392     }
2393#endif
2394     return;
2395}
2396// Writes out model (without names)
2397void 
2398ClpDynamicMatrix::writeMps(const char * name)
2399{
2400  int numberTotalRows = numberStaticRows_+numberSets_;
2401  int numberTotalColumns = firstDynamic_+numberGubColumns_;
2402  // over estimate
2403  int numberElements = getNumElements()+startColumn_[numberGubColumns_]
2404    + numberGubColumns_;
2405  double * columnLower = new double [numberTotalColumns];
2406  double * columnUpper = new double [numberTotalColumns];
2407  double * cost = new double [numberTotalColumns];
2408  double * rowLower = new double [numberTotalRows];
2409  double * rowUpper = new double [numberTotalRows];
2410  CoinBigIndex * start = new CoinBigIndex[numberTotalColumns+1];
2411  int * row = new int [numberElements];
2412  double * element = new double [numberElements];
2413  // Fill in
2414  const CoinBigIndex * startA = getVectorStarts();
2415  const int * lengthA = getVectorLengths();
2416  const int * rowA = getIndices();
2417  const double * elementA = getElements();
2418  const double * columnLowerA = model_->columnLower();
2419  const double * columnUpperA = model_->columnUpper();
2420  const double * costA = model_->objective();
2421  const double * rowLowerA = model_->rowLower();
2422  const double * rowUpperA = model_->rowUpper();
2423  start[0]=0;
2424  numberElements=0;
2425  for (int i=0;i<firstDynamic_;i++) {
2426    columnLower[i] = columnLowerA[i];
2427    columnUpper[i] = columnUpperA[i];
2428    cost[i] = costA[i];
2429    for (CoinBigIndex j = startA[i];j<startA[i]+lengthA[i];j++) {
2430      row[numberElements] = rowA[j];
2431      element[numberElements++]=elementA[j];
2432    }
2433    start[i+1]=numberElements;
2434  }
2435  for (int i=0;i<numberStaticRows_;i++) {
2436    rowLower[i] = rowLowerA[i];
2437    rowUpper[i] = rowUpperA[i];
2438  }
2439  int putC=firstDynamic_;
2440  int putR=numberStaticRows_;
2441  for (int i=0;i<numberSets_;i++) {
2442    rowLower[putR]=lowerSet_[i];
2443    rowUpper[putR]=upperSet_[i];
2444    for (CoinBigIndex k=startSet_[i];k<startSet_[i+1];k++) {
2445      columnLower[putC]=columnLower_[k];
2446      columnUpper[putC]=columnUpper_[k];
2447      cost[putC]=cost_[k];
2448      putC++;
2449      for (CoinBigIndex j = startColumn_[k];j<startColumn_[k+1];j++) {
2450        row[numberElements] = row_[j];
2451        element[numberElements++]=element_[j];
2452      }
2453      row[numberElements] = putR;
2454      element[numberElements++]=1.0;
2455      start[putC]=numberElements;
2456    }
2457    putR++;
2458  }
2459
2460  assert (putR==numberTotalRows);
2461  assert (putC==numberTotalColumns);
2462  ClpSimplex modelOut;
2463  modelOut.loadProblem(numberTotalColumns,numberTotalRows,
2464                       start,row,element,
2465                       columnLower,columnUpper,cost,
2466                       rowLower,rowUpper);
2467  modelOut.writeMps(name);
2468  delete [] columnLower;
2469  delete [] columnUpper;
2470  delete [] cost;
2471  delete [] rowLower;
2472  delete [] rowUpper;
2473  delete [] start;
2474  delete [] row; 
2475  delete [] element;
2476}
2477// Adds in a column to gub structure (called from descendant)
2478int
2479ClpDynamicMatrix::addColumn(int numberEntries, const int * row, const double * element,
2480                            double cost, double lower, double upper, int iSet,
2481                            DynamicStatus status)
2482{
2483     // check if already in
2484     int j = startSet_[iSet];
2485     while (j >= 0) {
2486          if (startColumn_[j+1] - startColumn_[j] == numberEntries) {
2487               const int * row2 = row_ + startColumn_[j];
2488               const double * element2 = element_ + startColumn_[j];
2489               bool same = true;
2490               for (int k = 0; k < numberEntries; k++) {
2491                    if (row[k] != row2[k] || element[k] != element2[k]) {
2492                         same = false;
2493                         break;
2494                    }
2495               }
2496               if (same) {
2497                    bool odd = false;
2498                    if (cost != cost_[j])
2499                         odd = true;
2500                    if (columnLower_ && lower != columnLower_[j])
2501                         odd = true;
2502                    if (columnUpper_ && upper != columnUpper_[j])
2503                         odd = true;
2504                    if (odd) {
2505                         printf("seems odd - same els but cost,lo,up are %g,%g,%g and %g,%g,%g\n",
2506                                cost, lower, upper, cost_[j],
2507                                columnLower_ ? columnLower_[j] : 0.0,
2508                                columnUpper_ ? columnUpper_[j] : 1.0e100);
2509                    } else {
2510                         setDynamicStatus(j, status);
2511                         return j;
2512                    }
2513               }
2514          }
2515          j = next_[j];
2516     }
2517
2518     if (numberGubColumns_ == maximumGubColumns_ ||
2519               startColumn_[numberGubColumns_] + numberEntries > maximumElements_) {
2520          CoinBigIndex j;
2521          int i;
2522          int put = 0;
2523          int numberElements = 0;
2524          CoinBigIndex start = 0;
2525          // compress - leave ones at ub and basic
2526          int * which = new int [numberGubColumns_];
2527          for (i = 0; i < numberGubColumns_; i++) {
2528               CoinBigIndex end = startColumn_[i+1];
2529               // what about ubs if column generation?
2530               if (getDynamicStatus(i) != atLowerBound) {
2531                    // keep in
2532                    for (j = start; j < end; j++) {
2533                         row_[numberElements] = row_[j];
2534                         element_[numberElements++] = element_[j];
2535                    }
2536                    startColumn_[put+1] = numberElements;
2537                    cost_[put] = cost_[i];
2538                    if (columnLower_)
2539                         columnLower_[put] = columnLower_[i];
2540                    if (columnUpper_)
2541                         columnUpper_[put] = columnUpper_[i];
2542                    dynamicStatus_[put] = dynamicStatus_[i];
2543                    id_[put] = id_[i];
2544                    which[i] = put;
2545                    put++;
2546               } else {
2547                    // out
2548                    which[i] = -1;
2549               }
2550               start = end;
2551          }
2552          // now redo startSet_ and next_
2553          int * newNext = new int [maximumGubColumns_];
2554          for (int jSet = 0; jSet < numberSets_; jSet++) {
2555               int sequence = startSet_[jSet];
2556               while (which[sequence] < 0) {
2557                    // out
2558                    assert (next_[sequence] >= 0);
2559                    sequence = next_[sequence];
2560               }
2561               startSet_[jSet] = which[sequence];
2562               int last = which[sequence];
2563               while (next_[sequence] >= 0) {
2564                    sequence = next_[sequence];
2565                    if(which[sequence] >= 0) {
2566                         // keep
2567                         newNext[last] = which[sequence];
2568                         last = which[sequence];
2569                    }
2570               }
2571               newNext[last] = -jSet - 1;
2572          }
2573          delete [] next_;
2574          next_ = newNext;
2575          delete [] which;
2576          abort();
2577     }
2578     CoinBigIndex start = startColumn_[numberGubColumns_];
2579     CoinMemcpyN(row, numberEntries, row_ + start);
2580     CoinMemcpyN(element, numberEntries, element_ + start);
2581     startColumn_[numberGubColumns_+1] = start + numberEntries;
2582     cost_[numberGubColumns_] = cost;
2583     if (columnLower_)
2584          columnLower_[numberGubColumns_] = lower;
2585     else
2586          assert (!lower);
2587     if (columnUpper_)
2588          columnUpper_[numberGubColumns_] = upper;
2589     else
2590          assert (upper > 1.0e20);
2591     setDynamicStatus(numberGubColumns_, status);
2592     // Do next_
2593     j = startSet_[iSet];
2594     startSet_[iSet] = numberGubColumns_;
2595     next_[numberGubColumns_] = j;
2596     numberGubColumns_++;
2597     return numberGubColumns_ - 1;
2598}
2599// Returns which set a variable is in
2600int
2601ClpDynamicMatrix::whichSet (int sequence) const
2602{
2603     while (next_[sequence] >= 0)
2604          sequence = next_[sequence];
2605     int iSet = - next_[sequence] - 1;
2606     return iSet;
2607}
Note: See TracBrowser for help on using the repository browser.