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

Last change on this file since 1677 was 1677, checked in by forrest, 9 years ago

gub

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