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

Last change on this file since 800 was 800, checked in by ladanyi, 13 years ago

finishing conversion to svn

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