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

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

For GUB

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 108.8 KB
Line 
1/* $Id: ClpDynamicMatrix.cpp 1676 2011-01-20 17:00:28Z 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_);
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_, CoinMax(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_);
262          assert (dynamicStatus);
263          dynamicStatus_ = ClpCopyOfArray(dynamicStatus, 2*numberGubColumns_);
264     } else {
265          status_ = new unsigned char [2*numberSets_];
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_);
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 (getStatus(iSet) != model->getStatus(sequenceOut))
868               printf("** set %d status %d, var status %d\n", iSet,
869                      getStatus(iSet), model->getStatus(sequenceOut));
870     }
871     ClpMatrixBase::updatePivot(model, oldInValue, oldOutValue);
872#ifdef CLP_DEBUG
873     char * inSmall = new char [numberGubColumns_];
874     memset(inSmall, 0, numberGubColumns_);
875     const double * solution = model->solutionRegion();
876     for (int i = 0; i < numberGubColumns_; i++)
877          if (getDynamicStatus(i) == ClpDynamicMatrix::inSmall)
878               inSmall[i] = 1;
879     for (int i = firstDynamic_; i < firstAvailable_; i++) {
880          int k = id_[i-firstDynamic_];
881          inSmall[k] = 0;
882          //if (k>=23289&&k<23357&&solution[i])
883          //printf("var %d (in small %d) has value %g\n",k,i,solution[i]);
884     }
885     for (int i = 0; i < numberGubColumns_; i++)
886          assert (!inSmall[i]);
887     delete [] inSmall;
888#ifndef NDEBUG
889     for (int i = 0; i < numberActiveSets_; i++) {
890          int iSet = fromIndex_[i];
891          assert (toIndex_[iSet] == i);
892          //if (getStatus(iSet)!=model->getRowStatus(i+numberStaticRows_))
893          //printf("*** set %d status %d, var status %d\n",iSet,
894          //     getStatus(iSet),model->getRowStatus(i+numberStaticRows_));
895          //assert (model->getRowStatus(i+numberStaticRows_)==getStatus(iSet));
896          //if (iSet==1035) {
897          //printf("rhs for set %d (%d) is %g %g - cost %g\n",iSet,i,model->lowerRegion(0)[i+numberStaticRows_],
898          //     model->upperRegion(0)[i+numberStaticRows_],model->costRegion(0)[i+numberStaticRows_]);
899          //}
900     }
901#endif
902#endif
903     return 0;
904}
905/*
906     utility dual function for dealing with dynamic constraints
907     mode=n see ClpGubMatrix.hpp for definition
908     Remember to update here when settled down
909*/
910void
911ClpDynamicMatrix::dualExpanded(ClpSimplex * model,
912                               CoinIndexedVector * /*array*/,
913                               double * /*other*/, int mode)
914{
915     switch (mode) {
916          // modify costs before transposeUpdate
917     case 0:
918          break;
919          // create duals for key variables (without check on dual infeasible)
920     case 1:
921          break;
922          // as 1 but check slacks and compute djs
923     case 2: {
924          // do pivots
925          int * pivotVariable = model->pivotVariable();
926          int numberRows = numberStaticRows_ + numberActiveSets_;
927          int numberColumns = model->numberColumns();
928          for (int iRow = 0; iRow < numberRows; iRow++) {
929               int iPivot = pivotVariable[iRow];
930               if (iPivot < numberColumns)
931                    backToPivotRow_[iPivot] = iRow;
932          }
933          if (noCheck_ >= 0) {
934               if (infeasibilityWeight_ != model_->infeasibilityCost()) {
935                    // don't bother checking
936                    sumDualInfeasibilities_ = 100.0;
937                    numberDualInfeasibilities_ = 1;
938                    sumOfRelaxedDualInfeasibilities_ = 100.0;
939                    return;
940               }
941          }
942          // In theory we should subtract out ones we have done but ....
943          // If key slack then dual 0.0
944          // If not then slack could be dual infeasible
945          // dj for key is zero so that defines dual on set
946          int i;
947          double * dual = model->dualRowSolution();
948          double dualTolerance = model->dualTolerance();
949          double relaxedTolerance = dualTolerance;
950          // we can't really trust infeasibilities if there is dual error
951          double error = CoinMin(1.0e-2, model->largestDualError());
952          // allow tolerance at least slightly bigger than standard
953          relaxedTolerance = relaxedTolerance +  error;
954          // but we will be using difference
955          relaxedTolerance -= dualTolerance;
956          sumDualInfeasibilities_ = 0.0;
957          numberDualInfeasibilities_ = 0;
958          sumOfRelaxedDualInfeasibilities_ = 0.0;
959          for (i = 0; i < numberSets_; i++) {
960               double value = 0.0;
961               int gubRow = toIndex_[i];
962               if (gubRow < 0) {
963                    int kColumn = keyVariable_[i];
964                    if (kColumn < maximumGubColumns_) {
965                         // dj without set
966                         value = cost_[kColumn];
967                         for (CoinBigIndex j = startColumn_[kColumn];
968                                   j < startColumn_[kColumn+1]; j++) {
969                              int iRow = row_[j];
970                              value -= dual[iRow] * element_[j];
971                         }
972                         double infeasibility = 0.0;
973                         if (getStatus(i) == ClpSimplex::atLowerBound) {
974                              if (-value > dualTolerance)
975                                   infeasibility = -value - dualTolerance;
976                         } else if (getStatus(i) == ClpSimplex::atUpperBound) {
977                              if (value > dualTolerance)
978                                   infeasibility = value - dualTolerance;
979                         }
980                         if (infeasibility > 0.0) {
981                              sumDualInfeasibilities_ += infeasibility;
982                              if (infeasibility > relaxedTolerance)
983                                   sumOfRelaxedDualInfeasibilities_ += infeasibility;
984                              numberDualInfeasibilities_ ++;
985                         }
986                    }
987               } else {
988                    value = dual[gubRow+numberStaticRows_];
989               }
990               // Now subtract out from all
991               int k = startSet_[i];
992               while (k >= 0) {
993                    if (getDynamicStatus(k) != inSmall) {
994                         double djValue = cost_[k] - value;
995                         for (CoinBigIndex j = startColumn_[k];
996                                   j < startColumn_[k+1]; j++) {
997                              int iRow = row_[j];
998                              djValue -= dual[iRow] * element_[j];
999                         }
1000                         double infeasibility = 0.0;
1001                         if (getDynamicStatus(k) == atLowerBound) {
1002                              if (djValue < -dualTolerance)
1003                                   infeasibility = -djValue - dualTolerance;
1004                         } else if (getDynamicStatus(k) == atUpperBound) {
1005                              // at upper bound
1006                              if (djValue > dualTolerance)
1007                                   infeasibility = djValue - dualTolerance;
1008                         }
1009                         if (infeasibility > 0.0) {
1010                              sumDualInfeasibilities_ += infeasibility;
1011                              if (infeasibility > relaxedTolerance)
1012                                   sumOfRelaxedDualInfeasibilities_ += infeasibility;
1013                              numberDualInfeasibilities_ ++;
1014                         }
1015                    }
1016                    k = next_[k]; //onto next in set
1017               }
1018          }
1019     }
1020     infeasibilityWeight_ = -1.0;
1021     break;
1022     // Report on infeasibilities of key variables
1023     case 3: {
1024          model->setSumDualInfeasibilities(model->sumDualInfeasibilities() +
1025                                           sumDualInfeasibilities_);
1026          model->setNumberDualInfeasibilities(model->numberDualInfeasibilities() +
1027                                              numberDualInfeasibilities_);
1028          model->setSumOfRelaxedDualInfeasibilities(model->sumOfRelaxedDualInfeasibilities() +
1029                    sumOfRelaxedDualInfeasibilities_);
1030     }
1031     break;
1032     // modify costs before transposeUpdate for partial pricing
1033     case 4:
1034          break;
1035     }
1036}
1037/*
1038     general utility function for dealing with dynamic constraints
1039     mode=n see ClpGubMatrix.hpp for definition
1040     Remember to update here when settled down
1041*/
1042int
1043ClpDynamicMatrix::generalExpanded(ClpSimplex * model, int mode, int &number)
1044{
1045     int returnCode = 0;
1046     switch (mode) {
1047          // Fill in pivotVariable
1048     case 0: {
1049          // If no effective rhs - form it
1050          if (!rhsOffset_) {
1051               rhsOffset_ = new double[model->numberRows()];
1052               rhsOffset(model, true);
1053          }
1054          int i;
1055          int numberBasic = number;
1056          int numberColumns = model->numberColumns();
1057          // Use different array so can build from true pivotVariable_
1058          //int * pivotVariable = model->pivotVariable();
1059          int * pivotVariable = model->rowArray(0)->getIndices();
1060          for (i = 0; i < numberColumns; i++) {
1061               if (model->getColumnStatus(i) == ClpSimplex::basic)
1062                    pivotVariable[numberBasic++] = i;
1063          }
1064          number = numberBasic;
1065     }
1066     break;
1067     // Do initial extra rows + maximum basic
1068     case 2: {
1069          number = model->numberRows();
1070     }
1071     break;
1072     // Before normal replaceColumn
1073     case 3: {
1074          if (numberActiveSets_ + numberStaticRows_ == model_->numberRows()) {
1075               // no space - re-factorize
1076               returnCode = 4;
1077               number = -1; // say no need for normal replaceColumn
1078          }
1079     }
1080     break;
1081     // To see if can dual or primal
1082     case 4: {
1083          returnCode = 1;
1084     }
1085     break;
1086     // save status
1087     case 5: {
1088       memcpy(status_+numberSets_,status_,numberSets_);
1089       memcpy(dynamicStatus_+maximumGubColumns_,
1090              dynamicStatus_,maximumGubColumns_);
1091     }
1092     break;
1093     // restore status
1094     case 6: {
1095       memcpy(status_,status_+numberSets_,numberSets_);
1096       memcpy(dynamicStatus_,dynamicStatus_+maximumGubColumns_,
1097              maximumGubColumns_);
1098       initialProblem();
1099     }
1100     break;
1101     // unflag all variables
1102     case 8: {
1103          for (int i = 0; i < numberGubColumns_; i++) {
1104               if (flagged(i)) {
1105                    unsetFlagged(i);
1106                    returnCode++;
1107               }
1108          }
1109     }
1110     break;
1111     // redo costs in primal
1112     case 9: {
1113          double * cost = model->costRegion();
1114          double * solution = model->solutionRegion();
1115          double * columnLower = model->lowerRegion();
1116          double * columnUpper = model->upperRegion();
1117          int i;
1118          bool doCosts = (number & 4) != 0;
1119          bool doBounds = (number & 1) != 0;
1120          for ( i = firstDynamic_; i < firstAvailable_; i++) {
1121               int jColumn = id_[i-firstDynamic_];
1122               if (doBounds) {
1123                    if (!columnLower_ && !columnUpper_) {
1124                         columnLower[i] = 0.0;
1125                         columnUpper[i] = COIN_DBL_MAX;
1126                    }  else {
1127                         if (columnLower_)
1128                              columnLower[i] = columnLower_[jColumn];
1129                         else
1130                              columnLower[i] = 0.0;
1131                         if (columnUpper_)
1132                              columnUpper[i] = columnUpper_[jColumn];
1133                         else
1134                              columnUpper[i] = COIN_DBL_MAX;
1135                    }
1136               }
1137               if (doCosts) {
1138                    cost[i] = cost_[jColumn];
1139                    // Original bounds
1140                    if (model->nonLinearCost())
1141                         model->nonLinearCost()->setOne(i, solution[i],
1142                                                        this->columnLower(jColumn),
1143                                                        this->columnUpper(jColumn), cost_[jColumn]);
1144               }
1145          }
1146          // and active sets
1147          for (i = 0; i < numberActiveSets_; i++ ) {
1148               int iSet = fromIndex_[i];
1149               int iSequence = lastDynamic_ + numberStaticRows_ + i;
1150               if (doBounds) {
1151                    if (lowerSet_[iSet] > -1.0e20)
1152                         columnLower[iSequence] = lowerSet_[iSet];
1153                    else
1154                         columnLower[iSequence] = -COIN_DBL_MAX;
1155                    if (upperSet_[iSet] < 1.0e20)
1156                         columnUpper[iSequence] = upperSet_[iSet];
1157                    else
1158                         columnUpper[iSequence] = COIN_DBL_MAX;
1159               }
1160               if (doCosts) {
1161                    if (model->nonLinearCost()) {
1162                         double trueLower;
1163                         if (lowerSet_[iSet] > -1.0e20)
1164                              trueLower = lowerSet_[iSet];
1165                         else
1166                              trueLower = -COIN_DBL_MAX;
1167                         double trueUpper;
1168                         if (upperSet_[iSet] < 1.0e20)
1169                              trueUpper = upperSet_[iSet];
1170                         else
1171                              trueUpper = COIN_DBL_MAX;
1172                         model->nonLinearCost()->setOne(iSequence, solution[iSequence],
1173                                                        trueLower, trueUpper, 0.0);
1174                    }
1175               }
1176          }
1177     }
1178     break;
1179     // return 1 if there may be changing bounds on variable (column generation)
1180     case 10: {
1181          // return 1 as bounds on rhs will change
1182          returnCode = 1;
1183     }
1184     break;
1185     // make sure set is clean
1186     case 7: {
1187          // first flag
1188          if (number >= firstDynamic_ && number < lastDynamic_) {
1189            int sequence = id_[number-firstDynamic_];
1190            setFlagged(sequence);
1191            //model->clearFlagged(number);
1192          } else if (number>=model_->numberColumns()+numberStaticRows_) {
1193            // slack
1194            int iSet = fromIndex_[number-model_->numberColumns()-
1195                                  numberStaticRows_];
1196            setFlaggedSlack(iSet);
1197            //model->clearFlagged(number);
1198          }
1199     }
1200     case 11: {
1201          //int sequenceIn = model->sequenceIn();
1202          if (number >= firstDynamic_ && number < lastDynamic_) {
1203            //assert (number == model->sequenceIn());
1204               // take out variable (but leave key)
1205               double * cost = model->costRegion();
1206               double * columnLower = model->lowerRegion();
1207               double * columnUpper = model->upperRegion();
1208               double * solution = model->solutionRegion();
1209               int * length = matrix_->getMutableVectorLengths();
1210#ifndef NDEBUG
1211               if (length[number]) {
1212                    int * row = matrix_->getMutableIndices();
1213                    CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
1214                    int iRow = row[startColumn[number] + length[number] - 1];
1215                    int which = iRow - numberStaticRows_;
1216                    assert (which >= 0);
1217                    int iSet = fromIndex_[which];
1218                    assert (toIndex_[iSet] == which);
1219               }
1220#endif
1221               // no need firstAvailable_--;
1222               solution[firstAvailable_] = 0.0;
1223               cost[firstAvailable_] = 0.0;
1224               length[firstAvailable_] = 0;
1225               model->nonLinearCost()->setOne(firstAvailable_, 0.0, 0.0, COIN_DBL_MAX, 0.0);
1226               model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
1227               columnLower[firstAvailable_] = 0.0;
1228               columnUpper[firstAvailable_] = COIN_DBL_MAX;
1229
1230               // not really in small problem
1231               int iBig = id_[number-firstDynamic_];
1232               if (model->getStatus(number) == ClpSimplex::atLowerBound) {
1233                    setDynamicStatus(iBig, atLowerBound);
1234                    if (columnLower_)
1235                         modifyOffset(number, columnLower_[iBig]);
1236               } else {
1237                    setDynamicStatus(iBig, atUpperBound);
1238                    modifyOffset(number, columnUpper_[iBig]);
1239               }
1240          } else if (number>=model_->numberColumns()+numberStaticRows_) {
1241            // slack
1242            int iSet = fromIndex_[number-model_->numberColumns()-
1243                                  numberStaticRows_];
1244            printf("what now - set %d\n",iSet);
1245          }
1246     }
1247     break;
1248     default:
1249          break;
1250     }
1251     return returnCode;
1252}
1253/* Purely for column generation and similar ideas.  Allows
1254   matrix and any bounds or costs to be updated (sensibly).
1255   Returns non-zero if any changes.
1256*/
1257int
1258ClpDynamicMatrix::refresh(ClpSimplex * model)
1259{
1260     // If at beginning then create initial problem
1261     if (firstDynamic_ == firstAvailable_) {
1262          initialProblem();
1263          return 1;
1264     } else if (!model->nonLinearCost()) {
1265          // will be same as last time
1266          return 1;
1267     }
1268     // lookup array
1269     int * active = new int [numberActiveSets_];
1270     CoinZeroN(active, numberActiveSets_);
1271     int iColumn;
1272     int numberColumns = model->numberColumns();
1273     double * element =  matrix_->getMutableElements();
1274     int * row = matrix_->getMutableIndices();
1275     CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
1276     int * length = matrix_->getMutableVectorLengths();
1277     double * cost = model->costRegion();
1278     double * columnLower = model->lowerRegion();
1279     double * columnUpper = model->upperRegion();
1280     CoinBigIndex numberElements = startColumn[firstDynamic_];
1281     // first just do lookup and basic stuff
1282     int currentNumber = firstAvailable_;
1283     firstAvailable_ = firstDynamic_;
1284     double objectiveChange = 0.0;
1285     double * solution = model->solutionRegion();
1286     int currentNumberActiveSets = numberActiveSets_;
1287     for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
1288          int iRow = row[startColumn[iColumn] + length[iColumn] - 1];
1289          int which = iRow - numberStaticRows_;
1290          assert (which >= 0);
1291          int iSet = fromIndex_[which];
1292          assert (toIndex_[iSet] == which);
1293          if (model->getStatus(iColumn) == ClpSimplex::basic) {
1294               active[which]++;
1295               // may as well make key
1296               keyVariable_[iSet] = id_[iColumn-firstDynamic_];
1297          }
1298     }
1299     int i;
1300     numberActiveSets_ = 0;
1301     int numberDeleted = 0;
1302     for (i = 0; i < currentNumberActiveSets; i++) {
1303          int iRow = i + numberStaticRows_;
1304          int numberActive = active[i];
1305          int iSet = fromIndex_[i];
1306          if (model->getRowStatus(iRow) == ClpSimplex::basic) {
1307               numberActive++;
1308               // may as well make key
1309               keyVariable_[iSet] = maximumGubColumns_ + iSet;
1310          }
1311          if (numberActive > 1) {
1312               // keep
1313               active[i] = numberActiveSets_;
1314               numberActiveSets_++;
1315          } else {
1316               active[i] = -1;
1317          }
1318     }
1319
1320     for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
1321          int iRow = row[startColumn[iColumn] + length[iColumn] - 1];
1322          int which = iRow - numberStaticRows_;
1323          int jColumn = id_[iColumn-firstDynamic_];
1324          if (model->getStatus(iColumn) == ClpSimplex::atLowerBound ||
1325                    model->getStatus(iColumn) == ClpSimplex::atUpperBound) {
1326               double value = solution[iColumn];
1327               bool toLowerBound = true;
1328               assert (jColumn>=0);assert (iColumn>=0);
1329               if (columnUpper_) {
1330                    if (!columnLower_) {
1331                         if (fabs(value - columnUpper_[jColumn]) < fabs(value))
1332                              toLowerBound = false;
1333                    } else if (fabs(value - columnUpper_[jColumn]) < fabs(value - columnLower_[jColumn])) {
1334                         toLowerBound = false;
1335                    }
1336               }
1337               if (toLowerBound) {
1338                    // throw out to lower bound
1339                    if (columnLower_) {
1340                         setDynamicStatus(jColumn, atLowerBound);
1341                         // treat solution as if exactly at a bound
1342                         double value = columnLower[iColumn];
1343                         objectiveChange += cost[iColumn] * value;
1344                    } else {
1345                         setDynamicStatus(jColumn, atLowerBound);
1346                    }
1347               } else {
1348                    // throw out to upper bound
1349                    setDynamicStatus(jColumn, atUpperBound);
1350                    double value = columnUpper[iColumn];
1351                    objectiveChange += cost[iColumn] * value;
1352               }
1353               numberDeleted++;
1354          } else {
1355               assert(model->getStatus(iColumn) != ClpSimplex::superBasic); // deal with later
1356               int iPut = active[which];
1357               if (iPut >= 0) {
1358                    // move
1359                    id_[firstAvailable_-firstDynamic_] = jColumn;
1360                    int numberThis = startColumn_[jColumn+1] - startColumn_[jColumn] + 1;
1361                    length[firstAvailable_] = numberThis;
1362                    cost[firstAvailable_] = cost[iColumn];
1363                    columnLower[firstAvailable_] = columnLower[iColumn];
1364                    columnUpper[firstAvailable_] = columnUpper[iColumn];
1365                    model->nonLinearCost()->setOne(firstAvailable_, solution[iColumn], 0.0, COIN_DBL_MAX,
1366                                                   cost_[jColumn]);
1367                    CoinBigIndex base = startColumn_[jColumn];
1368                    for (int j = 0; j < numberThis - 1; j++) {
1369                         row[numberElements] = row_[base+j];
1370                         element[numberElements++] = element_[base+j];
1371                    }
1372                    row[numberElements] = iPut + numberStaticRows_;
1373                    element[numberElements++] = 1.0;
1374                    model->setStatus(firstAvailable_, model->getStatus(iColumn));
1375                    solution[firstAvailable_] = solution[iColumn];
1376                    firstAvailable_++;
1377                    startColumn[firstAvailable_] = numberElements;
1378               }
1379          }
1380     }
1381     // zero solution
1382     CoinZeroN(solution + firstAvailable_, currentNumber - firstAvailable_);
1383     // zero cost
1384     CoinZeroN(cost + firstAvailable_, currentNumber - firstAvailable_);
1385     // zero lengths
1386     CoinZeroN(length + firstAvailable_, currentNumber - firstAvailable_);
1387     for ( iColumn = firstAvailable_; iColumn < currentNumber; iColumn++) {
1388          model->nonLinearCost()->setOne(iColumn, 0.0, 0.0, COIN_DBL_MAX, 0.0);
1389          model->setStatus(iColumn, ClpSimplex::atLowerBound);
1390          columnLower[iColumn] = 0.0;
1391          columnUpper[iColumn] = COIN_DBL_MAX;
1392     }
1393     // move rhs and set rest to infinite
1394     numberActiveSets_ = 0;
1395     for (i = 0; i < currentNumberActiveSets; i++) {
1396          int iSet = fromIndex_[i];
1397          assert (toIndex_[iSet] == i);
1398          int iRow = i + numberStaticRows_;
1399          int iPut = active[i];
1400          if (iPut >= 0) {
1401               int in = iRow + numberColumns;
1402               int out = iPut + numberColumns + numberStaticRows_;
1403               solution[out] = solution[in];
1404               columnLower[out] = columnLower[in];
1405               columnUpper[out] = columnUpper[in];
1406               cost[out] = cost[in];
1407               double trueLower;
1408               if (lowerSet_[iSet] > -1.0e20)
1409                    trueLower = lowerSet_[iSet];
1410               else
1411                    trueLower = -COIN_DBL_MAX;
1412               double trueUpper;
1413               if (upperSet_[iSet] < 1.0e20)
1414                    trueUpper = upperSet_[iSet];
1415               else
1416                    trueUpper = COIN_DBL_MAX;
1417               model->nonLinearCost()->setOne(out, solution[out],
1418                                              trueLower, trueUpper, 0.0);
1419               model->setStatus(out, model->getStatus(in));
1420               toIndex_[iSet] = numberActiveSets_;
1421               rhsOffset_[numberActiveSets_+numberStaticRows_] = rhsOffset_[i+numberStaticRows_];
1422               fromIndex_[numberActiveSets_++] = fromIndex_[i];
1423          } else {
1424               // adjust offset
1425               // put out as key
1426               int jColumn = keyVariable_[iSet];
1427               toIndex_[iSet] = -1;
1428               if (jColumn < maximumGubColumns_) {
1429                    setDynamicStatus(jColumn, soloKey);
1430                    double value = keyValue(iSet);
1431                    objectiveChange += cost_[jColumn] * value;
1432                    modifyOffset(jColumn, -value);
1433               }
1434          }
1435     }
1436     model->setObjectiveOffset(model->objectiveOffset() - objectiveChange);
1437     delete [] active;
1438     for (i = numberActiveSets_; i < currentNumberActiveSets; i++) {
1439          int iSequence = i + numberStaticRows_ + numberColumns;
1440          solution[iSequence] = 0.0;
1441          columnLower[iSequence] = -COIN_DBL_MAX;
1442          columnUpper[iSequence] = COIN_DBL_MAX;
1443          cost[iSequence] = 0.0;
1444          model->nonLinearCost()->setOne(iSequence, solution[iSequence],
1445                                         columnLower[iSequence],
1446                                         columnUpper[iSequence], 0.0);
1447          model->setStatus(iSequence, ClpSimplex::basic);
1448          rhsOffset_[i+numberStaticRows_] = 0.0;
1449     }
1450     if (currentNumberActiveSets != numberActiveSets_ || numberDeleted) {
1451          // deal with pivotVariable
1452          int * pivotVariable = model->pivotVariable();
1453          int sequence = firstDynamic_;
1454          int iRow = 0;
1455          int base1 = firstDynamic_;
1456          int base2 = lastDynamic_;
1457          int base3 = numberColumns + numberStaticRows_;
1458          int numberRows = numberStaticRows_ + currentNumberActiveSets;
1459          while (sequence < firstAvailable_) {
1460               int iPivot = pivotVariable[iRow];
1461               while (iPivot < base1 || (iPivot >= base2 && iPivot < base3)) {
1462                    iPivot = pivotVariable[++iRow];
1463               }
1464               pivotVariable[iRow++] = sequence;
1465               sequence++;
1466          }
1467          // move normal basic ones down
1468          int iPut = iRow;
1469          for (; iRow < numberRows; iRow++) {
1470               int iPivot = pivotVariable[iRow];
1471               if (iPivot < base1 || (iPivot >= base2 && iPivot < base3)) {
1472                    pivotVariable[iPut++] = iPivot;
1473               }
1474          }
1475          // and basic others
1476          for (i = 0; i < numberActiveSets_; i++) {
1477               if (model->getRowStatus(i + numberStaticRows_) == ClpSimplex::basic) {
1478                    pivotVariable[iPut++] = i + base3;
1479               }
1480          }
1481          for (i = numberActiveSets_; i < currentNumberActiveSets; i++) {
1482               pivotVariable[iPut++] = i + base3;
1483          }
1484          assert (iPut == numberRows);
1485     }
1486#ifdef CLP_DEBUG
1487#if 0
1488     printf("row for set 244 is %d, row status %d value %g ", toIndex_[244], status_[244],
1489            keyValue(244));
1490     int jj = startSet_[244];
1491     while (jj >= 0) {
1492          printf("var %d status %d ", jj, dynamicStatus_[jj]);
1493          jj = next_[jj];
1494     }
1495     printf("\n");
1496#endif
1497#ifdef CLP_DEBUG
1498     {
1499          // debug coding to analyze set
1500          int i;
1501          int kSet = -6;
1502          if (kSet >= 0) {
1503               int * back = new int [numberGubColumns_];
1504               for (i = 0; i < numberGubColumns_; i++)
1505                    back[i] = -1;
1506               for (i = firstDynamic_; i < firstAvailable_; i++)
1507                    back[id_[i-firstDynamic_]] = i;
1508               int sequence = startSet_[kSet];
1509               if (toIndex_[kSet] < 0) {
1510                    printf("Not in - Set %d - slack status %d, key %d\n", kSet, status_[kSet], keyVariable_[kSet]);
1511                    while (sequence >= 0) {
1512                         printf("( %d status %d ) ", sequence, getDynamicStatus(sequence));
1513                         sequence = next_[sequence];
1514                    }
1515               } else {
1516                    int iRow = numberStaticRows_ + toIndex_[kSet];
1517                    printf("In - Set %d - slack status %d, key %d offset %g slack %g\n", kSet, status_[kSet], keyVariable_[kSet],
1518                           rhsOffset_[iRow], model->solutionRegion(0)[iRow]);
1519                    while (sequence >= 0) {
1520                         int iBack = back[sequence];
1521                         printf("( %d status %d value %g) ", sequence, getDynamicStatus(sequence), model->solutionRegion()[iBack]);
1522                         sequence = next_[sequence];
1523                    }
1524               }
1525               printf("\n");
1526               delete [] back;
1527          }
1528     }
1529#endif
1530     int n = numberActiveSets_;
1531     for (i = 0; i < numberSets_; i++) {
1532          if (toIndex_[i] < 0) {
1533               //assert(keyValue(i)>=lowerSet_[i]&&keyValue(i)<=upperSet_[i]);
1534               n++;
1535          }
1536          int k=0;
1537          for (int j=startSet_[i];j<startSet_[i+1];j++) {
1538            if (getDynamicStatus(j)==soloKey)
1539              k++;
1540          }
1541          assert (k<2);
1542     }
1543     assert (n == numberSets_);
1544#endif
1545     return 1;
1546}
1547void
1548ClpDynamicMatrix::times(double scalar,
1549                        const double * x, double * y) const
1550{
1551     if (model_->specialOptions() != 16) {
1552          ClpPackedMatrix::times(scalar, x, y);
1553     } else {
1554          int iRow;
1555          const double * element =  matrix_->getElements();
1556          const int * row = matrix_->getIndices();
1557          const CoinBigIndex * startColumn = matrix_->getVectorStarts();
1558          const int * length = matrix_->getVectorLengths();
1559          int * pivotVariable = model_->pivotVariable();
1560          for (iRow = 0; iRow < numberStaticRows_ + numberActiveSets_; iRow++) {
1561               y[iRow] -= scalar * rhsOffset_[iRow];
1562               int iColumn = pivotVariable[iRow];
1563               if (iColumn < lastDynamic_) {
1564                    CoinBigIndex j;
1565                    double value = scalar * x[iColumn];
1566                    if (value) {
1567                         for (j = startColumn[iColumn];
1568                                   j < startColumn[iColumn] + length[iColumn]; j++) {
1569                              int jRow = row[j];
1570                              y[jRow] += value * element[j];
1571                         }
1572                    }
1573               }
1574          }
1575     }
1576}
1577// Modifies rhs offset
1578void
1579ClpDynamicMatrix::modifyOffset(int sequence, double amount)
1580{
1581     if (amount) {
1582          assert (rhsOffset_);
1583          CoinBigIndex j;
1584          for (j = startColumn_[sequence]; j < startColumn_[sequence+1]; j++) {
1585               int iRow = row_[j];
1586               rhsOffset_[iRow] += amount * element_[j];
1587          }
1588     }
1589}
1590// Gets key value when none in small
1591double
1592ClpDynamicMatrix::keyValue(int iSet) const
1593{
1594     double value = 0.0;
1595     if (toIndex_[iSet] < 0) {
1596          int key = keyVariable_[iSet];
1597          if (key < maximumGubColumns_) {
1598               if (getStatus(iSet) == ClpSimplex::atLowerBound)
1599                    value = lowerSet_[iSet];
1600               else
1601                    value = upperSet_[iSet];
1602               int numberKey = 0;
1603               int j = startSet_[iSet];
1604               while (j >= 0) {
1605                    DynamicStatus status = getDynamicStatus(j);
1606                    assert (status != inSmall);
1607                    if (status == soloKey) {
1608                         numberKey++;
1609                    } else if (status == atUpperBound) {
1610                         value -= columnUpper_[j];
1611                    } else if (columnLower_) {
1612                         value -= columnLower_[j];
1613                    }
1614                    j = next_[j]; //onto next in set
1615               }
1616               assert (numberKey == 1);
1617          } else {
1618               int j = startSet_[iSet];
1619               while (j >= 0) {
1620                    DynamicStatus status = getDynamicStatus(j);
1621                    assert (status != inSmall);
1622                    assert (status != soloKey);
1623                    if (status == atUpperBound) {
1624                         value += columnUpper_[j];
1625                    } else if (columnLower_) {
1626                         value += columnLower_[j];
1627                    }
1628                    j = next_[j]; //onto next in set
1629               }
1630#if 0
1631               // slack must be feasible
1632               double oldValue=value;
1633               value = CoinMax(value,lowerSet_[iSet]);
1634               value = CoinMin(value,upperSet_[iSet]);
1635               if (value!=oldValue)
1636                 printf("using %g (not %g) for slack on set %d (%g,%g)\n",
1637                        value,oldValue,iSet,lowerSet_[iSet],upperSet_[iSet]);
1638#endif
1639            }
1640     }
1641     return value;
1642}
1643// Switches off dj checking each factorization (for BIG models)
1644void
1645ClpDynamicMatrix::switchOffCheck()
1646{
1647     noCheck_ = 0;
1648     infeasibilityWeight_ = 0.0;
1649}
1650/* Creates a variable.  This is called after partial pricing and may modify matrix.
1651   May update bestSequence.
1652*/
1653void
1654ClpDynamicMatrix::createVariable(ClpSimplex * model, int & bestSequence)
1655{
1656     int numberRows = model->numberRows();
1657     int slackOffset = lastDynamic_ + numberRows;
1658     int structuralOffset = slackOffset + numberSets_;
1659     int bestSequence2 = savedBestSequence_ - structuralOffset;
1660     if (bestSequence >= slackOffset) {
1661          double * columnLower = model->lowerRegion();
1662          double * columnUpper = model->upperRegion();
1663          double * solution = model->solutionRegion();
1664          double * reducedCost = model->djRegion();
1665          const double * duals = model->dualRowSolution();
1666          if (toIndex_[savedBestSet_] < 0) {
1667               // need to put key into basis
1668               int newRow = numberActiveSets_ + numberStaticRows_;
1669               model->dualRowSolution()[newRow] = savedBestGubDual_;
1670               double valueOfKey = keyValue(savedBestSet_); // done before toIndex_ set
1671               toIndex_[savedBestSet_] = numberActiveSets_;
1672               fromIndex_[numberActiveSets_++] = savedBestSet_;
1673               int iSequence = lastDynamic_ + newRow;
1674               // we need to get lower and upper correct
1675               double shift = 0.0;
1676               int j = startSet_[savedBestSet_];
1677               while (j >= 0) {
1678                    if (getDynamicStatus(j) == atUpperBound)
1679                         shift += columnUpper_[j];
1680                    else if (getDynamicStatus(j) == atLowerBound && columnLower_)
1681                         shift += columnLower_[j];
1682                    j = next_[j]; //onto next in set
1683               }
1684               if (lowerSet_[savedBestSet_] > -1.0e20)
1685                    columnLower[iSequence] = lowerSet_[savedBestSet_];
1686               else
1687                    columnLower[iSequence] = -COIN_DBL_MAX;
1688               if (upperSet_[savedBestSet_] < 1.0e20)
1689                    columnUpper[iSequence] = upperSet_[savedBestSet_];
1690               else
1691                    columnUpper[iSequence] = COIN_DBL_MAX;
1692#ifdef CLP_DEBUG
1693               if (model_->logLevel() == 63) {
1694                    printf("%d in in set %d, key is %d rhs %g %g - keyvalue %g\n",
1695                           bestSequence2, savedBestSet_, keyVariable_[savedBestSet_],
1696                           columnLower[iSequence], columnUpper[iSequence], valueOfKey);
1697                    int j = startSet_[savedBestSet_];
1698                    while (j >= 0) {
1699                         if (getDynamicStatus(j) == atUpperBound)
1700                              printf("%d atup ", j);
1701                         else if (getDynamicStatus(j) == atLowerBound)
1702                              printf("%d atlo ", j);
1703                         else if (getDynamicStatus(j) == soloKey)
1704                              printf("%d solo ", j);
1705                         else
1706                              abort();
1707                         j = next_[j]; //onto next in set
1708                    }
1709                    printf("\n");
1710               }
1711#endif
1712               if (keyVariable_[savedBestSet_] < maximumGubColumns_) {
1713                    // slack not key
1714                    model_->pivotVariable()[newRow] = firstAvailable_;
1715                    backToPivotRow_[firstAvailable_] = newRow;
1716                    model->setStatus(iSequence, getStatus(savedBestSet_));
1717                    model->djRegion()[iSequence] = savedBestGubDual_;
1718                    solution[iSequence] = valueOfKey;
1719                    // create variable and pivot in
1720                    int key = keyVariable_[savedBestSet_];
1721                    setDynamicStatus(key, inSmall);
1722                    double * element =  matrix_->getMutableElements();
1723                    int * row = matrix_->getMutableIndices();
1724                    CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
1725                    int * length = matrix_->getMutableVectorLengths();
1726                    CoinBigIndex numberElements = startColumn[firstAvailable_];
1727                    int numberThis = startColumn_[key+1] - startColumn_[key] + 1;
1728                    if (numberElements + numberThis > numberElements_) {
1729                         // need to redo
1730                         numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
1731                         matrix_->reserve(lastDynamic_, numberElements_);
1732                         element =  matrix_->getMutableElements();
1733                         row = matrix_->getMutableIndices();
1734                         // these probably okay but be safe
1735                         startColumn = matrix_->getMutableVectorStarts();
1736                         length = matrix_->getMutableVectorLengths();
1737                    }
1738                    // already set startColumn[firstAvailable_]=numberElements;
1739                    length[firstAvailable_] = numberThis;
1740                    model->costRegion()[firstAvailable_] = cost_[key];
1741                    CoinBigIndex base = startColumn_[key];
1742                    for (int j = 0; j < numberThis - 1; j++) {
1743                         row[numberElements] = row_[base+j];
1744                         element[numberElements++] = element_[base+j];
1745                    }
1746                    row[numberElements] = newRow;
1747                    element[numberElements++] = 1.0;
1748                    id_[firstAvailable_-firstDynamic_] = key;
1749                    model->setObjectiveOffset(model->objectiveOffset() + cost_[key]*
1750                                              valueOfKey);
1751                    model->solutionRegion()[firstAvailable_] = valueOfKey;
1752                    model->setStatus(firstAvailable_, ClpSimplex::basic);
1753                    // ***** need to adjust effective rhs
1754                    if (!columnLower_ && !columnUpper_) {
1755                         columnLower[firstAvailable_] = 0.0;
1756                         columnUpper[firstAvailable_] = COIN_DBL_MAX;
1757                    }  else {
1758                         if (columnLower_)
1759                              columnLower[firstAvailable_] = columnLower_[key];
1760                         else
1761                              columnLower[firstAvailable_] = 0.0;
1762                         if (columnUpper_)
1763                              columnUpper[firstAvailable_] = columnUpper_[key];
1764                         else
1765                              columnUpper[firstAvailable_] = COIN_DBL_MAX;
1766                    }
1767                    model->nonLinearCost()->setOne(firstAvailable_, solution[firstAvailable_],
1768                                                   columnLower[firstAvailable_],
1769                                                   columnUpper[firstAvailable_], cost_[key]);
1770                    startColumn[firstAvailable_+1] = numberElements;
1771                    reducedCost[firstAvailable_] = 0.0;
1772                    modifyOffset(key, valueOfKey);
1773                    rhsOffset_[newRow] = -shift; // sign?
1774#ifdef CLP_DEBUG
1775                    model->rowArray(1)->checkClear();
1776#endif
1777                    // now pivot in
1778                    unpack(model, model->rowArray(1), firstAvailable_);
1779                    model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(1));
1780                    double alpha = model->rowArray(1)->denseVector()[newRow];
1781                    int updateStatus =
1782                         model->factorization()->replaceColumn(model,
1783                                   model->rowArray(2),
1784                                   model->rowArray(1),
1785                                   newRow, alpha);
1786                    model->rowArray(1)->clear();
1787                    if (updateStatus) {
1788                         if (updateStatus == 3) {
1789                              // out of memory
1790                              // increase space if not many iterations
1791                              if (model->factorization()->pivots() <
1792                                        0.5 * model->factorization()->maximumPivots() &&
1793                                        model->factorization()->pivots() < 400)
1794                                   model->factorization()->areaFactor(
1795                                        model->factorization()->areaFactor() * 1.1);
1796                         } else {
1797                              printf("Bad returncode %d from replaceColumn\n", updateStatus);
1798                         }
1799                         bestSequence = -1;
1800                         return;
1801                    }
1802                    // firstAvailable_ only finally updated if good pivot (in updatePivot)
1803                    // otherwise it reverts to firstAvailableBefore_
1804                    firstAvailable_++;
1805               } else {
1806                    // slack key
1807                    model->setStatus(iSequence, ClpSimplex::basic);
1808                    model->djRegion()[iSequence] = 0.0;
1809                    solution[iSequence] = valueOfKey+shift;
1810                    rhsOffset_[newRow] = -shift; // sign?
1811               }
1812               // correct slack
1813               model->costRegion()[iSequence] = 0.0;
1814               model->nonLinearCost()->setOne(iSequence, solution[iSequence], columnLower[iSequence],
1815                                              columnUpper[iSequence], 0.0);
1816          }
1817          if (savedBestSequence_ >= structuralOffset) {
1818               // recompute dj and create
1819               double value = cost_[bestSequence2] - savedBestGubDual_;
1820               for (CoinBigIndex jBigIndex = startColumn_[bestSequence2];
1821                         jBigIndex < startColumn_[bestSequence2+1]; jBigIndex++) {
1822                    int jRow = row_[jBigIndex];
1823                    value -= duals[jRow] * element_[jBigIndex];
1824               }
1825               int gubRow = toIndex_[savedBestSet_] + numberStaticRows_;
1826               double * element =  matrix_->getMutableElements();
1827               int * row = matrix_->getMutableIndices();
1828               CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
1829               int * length = matrix_->getMutableVectorLengths();
1830               CoinBigIndex numberElements = startColumn[firstAvailable_];
1831               int numberThis = startColumn_[bestSequence2+1] - startColumn_[bestSequence2] + 1;
1832               if (numberElements + numberThis > numberElements_) {
1833                    // need to redo
1834                    numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
1835                    matrix_->reserve(lastDynamic_, numberElements_);
1836                    element =  matrix_->getMutableElements();
1837                    row = matrix_->getMutableIndices();
1838                    // these probably okay but be safe
1839                    startColumn = matrix_->getMutableVectorStarts();
1840                    length = matrix_->getMutableVectorLengths();
1841               }
1842               // already set startColumn[firstAvailable_]=numberElements;
1843               length[firstAvailable_] = numberThis;
1844               model->costRegion()[firstAvailable_] = cost_[bestSequence2];
1845               CoinBigIndex base = startColumn_[bestSequence2];
1846               for (int j = 0; j < numberThis - 1; j++) {
1847                    row[numberElements] = row_[base+j];
1848                    element[numberElements++] = element_[base+j];
1849               }
1850               row[numberElements] = gubRow;
1851               element[numberElements++] = 1.0;
1852               id_[firstAvailable_-firstDynamic_] = bestSequence2;
1853               //printf("best %d\n",bestSequence2);
1854               model->solutionRegion()[firstAvailable_] = 0.0;
1855               model->clearFlagged(firstAvailable_);
1856               if (!columnLower_ && !columnUpper_) {
1857                    model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
1858                    columnLower[firstAvailable_] = 0.0;
1859                    columnUpper[firstAvailable_] = COIN_DBL_MAX;
1860               }  else {
1861                    DynamicStatus status = getDynamicStatus(bestSequence2);
1862                    if (columnLower_)
1863                         columnLower[firstAvailable_] = columnLower_[bestSequence2];
1864                    else
1865                         columnLower[firstAvailable_] = 0.0;
1866                    if (columnUpper_)
1867                         columnUpper[firstAvailable_] = columnUpper_[bestSequence2];
1868                    else
1869                         columnUpper[firstAvailable_] = COIN_DBL_MAX;
1870                    if (status == atLowerBound) {
1871                         solution[firstAvailable_] = columnLower[firstAvailable_];
1872                         model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
1873                    } else {
1874                         solution[firstAvailable_] = columnUpper[firstAvailable_];
1875                         model->setStatus(firstAvailable_, ClpSimplex::atUpperBound);
1876                    }
1877               }
1878               model->setObjectiveOffset(model->objectiveOffset() + cost_[bestSequence2]*
1879                                         solution[firstAvailable_]);
1880               model->nonLinearCost()->setOne(firstAvailable_, solution[firstAvailable_],
1881                                              columnLower[firstAvailable_],
1882                                              columnUpper[firstAvailable_], cost_[bestSequence2]);
1883               bestSequence = firstAvailable_;
1884               // firstAvailable_ only updated if good pivot (in updatePivot)
1885               startColumn[firstAvailable_+1] = numberElements;
1886               //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,savedBestGubDual_);
1887               reducedCost[bestSequence] = value;
1888          } else {
1889               // slack - row must just have been created
1890               assert (toIndex_[savedBestSet_] == numberActiveSets_ - 1);
1891               int newRow = numberStaticRows_ + numberActiveSets_ - 1;
1892               bestSequence = lastDynamic_ + newRow;
1893               reducedCost[bestSequence] = savedBestGubDual_;
1894          }
1895     }
1896     // clear for next iteration
1897     savedBestSequence_ = -1;
1898}
1899// Returns reduced cost of a variable
1900double
1901ClpDynamicMatrix::reducedCost(ClpSimplex * model, int sequence) const
1902{
1903     int numberRows = model->numberRows();
1904     int slackOffset = lastDynamic_ + numberRows;
1905     if (sequence < slackOffset)
1906          return model->djRegion()[sequence];
1907     else
1908          return savedBestDj_;
1909}
1910// Does gub crash
1911void
1912ClpDynamicMatrix::gubCrash()
1913{
1914     // Do basis - cheapest or slack if feasible
1915     int longestSet = 0;
1916     int iSet;
1917     for (iSet = 0; iSet < numberSets_; iSet++) {
1918          int n = 0;
1919          int j = startSet_[iSet];
1920          while (j >= 0) {
1921               n++;
1922               j = next_[j];
1923          }
1924          longestSet = CoinMax(longestSet, n);
1925     }
1926     double * upper = new double[longestSet+1];
1927     double * cost = new double[longestSet+1];
1928     double * lower = new double[longestSet+1];
1929     double * solution = new double[longestSet+1];
1930     int * back = new int[longestSet+1];
1931     double tolerance = model_->primalTolerance();
1932     double objectiveOffset = 0.0;
1933     for (iSet = 0; iSet < numberSets_; iSet++) {
1934          int iBasic = -1;
1935          double value = 0.0;
1936          // find cheapest
1937          int numberInSet = 0;
1938          int j = startSet_[iSet];
1939          while (j >= 0) {
1940               if (!columnLower_)
1941                    lower[numberInSet] = 0.0;
1942               else
1943                    lower[numberInSet] = columnLower_[j];
1944               if (!columnUpper_)
1945                    upper[numberInSet] = COIN_DBL_MAX;
1946               else
1947                    upper[numberInSet] = columnUpper_[j];
1948               back[numberInSet++] = j;
1949               j = next_[j];
1950          }
1951          CoinFillN(solution, numberInSet, 0.0);
1952          // and slack
1953          iBasic = numberInSet;
1954          solution[iBasic] = -value;
1955          lower[iBasic] = -upperSet_[iSet];
1956          upper[iBasic] = -lowerSet_[iSet];
1957          int kphase;
1958          if (value < lowerSet_[iSet] - tolerance || value > upperSet_[iSet] + tolerance) {
1959               // infeasible
1960               kphase = 0;
1961               // remember bounds are flipped so opposite to natural
1962               if (value < lowerSet_[iSet] - tolerance)
1963                    cost[iBasic] = 1.0;
1964               else
1965                    cost[iBasic] = -1.0;
1966               CoinZeroN(cost, numberInSet);
1967               double dualTolerance = model_->dualTolerance();
1968               for (int iphase = kphase; iphase < 2; iphase++) {
1969                    if (iphase) {
1970                         cost[numberInSet] = 0.0;
1971                         for (int j = 0; j < numberInSet; j++)
1972                              cost[j] = cost_[back[j]];
1973                    }
1974                    // now do one row lp
1975                    bool improve = true;
1976                    while (improve) {
1977                         improve = false;
1978                         double dual = cost[iBasic];
1979                         int chosen = -1;
1980                         double best = dualTolerance;
1981                         int way = 0;
1982                         for (int i = 0; i <= numberInSet; i++) {
1983                              double dj = cost[i] - dual;
1984                              double improvement = 0.0;
1985                              double distance = 0.0;
1986                              if (iphase || i < numberInSet)
1987                                   assert (solution[i] >= lower[i] && solution[i] <= upper[i]);
1988                              if (dj > dualTolerance)
1989                                   improvement = dj * (solution[i] - lower[i]);
1990                              else if (dj < -dualTolerance)
1991                                   improvement = dj * (solution[i] - upper[i]);
1992                              if (improvement > best) {
1993                                   best = improvement;
1994                                   chosen = i;
1995                                   if (dj < 0.0) {
1996                                        way = 1;
1997                                        distance = upper[i] - solution[i];
1998                                   } else {
1999                                        way = -1;
2000                                        distance = solution[i] - lower[i];
2001                                   }
2002                              }
2003                         }
2004                         if (chosen >= 0) {
2005                              improve = true;
2006                              // now see how far
2007                              if (way > 0) {
2008                                   // incoming increasing so basic decreasing
2009                                   // if phase 0 then go to nearest bound
2010                                   double distance = upper[chosen] - solution[chosen];
2011                                   double basicDistance;
2012                                   if (!iphase) {
2013                                        assert (iBasic == numberInSet);
2014                                        assert (solution[iBasic] > upper[iBasic]);
2015                                        basicDistance = solution[iBasic] - upper[iBasic];
2016                                   } else {
2017                                        basicDistance = solution[iBasic] - lower[iBasic];
2018                                   }
2019                                   // need extra coding for unbounded
2020                                   assert (CoinMin(distance, basicDistance) < 1.0e20);
2021                                   if (distance > basicDistance) {
2022                                        // incoming becomes basic
2023                                        solution[chosen] += basicDistance;
2024                                        if (!iphase)
2025                                             solution[iBasic] = upper[iBasic];
2026                                        else
2027                                             solution[iBasic] = lower[iBasic];
2028                                        iBasic = chosen;
2029                                   } else {
2030                                        // flip
2031                                        solution[chosen] = upper[chosen];
2032                                        solution[iBasic] -= distance;
2033                                   }
2034                              } else {
2035                                   // incoming decreasing so basic increasing
2036                                   // if phase 0 then go to nearest bound
2037                                   double distance = solution[chosen] - lower[chosen];
2038                                   double basicDistance;
2039                                   if (!iphase) {
2040                                        assert (iBasic == numberInSet);
2041                                        assert (solution[iBasic] < lower[iBasic]);
2042                                        basicDistance = lower[iBasic] - solution[iBasic];
2043                                   } else {
2044                                        basicDistance = upper[iBasic] - solution[iBasic];
2045                                   }
2046                                   // need extra coding for unbounded - for now just exit
2047                                   if (CoinMin(distance, basicDistance) > 1.0e20) {
2048                                        printf("unbounded on set %d\n", iSet);
2049                                        iphase = 1;
2050                                        iBasic = numberInSet;
2051                                        break;
2052                                   }
2053                                   if (distance > basicDistance) {
2054                                        // incoming becomes basic
2055                                        solution[chosen] -= basicDistance;
2056                                        if (!iphase)
2057                                             solution[iBasic] = lower[iBasic];
2058                                        else
2059                                             solution[iBasic] = upper[iBasic];
2060                                        iBasic = chosen;
2061                                   } else {
2062                                        // flip
2063                                        solution[chosen] = lower[chosen];
2064                                        solution[iBasic] += distance;
2065                                   }
2066                              }
2067                              if (!iphase) {
2068                                   if(iBasic < numberInSet)
2069                                        break; // feasible
2070                                   else if (solution[iBasic] >= lower[iBasic] &&
2071                                             solution[iBasic] <= upper[iBasic])
2072                                        break; // feasible (on flip)
2073                              }
2074                         }
2075                    }
2076               }
2077          }
2078          // do solution i.e. bounds
2079          if (columnLower_ || columnUpper_) {
2080               for (int j = 0; j < numberInSet; j++) {
2081                    if (j != iBasic) {
2082                         objectiveOffset += solution[j] * cost[j];
2083                         if (columnLower_ && columnUpper_) {
2084                              if (fabs(solution[j] - columnLower_[back[j]]) >
2085                                        fabs(solution[j] - columnUpper_[back[j]]))
2086                                   setDynamicStatus(back[j], atUpperBound);
2087                         } else if (columnUpper_ && solution[j] > 0.0) {
2088                              setDynamicStatus(back[j], atUpperBound);
2089                         } else {
2090                              setDynamicStatus(back[j], atLowerBound);
2091                              assert(!solution[j]);
2092                         }
2093                    }
2094               }
2095          }
2096          // convert iBasic back and do bounds
2097          if (iBasic == numberInSet) {
2098               // slack basic
2099               setStatus(iSet, ClpSimplex::basic);
2100               iBasic = iSet + maximumGubColumns_;
2101          } else {
2102               iBasic = back[iBasic];
2103               setDynamicStatus(iBasic, soloKey);
2104               // remember bounds flipped
2105               if (upper[numberInSet] == lower[numberInSet])
2106                    setStatus(iSet, ClpSimplex::isFixed);
2107               else if (solution[numberInSet] == upper[numberInSet])
2108                    setStatus(iSet, ClpSimplex::atLowerBound);
2109               else if (solution[numberInSet] == lower[numberInSet])
2110                    setStatus(iSet, ClpSimplex::atUpperBound);
2111               else
2112                    abort();
2113          }
2114          keyVariable_[iSet] = iBasic;
2115     }
2116     model_->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
2117     delete [] lower;
2118     delete [] solution;
2119     delete [] upper;
2120     delete [] cost;
2121     delete [] back;
2122     // make sure matrix is in good shape
2123     matrix_->orderMatrix();
2124}
2125// Populates initial matrix from dynamic status
2126void
2127ClpDynamicMatrix::initialProblem()
2128{
2129     int iSet;
2130     double * element =  matrix_->getMutableElements();
2131     int * row = matrix_->getMutableIndices();
2132     CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
2133     int * length = matrix_->getMutableVectorLengths();
2134     double * cost = model_->objective();
2135     double * solution = model_->primalColumnSolution();
2136     double * columnLower = model_->columnLower();
2137     double * columnUpper = model_->columnUpper();
2138     double * rowSolution = model_->primalRowSolution();
2139     double * rowLower = model_->rowLower();
2140     double * rowUpper = model_->rowUpper();
2141     CoinBigIndex numberElements = startColumn[firstDynamic_];
2142
2143     firstAvailable_ = firstDynamic_;
2144     numberActiveSets_ = 0;
2145     for (iSet = 0; iSet < numberSets_; iSet++) {
2146          toIndex_[iSet] = -1;
2147          int numberActive = 0;
2148          int whichKey = -1;
2149          if (getStatus(iSet) == ClpSimplex::basic) {
2150               whichKey = maximumGubColumns_ + iSet;
2151               numberActive = 1;
2152          } else {
2153               whichKey = -1;
2154          }
2155          int j = startSet_[iSet];
2156          while (j >= 0) {
2157               assert (getDynamicStatus(j) != soloKey || whichKey == -1);
2158               if (getDynamicStatus(j) == inSmall) {
2159                    numberActive++;
2160               } else if (getDynamicStatus(j) == soloKey) {
2161                    whichKey = j;
2162                    numberActive++;
2163               }
2164               j = next_[j]; //onto next in set
2165          }
2166          if (numberActive > 1) {
2167               int iRow = numberActiveSets_ + numberStaticRows_;
2168               rowSolution[iRow] = 0.0;
2169               double lowerValue;
2170               if (lowerSet_[iSet] > -1.0e20)
2171                    lowerValue = lowerSet_[iSet];
2172               else
2173                    lowerValue = -COIN_DBL_MAX;
2174               double upperValue;
2175               if (upperSet_[iSet] < 1.0e20)
2176                    upperValue = upperSet_[iSet];
2177               else
2178                    upperValue = COIN_DBL_MAX;
2179               rowLower[iRow] = lowerValue;
2180               rowUpper[iRow] = upperValue;
2181               if (getStatus(iSet) == ClpSimplex::basic) {
2182                    model_->setRowStatus(iRow, ClpSimplex::basic);
2183                    rowSolution[iRow] = 0.0;
2184               } else if (getStatus(iSet) == ClpSimplex::atLowerBound) {
2185                    model_->setRowStatus(iRow, ClpSimplex::atLowerBound);
2186                    rowSolution[iRow] = lowerValue;
2187               } else {
2188                    model_->setRowStatus(iRow, ClpSimplex::atUpperBound);
2189                    rowSolution[iRow] = upperValue;
2190               }
2191               j = startSet_[iSet];
2192               while (j >= 0) {
2193                    DynamicStatus status = getDynamicStatus(j);
2194                    if (status == inSmall) {
2195                         int numberThis = startColumn_[j+1] - startColumn_[j] + 1;
2196                         if (numberElements + numberThis > numberElements_) {
2197                              // need to redo
2198                              numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
2199                              matrix_->reserve(lastDynamic_, numberElements_);
2200                              element =  matrix_->getMutableElements();
2201                              row = matrix_->getMutableIndices();
2202                              // these probably okay but be safe
2203                              startColumn = matrix_->getMutableVectorStarts();
2204                              length = matrix_->getMutableVectorLengths();
2205                         }
2206                         length[firstAvailable_] = numberThis;
2207                         cost[firstAvailable_] = cost_[j];
2208                         CoinBigIndex base = startColumn_[j];
2209                         for (int k = 0; k < numberThis - 1; k++) {
2210                              row[numberElements] = row_[base+k];
2211                              element[numberElements++] = element_[base+k];
2212                         }
2213                         row[numberElements] = iRow;
2214                         element[numberElements++] = 1.0;
2215                         id_[firstAvailable_-firstDynamic_] = j;
2216                         solution[firstAvailable_] = 0.0;
2217                         model_->setStatus(firstAvailable_, ClpSimplex::basic);
2218                         if (!columnLower_ && !columnUpper_) {
2219                              columnLower[firstAvailable_] = 0.0;
2220                              columnUpper[firstAvailable_] = COIN_DBL_MAX;
2221                         }  else {
2222                              if (columnLower_)
2223                                   columnLower[firstAvailable_] = columnLower_[j];
2224                              else
2225                                   columnLower[firstAvailable_] = 0.0;
2226                              if (columnUpper_)
2227                                   columnUpper[firstAvailable_] = columnUpper_[j];
2228                              else
2229                                   columnUpper[firstAvailable_] = COIN_DBL_MAX;
2230                              if (status != atUpperBound) {
2231                                   solution[firstAvailable_] = columnLower[firstAvailable_];
2232                              } else {
2233                                   solution[firstAvailable_] = columnUpper[firstAvailable_];
2234                              }
2235                         }
2236                         firstAvailable_++;
2237                         startColumn[firstAvailable_] = numberElements;
2238                    }
2239                    j = next_[j]; //onto next in set
2240               }
2241               model_->setRowStatus(numberActiveSets_ + numberStaticRows_, getStatus(iSet));
2242               toIndex_[iSet] = numberActiveSets_;
2243               fromIndex_[numberActiveSets_++] = iSet;
2244          } else {
2245            // solo key
2246            bool needKey;
2247            if (numberActive) {
2248              if (whichKey<maximumGubColumns_) {
2249                // structural - assume ok
2250                needKey = false;
2251              } else {
2252                // slack
2253                keyVariable_[iSet] = maximumGubColumns_ + iSet;
2254                double value = keyValue(iSet);
2255                if (value<lowerSet_[iSet]-1.0e-8||
2256                    value>upperSet_[iSet]+1.0e-8)
2257                  needKey=true;
2258              }
2259            } else {
2260              needKey = true;
2261            }
2262            if (needKey) {
2263              // all to lb then up some (slack/null if possible)
2264              int length=99999999;
2265              int which=-1;
2266              double sum=0.0;
2267              for (int iColumn=startSet_[iSet];iColumn<startSet_[iSet+1];iColumn++) {
2268                setDynamicStatus(iColumn,atLowerBound);
2269                sum += columnLower_[iColumn];
2270                if (length>startColumn_[iColumn+1]-startColumn_[iColumn]) {
2271                  which=iColumn;
2272                  length=startColumn_[iColumn+1]-startColumn_[iColumn];
2273                }
2274              }
2275              if (sum>lowerSet_[iSet]-1.0e-8) {
2276                // slack can be basic
2277                setStatus(iSet,ClpSimplex::basic);
2278                keyVariable_[iSet] = maximumGubColumns_ + iSet;
2279              } else {
2280                // use shortest
2281                setDynamicStatus(which,soloKey);
2282                keyVariable_[iSet] = which;
2283                setStatus(iSet,ClpSimplex::atLowerBound);
2284              }
2285            }
2286          }
2287          assert (toIndex_[iSet] >= 0 || whichKey >= 0);
2288          keyVariable_[iSet] = whichKey;
2289     }
2290     numberActiveColumns_ = firstAvailable_;
2291     return;
2292}
2293// Writes out model (without names)
2294void 
2295ClpDynamicMatrix::writeMps(const char * name)
2296{
2297  int numberTotalRows = numberStaticRows_+numberSets_;
2298  int numberTotalColumns = firstDynamic_+numberGubColumns_;
2299  // over estimate
2300  int numberElements = getNumElements()+startColumn_[numberGubColumns_]
2301    + numberGubColumns_;
2302  double * columnLower = new double [numberTotalColumns];
2303  double * columnUpper = new double [numberTotalColumns];
2304  double * cost = new double [numberTotalColumns];
2305  double * rowLower = new double [numberTotalRows];
2306  double * rowUpper = new double [numberTotalRows];
2307  CoinBigIndex * start = new CoinBigIndex[numberTotalColumns+1];
2308  int * row = new int [numberElements];
2309  double * element = new double [numberElements];
2310  // Fill in
2311  const CoinBigIndex * startA = getVectorStarts();
2312  const int * lengthA = getVectorLengths();
2313  const int * rowA = getIndices();
2314  const double * elementA = getElements();
2315  const double * columnLowerA = model_->columnLower();
2316  const double * columnUpperA = model_->columnUpper();
2317  const double * costA = model_->objective();
2318  const double * rowLowerA = model_->rowLower();
2319  const double * rowUpperA = model_->rowUpper();
2320  start[0]=0;
2321  numberElements=0;
2322  for (int i=0;i<firstDynamic_;i++) {
2323    columnLower[i] = columnLowerA[i];
2324    columnUpper[i] = columnUpperA[i];
2325    cost[i] = costA[i];
2326    for (CoinBigIndex j = startA[i];j<startA[i]+lengthA[i];j++) {
2327      row[numberElements] = rowA[j];
2328      element[numberElements++]=elementA[j];
2329    }
2330    start[i+1]=numberElements;
2331  }
2332  for (int i=0;i<numberStaticRows_;i++) {
2333    rowLower[i] = rowLowerA[i];
2334    rowUpper[i] = rowUpperA[i];
2335  }
2336  int putC=firstDynamic_;
2337  int putR=numberStaticRows_;
2338  for (int i=0;i<numberSets_;i++) {
2339    rowLower[putR]=lowerSet_[i];
2340    rowUpper[putR]=upperSet_[i];
2341    for (CoinBigIndex k=startSet_[i];k<startSet_[i+1];k++) {
2342      columnLower[putC]=columnLower_[k];
2343      columnUpper[putC]=columnUpper_[k];
2344      cost[putC]=cost_[k];
2345      putC++;
2346      for (CoinBigIndex j = startColumn_[k];j<startColumn_[k+1];j++) {
2347        row[numberElements] = row_[j];
2348        element[numberElements++]=element_[j];
2349      }
2350      row[numberElements] = putR;
2351      element[numberElements++]=1.0;
2352      start[putC]=numberElements;
2353    }
2354    putR++;
2355  }
2356
2357  assert (putR==numberTotalRows);
2358  assert (putC==numberTotalColumns);
2359  ClpSimplex modelOut;
2360  modelOut.loadProblem(numberTotalColumns,numberTotalRows,
2361                       start,row,element,
2362                       columnLower,columnUpper,cost,
2363                       rowLower,rowUpper);
2364  modelOut.writeMps(name);
2365  delete [] columnLower;
2366  delete [] columnUpper;
2367  delete [] cost;
2368  delete [] rowLower;
2369  delete [] rowUpper;
2370  delete [] start;
2371  delete [] row; 
2372  delete [] element;
2373}
2374// Adds in a column to gub structure (called from descendant)
2375int
2376ClpDynamicMatrix::addColumn(int numberEntries, const int * row, const double * element,
2377                            double cost, double lower, double upper, int iSet,
2378                            DynamicStatus status)
2379{
2380     // check if already in
2381     int j = startSet_[iSet];
2382     while (j >= 0) {
2383          if (startColumn_[j+1] - startColumn_[j] == numberEntries) {
2384               const int * row2 = row_ + startColumn_[j];
2385               const double * element2 = element_ + startColumn_[j];
2386               bool same = true;
2387               for (int k = 0; k < numberEntries; k++) {
2388                    if (row[k] != row2[k] || element[k] != element2[k]) {
2389                         same = false;
2390                         break;
2391                    }
2392               }
2393               if (same) {
2394                    bool odd = false;
2395                    if (cost != cost_[j])
2396                         odd = true;
2397                    if (columnLower_ && lower != columnLower_[j])
2398                         odd = true;
2399                    if (columnUpper_ && upper != columnUpper_[j])
2400                         odd = true;
2401                    if (odd) {
2402                         printf("seems odd - same els but cost,lo,up are %g,%g,%g and %g,%g,%g\n",
2403                                cost, lower, upper, cost_[j],
2404                                columnLower_ ? columnLower_[j] : 0.0,
2405                                columnUpper_ ? columnUpper_[j] : 1.0e100);
2406                    } else {
2407                         setDynamicStatus(j, status);
2408                         return j;
2409                    }
2410               }
2411          }
2412          j = next_[j];
2413     }
2414
2415     if (numberGubColumns_ == maximumGubColumns_ ||
2416               startColumn_[numberGubColumns_] + numberEntries > maximumElements_) {
2417          CoinBigIndex j;
2418          int i;
2419          int put = 0;
2420          int numberElements = 0;
2421          CoinBigIndex start = 0;
2422          // compress - leave ones at ub and basic
2423          int * which = new int [numberGubColumns_];
2424          for (i = 0; i < numberGubColumns_; i++) {
2425               CoinBigIndex end = startColumn_[i+1];
2426               // what about ubs if column generation?
2427               if (getDynamicStatus(i) != atLowerBound) {
2428                    // keep in
2429                    for (j = start; j < end; j++) {
2430                         row_[numberElements] = row_[j];
2431                         element_[numberElements++] = element_[j];
2432                    }
2433                    startColumn_[put+1] = numberElements;
2434                    cost_[put] = cost_[i];
2435                    if (columnLower_)
2436                         columnLower_[put] = columnLower_[i];
2437                    if (columnUpper_)
2438                         columnUpper_[put] = columnUpper_[i];
2439                    dynamicStatus_[put] = dynamicStatus_[i];
2440                    id_[put] = id_[i];
2441                    which[i] = put;
2442                    put++;
2443               } else {
2444                    // out
2445                    which[i] = -1;
2446               }
2447               start = end;
2448          }
2449          // now redo startSet_ and next_
2450          int * newNext = new int [maximumGubColumns_];
2451          for (int jSet = 0; jSet < numberSets_; jSet++) {
2452               int sequence = startSet_[jSet];
2453               while (which[sequence] < 0) {
2454                    // out
2455                    assert (next_[sequence] >= 0);
2456                    sequence = next_[sequence];
2457               }
2458               startSet_[jSet] = which[sequence];
2459               int last = which[sequence];
2460               while (next_[sequence] >= 0) {
2461                    sequence = next_[sequence];
2462                    if(which[sequence] >= 0) {
2463                         // keep
2464                         newNext[last] = which[sequence];
2465                         last = which[sequence];
2466                    }
2467               }
2468               newNext[last] = -jSet - 1;
2469          }
2470          delete [] next_;
2471          next_ = newNext;
2472          delete [] which;
2473          abort();
2474     }
2475     CoinBigIndex start = startColumn_[numberGubColumns_];
2476     CoinMemcpyN(row, numberEntries, row_ + start);
2477     CoinMemcpyN(element, numberEntries, element_ + start);
2478     startColumn_[numberGubColumns_+1] = start + numberEntries;
2479     cost_[numberGubColumns_] = cost;
2480     if (columnLower_)
2481          columnLower_[numberGubColumns_] = lower;
2482     else
2483          assert (!lower);
2484     if (columnUpper_)
2485          columnUpper_[numberGubColumns_] = upper;
2486     else
2487          assert (upper > 1.0e20);
2488     setDynamicStatus(numberGubColumns_, status);
2489     // Do next_
2490     j = startSet_[iSet];
2491     startSet_[iSet] = numberGubColumns_;
2492     next_[numberGubColumns_] = j;
2493     numberGubColumns_++;
2494     return numberGubColumns_ - 1;
2495}
2496// Returns which set a variable is in
2497int
2498ClpDynamicMatrix::whichSet (int sequence) const
2499{
2500     while (next_[sequence] >= 0)
2501          sequence = next_[sequence];
2502     int iSet = - next_[sequence] - 1;
2503     return iSet;
2504}
Note: See TracBrowser for help on using the repository browser.