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

Last change on this file since 754 was 754, checked in by andreasw, 14 years ago

first version

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