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

Last change on this file since 1370 was 1370, checked in by forrest, 10 years ago

add ids

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