source: stable/1.15/Clp/src/ClpDynamicMatrix.cpp @ 1941

Last change on this file since 1941 was 1941, checked in by stefan, 6 years ago

sync with trunk rev 1940

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