source: trunk/ClpGubDynamicMatrix.cpp @ 437

Last change on this file since 437 was 401, checked in by forrest, 16 years ago

quadratic

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 64.3 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 j=startColumn_[bestSequence];
464             j<startColumn_[bestSequence+1];j++) {
465          int jRow=row_[j];
466          value -= duals[jRow]*element_[j];
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 iColumn=0;iColumn<firstAvailable_;iColumn++) 
1412    next_[iColumn]=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        for (int iSet=0;iSet<numberSets_;iSet++) {
1707          iColumn = keyVariable_[iSet];
1708          if (iColumn<numberColumns) {
1709            int iSequence = id_[iColumn-firstDynamic_];
1710            solution[iSequence]=0.0;
1711            double b=0.0;
1712            // key is structural - where is slack
1713            iStatus = getStatus(iSet);
1714            assert (iStatus!=ClpSimplex::basic);
1715            if (iStatus==ClpSimplex::atLowerBound)
1716              b=lowerSet_[iSet];
1717            else
1718              b=upperSet_[iSet];
1719            // subtract out others at bounds
1720            for (int j=fullStart_[iSet];j<fullStart_[iSet+1];j++) 
1721              b -= solution[j];
1722            solution[iSequence]=b;
1723          }
1724        }
1725        for (iColumn=0;iColumn<numberGubColumns_;iColumn++) {
1726          double value = solution[iColumn];
1727          if (value) {
1728            for (CoinBigIndex j= startColumn_[iColumn];j<startColumn_[iColumn+1];j++) {
1729              int iRow = row_[j];
1730              rhsOffset_[iRow] -= element_[j]*value;
1731            }
1732          }
1733        }
1734        // now do lower and upper bounds on sets
1735        // and offset
1736        double objectiveOffset = 0.0;
1737        for (int iSet=0;iSet<numberSets_;iSet++) {
1738          iColumn = keyVariable_[iSet];
1739          double shift=0.0;
1740          for (CoinBigIndex j=fullStart_[iSet];j<fullStart_[iSet+1];j++) {
1741            if (getDynamicStatus(j)!=inSmall) {
1742              double value=0.0;
1743              if (getDynamicStatus(j)==atLowerBound) {
1744                if (lowerColumn_) 
1745                  value= lowerColumn_[j];
1746              } else { 
1747                value= upperColumn_[j];
1748              }
1749              if (j!=iColumn) 
1750                shift += value;
1751              objectiveOffset += value*cost_[j];
1752            }
1753          }
1754          if (lowerSet_[iSet]>-1.0e20)
1755            lower_[iSet] = lowerSet_[iSet]-shift;
1756          if (upperSet_[iSet]<1.0e20)
1757            upper_[iSet] = upperSet_[iSet]-shift;
1758        }
1759        delete [] solution;
1760        model->setObjectiveOffset(objectiveOffset_-objectiveOffset);
1761      } else {
1762        // no bounds
1763        ClpSimplex::Status iStatus;
1764        for (int iSet=0;iSet<numberSets_;iSet++) {
1765          int iColumn = keyVariable_[iSet];
1766          if (iColumn<numberColumns) {
1767            int iSequence = id_[iColumn-firstDynamic_];
1768            double b=0.0;
1769            // key is structural - where is slack
1770            iStatus = getStatus(iSet);
1771            assert (iStatus!=ClpSimplex::basic);
1772            if (iStatus==ClpSimplex::atLowerBound)
1773              b=lower_[iSet];
1774            else
1775              b=upper_[iSet];
1776            if (b) {
1777              for (CoinBigIndex j= startColumn_[iSequence];j<startColumn_[iSequence+1];j++) {
1778                int iRow = row_[j];
1779                rhsOffset_[iRow] -= element_[j]*b;
1780              }
1781            }
1782          }
1783        }
1784      }
1785#ifdef CLP_DEBUG
1786      if (saveE) {
1787        for (iRow=0;iRow<numberRows;iRow++) {
1788          if (fabs(saveE[iRow]-rhsOffset_[iRow])>1.0e-3)
1789            printf("** %d - old eff %g new %g\n",iRow,saveE[iRow],rhsOffset_[iRow]);
1790        }
1791        delete [] saveE;
1792      }
1793#endif
1794      lastRefresh_ = model->numberIterations();
1795    }
1796  }
1797  return rhsOffset_;
1798}
1799/*
1800  update information for a pivot (and effective rhs)
1801*/
1802int 
1803ClpGubDynamicMatrix::updatePivot(ClpSimplex * model,double oldInValue, double oldOutValue)
1804{
1805 
1806  // now update working model
1807  int sequenceIn = model->sequenceIn();
1808  int sequenceOut = model->sequenceOut();
1809  bool doPrinting = (model->messageHandler()->logLevel()==63);
1810  bool print=false;
1811  int iSet;
1812  int trueIn=-1;
1813  int trueOut=-1;
1814  int numberRows = model->numberRows();
1815  int numberColumns = model->numberColumns();
1816  if (sequenceIn==firstAvailable_) {
1817    if (doPrinting)
1818      printf("New variable ");
1819    if (sequenceIn!=sequenceOut) {
1820      insertNonBasic(firstAvailable_,backward_[firstAvailable_]);
1821      setDynamicStatus(id_[sequenceIn-firstDynamic_],inSmall);
1822      firstAvailable_++;
1823    } else {
1824      int bigSequence = id_[sequenceIn-firstDynamic_];
1825      if (model->getStatus(sequenceIn)==ClpSimplex::atUpperBound) 
1826        setDynamicStatus(bigSequence,atUpperBound);
1827      else
1828        setDynamicStatus(bigSequence,atLowerBound);
1829    }
1830    synchronize(model,8);
1831  }
1832  if (sequenceIn<lastDynamic_) {
1833    iSet = backward_[sequenceIn];
1834    if (iSet>=0) {
1835      int bigSequence = id_[sequenceIn-firstDynamic_];
1836      trueIn=bigSequence+numberRows+numberColumns+numberSets_;
1837      if (doPrinting)
1838        printf(" incoming set %d big seq %d",iSet,bigSequence);
1839      print=true;
1840    }
1841  } else if (sequenceIn>=numberRows+numberColumns) {
1842    trueIn = numberRows+numberColumns+gubSlackIn_;
1843  }
1844  if (sequenceOut<lastDynamic_) {
1845    iSet = backward_[sequenceOut];
1846    if (iSet>=0) {
1847      int bigSequence = id_[sequenceOut-firstDynamic_];
1848      trueOut=bigSequence+firstDynamic_;
1849      if (model->getStatus(sequenceOut)==ClpSimplex::atUpperBound) 
1850        setDynamicStatus(bigSequence,atUpperBound);
1851      else
1852        setDynamicStatus(bigSequence,atLowerBound);
1853      if (doPrinting)
1854        printf(" ,outgoing set %d big seq %d,",iSet,bigSequence);
1855      print=true;
1856      model->setSequenceIn(sequenceOut);
1857      synchronize(model,8);
1858      model->setSequenceIn(sequenceIn);
1859    }
1860  }
1861  if (print&&doPrinting)
1862    printf("\n");
1863  ClpGubMatrix::updatePivot(model,oldInValue,oldOutValue);
1864  // Redo true in and out
1865  if (trueIn>=0)
1866    trueSequenceIn_=trueIn;
1867  if (trueOut>=0)
1868    trueSequenceOut_=trueOut;
1869  if (doPrinting&&0) {
1870    for (int i=0;i<numberSets_;i++) {
1871      printf("set %d key %d lower %g upper %g\n",i,keyVariable_[i],lower_[i],upper_[i]);
1872      for (int j=fullStart_[i];j<fullStart_[i+1];j++) 
1873        if (getDynamicStatus(j)==atUpperBound) {
1874          bool print=true;
1875          for (int k=firstDynamic_;k<firstAvailable_;k++) {
1876            if (id_[k-firstDynamic_]==j)
1877              print=false;
1878            if (id_[k-firstDynamic_]==j)
1879              assert(getDynamicStatus(j)==inSmall);
1880          }
1881          if (print)
1882            printf("variable %d at ub\n",j);
1883        }
1884    }
1885  }
1886#ifdef CLP_DEBUG
1887  char * inSmall = new char [numberGubColumns_];
1888  memset(inSmall,0,numberGubColumns_);
1889  for (int i=0;i<numberGubColumns_;i++)
1890    if (getDynamicStatus(i)==ClpGubDynamicMatrix::inSmall)
1891      inSmall[i]=1;
1892  for (int i=firstDynamic_;i<firstAvailable_;i++) {
1893    int k=id_[i-firstDynamic_];
1894    inSmall[k]=0;
1895  }
1896  for (int i=0;i<numberGubColumns_;i++)
1897    assert (!inSmall[i]);
1898  delete [] inSmall;
1899#endif
1900  return 0;
1901}
1902void 
1903ClpGubDynamicMatrix::times(double scalar,
1904                           const double * x, double * y) const
1905{
1906  if (model_->specialOptions()!=16) {
1907    ClpPackedMatrix::times(scalar,x,y);
1908  } else {
1909    int iRow;
1910    int numberColumns = model_->numberColumns();
1911    int numberRows = model_->numberRows();
1912    const double * element =  matrix_->getElements();
1913    const int * row = matrix_->getIndices();
1914    const CoinBigIndex * startColumn = matrix_->getVectorStarts();
1915    const int * length = matrix_->getVectorLengths();
1916    int * pivotVariable = model_->pivotVariable();
1917    int numberToDo=0;
1918    for (iRow=0;iRow<numberRows;iRow++) {
1919      y[iRow] -= scalar*rhsOffset_[iRow];
1920      int iColumn = pivotVariable[iRow];
1921      if (iColumn<numberColumns) {
1922        int iSet = backward_[iColumn];
1923        if (iSet>=0&&toIndex_[iSet]<0) {
1924          toIndex_[iSet]=0;
1925          fromIndex_[numberToDo++]=iSet;
1926        }
1927        CoinBigIndex j;
1928        double value = scalar*x[iColumn];
1929        if (value) {
1930          for (j=startColumn[iColumn];
1931               j<startColumn[iColumn]+length[iColumn];j++) {
1932            int jRow=row[j];
1933            y[jRow] += value*element[j];
1934          }
1935        }
1936      }
1937    }
1938    // and gubs which are interacting
1939    for (int jSet=0;jSet<numberToDo;jSet++) {
1940      int iSet = fromIndex_[jSet];
1941      toIndex_[iSet]=-1;
1942      int iKey=keyVariable_[iSet];
1943      if (iKey<numberColumns) {
1944        double valueKey;
1945        if (getStatus(iSet)==ClpSimplex::atLowerBound) 
1946          valueKey = lower_[iSet];
1947        else
1948          valueKey = upper_[iSet];
1949        double value = scalar*(x[iKey]-valueKey);
1950        if (value) {
1951          for (CoinBigIndex j=startColumn[iKey];
1952               j<startColumn[iKey]+length[iKey];j++) {
1953            int jRow=row[j];
1954            y[jRow] += value*element[j];
1955          }
1956        }
1957      }
1958    }
1959  }
1960}
1961/* Just for debug - may be extended to other matrix types later.
1962   Returns number of primal infeasibilities.
1963*/
1964int 
1965ClpGubDynamicMatrix::checkFeasible(ClpSimplex * model) const
1966{
1967  int numberRows = model_->numberRows();
1968  double * rhs = new double[numberRows];
1969  int numberColumns = model_->numberColumns();
1970  int iRow;
1971  CoinZeroN(rhs,numberRows);
1972  // do ones at bounds before gub
1973  const double * smallSolution = model_->solutionRegion();
1974  const double * element =matrix_->getElements();
1975  const int * row = matrix_->getIndices();
1976  const CoinBigIndex * startColumn = matrix_->getVectorStarts();
1977  const int * length = matrix_->getVectorLengths();
1978  int iColumn;
1979  int numberInfeasible=0;
1980  const double * rowLower = model_->rowLower();
1981  const double * rowUpper = model_->rowUpper();
1982  for (iRow=0;iRow<numberRows;iRow++) {
1983    double value=smallSolution[numberColumns+iRow];
1984    if (value<rowLower[iRow]-1.0e-5||
1985        value>rowUpper[iRow]+1.0e-5) {
1986      //printf("row %d %g %g %g\n",
1987      //     iRow,rowLower[iRow],value,rowUpper[iRow]);
1988      numberInfeasible++;
1989    }
1990    rhs[iRow]=value;
1991  }
1992  const double * columnLower = model_->columnLower();
1993  const double * columnUpper = model_->columnUpper();
1994  for (iColumn=0;iColumn<firstDynamic_;iColumn++) {
1995    double value=smallSolution[iColumn];
1996    if (value<columnLower[iColumn]-1.0e-5||
1997        value>columnUpper[iColumn]+1.0e-5) {
1998      //printf("column %d %g %g %g\n",
1999      //     iColumn,columnLower[iColumn],value,columnUpper[iColumn]);
2000      numberInfeasible++;
2001    }
2002    for (CoinBigIndex j=startColumn[iColumn];
2003         j<startColumn[iColumn]+length[iColumn];j++) {
2004      int jRow=row[j];
2005      rhs[jRow] -= value*element[j];
2006    }
2007  }
2008  double * solution = new double [numberGubColumns_];
2009  for (iColumn=0;iColumn<numberGubColumns_;iColumn++) {
2010    double value=0.0;
2011    if(getDynamicStatus(iColumn)==atUpperBound)
2012      value = upperColumn_[iColumn];
2013    else if (lowerColumn_)
2014      value = lowerColumn_[iColumn];
2015    solution[iColumn]=value;
2016  }
2017  // ones in small and gub
2018  for (iColumn=firstDynamic_;iColumn<firstAvailable_;iColumn++) {
2019    int jFull = id_[iColumn-firstDynamic_];
2020    solution[jFull]=smallSolution[iColumn];
2021  }
2022  // fill in all basic in small model
2023  int * pivotVariable = model_->pivotVariable();
2024  for (iRow=0;iRow<numberRows;iRow++) {
2025    int iColumn = pivotVariable[iRow];
2026    if (iColumn>=firstDynamic_&&iColumn<lastDynamic_) {
2027      int iSequence = id_[iColumn-firstDynamic_];
2028      solution[iSequence]=smallSolution[iColumn];
2029    }
2030  }
2031  // and now compute value to use for key
2032  ClpSimplex::Status iStatus;
2033  for (int iSet=0;iSet<numberSets_;iSet++) {
2034    iColumn = keyVariable_[iSet];
2035    if (iColumn<numberColumns) {
2036      int iSequence = id_[iColumn-firstDynamic_];
2037      solution[iSequence]=0.0;
2038      double b=0.0;
2039      // key is structural - where is slack
2040      iStatus = getStatus(iSet);
2041      assert (iStatus!=ClpSimplex::basic);
2042      if (iStatus==ClpSimplex::atLowerBound)
2043        b=lower_[iSet];
2044      else
2045        b=upper_[iSet];
2046      // subtract out others at bounds
2047      for (int j=fullStart_[iSet];j<fullStart_[iSet+1];j++) 
2048        b -= solution[j];
2049      solution[iSequence]=b;
2050    }
2051  }
2052  for (iColumn=0;iColumn<numberGubColumns_;iColumn++) {
2053    double value = solution[iColumn];
2054    if ((lowerColumn_&&value<lowerColumn_[iColumn]-1.0e-5)||
2055        (!lowerColumn_&&value<-1.0e-5)||
2056        (upperColumn_&&value>upperColumn_[iColumn]+1.0e-5)) {
2057      //printf("column %d %g %g %g\n",
2058      //     iColumn,lowerColumn_[iColumn],value,upperColumn_[iColumn]);
2059      numberInfeasible++;
2060    }
2061    if (value) {
2062      for (CoinBigIndex j= startColumn_[iColumn];j<startColumn_[iColumn+1];j++) {
2063        int iRow = row_[j];
2064        rhs[iRow] -= element_[j]*value;
2065      }
2066    }
2067  }
2068  for (iRow=0;iRow<numberRows;iRow++) {
2069    if (fabs(rhs[iRow])>1.0e-5)
2070      printf("rhs mismatch %d %g\n",iRow,rhs[iRow]);
2071  }
2072  delete [] solution;
2073  delete [] rhs;
2074  return numberInfeasible;
2075}
Note: See TracBrowser for help on using the repository browser.