source: trunk/Clp/src/ClpGubDynamicMatrix.cpp @ 1034

Last change on this file since 1034 was 1034, checked in by forrest, 13 years ago

moving branches/devel to trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 67.4 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5#include <cstdio>
6
7#include "CoinPragma.hpp"
8#include "CoinIndexedVector.hpp"
9#include "CoinHelperFunctions.hpp"
10
11#include "ClpSimplex.hpp"
12#include "ClpFactorization.hpp"
13#include "ClpQuadraticObjective.hpp"
14#include "ClpNonLinearCost.hpp"
15// at end to get min/max!
16#include "ClpGubDynamicMatrix.hpp"
17#include "ClpMessage.hpp"
18//#define CLP_DEBUG
19//#define CLP_DEBUG_PRINT
20//#############################################################################
21// Constructors / Destructor / Assignment
22//#############################################################################
23
24//-------------------------------------------------------------------
25// Default Constructor
26//-------------------------------------------------------------------
27ClpGubDynamicMatrix::ClpGubDynamicMatrix () 
28  : ClpGubMatrix(),
29    objectiveOffset_(0.0),
30    startColumn_(NULL),
31    row_(NULL),
32    element_(NULL),
33    cost_(NULL),
34    fullStart_(NULL),
35    id_(NULL),
36    dynamicStatus_(NULL),
37    lowerColumn_(NULL),
38    upperColumn_(NULL),
39    lowerSet_(NULL),
40    upperSet_(NULL),
41    numberGubColumns_(0),
42    firstAvailable_(0),
43    savedFirstAvailable_(0),
44    firstDynamic_(0),
45    lastDynamic_(0),
46    numberElements_(0)
47{
48  setType(13);
49}
50
51//-------------------------------------------------------------------
52// Copy constructor
53//-------------------------------------------------------------------
54ClpGubDynamicMatrix::ClpGubDynamicMatrix (const ClpGubDynamicMatrix & rhs) 
55: ClpGubMatrix(rhs)
56{ 
57  objectiveOffset_ = rhs.objectiveOffset_;
58  numberGubColumns_ = rhs.numberGubColumns_;
59  firstAvailable_ = rhs.firstAvailable_;
60  savedFirstAvailable_ = rhs.savedFirstAvailable_;
61  firstDynamic_ = rhs.firstDynamic_;
62  lastDynamic_ = rhs.lastDynamic_;
63  numberElements_ = rhs.numberElements_;
64  startColumn_ = ClpCopyOfArray(rhs.startColumn_,numberGubColumns_+1);
65  CoinBigIndex numberElements = startColumn_[numberGubColumns_];
66  row_ = ClpCopyOfArray(rhs.row_,numberElements);;
67  element_ = ClpCopyOfArray(rhs.element_,numberElements);;
68  cost_ = ClpCopyOfArray(rhs.cost_,numberGubColumns_);
69  fullStart_ = ClpCopyOfArray(rhs.fullStart_,numberSets_+1);
70  id_ = ClpCopyOfArray(rhs.id_,lastDynamic_-firstDynamic_);
71  lowerColumn_ = ClpCopyOfArray(rhs.lowerColumn_,numberGubColumns_);
72  upperColumn_ = ClpCopyOfArray(rhs.upperColumn_,numberGubColumns_);
73  dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_,numberGubColumns_);
74  lowerSet_ = ClpCopyOfArray(rhs.lowerSet_,numberSets_);
75  upperSet_ = ClpCopyOfArray(rhs.upperSet_,numberSets_);
76}
77
78/* This is the real constructor*/
79ClpGubDynamicMatrix::ClpGubDynamicMatrix(ClpSimplex * model, int numberSets,
80                                         int numberGubColumns, const int * starts,
81                                         const double * lower, const double * upper,
82                                         const CoinBigIndex * startColumn, const int * row,
83                                         const double * element, const double * cost,
84                                         const double * lowerColumn, const double * upperColumn,
85                                         const unsigned char * status)
86  : ClpGubMatrix()
87{
88  objectiveOffset_ = model->objectiveOffset();
89  model_ = model;
90  numberSets_ = numberSets;
91  numberGubColumns_ = numberGubColumns;
92  fullStart_ = ClpCopyOfArray(starts,numberSets_+1);
93  lower_ = ClpCopyOfArray(lower,numberSets_);
94  upper_ = ClpCopyOfArray(upper,numberSets_);
95  int numberColumns = model->numberColumns();
96  int numberRows = model->numberRows();
97  // Number of columns needed
98  int numberGubInSmall = numberSets_+numberRows + 2*model->factorizationFrequency()+2;
99  // for small problems this could be too big
100  //numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_);
101  int numberNeeded = numberGubInSmall + numberColumns;
102  firstAvailable_ = numberColumns;
103  savedFirstAvailable_ = numberColumns;
104  firstDynamic_ = numberColumns;
105  lastDynamic_ = numberNeeded;
106  startColumn_ = ClpCopyOfArray(startColumn,numberGubColumns_+1);
107  CoinBigIndex numberElements = startColumn_[numberGubColumns_];
108  row_ = ClpCopyOfArray(row,numberElements);
109  element_ = new float[numberElements];
110  CoinBigIndex i;
111  for (i=0;i<numberElements;i++)
112    element_[i]=element[i];
113  cost_ = new float[numberGubColumns_];
114  for (i=0;i<numberGubColumns_;i++) {
115    cost_[i]=cost[i];
116    // need sorted
117    CoinSort_2(row_+startColumn_[i],row_+startColumn_[i+1],element_+startColumn_[i]);
118  }
119  if (lowerColumn) {
120    lowerColumn_ = new float[numberGubColumns_];
121    for (i=0;i<numberGubColumns_;i++) 
122      lowerColumn_[i]=lowerColumn[i];
123  } else {
124    lowerColumn_=NULL;
125  }
126  if (upperColumn) {
127    upperColumn_ = new float[numberGubColumns_];
128    for (i=0;i<numberGubColumns_;i++) 
129      upperColumn_[i]=upperColumn[i];
130  } else {
131    upperColumn_=NULL;
132  }
133  if (upperColumn||lowerColumn) {
134    lowerSet_ = new float[numberSets_];
135    for (i=0;i<numberSets_;i++) {
136      if (lower[i]>-1.0e20)
137        lowerSet_[i]=lower[i];
138      else
139        lowerSet_[i] = -1.0e30;
140    }
141    upperSet_ = new float[numberSets_];
142    for (i=0;i<numberSets_;i++) {
143      if (upper[i]<1.0e20)
144        upperSet_[i]=upper[i];
145      else
146        upperSet_[i] = 1.0e30;
147    }
148  } else {
149    lowerSet_=NULL;
150    upperSet_=NULL;
151  }
152  start_=NULL;
153  end_=NULL;
154  dynamicStatus_=NULL;
155  id_ = new int[numberGubInSmall]; 
156  for (i=0;i<numberGubInSmall;i++)
157    id_[i]=-1;
158  ClpPackedMatrix* originalMatrixA =
159    dynamic_cast< ClpPackedMatrix*>(model->clpMatrix());
160  assert (originalMatrixA);
161  CoinPackedMatrix * originalMatrix = originalMatrixA->getPackedMatrix();
162  originalMatrixA->setMatrixNull(); // so can be deleted safely
163  // guess how much space needed
164  double guess = originalMatrix->getNumElements()+10;
165  guess /= (double) numberColumns;
166  guess *= 2*numberGubColumns_;
167  numberElements_ = (int) CoinMin(guess,10000000.0);
168  numberElements_ = CoinMin(numberElements_,numberElements)+originalMatrix->getNumElements();
169  matrix_ = originalMatrix;
170  flags_ &= ~1;
171  // resize model (matrix stays same)
172  model->resize(numberRows,numberNeeded);
173  if (upperColumn_) {
174    // set all upper bounds so we have enough space
175    double * columnUpper = model->columnUpper();
176    for(i=firstDynamic_;i<lastDynamic_;i++)
177      columnUpper[i]=1.0e10;
178  }
179  // resize matrix
180  // extra 1 is so can keep number of elements handy
181  originalMatrix->reserve(numberNeeded,numberElements_,true);
182  originalMatrix->reserve(numberNeeded+1,numberElements_,false);
183  originalMatrix->getMutableVectorStarts()[numberColumns]=originalMatrix->getNumElements();
184  // redo number of columns
185  numberColumns = matrix_->getNumCols();
186  backward_ = new int[numberNeeded];
187  backToPivotRow_ = new int[numberNeeded];
188  // We know a bit better
189  delete [] changeCost_;
190  changeCost_ = new double [numberRows+numberSets_];
191  keyVariable_ = new int[numberSets_];
192  // signal to need new ordering
193  next_ = NULL;
194  for (int iColumn=0;iColumn<numberNeeded;iColumn++) 
195    backward_[iColumn]=-1;
196
197  firstGub_=firstDynamic_;
198  lastGub_=lastDynamic_;
199  if (!lowerColumn_&&!upperColumn_)
200    gubType_=8;
201  if (status) {
202    status_ = ClpCopyOfArray(status,numberSets_);
203  } else {
204    status_= new unsigned char [numberSets_];
205    memset(status_,0,numberSets_);
206    int i;
207    for (i=0;i<numberSets_;i++) {
208      // make slack key
209      setStatus(i,ClpSimplex::basic);
210    }
211  }
212  saveStatus_= new unsigned char [numberSets_];
213  memset(saveStatus_,0,numberSets_);
214  savedKeyVariable_= new int [numberSets_];
215  memset(savedKeyVariable_,0,numberSets_*sizeof(int));
216}
217
218//-------------------------------------------------------------------
219// Destructor
220//-------------------------------------------------------------------
221ClpGubDynamicMatrix::~ClpGubDynamicMatrix ()
222{
223  delete [] startColumn_;
224  delete [] row_;
225  delete [] element_;
226  delete [] cost_;
227  delete [] fullStart_;
228  delete [] id_;
229  delete [] dynamicStatus_;
230  delete [] lowerColumn_;
231  delete [] upperColumn_;
232  delete [] lowerSet_;
233  delete [] upperSet_;
234}
235
236//----------------------------------------------------------------
237// Assignment operator
238//-------------------------------------------------------------------
239ClpGubDynamicMatrix &
240ClpGubDynamicMatrix::operator=(const ClpGubDynamicMatrix& rhs)
241{
242  if (this != &rhs) {
243    ClpGubMatrix::operator=(rhs);
244    delete [] startColumn_;
245    delete [] row_;
246    delete [] element_;
247    delete [] cost_;
248    delete [] fullStart_;
249    delete [] id_;
250    delete [] dynamicStatus_;
251    delete [] lowerColumn_;
252    delete [] upperColumn_;
253    delete [] lowerSet_;
254    delete [] upperSet_;
255    objectiveOffset_ = rhs.objectiveOffset_;
256    numberGubColumns_ = rhs.numberGubColumns_;
257    firstAvailable_ = rhs.firstAvailable_;
258    savedFirstAvailable_ = rhs.savedFirstAvailable_;
259    firstDynamic_ = rhs.firstDynamic_;
260    lastDynamic_ = rhs.lastDynamic_;
261    numberElements_ = rhs.numberElements_;
262    startColumn_ = ClpCopyOfArray(rhs.startColumn_,numberGubColumns_+1);
263    int numberElements = startColumn_[numberGubColumns_];
264    row_ = ClpCopyOfArray(rhs.row_,numberElements);;
265    element_ = ClpCopyOfArray(rhs.element_,numberElements);;
266    cost_ = ClpCopyOfArray(rhs.cost_,numberGubColumns_);
267    fullStart_ = ClpCopyOfArray(rhs.fullStart_,numberSets_+1);
268    id_ = ClpCopyOfArray(rhs.id_,lastDynamic_-firstDynamic_);
269    lowerColumn_ = ClpCopyOfArray(rhs.lowerColumn_,numberGubColumns_);
270    upperColumn_ = ClpCopyOfArray(rhs.upperColumn_,numberGubColumns_);
271    dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_,numberGubColumns_);
272    lowerSet_ = ClpCopyOfArray(rhs.lowerSet_,numberSets_);
273    upperSet_ = ClpCopyOfArray(rhs.upperSet_,numberSets_);
274  }
275  return *this;
276}
277//-------------------------------------------------------------------
278// Clone
279//-------------------------------------------------------------------
280ClpMatrixBase * ClpGubDynamicMatrix::clone() const
281{
282  return new ClpGubDynamicMatrix(*this);
283}
284// Partial pricing
285void 
286ClpGubDynamicMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
287                              int & bestSequence, int & numberWanted)
288{
289  assert(!model->rowScale());
290  numberWanted=currentWanted_;
291  if (!numberSets_) {
292    // no gub
293    ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
294    return;
295  } else {
296    // and do some proportion of full set
297    int startG2 = (int) (startFraction*numberSets_);
298    int endG2 = (int) (endFraction*numberSets_+0.1);
299    endG2 = CoinMin(endG2,numberSets_);
300    //printf("gub price - set start %d end %d\n",
301    //   startG2,endG2);
302    double tolerance=model->currentDualTolerance();
303    double * reducedCost = model->djRegion();
304    const double * duals = model->dualRowSolution();
305    double * cost = model->costRegion();
306    double bestDj;
307    int numberRows = model->numberRows();
308    int numberColumns = lastDynamic_;
309    // If nothing found yet can go all the way to end
310    int endAll = endG2;
311    if (bestSequence<0&&!startG2)
312      endAll = numberSets_;
313    if (bestSequence>=0)
314      bestDj = fabs(reducedCost[bestSequence]);
315    else
316      bestDj=tolerance;
317    int saveSequence = bestSequence;
318    double djMod=0.0;
319    double infeasibilityCost = model->infeasibilityCost();
320    double bestDjMod=0.0;
321    //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
322    //     startG2,endG2,numberWanted);
323    int bestType=-1;
324    int bestSet=-1;
325    const double * element =matrix_->getElements();
326    const int * row = matrix_->getIndices();
327    const CoinBigIndex * startColumn = matrix_->getVectorStarts();
328    int * length = matrix_->getMutableVectorLengths();
329#if 0   
330    // make sure first available is clean (in case last iteration rejected)
331    cost[firstAvailable_]=0.0;
332    length[firstAvailable_]=0;
333    model->nonLinearCost()->setOne(firstAvailable_,0.0,0.0,COIN_DBL_MAX,0.0);
334    model->setStatus(firstAvailable_,ClpSimplex::atLowerBound);
335    {
336      for (int i=firstAvailable_;i<lastDynamic_;i++)
337        assert(!cost[i]);
338    }
339#endif
340#ifdef CLP_DEBUG
341 {
342   for (int i=firstDynamic_;i<firstAvailable_;i++) {
343     assert (getDynamicStatus(id_[i-firstDynamic_])==inSmall);
344   }
345 }
346#endif
347    int minSet = minimumObjectsScan_<0 ? 5 : minimumObjectsScan_; 
348    int minNeg = minimumGoodReducedCosts_<0 ? 5 : minimumGoodReducedCosts_;
349    for (int iSet = startG2;iSet<endAll;iSet++) {
350      if (numberWanted+minNeg<originalWanted_&&iSet>startG2+minSet) {
351        // give up
352        numberWanted=0;
353        break;
354      } else if (iSet==endG2&&bestSequence>=0) {
355        break;
356      }
357      CoinBigIndex j;
358      int iBasic = keyVariable_[iSet];
359      if (iBasic>=numberColumns) {
360        djMod = - weight(iSet)*infeasibilityCost;
361      } else {
362        // get dj without
363        assert (model->getStatus(iBasic)==ClpSimplex::basic);
364        djMod=0.0;
365       
366        for (j=startColumn[iBasic];
367             j<startColumn[iBasic]+length[iBasic];j++) {
368          int jRow=row[j];
369          djMod -= duals[jRow]*element[j];
370        }
371        djMod += cost[iBasic];
372        // See if gub slack possible - dj is djMod
373        if (getStatus(iSet)==ClpSimplex::atLowerBound) {
374          double value = -djMod;
375          if (value>tolerance) {
376            numberWanted--;
377            if (value>bestDj) {
378              // check flagged variable and correct dj
379              if (!flagged(iSet)) {
380                bestDj=value;
381                bestSequence = numberRows+numberColumns+iSet;
382                bestDjMod = djMod;
383                bestType=0;
384                bestSet = iSet;
385              } else {
386                // just to make sure we don't exit before got something
387                numberWanted++;
388                abort();
389              }
390            }
391          }
392        } else if (getStatus(iSet)==ClpSimplex::atUpperBound) {
393          double value = djMod;
394          if (value>tolerance) {
395            numberWanted--;
396            if (value>bestDj) {
397              // check flagged variable and correct dj
398              if (!flagged(iSet)) {
399                bestDj=value;
400                bestSequence = numberRows+numberColumns+iSet;
401                bestDjMod = djMod;
402                bestType=0;
403                bestSet = iSet;
404              } else {
405                // just to make sure we don't exit before got something
406                numberWanted++;
407                abort();
408              }
409            }
410          }
411        }
412      }
413      for (int iSequence=fullStart_[iSet];iSequence<fullStart_[iSet+1];iSequence++) {
414        DynamicStatus status = getDynamicStatus(iSequence);
415        if (status!=inSmall) {
416          double value=cost_[iSequence]-djMod;
417          for (j=startColumn_[iSequence];
418               j<startColumn_[iSequence+1];j++) {
419            int jRow=row_[j];
420            value -= duals[jRow]*element_[j];
421          }
422          // change sign if at lower bound
423          if (status==atLowerBound)
424            value = -value;
425          if (value>tolerance) {
426            numberWanted--;
427            if (value>bestDj) {
428              // check flagged variable and correct dj
429              if (!flagged(iSequence)) {
430                bestDj=value;
431                bestSequence = iSequence;
432                bestDjMod = djMod;
433                bestType=1;
434                bestSet = iSet;
435              } else {
436                // just to make sure we don't exit before got something
437                numberWanted++;
438              }
439            }
440          }
441        }
442      }
443      if (numberWanted<=0) {
444        numberWanted=0;
445        break;
446      }
447    }
448    // Do packed part before gub and small gub - but lightly
449    int saveMinNeg=minimumGoodReducedCosts_;
450    int saveSequence2 = bestSequence;
451    if (bestSequence>=0)
452      minimumGoodReducedCosts_=-2;
453    int saveLast=lastGub_;
454    lastGub_=firstAvailable_;
455    currentWanted_=numberWanted;
456    ClpGubMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
457    minimumGoodReducedCosts_=saveMinNeg;
458    lastGub_=saveLast;
459    if (bestSequence!=saveSequence2) {
460      bestType=-1; // in normal or small gub part
461      saveSequence=bestSequence;
462    }
463    if (bestSequence!=saveSequence||bestType>=0) {
464      double * lowerColumn = model->lowerRegion();
465      double * upperColumn = model->upperRegion();
466      double * solution = model->solutionRegion();
467      if (bestType>0) {
468        // recompute dj and create
469        double value=cost_[bestSequence]-bestDjMod;
470        for (CoinBigIndex jBigIndex=startColumn_[bestSequence];
471             jBigIndex<startColumn_[bestSequence+1];jBigIndex++) {
472          int jRow=row_[jBigIndex];
473          value -= duals[jRow]*element_[jBigIndex];
474        }
475        double * element =  matrix_->getMutableElements();
476        int * row = matrix_->getMutableIndices();
477        CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
478        int * length = matrix_->getMutableVectorLengths();
479        CoinBigIndex numberElements = startColumn[firstAvailable_];
480        int numberThis = startColumn_[bestSequence+1]-startColumn_[bestSequence];
481        if (numberElements+numberThis>numberElements_) {
482          // need to redo
483          numberElements_ = CoinMax(3*numberElements_/2,numberElements+numberThis);
484          matrix_->reserve(numberColumns,numberElements_);
485          element =  matrix_->getMutableElements();
486          row = matrix_->getMutableIndices();
487          // these probably okay but be safe
488          startColumn = matrix_->getMutableVectorStarts();
489          length = matrix_->getMutableVectorLengths();
490        }
491        // already set startColumn[firstAvailable_]=numberElements;
492        length[firstAvailable_]=numberThis;
493        model->costRegion()[firstAvailable_]=cost_[bestSequence];
494        CoinBigIndex base = startColumn_[bestSequence];
495        for (int j=0;j<numberThis;j++) {
496          row[numberElements]=row_[base+j];
497          element[numberElements++]=element_[base+j];
498        }
499        id_[firstAvailable_-firstDynamic_]=bestSequence;
500        //printf("best %d\n",bestSequence);
501        backward_[firstAvailable_]=bestSet;
502        model->solutionRegion()[firstAvailable_]=0.0;
503        if (!lowerColumn_&&!upperColumn_) {
504          model->setStatus(firstAvailable_,ClpSimplex::atLowerBound);
505          lowerColumn[firstAvailable_] = 0.0;
506          upperColumn[firstAvailable_] = COIN_DBL_MAX;
507        }  else {
508          DynamicStatus status = getDynamicStatus(bestSequence);
509          if (lowerColumn_) 
510            lowerColumn[firstAvailable_] = lowerColumn_[bestSequence];
511          else
512            lowerColumn[firstAvailable_] = 0.0;
513          if (upperColumn_)
514            upperColumn[firstAvailable_] = upperColumn_[bestSequence];
515          else
516            upperColumn[firstAvailable_] = COIN_DBL_MAX;
517          if (status==atLowerBound) {
518            solution[firstAvailable_]=lowerColumn[firstAvailable_];
519            model->setStatus(firstAvailable_,ClpSimplex::atLowerBound);
520          } else {
521            solution[firstAvailable_]=upperColumn[firstAvailable_];
522            model->setStatus(firstAvailable_,ClpSimplex::atUpperBound);
523          }
524        }
525        model->nonLinearCost()->setOne(firstAvailable_,solution[firstAvailable_],
526                                       lowerColumn[firstAvailable_],
527                                       upperColumn[firstAvailable_],cost_[bestSequence]);
528        bestSequence=firstAvailable_;
529        // firstAvailable_ only updated if good pivot (in updatePivot)
530        startColumn[firstAvailable_+1]=numberElements;
531        //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,bestDjMod);
532        reducedCost[bestSequence] = value;
533        gubSlackIn_=-1;
534      } else {
535        // slack - make last column
536        gubSlackIn_= bestSequence - numberRows - numberColumns;
537        bestSequence = numberColumns + 2*numberRows;
538        reducedCost[bestSequence] = bestDjMod;
539        //printf("price slack %d - gubpi %g\n",gubSlackIn_,bestDjMod);
540        model->setStatus(bestSequence,getStatus(gubSlackIn_));
541        if (getStatus(gubSlackIn_)==ClpSimplex::atUpperBound)
542          solution[bestSequence] = upper_[gubSlackIn_];
543        else
544          solution[bestSequence] = lower_[gubSlackIn_];
545        lowerColumn[bestSequence] = lower_[gubSlackIn_];
546        upperColumn[bestSequence] = upper_[gubSlackIn_];
547        model->costRegion()[bestSequence] = 0.0;
548        model->nonLinearCost()->setOne(bestSequence,solution[bestSequence],lowerColumn[bestSequence],
549                                       upperColumn[bestSequence],0.0);
550      }
551      savedBestSequence_ = bestSequence;
552      savedBestDj_ = reducedCost[savedBestSequence_];
553    }
554    // See if may be finished
555    if (!startG2&&bestSequence<0)
556      infeasibilityWeight_=model_->infeasibilityCost();
557    else if (bestSequence>=0) 
558      infeasibilityWeight_=-1.0;
559  }
560  currentWanted_=numberWanted;
561}
562// This is local to Gub to allow synchronization when status is good
563int 
564ClpGubDynamicMatrix::synchronize(ClpSimplex * model,int mode)
565{
566  int returnNumber=0;
567  switch (mode) {
568  case 0:
569    {
570#ifdef CLP_DEBUG
571      {
572        for (int i=0;i<numberSets_;i++)
573          assert(toIndex_[i]==-1);
574      }
575#endif
576      // lookup array
577      int * lookup = new int[lastDynamic_];
578      int iColumn;
579      int numberColumns = model->numberColumns();
580      double * element =  matrix_->getMutableElements();
581      int * row = matrix_->getMutableIndices();
582      CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
583      int * length = matrix_->getMutableVectorLengths();
584      double * cost = model->costRegion();
585      double * lowerColumn = model->lowerRegion();
586      double * upperColumn = model->upperRegion();
587      int * pivotVariable = model->pivotVariable();
588      CoinBigIndex numberElements=startColumn[firstDynamic_];
589      // first just do lookup and basic stuff
590      int currentNumber = firstAvailable_;
591      firstAvailable_ = firstDynamic_;
592      int numberToDo=0;
593      double objectiveChange=0.0;
594      double * solution = model->solutionRegion();
595      for (iColumn=firstDynamic_;iColumn<currentNumber;iColumn++) {
596        int iSet = backward_[iColumn];
597        if (toIndex_[iSet]<0) {
598          toIndex_[iSet]=0;
599          fromIndex_[numberToDo++]=iSet;
600        }
601        if (model->getStatus(iColumn)==ClpSimplex::basic||iColumn==keyVariable_[iSet]) {
602          lookup[iColumn]=firstAvailable_;
603          if (iColumn!=keyVariable_[iSet]) {
604            int iPivot = backToPivotRow_[iColumn];
605            backToPivotRow_[firstAvailable_]=iPivot;
606            pivotVariable[iPivot]=firstAvailable_;
607          }
608          firstAvailable_++;
609        } else {
610          int jColumn = id_[iColumn-firstDynamic_];
611          setDynamicStatus(jColumn,atLowerBound);
612          if (lowerColumn_||upperColumn_) {
613            if (model->getStatus(iColumn)==ClpSimplex::atUpperBound)
614              setDynamicStatus(jColumn,atUpperBound);
615            // treat solution as if exactly at a bound
616            double value = solution[iColumn];
617            if (fabs(value-lowerColumn[iColumn])<fabs(value-upperColumn[iColumn]))
618              value = lowerColumn[iColumn];
619            else
620              value = upperColumn[iColumn];
621            objectiveChange += cost[iColumn]*value;
622            // redo lower and upper on sets
623            double shift=value;
624            if (lowerSet_[iSet]>-1.0e20)
625              lower_[iSet] = lowerSet_[iSet]-shift;
626            if (upperSet_[iSet]<1.0e20)
627              upper_[iSet] = upperSet_[iSet]-shift;
628          }
629          lookup[iColumn]=-1;
630        }
631      }
632      model->setObjectiveOffset(model->objectiveOffset()+objectiveChange);
633      firstAvailable_ = firstDynamic_;
634      for (iColumn=firstDynamic_;iColumn<currentNumber;iColumn++) {
635        if (lookup[iColumn]>=0) {
636          // move
637          int jColumn = id_[iColumn-firstDynamic_];
638          id_[firstAvailable_-firstDynamic_] = jColumn;
639          int numberThis = startColumn_[jColumn+1]-startColumn_[jColumn];
640          length[firstAvailable_]=numberThis;
641          cost[firstAvailable_]=cost[iColumn];
642          lowerColumn[firstAvailable_]=lowerColumn[iColumn];
643          upperColumn[firstAvailable_]=upperColumn[iColumn];
644          double originalLower = lowerColumn_ ? lowerColumn_[jColumn] : 0.0;
645          double originalUpper = upperColumn_ ? upperColumn_[jColumn] : COIN_DBL_MAX;
646          if (originalUpper>1.0e30)
647            originalUpper = COIN_DBL_MAX;
648          model->nonLinearCost()->setOne(firstAvailable_,solution[iColumn],
649                                         originalLower,originalUpper,
650                                         cost_[jColumn]);
651          CoinBigIndex base = startColumn_[jColumn];
652          for (int j=0;j<numberThis;j++) {
653            row[numberElements]=row_[base+j];
654            element[numberElements++]=element_[base+j];
655          }
656          model->setStatus(firstAvailable_,model->getStatus(iColumn));
657          backward_[firstAvailable_] = backward_[iColumn];
658          solution[firstAvailable_] = solution[iColumn];
659          firstAvailable_++;
660          startColumn[firstAvailable_]=numberElements;
661        }
662      }
663      // clean up next_
664      int * temp = new int [firstAvailable_];
665      for (int jSet=0;jSet<numberToDo;jSet++) {
666        int iSet = fromIndex_[jSet];
667        toIndex_[iSet]=-1;
668        int last = keyVariable_[iSet];
669        int j = next_[last];
670        bool setTemp=true;
671        if (last<lastDynamic_) {
672          last = lookup[last];
673          assert (last>=0);
674          keyVariable_[iSet]=last;
675        } else if (j>=0) {
676          int newJ = lookup[j];
677          assert (newJ>=0);
678          j=next_[j];
679          next_[last]=newJ;
680          last = newJ;
681        } else {
682          next_[last] = -(iSet+numberColumns+1);
683          setTemp=false;
684        }
685        while (j>=0) {
686          int newJ = lookup[j];
687          assert (newJ>=0);
688          temp[last]=newJ;
689          last = newJ;
690          j=next_[j];
691        }
692        if (setTemp)
693          temp[last]=-(keyVariable_[iSet]+1);
694        if (lowerSet_) {
695          // we only need to get lower_ and upper_ correct
696          double shift=0.0;
697          for (int j=fullStart_[iSet];j<fullStart_[iSet+1];j++) 
698            if (getDynamicStatus(j)==atUpperBound) 
699              shift += upperColumn_[j];
700            else if (getDynamicStatus(j)==atLowerBound&&lowerColumn_)
701              shift += lowerColumn_[j];
702          if (lowerSet_[iSet]>-1.0e20)
703            lower_[iSet] = lowerSet_[iSet]-shift;
704          if (upperSet_[iSet]<1.0e20)
705            upper_[iSet] = upperSet_[iSet]-shift;
706        }
707      }
708      // move to next_
709      memcpy(next_+firstDynamic_,temp+firstDynamic_,(firstAvailable_-firstDynamic_)*sizeof(int));
710      // if odd iterations may be one out so adjust currentNumber
711      currentNumber=CoinMin(currentNumber+1,lastDynamic_);
712      // zero solution
713      CoinZeroN(solution+firstAvailable_,currentNumber-firstAvailable_);
714      // zero cost
715      CoinZeroN(cost+firstAvailable_,currentNumber-firstAvailable_);
716      // zero lengths
717      CoinZeroN(length+firstAvailable_,currentNumber-firstAvailable_);
718      for ( iColumn=firstAvailable_;iColumn<currentNumber;iColumn++) {
719        model->nonLinearCost()->setOne(iColumn,0.0,0.0,COIN_DBL_MAX,0.0);
720        model->setStatus(iColumn,ClpSimplex::atLowerBound);
721        backward_[iColumn]=-1;
722      }
723      delete [] lookup;
724      delete [] temp;
725      // make sure fromIndex clean
726      fromIndex_[0]=-1;
727      //#define CLP_DEBUG
728#ifdef CLP_DEBUG
729      // debug
730      {
731        int i;
732        int numberRows=model->numberRows();
733        char * xxxx = new char[numberColumns];
734        memset(xxxx,0,numberColumns);
735        for (i=0;i<numberRows;i++) {
736          int iPivot = pivotVariable[i];
737          assert (model->getStatus(iPivot)==ClpSimplex::basic);
738          if (iPivot<numberColumns&&backward_[iPivot]>=0)
739            xxxx[iPivot]=1;
740        }
741        for (i=0;i<numberSets_;i++) {
742          int key=keyVariable_[i];
743          int iColumn =next_[key];
744          int k=0;
745          while(iColumn>=0) {
746            k++;
747            assert (k<100);
748            assert (backward_[iColumn]==i);
749            iColumn=next_[iColumn];
750          }
751          int stop = -(key+1);
752          while (iColumn!=stop) {
753            assert (iColumn<0);
754            iColumn = -iColumn-1;
755            k++;
756            assert (k<100);
757            assert (backward_[iColumn]==i);
758            iColumn=next_[iColumn];
759          }
760          iColumn =next_[key];
761          while (iColumn>=0) {
762            assert (xxxx[iColumn]);
763            xxxx[iColumn]=0;
764            iColumn=next_[iColumn];
765          }
766        }
767        for (i=0;i<numberColumns;i++) {
768          if (i<numberColumns&&backward_[i]>=0) {
769            assert (!xxxx[i]||i==keyVariable_[backward_[i]]);
770          }
771        }
772        delete [] xxxx;
773      }
774      {
775        for (int i=0;i<numberSets_;i++)
776          assert(toIndex_[i]==-1);
777      }
778#endif
779      savedFirstAvailable_ = firstAvailable_;
780    }
781    break;
782    // flag a variable
783  case 1:
784    {
785      // id will be sitting at firstAvailable
786      int sequence = id_[firstAvailable_-firstDynamic_];
787      assert (!flagged(sequence));
788      setFlagged(sequence);
789      model->clearFlagged(firstAvailable_);
790    }
791    break;
792    // unflag all variables
793  case 2:
794    {
795      for (int i=0;i<numberGubColumns_;i++) {
796        if (flagged(i)) {
797          unsetFlagged(i);
798          returnNumber++;
799        }
800      }
801    }
802    break;
803    //  just reset costs and bounds (primal)
804  case 3:
805    {
806      double * cost = model->costRegion();
807      double * solution = model->solutionRegion();
808      double * lowerColumn = model->columnLower();
809      double * upperColumn = model->columnUpper();
810      for (int i=firstDynamic_;i<firstAvailable_;i++) {
811        int jColumn = id_[i-firstDynamic_];
812        cost[i]=cost_[jColumn];
813        if (!lowerColumn_&&!upperColumn_) {
814          lowerColumn[i] = 0.0;
815          upperColumn[i] = COIN_DBL_MAX;
816        }  else {
817          if (lowerColumn_) 
818            lowerColumn[i] = lowerColumn_[jColumn];
819          else
820            lowerColumn[i] = 0.0;
821          if (upperColumn_)
822            upperColumn[i] = upperColumn_[jColumn];
823          else
824            upperColumn[i] = COIN_DBL_MAX;
825        }
826        if (model->nonLinearCost())
827          model->nonLinearCost()->setOne(i,solution[i],
828                                         lowerColumn[i],
829                                         upperColumn[i],cost_[jColumn]);
830      }
831      if (!model->numberIterations()&&rhsOffset_) {
832        lastRefresh_ = - refreshFrequency_; // force refresh
833      }
834    }
835    break;
836    // and get statistics for column generation
837  case 4:
838    {
839      // In theory we should subtract out ones we have done but ....
840      // If key slack then dual 0.0
841      // If not then slack could be dual infeasible
842      // dj for key is zero so that defines dual on set
843      int i;
844      int numberColumns = model->numberColumns();
845      double * dual = model->dualRowSolution();
846      double infeasibilityCost = model->infeasibilityCost();
847      double dualTolerance = model->dualTolerance();
848      double relaxedTolerance=dualTolerance;
849      // we can't really trust infeasibilities if there is dual error
850      double error = CoinMin(1.0e-2,model->largestDualError());
851      // allow tolerance at least slightly bigger than standard
852      relaxedTolerance = relaxedTolerance +  error;
853      // but we will be using difference
854      relaxedTolerance -= dualTolerance;
855      double objectiveOffset = 0.0;
856      for (i=0;i<numberSets_;i++) {
857        int kColumn = keyVariable_[i];
858        double value=0.0;
859        if (kColumn<numberColumns) {
860          kColumn = id_[kColumn-firstDynamic_];
861          // dj without set
862          value = cost_[kColumn];
863          for (CoinBigIndex j=startColumn_[kColumn];
864               j<startColumn_[kColumn+1];j++) {
865            int iRow = row_[j];
866            value -= dual[iRow]*element_[j];
867          }
868          double infeasibility = 0.0;
869          if (getStatus(i)==ClpSimplex::atLowerBound) {
870            if (-value>dualTolerance) 
871              infeasibility=-value-dualTolerance;
872          } else if (getStatus(i)==ClpSimplex::atUpperBound) {
873            if (value>dualTolerance) 
874              infeasibility=value-dualTolerance;
875          }
876          if (infeasibility>0.0) {
877            sumDualInfeasibilities_ += infeasibility;
878            if (infeasibility>relaxedTolerance) 
879              sumOfRelaxedDualInfeasibilities_ += infeasibility;
880            numberDualInfeasibilities_ ++;
881          }
882        } else {
883          // slack key - may not be feasible
884          assert (getStatus(i)==ClpSimplex::basic);
885          // negative as -1.0 for slack
886          value=-weight(i)*infeasibilityCost;
887        }
888        // Now subtract out from all
889        for (CoinBigIndex k= fullStart_[i];k<fullStart_[i+1];k++) {
890          if (getDynamicStatus(k)!=inSmall) {
891            double djValue = cost_[k]-value;
892            for (CoinBigIndex j=startColumn_[k];
893                 j<startColumn_[k+1];j++) {
894              int iRow = row_[j];
895              djValue -= dual[iRow]*element_[j];
896            }
897            double infeasibility=0.0;
898            double shift=0.0;
899            if (getDynamicStatus(k)==atLowerBound) {
900              if (lowerColumn_) 
901                shift= lowerColumn_[k];
902              if (djValue<-dualTolerance) 
903                infeasibility=-djValue-dualTolerance;
904            } else {
905              // at upper bound
906              shift= upperColumn_[k];
907              if (djValue>dualTolerance) 
908                infeasibility=djValue-dualTolerance;
909            }
910            objectiveOffset += shift*cost_[k];
911            if (infeasibility>0.0) {
912              sumDualInfeasibilities_ += infeasibility;
913              if (infeasibility>relaxedTolerance) 
914                sumOfRelaxedDualInfeasibilities_ += infeasibility;
915              numberDualInfeasibilities_ ++;
916            }
917          }
918        }
919      }
920      model->setObjectiveOffset(objectiveOffset_-objectiveOffset);
921    }
922    break;
923    // see if time to re-factorize
924  case 5:
925    {
926      if (firstAvailable_>numberSets_+model->numberRows()+ model->factorizationFrequency())
927        returnNumber=4;
928    }
929    break;
930    // return 1 if there may be changing bounds on variable (column generation)
931  case 6:
932    {
933      returnNumber = (lowerColumn_!=NULL||upperColumn_!=NULL) ? 1 : 0;
934#if 0
935      if (!returnNumber) {
936        // may be gub slacks
937        for (int i=0;i<numberSets_;i++) {
938          if (upper_[i]>lower_[i]) {
939            returnNumber=1;
940            break;
941          }
942        }
943      }
944#endif
945    }
946    break;
947    // restore firstAvailable_
948  case 7:
949    {
950      int iColumn;
951      int * length = matrix_->getMutableVectorLengths();
952      double * cost = model->costRegion();
953      double * solution = model->solutionRegion();
954      int currentNumber=firstAvailable_;
955      firstAvailable_ = savedFirstAvailable_;
956      // zero solution
957      CoinZeroN(solution+firstAvailable_,currentNumber-firstAvailable_);
958      // zero cost
959      CoinZeroN(cost+firstAvailable_,currentNumber-firstAvailable_);
960      // zero lengths
961      CoinZeroN(length+firstAvailable_,currentNumber-firstAvailable_);
962      for ( iColumn=firstAvailable_;iColumn<currentNumber;iColumn++) {
963        model->nonLinearCost()->setOne(iColumn,0.0,0.0,COIN_DBL_MAX,0.0);
964        model->setStatus(iColumn,ClpSimplex::atLowerBound);
965        backward_[iColumn]=-1;
966      }
967    }
968    break;
969    // make sure set is clean
970  case 8:
971    {
972      int sequenceIn = model->sequenceIn();
973      if (sequenceIn<model->numberColumns()) {
974        int iSet = backward_[sequenceIn];
975        if (iSet>=0&&lowerSet_) {
976          // we only need to get lower_ and upper_ correct
977          double shift=0.0;
978          for (int j=fullStart_[iSet];j<fullStart_[iSet+1];j++) 
979            if (getDynamicStatus(j)==atUpperBound) 
980              shift += upperColumn_[j];
981            else if (getDynamicStatus(j)==atLowerBound&&lowerColumn_)
982              shift += lowerColumn_[j];
983          if (lowerSet_[iSet]>-1.0e20)
984            lower_[iSet] = lowerSet_[iSet]-shift;
985          if (upperSet_[iSet]<1.0e20)
986            upper_[iSet] = upperSet_[iSet]-shift;
987        }
988        if (sequenceIn==firstAvailable_) {
989          // not really in small problem
990          int iBig=id_[sequenceIn-firstDynamic_];
991          if (model->getStatus(sequenceIn)==ClpSimplex::atLowerBound)
992            setDynamicStatus(iBig,atLowerBound);
993          else
994            setDynamicStatus(iBig,atUpperBound);
995        }
996      }
997    }
998    break;
999    // adjust lower,upper
1000  case 9:
1001    {
1002      int sequenceIn = model->sequenceIn();
1003      if (sequenceIn>=firstDynamic_&&sequenceIn<lastDynamic_&&lowerSet_) {
1004        int iSet = backward_[sequenceIn];
1005        assert (iSet>=0);
1006        int inBig = id_[sequenceIn-firstDynamic_];
1007        const double * solution = model->solutionRegion();
1008        setDynamicStatus(inBig,inSmall);
1009        if (lowerSet_[iSet]>-1.0e20)
1010          lower_[iSet] += solution[sequenceIn];
1011        if (upperSet_[iSet]<1.0e20)
1012          upper_[iSet] += solution[sequenceIn];
1013        model->setObjectiveOffset(model->objectiveOffset()-
1014                                  solution[sequenceIn]*cost_[inBig]);
1015      }
1016    }
1017  }
1018  return returnNumber;
1019}
1020// Add a new variable to a set
1021void 
1022ClpGubDynamicMatrix::insertNonBasic(int sequence, int iSet)
1023{
1024  int last = keyVariable_[iSet];
1025  int j=next_[last];
1026  while (j>=0) {
1027    last=j;
1028    j=next_[j];
1029  }
1030  next_[last]=-(sequence+1);
1031  next_[sequence]=j;
1032}
1033// Sets up an effective RHS and does gub crash if needed
1034void 
1035ClpGubDynamicMatrix::useEffectiveRhs(ClpSimplex * model, bool cheapest)
1036{
1037  // Do basis - cheapest or slack if feasible (unless cheapest set)
1038  int longestSet=0;
1039  int iSet;
1040  for (iSet=0;iSet<numberSets_;iSet++) 
1041    longestSet = CoinMax(longestSet,fullStart_[iSet+1]-fullStart_[iSet]);
1042   
1043  double * upper = new double[longestSet+1];
1044  double * cost = new double[longestSet+1];
1045  double * lower = new double[longestSet+1];
1046  double * solution = new double[longestSet+1];
1047  assert (!next_);
1048  delete [] next_;
1049  int numberColumns = model->numberColumns();
1050  next_ = new int[numberColumns+numberSets_+CoinMax(2*longestSet,lastDynamic_-firstDynamic_)];
1051  char * mark = new char[numberColumns];
1052  memset(mark,0,numberColumns);
1053  for (int iColumn=0;iColumn<numberColumns;iColumn++) 
1054    next_[iColumn]=INT_MAX;
1055  int i;
1056  int * keys = new int[numberSets_];
1057  int * back = new int[numberGubColumns_];
1058  CoinFillN(back,numberGubColumns_,-1);
1059  for (i=0;i<numberSets_;i++) 
1060    keys[i]=INT_MAX;
1061  delete [] dynamicStatus_;
1062  dynamicStatus_ = new unsigned char [numberGubColumns_];
1063  memset(dynamicStatus_,0,numberGubColumns_); // for clarity
1064  for (i=0;i<numberGubColumns_;i++)
1065    setDynamicStatus(i,atLowerBound);
1066  // set up chains
1067  for (i=firstDynamic_;i<lastDynamic_;i++){
1068    if (id_[i-firstDynamic_]>=0) {
1069      if (model->getStatus(i)==ClpSimplex::basic) 
1070        mark[i]=1;
1071      int iSet = backward_[i];
1072      assert (iSet>=0);
1073      int iNext = keys[iSet];
1074      next_[i]=iNext;
1075      keys[iSet]=i;
1076      back[id_[i-firstDynamic_]]=i;
1077    } else {
1078      model->setStatus(i,ClpSimplex::atLowerBound); 
1079      backward_[i]=-1;
1080    }
1081  }
1082  double * columnSolution = model->solutionRegion();
1083  int numberRows = getNumRows();
1084  toIndex_ = new int[numberSets_];
1085  for (iSet=0;iSet<numberSets_;iSet++) 
1086    toIndex_[iSet]=-1;
1087  fromIndex_ = new int [numberRows+numberSets_];
1088  double tolerance = model->primalTolerance();
1089  double * element =  matrix_->getMutableElements();
1090  int * row = matrix_->getMutableIndices();
1091  CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
1092  int * length = matrix_->getMutableVectorLengths();
1093  double objectiveOffset = 0.0;
1094  for (iSet=0;iSet<numberSets_;iSet++) {
1095    int j;
1096    int numberBasic=0;
1097    int iBasic=-1;
1098    int iStart = fullStart_[iSet];
1099    int iEnd=fullStart_[iSet+1];
1100    // find one with smallest length
1101    int smallest=numberRows+1;
1102    double value=0.0;
1103    j=keys[iSet];
1104    while (j!=INT_MAX) {
1105      if (model->getStatus(j)== ClpSimplex::basic) {
1106        if (length[j]<smallest) {
1107          smallest=length[j];
1108          iBasic=j;
1109        }
1110        numberBasic++;
1111      }
1112      value += columnSolution[j];
1113      j = next_[j];
1114    }
1115    bool done=false;
1116    if (numberBasic>1||(numberBasic==1&&getStatus(iSet)==ClpSimplex::basic)) {
1117      if (getStatus(iSet)==ClpSimplex::basic) 
1118        iBasic = iSet+numberColumns;// slack key - use
1119      done=true;
1120    } else if (numberBasic==1) {
1121      // see if can be key
1122      double thisSolution = columnSolution[iBasic];
1123      if (thisSolution<0.0) {
1124        value -= thisSolution;
1125        thisSolution = 0.0;
1126        columnSolution[iBasic]=thisSolution;
1127      }
1128      // try setting slack to a bound
1129      assert (upper_[iSet]<1.0e20||lower_[iSet]>-1.0e20);
1130      double cost1 = COIN_DBL_MAX;
1131      int whichBound=-1;
1132      if (upper_[iSet]<1.0e20) {
1133        // try slack at ub
1134        double newBasic = thisSolution +upper_[iSet]-value;
1135        if (newBasic>=-tolerance) {
1136          // can go
1137          whichBound=1;
1138          cost1 = newBasic*cost_[iBasic];
1139          // But if exact then may be good solution
1140          if (fabs(upper_[iSet]-value)<tolerance)
1141            cost1=-COIN_DBL_MAX;
1142        }
1143      }
1144      if (lower_[iSet]>-1.0e20) {
1145        // try slack at lb
1146        double newBasic = thisSolution +lower_[iSet]-value;
1147        if (newBasic>=-tolerance) {
1148          // can go but is it cheaper
1149          double cost2 = newBasic*cost_[iBasic];
1150          // But if exact then may be good solution
1151          if (fabs(lower_[iSet]-value)<tolerance)
1152            cost2=-COIN_DBL_MAX;
1153          if (cost2<cost1)
1154            whichBound=0;
1155        }
1156      }
1157      if (whichBound!=-1) {
1158        // key
1159        done=true;
1160        if (whichBound) {
1161          // slack to upper
1162          columnSolution[iBasic]=thisSolution + upper_[iSet]-value;
1163          setStatus(iSet,ClpSimplex::atUpperBound);
1164        } else {
1165          // slack to lower
1166          columnSolution[iBasic]=thisSolution + lower_[iSet]-value;
1167          setStatus(iSet,ClpSimplex::atLowerBound);
1168        }
1169      }
1170    }
1171    if (!done) {
1172      if (!cheapest) {
1173        // see if slack can be key
1174        if (value>=lower_[iSet]-tolerance&&value<=upper_[iSet]+tolerance) {
1175          done=true;
1176          setStatus(iSet,ClpSimplex::basic);
1177          iBasic=iSet+numberColumns;
1178        }
1179      }
1180      if (!done) {
1181        // set non basic if there was one
1182        if (iBasic>=0)
1183          model->setStatus(iBasic,ClpSimplex::atLowerBound);
1184        // find cheapest
1185        int numberInSet = iEnd-iStart;
1186        if (!lowerColumn_) {
1187          CoinZeroN(lower,numberInSet);
1188        } else {
1189          for (int j=0;j<numberInSet;j++)
1190            lower[j]= lowerColumn_[j+iStart];
1191        }
1192        if (!upperColumn_) {
1193          CoinFillN(upper,numberInSet,COIN_DBL_MAX);
1194        } else {
1195          for (int j=0;j<numberInSet;j++)
1196            upper[j]= upperColumn_[j+iStart];
1197        }
1198        CoinFillN(solution,numberInSet,0.0);
1199        // and slack
1200        iBasic=numberInSet;
1201        solution[iBasic]=-value;
1202        lower[iBasic]=-upper_[iSet];
1203        upper[iBasic]=-lower_[iSet];
1204        int kphase;
1205        if (value>=lower_[iSet]-tolerance&&value<=upper_[iSet]+tolerance) {
1206          // feasible
1207          kphase=1;
1208          cost[iBasic]=0.0;
1209          for (int j=0;j<numberInSet;j++)
1210            cost[j]= cost_[j+iStart];
1211        } else {
1212          // infeasible
1213          kphase=0;
1214          // remember bounds are flipped so opposite to natural
1215          if (value<lower_[iSet]-tolerance)
1216            cost[iBasic]=1.0;
1217          else
1218            cost[iBasic]=-1.0;
1219          CoinZeroN(cost,numberInSet);
1220        }
1221        double dualTolerance =model->dualTolerance();
1222        for (int iphase =kphase;iphase<2;iphase++) {
1223          if (iphase) {
1224            cost[numberInSet]=0.0;
1225            for (int j=0;j<numberInSet;j++)
1226              cost[j]= cost_[j+iStart];
1227          }
1228          // now do one row lp
1229          bool improve=true;
1230          while (improve) {
1231            improve=false;
1232            double dual = cost[iBasic];
1233            int chosen =-1;
1234            double best=dualTolerance;
1235            int way=0;
1236            for (int i=0;i<=numberInSet;i++) {
1237              double dj = cost[i]-dual;
1238              double improvement =0.0;
1239              double distance=0.0;
1240              if (iphase||i<numberInSet)
1241                assert (solution[i]>=lower[i]&&solution[i]<=upper[i]);
1242              if (dj>dualTolerance)
1243                improvement = dj*(solution[i]-lower[i]);
1244              else if (dj<-dualTolerance)
1245                improvement = dj*(solution[i]-upper[i]);
1246              if (improvement>best) {
1247                best=improvement;
1248                chosen=i;
1249                if (dj<0.0) {
1250                  way = 1;
1251                  distance = upper[i]-solution[i];
1252                } else {
1253                  way = -1;
1254                  distance = solution[i]-lower[i];
1255                }
1256              }
1257            }
1258            if (chosen>=0) {
1259              improve=true;
1260              // now see how far
1261              if (way>0) {
1262                // incoming increasing so basic decreasing
1263                // if phase 0 then go to nearest bound
1264                double distance=upper[chosen]-solution[chosen];
1265                double basicDistance;
1266                if (!iphase) {
1267                  assert (iBasic==numberInSet);
1268                  assert (solution[iBasic]>upper[iBasic]);
1269                  basicDistance = solution[iBasic]-upper[iBasic];
1270                } else {
1271                  basicDistance = solution[iBasic]-lower[iBasic];
1272                }
1273                // need extra coding for unbounded
1274                assert (CoinMin(distance,basicDistance)<1.0e20);
1275                if (distance>basicDistance) {
1276                  // incoming becomes basic
1277                  solution[chosen] += basicDistance;
1278                  if (!iphase) 
1279                    solution[iBasic]=upper[iBasic];
1280                  else 
1281                    solution[iBasic]=lower[iBasic];
1282                  iBasic = chosen;
1283                } else {
1284                  // flip
1285                  solution[chosen]=upper[chosen];
1286                  solution[iBasic] -= distance;
1287                }
1288              } else {
1289                // incoming decreasing so basic increasing
1290                // if phase 0 then go to nearest bound
1291                double distance=solution[chosen]-lower[chosen];
1292                double basicDistance;
1293                if (!iphase) {
1294                  assert (iBasic==numberInSet);
1295                  assert (solution[iBasic]<lower[iBasic]);
1296                  basicDistance = lower[iBasic]-solution[iBasic];
1297                } else {
1298                  basicDistance = upper[iBasic]-solution[iBasic];
1299                }
1300                // need extra coding for unbounded - for now just exit
1301                if (CoinMin(distance,basicDistance)>1.0e20) {
1302                  printf("unbounded on set %d\n",iSet);
1303                  iphase=1;
1304                  iBasic=numberInSet;
1305                  break;
1306                }
1307                if (distance>basicDistance) {
1308                  // incoming becomes basic
1309                  solution[chosen] -= basicDistance;
1310                  if (!iphase) 
1311                    solution[iBasic]=lower[iBasic];
1312                  else 
1313                    solution[iBasic]=upper[iBasic];
1314                  iBasic = chosen;
1315                } else {
1316                  // flip
1317                  solution[chosen]=lower[chosen];
1318                  solution[iBasic] += distance;
1319                }
1320              }
1321              if (!iphase) {
1322                if(iBasic<numberInSet)
1323                  break; // feasible
1324                else if (solution[iBasic]>=lower[iBasic]&&
1325                         solution[iBasic]<=upper[iBasic])
1326                  break; // feasible (on flip)
1327              }
1328            }
1329          }
1330        }
1331        // do solution i.e. bounds
1332        if (lowerColumn_||upperColumn_) {
1333          for (int j=0;j<numberInSet;j++) {
1334            if (j!=iBasic) {
1335              objectiveOffset += solution[j]*cost[j];
1336              if (lowerColumn_&&upperColumn_) {
1337                if (fabs(solution[j]-lowerColumn_[j+iStart])>
1338                    fabs(solution[j]-upperColumn_[j+iStart]))
1339                  setDynamicStatus(j+iStart,atUpperBound);
1340              } else if (upperColumn_&&solution[j]>0.0) {
1341                setDynamicStatus(j+iStart,atUpperBound);
1342              } else {
1343                setDynamicStatus(j+iStart,atLowerBound);
1344              }
1345            }
1346          }
1347        }
1348        // convert iBasic back and do bounds
1349        if (iBasic==numberInSet) {
1350          // slack basic
1351          setStatus(iSet,ClpSimplex::basic);
1352          iBasic=iSet+numberColumns;
1353        } else {
1354          iBasic += fullStart_[iSet];
1355          if (back[iBasic]>=0) {
1356            // exists
1357            iBasic = back[iBasic];
1358          } else {
1359            // create
1360            CoinBigIndex numberElements = startColumn[firstAvailable_];
1361            int numberThis = startColumn_[iBasic+1]-startColumn_[iBasic];
1362            if (numberElements+numberThis>numberElements_) {
1363              // need to redo
1364              numberElements_ = CoinMax(3*numberElements_/2,numberElements+numberThis);
1365              matrix_->reserve(numberColumns,numberElements_);
1366              element =  matrix_->getMutableElements();
1367              row = matrix_->getMutableIndices();
1368              // these probably okay but be safe
1369              startColumn = matrix_->getMutableVectorStarts();
1370              length = matrix_->getMutableVectorLengths();
1371            }
1372            length[firstAvailable_]=numberThis;
1373            model->costRegion()[firstAvailable_]=cost_[iBasic];
1374            if (lowerColumn_)
1375              model->lowerRegion()[firstAvailable_] = lowerColumn_[iBasic];
1376            else
1377              model->lowerRegion()[firstAvailable_] = 0.0;
1378            if (upperColumn_)
1379              model->upperRegion()[firstAvailable_] = upperColumn_[iBasic];
1380            else
1381              model->upperRegion()[firstAvailable_] = COIN_DBL_MAX;
1382            columnSolution[firstAvailable_]=solution[iBasic-fullStart_[iSet]];
1383            CoinBigIndex base = startColumn_[iBasic];
1384            for (int j=0;j<numberThis;j++) {
1385              row[numberElements]=row_[base+j];
1386              element[numberElements++]=element_[base+j];
1387            }
1388            // already set startColumn[firstAvailable_]=numberElements;
1389            id_[firstAvailable_-firstDynamic_]=iBasic;
1390            setDynamicStatus(iBasic,inSmall);
1391            backward_[firstAvailable_]=iSet;
1392            iBasic=firstAvailable_;
1393            firstAvailable_++;
1394            startColumn[firstAvailable_]=numberElements;
1395          }
1396          model->setStatus(iBasic,ClpSimplex::basic);
1397          // remember bounds flipped
1398          if (upper[numberInSet]==lower[numberInSet]) 
1399            setStatus(iSet,ClpSimplex::isFixed);
1400          else if (solution[numberInSet]==upper[numberInSet])
1401            setStatus(iSet,ClpSimplex::atLowerBound);
1402          else if (solution[numberInSet]==lower[numberInSet])
1403            setStatus(iSet,ClpSimplex::atUpperBound);
1404          else 
1405            abort();
1406        }
1407        for (j=iStart;j<iEnd;j++) {
1408          int iBack = back[j];
1409          if (iBack>=0) {
1410            if (model->getStatus(iBack)!=ClpSimplex::basic) {
1411              int inSet=j-iStart;
1412              columnSolution[iBack]=solution[inSet];
1413              if (upper[inSet]==lower[inSet]) 
1414                model->setStatus(iBack,ClpSimplex::isFixed);
1415              else if (solution[inSet]==upper[inSet])
1416                model->setStatus(iBack,ClpSimplex::atUpperBound);
1417              else if (solution[inSet]==lower[inSet])
1418                model->setStatus(iBack,ClpSimplex::atLowerBound);
1419            }
1420          }
1421        }
1422      }
1423    } 
1424    keyVariable_[iSet]=iBasic;
1425  }
1426  model->setObjectiveOffset(objectiveOffset_-objectiveOffset);
1427  delete [] lower;
1428  delete [] solution;
1429  delete [] upper;
1430  delete [] cost;
1431  // make sure matrix is in good shape
1432  matrix_->orderMatrix();
1433  // create effective rhs
1434  delete [] rhsOffset_;
1435  rhsOffset_ = new double[numberRows];
1436  // and redo chains
1437  memset(mark,0,numberColumns);
1438  for (int iColumnX=0;iColumnX<firstAvailable_;iColumnX++) 
1439    next_[iColumnX]=INT_MAX;
1440  for (i=0;i<numberSets_;i++) {
1441    keys[i]=INT_MAX;
1442    int iKey = keyVariable_[i];
1443    if (iKey<numberColumns)
1444      model->setStatus(iKey,ClpSimplex::basic);
1445  }
1446  // set up chains
1447  for (i=0;i<firstAvailable_;i++){
1448    if (model->getStatus(i)==ClpSimplex::basic) 
1449      mark[i]=1;
1450    int iSet = backward_[i];
1451    if (iSet>=0) {
1452      int iNext = keys[iSet];
1453      next_[i]=iNext;
1454      keys[iSet]=i;
1455    }
1456  }
1457  for (i=0;i<numberSets_;i++) {
1458    if (keys[i]!=INT_MAX) {
1459      // something in set
1460      int j;
1461      if (getStatus(i)!=ClpSimplex::basic) {
1462        // make sure fixed if it is
1463        if (upper_[i]==lower_[i])
1464          setStatus(i,ClpSimplex::isFixed);
1465        // slack not key - choose one with smallest length
1466        int smallest=numberRows+1;
1467        int key=-1;
1468        j = keys[i];
1469        while (1) {
1470          if (mark[j]&&length[j]<smallest) {
1471            key=j;
1472            smallest=length[j];
1473          }
1474          if (next_[j]!=INT_MAX) {
1475            j = next_[j];
1476          } else {
1477            // correct end
1478            next_[j]=-(keys[i]+1);
1479            break;
1480          }
1481        }
1482        if (key>=0) {
1483          keyVariable_[i]=key;
1484        } else {
1485          // nothing basic - make slack key
1486          //((ClpGubMatrix *)this)->setStatus(i,ClpSimplex::basic);
1487          // fudge to avoid const problem
1488          status_[i]=1;
1489        }
1490      } else {
1491        // slack key
1492        keyVariable_[i]=numberColumns+i;
1493        int j;
1494        double sol=0.0;
1495        j = keys[i];
1496        while (1) {
1497          sol += columnSolution[j];
1498          if (next_[j]!=INT_MAX) {
1499            j = next_[j];
1500          } else {
1501            // correct end
1502            next_[j]=-(keys[i]+1);
1503            break;
1504          }
1505        }
1506        if (sol>upper_[i]+tolerance) {
1507          setAbove(i);
1508        } else if (sol<lower_[i]-tolerance) {
1509          setBelow(i);
1510        } else {
1511          setFeasible(i);
1512        }
1513      }
1514      // Create next_
1515      int key = keyVariable_[i];
1516      redoSet(model,key,keys[i],i);
1517    } else {
1518      // nothing in set!
1519      next_[i+numberColumns]=-(i+numberColumns+1);
1520      keyVariable_[i]=numberColumns+i;
1521      double sol=0.0;
1522      if (sol>upper_[i]+tolerance) {
1523        setAbove(i);
1524      } else if (sol<lower_[i]-tolerance) {
1525        setBelow(i);
1526      } else {
1527        setFeasible(i);
1528      }
1529    }
1530  }
1531  delete [] keys;
1532  delete [] mark;
1533  delete [] back;
1534  rhsOffset(model,true);
1535}
1536/* Returns effective RHS if it is being used.  This is used for long problems
1537   or big gub or anywhere where going through full columns is
1538   expensive.  This may re-compute */
1539double * 
1540ClpGubDynamicMatrix::rhsOffset(ClpSimplex * model,bool forceRefresh,
1541                      bool check)
1542{
1543  //forceRefresh=true;
1544  //check=false;
1545#ifdef CLP_DEBUG
1546  double * saveE=NULL;
1547  if (rhsOffset_&&check) {
1548    int numberRows = model->numberRows();
1549    saveE = new double[numberRows];
1550  }
1551#endif
1552  if (rhsOffset_) {
1553#ifdef CLP_DEBUG
1554    if (check) {
1555      // no need - but check anyway
1556      int numberRows = model->numberRows();
1557      double * rhs = new double[numberRows];
1558      int numberColumns = model->numberColumns();
1559      int iRow;
1560      CoinZeroN(rhs,numberRows);
1561      // do ones at bounds before gub
1562      const double * smallSolution = model->solutionRegion();
1563      const double * element =matrix_->getElements();
1564      const int * row = matrix_->getIndices();
1565      const CoinBigIndex * startColumn = matrix_->getVectorStarts();
1566      const int * length = matrix_->getVectorLengths();
1567      int iColumn;
1568      for (iColumn=0;iColumn<firstDynamic_;iColumn++) {
1569        if (model->getStatus(iColumn)!=ClpSimplex::basic) {
1570          double value = smallSolution[iColumn];
1571          for (CoinBigIndex j=startColumn[iColumn];
1572               j<startColumn[iColumn]+length[iColumn];j++) {
1573            int jRow=row[j];
1574            rhs[jRow] -= value*element[j];
1575          }
1576        }
1577      }
1578      if (lowerColumn_||upperColumn_) {
1579        double * solution = new double [numberGubColumns_];
1580        for (iColumn=0;iColumn<numberGubColumns_;iColumn++) {
1581          double value=0.0;
1582          if(getDynamicStatus(iColumn)==atUpperBound)
1583            value = upperColumn_[iColumn];
1584          else if (lowerColumn_)
1585            value = lowerColumn_[iColumn];
1586          solution[iColumn]=value;
1587        }
1588        // ones at bounds in small and gub
1589        for (iColumn=firstDynamic_;iColumn<firstAvailable_;iColumn++) {
1590          int jFull = id_[iColumn-firstDynamic_];
1591          solution[jFull]=smallSolution[iColumn];
1592        }
1593        // zero all basic in small model
1594        int * pivotVariable = model->pivotVariable();
1595        for (iRow=0;iRow<numberRows;iRow++) {
1596          int iColumn = pivotVariable[iRow];
1597          if (iColumn>=firstDynamic_&&iColumn<lastDynamic_) {
1598            int iSequence = id_[iColumn-firstDynamic_];
1599            solution[iSequence]=0.0;
1600          }
1601        }
1602        // and now compute value to use for key
1603        ClpSimplex::Status iStatus;
1604        for (int iSet=0;iSet<numberSets_;iSet++) {
1605          iColumn = keyVariable_[iSet];
1606          if (iColumn<numberColumns) {
1607            int iSequence = id_[iColumn-firstDynamic_];
1608            solution[iSequence]=0.0;
1609            double b=0.0;
1610            // key is structural - where is slack
1611            iStatus = getStatus(iSet);
1612            assert (iStatus!=ClpSimplex::basic);
1613            if (iStatus==ClpSimplex::atLowerBound)
1614              b=lowerSet_[iSet];
1615            else
1616              b=upperSet_[iSet];
1617            // subtract out others at bounds
1618            for (int j=fullStart_[iSet];j<fullStart_[iSet+1];j++) 
1619              b -= solution[j];
1620            solution[iSequence]=b;
1621          }
1622        }
1623        for (iColumn=0;iColumn<numberGubColumns_;iColumn++) {
1624          double value = solution[iColumn];
1625          if (value) {
1626            for (CoinBigIndex j= startColumn_[iColumn];j<startColumn_[iColumn+1];j++) {
1627              int iRow = row_[j];
1628              rhs[iRow] -= element_[j]*value;
1629            }
1630          }
1631        }
1632        // now do lower and upper bounds on sets
1633        for (int iSet=0;iSet<numberSets_;iSet++) {
1634          iColumn = keyVariable_[iSet];
1635          double shift=0.0;
1636          for (int j=fullStart_[iSet];j<fullStart_[iSet+1];j++) {
1637            if (getDynamicStatus(j)!=inSmall&&j!=iColumn) {
1638              if (getDynamicStatus(j)==atLowerBound) {
1639                if (lowerColumn_)
1640                  shift += lowerColumn_[j];
1641              } else {
1642                shift += upperColumn_[j];
1643              }
1644            }
1645          }
1646          if (lowerSet_[iSet]>-1.0e20) 
1647            assert(fabs(lower_[iSet] - (lowerSet_[iSet]-shift))<1.0e-3);
1648          if (upperSet_[iSet]<1.0e20)
1649            assert(fabs(upper_[iSet] -( upperSet_[iSet]-shift))<1.0e-3);
1650        }
1651        delete [] solution;
1652      } else {
1653        // no bounds
1654        ClpSimplex::Status iStatus;
1655        for (int iSet=0;iSet<numberSets_;iSet++) {
1656          int iColumn = keyVariable_[iSet];
1657          if (iColumn<numberColumns) {
1658            int iSequence = id_[iColumn-firstDynamic_];
1659            double b=0.0;
1660            // key is structural - where is slack
1661            iStatus = getStatus(iSet);
1662            assert (iStatus!=ClpSimplex::basic);
1663            if (iStatus==ClpSimplex::atLowerBound)
1664              b=lower_[iSet];
1665            else
1666              b=upper_[iSet];
1667            if (b) {
1668              for (CoinBigIndex j= startColumn_[iSequence];j<startColumn_[iSequence+1];j++) {
1669                int iRow = row_[j];
1670                rhs[iRow] -= element_[j]*b;
1671              }
1672            }
1673          }
1674        }
1675      }
1676      for (iRow=0;iRow<numberRows;iRow++) {
1677        if (fabs(rhs[iRow]-rhsOffset_[iRow])>1.0e-3)
1678          printf("** bad effective %d - true %g old %g\n",iRow,rhs[iRow],rhsOffset_[iRow]);
1679      }
1680      memcpy(saveE,rhs,numberRows*sizeof(double));
1681      delete [] rhs;
1682    }
1683#endif
1684    if (forceRefresh||(refreshFrequency_&&model->numberIterations()>=
1685                       lastRefresh_+refreshFrequency_)) {
1686      int numberRows = model->numberRows();
1687      int numberColumns = model->numberColumns();
1688      int iRow;
1689      CoinZeroN(rhsOffset_,numberRows);
1690      // do ones at bounds before gub
1691      const double * smallSolution = model->solutionRegion();
1692      const double * element =matrix_->getElements();
1693      const int * row = matrix_->getIndices();
1694      const CoinBigIndex * startColumn = matrix_->getVectorStarts();
1695      const int * length = matrix_->getVectorLengths();
1696      int iColumn;
1697      for (iColumn=0;iColumn<firstDynamic_;iColumn++) {
1698        if (model->getStatus(iColumn)!=ClpSimplex::basic) {
1699          double value = smallSolution[iColumn];
1700          for (CoinBigIndex j=startColumn[iColumn];
1701               j<startColumn[iColumn]+length[iColumn];j++) {
1702            int jRow=row[j];
1703            rhsOffset_[jRow] -= value*element[j];
1704          }
1705        }
1706      }
1707      if (lowerColumn_||upperColumn_) {
1708        double * solution = new double [numberGubColumns_];
1709        for (iColumn=0;iColumn<numberGubColumns_;iColumn++) {
1710          double value=0.0;
1711          if(getDynamicStatus(iColumn)==atUpperBound)
1712            value = upperColumn_[iColumn];
1713          else if (lowerColumn_)
1714            value = lowerColumn_[iColumn];
1715          solution[iColumn]=value;
1716        }
1717        // ones in gub and in small problem
1718        for (iColumn=firstDynamic_;iColumn<firstAvailable_;iColumn++) {
1719          int jFull = id_[iColumn-firstDynamic_];
1720          solution[jFull]=smallSolution[iColumn];
1721        }
1722        // zero all basic in small model
1723        int * pivotVariable = model->pivotVariable();
1724        for (iRow=0;iRow<numberRows;iRow++) {
1725          int iColumn = pivotVariable[iRow];
1726          if (iColumn>=firstDynamic_&&iColumn<lastDynamic_) {
1727            int iSequence = id_[iColumn-firstDynamic_];
1728            solution[iSequence]=0.0;
1729          }
1730        }
1731        // and now compute value to use for key
1732        ClpSimplex::Status iStatus;
1733        int iSet;
1734        for ( iSet=0;iSet<numberSets_;iSet++) {
1735          iColumn = keyVariable_[iSet];
1736          if (iColumn<numberColumns) {
1737            int iSequence = id_[iColumn-firstDynamic_];
1738            solution[iSequence]=0.0;
1739            double b=0.0;
1740            // key is structural - where is slack
1741            iStatus = getStatus(iSet);
1742            assert (iStatus!=ClpSimplex::basic);
1743            if (iStatus==ClpSimplex::atLowerBound)
1744              b=lowerSet_[iSet];
1745            else
1746              b=upperSet_[iSet];
1747            // subtract out others at bounds
1748            for (int j=fullStart_[iSet];j<fullStart_[iSet+1];j++) 
1749              b -= solution[j];
1750            solution[iSequence]=b;
1751          }
1752        }
1753        for (iColumn=0;iColumn<numberGubColumns_;iColumn++) {
1754          double value = solution[iColumn];
1755          if (value) {
1756            for (CoinBigIndex j= startColumn_[iColumn];j<startColumn_[iColumn+1];j++) {
1757              int iRow = row_[j];
1758              rhsOffset_[iRow] -= element_[j]*value;
1759            }
1760          }
1761        }
1762        // now do lower and upper bounds on sets
1763        // and offset
1764        double objectiveOffset = 0.0;
1765        for ( iSet=0;iSet<numberSets_;iSet++) {
1766          iColumn = keyVariable_[iSet];
1767          double shift=0.0;
1768          for (CoinBigIndex j=fullStart_[iSet];j<fullStart_[iSet+1];j++) {
1769            if (getDynamicStatus(j)!=inSmall) {
1770              double value=0.0;
1771              if (getDynamicStatus(j)==atLowerBound) {
1772                if (lowerColumn_) 
1773                  value= lowerColumn_[j];
1774              } else { 
1775                value= upperColumn_[j];
1776              }
1777              if (j!=iColumn) 
1778                shift += value;
1779              objectiveOffset += value*cost_[j];
1780            }
1781          }
1782          if (lowerSet_[iSet]>-1.0e20)
1783            lower_[iSet] = lowerSet_[iSet]-shift;
1784          if (upperSet_[iSet]<1.0e20)
1785            upper_[iSet] = upperSet_[iSet]-shift;
1786        }
1787        delete [] solution;
1788        model->setObjectiveOffset(objectiveOffset_-objectiveOffset);
1789      } else {
1790        // no bounds
1791        ClpSimplex::Status iStatus;
1792        for (int iSet=0;iSet<numberSets_;iSet++) {
1793          int iColumn = keyVariable_[iSet];
1794          if (iColumn<numberColumns) {
1795            int iSequence = id_[iColumn-firstDynamic_];
1796            double b=0.0;
1797            // key is structural - where is slack
1798            iStatus = getStatus(iSet);
1799            assert (iStatus!=ClpSimplex::basic);
1800            if (iStatus==ClpSimplex::atLowerBound)
1801              b=lower_[iSet];
1802            else
1803              b=upper_[iSet];
1804            if (b) {
1805              for (CoinBigIndex j= startColumn_[iSequence];j<startColumn_[iSequence+1];j++) {
1806                int iRow = row_[j];
1807                rhsOffset_[iRow] -= element_[j]*b;
1808              }
1809            }
1810          }
1811        }
1812      }
1813#ifdef CLP_DEBUG
1814      if (saveE) {
1815        for (iRow=0;iRow<numberRows;iRow++) {
1816          if (fabs(saveE[iRow]-rhsOffset_[iRow])>1.0e-3)
1817            printf("** %d - old eff %g new %g\n",iRow,saveE[iRow],rhsOffset_[iRow]);
1818        }
1819        delete [] saveE;
1820      }
1821#endif
1822      lastRefresh_ = model->numberIterations();
1823    }
1824  }
1825  return rhsOffset_;
1826}
1827/*
1828  update information for a pivot (and effective rhs)
1829*/
1830int 
1831ClpGubDynamicMatrix::updatePivot(ClpSimplex * model,double oldInValue, double oldOutValue)
1832{
1833 
1834  // now update working model
1835  int sequenceIn = model->sequenceIn();
1836  int sequenceOut = model->sequenceOut();
1837  bool doPrinting = (model->messageHandler()->logLevel()==63);
1838  bool print=false;
1839  int iSet;
1840  int trueIn=-1;
1841  int trueOut=-1;
1842  int numberRows = model->numberRows();
1843  int numberColumns = model->numberColumns();
1844  if (sequenceIn==firstAvailable_) {
1845    if (doPrinting)
1846      printf("New variable ");
1847    if (sequenceIn!=sequenceOut) {
1848      insertNonBasic(firstAvailable_,backward_[firstAvailable_]);
1849      setDynamicStatus(id_[sequenceIn-firstDynamic_],inSmall);
1850      firstAvailable_++;
1851    } else {
1852      int bigSequence = id_[sequenceIn-firstDynamic_];
1853      if (model->getStatus(sequenceIn)==ClpSimplex::atUpperBound) 
1854        setDynamicStatus(bigSequence,atUpperBound);
1855      else
1856        setDynamicStatus(bigSequence,atLowerBound);
1857    }
1858    synchronize(model,8);
1859  }
1860  if (sequenceIn<lastDynamic_) {
1861    iSet = backward_[sequenceIn];
1862    if (iSet>=0) {
1863      int bigSequence = id_[sequenceIn-firstDynamic_];
1864      trueIn=bigSequence+numberRows+numberColumns+numberSets_;
1865      if (doPrinting)
1866        printf(" incoming set %d big seq %d",iSet,bigSequence);
1867      print=true;
1868    }
1869  } else if (sequenceIn>=numberRows+numberColumns) {
1870    trueIn = numberRows+numberColumns+gubSlackIn_;
1871  }
1872  if (sequenceOut<lastDynamic_) {
1873    iSet = backward_[sequenceOut];
1874    if (iSet>=0) {
1875      int bigSequence = id_[sequenceOut-firstDynamic_];
1876      trueOut=bigSequence+firstDynamic_;
1877      if (getDynamicStatus(bigSequence)!=inSmall) {
1878        if (model->getStatus(sequenceOut)==ClpSimplex::atUpperBound) 
1879          setDynamicStatus(bigSequence,atUpperBound);
1880        else
1881          setDynamicStatus(bigSequence,atLowerBound);
1882      }
1883      if (doPrinting)
1884        printf(" ,outgoing set %d big seq %d,",iSet,bigSequence);
1885      print=true;
1886      model->setSequenceIn(sequenceOut);
1887      synchronize(model,8);
1888      model->setSequenceIn(sequenceIn);
1889    }
1890  }
1891  if (print&&doPrinting)
1892    printf("\n");
1893  ClpGubMatrix::updatePivot(model,oldInValue,oldOutValue);
1894  // Redo true in and out
1895  if (trueIn>=0)
1896    trueSequenceIn_=trueIn;
1897  if (trueOut>=0)
1898    trueSequenceOut_=trueOut;
1899  if (doPrinting&&0) {
1900    for (int i=0;i<numberSets_;i++) {
1901      printf("set %d key %d lower %g upper %g\n",i,keyVariable_[i],lower_[i],upper_[i]);
1902      for (int j=fullStart_[i];j<fullStart_[i+1];j++) 
1903        if (getDynamicStatus(j)==atUpperBound) {
1904          bool print=true;
1905          for (int k=firstDynamic_;k<firstAvailable_;k++) {
1906            if (id_[k-firstDynamic_]==j)
1907              print=false;
1908            if (id_[k-firstDynamic_]==j)
1909              assert(getDynamicStatus(j)==inSmall);
1910          }
1911          if (print)
1912            printf("variable %d at ub\n",j);
1913        }
1914    }
1915  }
1916#ifdef CLP_DEBUG
1917  char * inSmall = new char [numberGubColumns_];
1918  memset(inSmall,0,numberGubColumns_);
1919  for (int i=0;i<numberGubColumns_;i++)
1920    if (getDynamicStatus(i)==ClpGubDynamicMatrix::inSmall)
1921      inSmall[i]=1;
1922  for (int i=firstDynamic_;i<firstAvailable_;i++) {
1923    int k=id_[i-firstDynamic_];
1924    inSmall[k]=0;
1925  }
1926  for (int i=0;i<numberGubColumns_;i++)
1927    assert (!inSmall[i]);
1928  delete [] inSmall;
1929#endif
1930  return 0;
1931}
1932void 
1933ClpGubDynamicMatrix::times(double scalar,
1934                           const double * x, double * y) const
1935{
1936  if (model_->specialOptions()!=16) {
1937    ClpPackedMatrix::times(scalar,x,y);
1938  } else {
1939    int iRow;
1940    int numberColumns = model_->numberColumns();
1941    int numberRows = model_->numberRows();
1942    const double * element =  matrix_->getElements();
1943    const int * row = matrix_->getIndices();
1944    const CoinBigIndex * startColumn = matrix_->getVectorStarts();
1945    const int * length = matrix_->getVectorLengths();
1946    int * pivotVariable = model_->pivotVariable();
1947    int numberToDo=0;
1948    for (iRow=0;iRow<numberRows;iRow++) {
1949      y[iRow] -= scalar*rhsOffset_[iRow];
1950      int iColumn = pivotVariable[iRow];
1951      if (iColumn<numberColumns) {
1952        int iSet = backward_[iColumn];
1953        if (iSet>=0&&toIndex_[iSet]<0) {
1954          toIndex_[iSet]=0;
1955          fromIndex_[numberToDo++]=iSet;
1956        }
1957        CoinBigIndex j;
1958        double value = scalar*x[iColumn];
1959        if (value) {
1960          for (j=startColumn[iColumn];
1961               j<startColumn[iColumn]+length[iColumn];j++) {
1962            int jRow=row[j];
1963            y[jRow] += value*element[j];
1964          }
1965        }
1966      }
1967    }
1968    // and gubs which are interacting
1969    for (int jSet=0;jSet<numberToDo;jSet++) {
1970      int iSet = fromIndex_[jSet];
1971      toIndex_[iSet]=-1;
1972      int iKey=keyVariable_[iSet];
1973      if (iKey<numberColumns) {
1974        double valueKey;
1975        if (getStatus(iSet)==ClpSimplex::atLowerBound) 
1976          valueKey = lower_[iSet];
1977        else
1978          valueKey = upper_[iSet];
1979        double value = scalar*(x[iKey]-valueKey);
1980        if (value) {
1981          for (CoinBigIndex j=startColumn[iKey];
1982               j<startColumn[iKey]+length[iKey];j++) {
1983            int jRow=row[j];
1984            y[jRow] += value*element[j];
1985          }
1986        }
1987      }
1988    }
1989  }
1990}
1991/* Just for debug - may be extended to other matrix types later.
1992   Returns number and sum of primal infeasibilities.
1993*/
1994int 
1995ClpGubDynamicMatrix::checkFeasible(ClpSimplex * model,double & sum) const
1996{
1997  int numberRows = model_->numberRows();
1998  double * rhs = new double[numberRows];
1999  int numberColumns = model_->numberColumns();
2000  int iRow;
2001  CoinZeroN(rhs,numberRows);
2002  // do ones at bounds before gub
2003  const double * smallSolution = model_->solutionRegion();
2004  const double * element =matrix_->getElements();
2005  const int * row = matrix_->getIndices();
2006  const CoinBigIndex * startColumn = matrix_->getVectorStarts();
2007  const int * length = matrix_->getVectorLengths();
2008  int iColumn;
2009  int numberInfeasible=0;
2010  const double * rowLower = model_->rowLower();
2011  const double * rowUpper = model_->rowUpper();
2012  sum=0.0;
2013  for (iRow=0;iRow<numberRows;iRow++) {
2014    double value=smallSolution[numberColumns+iRow];
2015    if (value<rowLower[iRow]-1.0e-5||
2016        value>rowUpper[iRow]+1.0e-5) {
2017      //printf("row %d %g %g %g\n",
2018      //     iRow,rowLower[iRow],value,rowUpper[iRow]);
2019      numberInfeasible++;
2020      sum += CoinMax(rowLower[iRow]-value,value-rowUpper[iRow]);
2021    }
2022    rhs[iRow]=value;
2023  }
2024  const double * columnLower = model_->columnLower();
2025  const double * columnUpper = model_->columnUpper();
2026  for (iColumn=0;iColumn<firstDynamic_;iColumn++) {
2027    double value=smallSolution[iColumn];
2028    if (value<columnLower[iColumn]-1.0e-5||
2029        value>columnUpper[iColumn]+1.0e-5) {
2030      //printf("column %d %g %g %g\n",
2031      //     iColumn,columnLower[iColumn],value,columnUpper[iColumn]);
2032      numberInfeasible++;
2033      sum += CoinMax(columnLower[iColumn]-value,value-columnUpper[iColumn]);
2034    }
2035    for (CoinBigIndex j=startColumn[iColumn];
2036         j<startColumn[iColumn]+length[iColumn];j++) {
2037      int jRow=row[j];
2038      rhs[jRow] -= value*element[j];
2039    }
2040  }
2041  double * solution = new double [numberGubColumns_];
2042  for (iColumn=0;iColumn<numberGubColumns_;iColumn++) {
2043    double value=0.0;
2044    if(getDynamicStatus(iColumn)==atUpperBound)
2045      value = upperColumn_[iColumn];
2046    else if (lowerColumn_)
2047      value = lowerColumn_[iColumn];
2048    solution[iColumn]=value;
2049  }
2050  // ones in small and gub
2051  for (iColumn=firstDynamic_;iColumn<firstAvailable_;iColumn++) {
2052    int jFull = id_[iColumn-firstDynamic_];
2053    solution[jFull]=smallSolution[iColumn];
2054  }
2055  // fill in all basic in small model
2056  int * pivotVariable = model_->pivotVariable();
2057  for (iRow=0;iRow<numberRows;iRow++) {
2058    int iColumn = pivotVariable[iRow];
2059    if (iColumn>=firstDynamic_&&iColumn<lastDynamic_) {
2060      int iSequence = id_[iColumn-firstDynamic_];
2061      solution[iSequence]=smallSolution[iColumn];
2062    }
2063  }
2064  // and now compute value to use for key
2065  ClpSimplex::Status iStatus;
2066  for (int iSet=0;iSet<numberSets_;iSet++) {
2067    iColumn = keyVariable_[iSet];
2068    if (iColumn<numberColumns) {
2069      int iSequence = id_[iColumn-firstDynamic_];
2070      solution[iSequence]=0.0;
2071      double b=0.0;
2072      // key is structural - where is slack
2073      iStatus = getStatus(iSet);
2074      assert (iStatus!=ClpSimplex::basic);
2075      if (iStatus==ClpSimplex::atLowerBound)
2076        b=lower_[iSet];
2077      else
2078        b=upper_[iSet];
2079      // subtract out others at bounds
2080      for (int j=fullStart_[iSet];j<fullStart_[iSet+1];j++) 
2081        b -= solution[j];
2082      solution[iSequence]=b;
2083    }
2084  }
2085  for (iColumn=0;iColumn<numberGubColumns_;iColumn++) {
2086    double value = solution[iColumn];
2087    if ((lowerColumn_&&value<lowerColumn_[iColumn]-1.0e-5)||
2088        (!lowerColumn_&&value<-1.0e-5)||
2089        (upperColumn_&&value>upperColumn_[iColumn]+1.0e-5)) {
2090      //printf("column %d %g %g %g\n",
2091      //     iColumn,lowerColumn_[iColumn],value,upperColumn_[iColumn]);
2092      numberInfeasible++;
2093    }
2094    if (value) {
2095      for (CoinBigIndex j= startColumn_[iColumn];j<startColumn_[iColumn+1];j++) {
2096        int iRow = row_[j];
2097        rhs[iRow] -= element_[j]*value;
2098      }
2099    }
2100  }
2101  for (iRow=0;iRow<numberRows;iRow++) {
2102    if (fabs(rhs[iRow])>1.0e-5)
2103      printf("rhs mismatch %d %g\n",iRow,rhs[iRow]);
2104  }
2105  delete [] solution;
2106  delete [] rhs;
2107  return numberInfeasible;
2108}
2109// Cleans data after setWarmStart
2110void 
2111ClpGubDynamicMatrix::cleanData(ClpSimplex * model)
2112{
2113  // and redo chains
2114  int numberColumns = model->numberColumns();
2115  int iColumn;
2116  // do backward
2117  int * mark = new int [numberGubColumns_];
2118  for (iColumn=0;iColumn<numberGubColumns_;iColumn++) 
2119    mark[iColumn]=-1;
2120  int i;
2121  for (i=0;i<firstDynamic_;i++) {
2122    assert (backward_[i]==-1);
2123    next_[i]=-1;
2124  }
2125  for (i=firstDynamic_;i<firstAvailable_;i++) {
2126    iColumn = id_[i-firstDynamic_];
2127    mark[iColumn]=i;
2128  }
2129  for (i=0;i<numberSets_;i++) {
2130    int iKey=keyVariable_[i];
2131    int lastNext = -1;
2132    int firstNext = -1;
2133    for (CoinBigIndex k= fullStart_[i];k<fullStart_[i+1];k++) {
2134      iColumn = mark[k];
2135      if (iColumn>=0) {
2136        if (iColumn!=iKey) {
2137          if (lastNext>=0)
2138            next_[lastNext]=iColumn;
2139          else
2140            firstNext = iColumn;
2141          lastNext=iColumn;
2142        }
2143        backward_[iColumn]=i;
2144      }
2145    }
2146    setFeasible(i);
2147    if (firstNext>=0) {
2148      // others
2149      next_[iKey]=firstNext;
2150      next_[lastNext]=-(iKey+1);
2151    } else if (iKey<numberColumns) {
2152      next_[iKey]=-(iKey+1);
2153    }
2154  }
2155  delete [] mark;
2156  // fill matrix
2157  double * element =  matrix_->getMutableElements();
2158  int * row = matrix_->getMutableIndices();
2159  CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
2160  int * length = matrix_->getMutableVectorLengths();
2161  CoinBigIndex numberElements = startColumn[firstDynamic_];
2162  for (i=firstDynamic_;i<firstAvailable_;i++) {
2163    int iColumn = id_[i-firstDynamic_];
2164    int numberThis = startColumn_[iColumn+1]-startColumn_[iColumn];
2165    length[i]=numberThis;
2166    for (CoinBigIndex jBigIndex=startColumn_[iColumn];
2167         jBigIndex<startColumn_[iColumn+1];jBigIndex++) {
2168      row[numberElements] = row_[jBigIndex];
2169      element[numberElements++] = element_[jBigIndex];
2170    }
2171    startColumn[i+1]=numberElements;
2172  }
2173}
Note: See TracBrowser for help on using the repository browser.