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

Last change on this file since 2030 was 1732, checked in by forrest, 8 years ago

various fixes plus slightly weighted pricing plus lagomory switches

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