source: trunk/Clp/src/ClpGubMatrix.cpp @ 2385

Last change on this file since 2385 was 2385, checked in by unxusr, 4 months ago

formatting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 138.5 KB
Line 
1/* $Id: ClpGubMatrix.cpp 2385 2019-01-06 19:43:06Z unxusr $ */
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#include <cstdio>
7
8#include "CoinPragma.hpp"
9#include "CoinIndexedVector.hpp"
10#include "CoinHelperFunctions.hpp"
11
12#include "ClpSimplex.hpp"
13#include "ClpFactorization.hpp"
14#include "ClpQuadraticObjective.hpp"
15#include "ClpNonLinearCost.hpp"
16// at end to get min/max!
17#include "ClpGubMatrix.hpp"
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//-------------------------------------------------------------------
29ClpGubMatrix::ClpGubMatrix()
30  : ClpPackedMatrix()
31  , sumDualInfeasibilities_(0.0)
32  , sumPrimalInfeasibilities_(0.0)
33  , sumOfRelaxedDualInfeasibilities_(0.0)
34  , sumOfRelaxedPrimalInfeasibilities_(0.0)
35  , infeasibilityWeight_(0.0)
36  , start_(NULL)
37  , end_(NULL)
38  , lower_(NULL)
39  , upper_(NULL)
40  , status_(NULL)
41  , saveStatus_(NULL)
42  , savedKeyVariable_(NULL)
43  , backward_(NULL)
44  , backToPivotRow_(NULL)
45  , changeCost_(NULL)
46  , keyVariable_(NULL)
47  , next_(NULL)
48  , toIndex_(NULL)
49  , fromIndex_(NULL)
50  , model_(NULL)
51  , numberDualInfeasibilities_(0)
52  , numberPrimalInfeasibilities_(0)
53  , noCheck_(-1)
54  , numberSets_(0)
55  , saveNumber_(0)
56  , possiblePivotKey_(0)
57  , gubSlackIn_(-1)
58  , firstGub_(0)
59  , lastGub_(0)
60  , gubType_(0)
61{
62  setType(16);
63}
64
65//-------------------------------------------------------------------
66// Copy constructor
67//-------------------------------------------------------------------
68ClpGubMatrix::ClpGubMatrix(const ClpGubMatrix &rhs)
69  : ClpPackedMatrix(rhs)
70{
71  numberSets_ = rhs.numberSets_;
72  saveNumber_ = rhs.saveNumber_;
73  possiblePivotKey_ = rhs.possiblePivotKey_;
74  gubSlackIn_ = rhs.gubSlackIn_;
75  start_ = ClpCopyOfArray(rhs.start_, numberSets_);
76  end_ = ClpCopyOfArray(rhs.end_, numberSets_);
77  lower_ = ClpCopyOfArray(rhs.lower_, numberSets_);
78  upper_ = ClpCopyOfArray(rhs.upper_, numberSets_);
79  status_ = ClpCopyOfArray(rhs.status_, numberSets_);
80  saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberSets_);
81  savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_, numberSets_);
82  int numberColumns = getNumCols();
83  backward_ = ClpCopyOfArray(rhs.backward_, numberColumns);
84  backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, numberColumns);
85  changeCost_ = ClpCopyOfArray(rhs.changeCost_, getNumRows() + numberSets_);
86  fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + numberSets_ + 1);
87  keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_);
88  // find longest set
89  int *longest = new int[numberSets_];
90  CoinZeroN(longest, numberSets_);
91  int j;
92  for (j = 0; j < numberColumns; j++) {
93    int iSet = backward_[j];
94    if (iSet >= 0)
95      longest[iSet]++;
96  }
97  int length = 0;
98  for (j = 0; j < numberSets_; j++)
99    length = CoinMax(length, longest[j]);
100  next_ = ClpCopyOfArray(rhs.next_, numberColumns + numberSets_ + 2 * length);
101  toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_);
102  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
103  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
104  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
105  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
106  infeasibilityWeight_ = rhs.infeasibilityWeight_;
107  numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
108  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
109  noCheck_ = rhs.noCheck_;
110  firstGub_ = rhs.firstGub_;
111  lastGub_ = rhs.lastGub_;
112  gubType_ = rhs.gubType_;
113  model_ = rhs.model_;
114}
115
116//-------------------------------------------------------------------
117// assign matrix (for space reasons)
118//-------------------------------------------------------------------
119ClpGubMatrix::ClpGubMatrix(CoinPackedMatrix *rhs)
120  : ClpPackedMatrix(rhs)
121  , sumDualInfeasibilities_(0.0)
122  , sumPrimalInfeasibilities_(0.0)
123  , sumOfRelaxedDualInfeasibilities_(0.0)
124  , sumOfRelaxedPrimalInfeasibilities_(0.0)
125  , infeasibilityWeight_(0.0)
126  , start_(NULL)
127  , end_(NULL)
128  , lower_(NULL)
129  , upper_(NULL)
130  , status_(NULL)
131  , saveStatus_(NULL)
132  , savedKeyVariable_(NULL)
133  , backward_(NULL)
134  , backToPivotRow_(NULL)
135  , changeCost_(NULL)
136  , keyVariable_(NULL)
137  , next_(NULL)
138  , toIndex_(NULL)
139  , fromIndex_(NULL)
140  , model_(NULL)
141  , numberDualInfeasibilities_(0)
142  , numberPrimalInfeasibilities_(0)
143  , noCheck_(-1)
144  , numberSets_(0)
145  , saveNumber_(0)
146  , possiblePivotKey_(0)
147  , gubSlackIn_(-1)
148  , firstGub_(0)
149  , lastGub_(0)
150  , gubType_(0)
151{
152  setType(16);
153}
154
155/* This takes over ownership (for space reasons) and is the
156   real constructor*/
157ClpGubMatrix::ClpGubMatrix(ClpPackedMatrix *matrix, int numberSets,
158  const int *start, const int *end,
159  const double *lower, const double *upper,
160  const unsigned char *status)
161  : ClpPackedMatrix(matrix->matrix())
162  , sumDualInfeasibilities_(0.0)
163  , sumPrimalInfeasibilities_(0.0)
164  , sumOfRelaxedDualInfeasibilities_(0.0)
165  , sumOfRelaxedPrimalInfeasibilities_(0.0)
166  , numberDualInfeasibilities_(0)
167  , numberPrimalInfeasibilities_(0)
168  , saveNumber_(0)
169  , possiblePivotKey_(0)
170  , gubSlackIn_(-1)
171{
172  model_ = NULL;
173  numberSets_ = numberSets;
174  start_ = ClpCopyOfArray(start, numberSets_);
175  end_ = ClpCopyOfArray(end, numberSets_);
176  lower_ = ClpCopyOfArray(lower, numberSets_);
177  upper_ = ClpCopyOfArray(upper, numberSets_);
178  // Check valid and ordered
179  int last = -1;
180  int numberColumns = matrix_->getNumCols();
181  int numberRows = matrix_->getNumRows();
182  backward_ = new int[numberColumns];
183  backToPivotRow_ = new int[numberColumns];
184  changeCost_ = new double[numberRows + numberSets_];
185  keyVariable_ = new int[numberSets_];
186  // signal to need new ordering
187  next_ = NULL;
188  for (int iColumn = 0; iColumn < numberColumns; iColumn++)
189    backward_[iColumn] = -1;
190
191  int iSet;
192  for (iSet = 0; iSet < numberSets_; iSet++) {
193    // set key variable as slack
194    keyVariable_[iSet] = iSet + numberColumns;
195    if (start_[iSet] < 0 || start_[iSet] >= numberColumns)
196      throw CoinError("Index out of range", "constructor", "ClpGubMatrix");
197    if (end_[iSet] < 0 || end_[iSet] > numberColumns)
198      throw CoinError("Index out of range", "constructor", "ClpGubMatrix");
199    if (end_[iSet] <= start_[iSet])
200      throw CoinError("Empty or negative set", "constructor", "ClpGubMatrix");
201    if (start_[iSet] < last)
202      throw CoinError("overlapping or non-monotonic sets", "constructor", "ClpGubMatrix");
203    last = end_[iSet];
204    int j;
205    for (j = start_[iSet]; j < end_[iSet]; j++)
206      backward_[j] = iSet;
207  }
208  // Find type of gub
209  firstGub_ = numberColumns + 1;
210  lastGub_ = -1;
211  int i;
212  for (i = 0; i < numberColumns; i++) {
213    if (backward_[i] >= 0) {
214      firstGub_ = CoinMin(firstGub_, i);
215      lastGub_ = CoinMax(lastGub_, i);
216    }
217  }
218  gubType_ = 0;
219  // adjust lastGub_
220  if (lastGub_ > 0)
221    lastGub_++;
222  for (i = firstGub_; i < lastGub_; i++) {
223    if (backward_[i] < 0) {
224      gubType_ = 1;
225      printf("interior non gub %d\n", i);
226      break;
227    }
228  }
229  if (status) {
230    status_ = ClpCopyOfArray(status, numberSets_);
231  } else {
232    status_ = new unsigned char[numberSets_];
233    memset(status_, 0, numberSets_);
234    int i;
235    for (i = 0; i < numberSets_; i++) {
236      // make slack key
237      setStatus(i, ClpSimplex::basic);
238    }
239  }
240  saveStatus_ = new unsigned char[numberSets_];
241  memset(saveStatus_, 0, numberSets_);
242  savedKeyVariable_ = new int[numberSets_];
243  memset(savedKeyVariable_, 0, numberSets_ * sizeof(int));
244  noCheck_ = -1;
245  infeasibilityWeight_ = 0.0;
246}
247
248ClpGubMatrix::ClpGubMatrix(const CoinPackedMatrix &rhs)
249  : ClpPackedMatrix(rhs)
250  , sumDualInfeasibilities_(0.0)
251  , sumPrimalInfeasibilities_(0.0)
252  , sumOfRelaxedDualInfeasibilities_(0.0)
253  , sumOfRelaxedPrimalInfeasibilities_(0.0)
254  , infeasibilityWeight_(0.0)
255  , start_(NULL)
256  , end_(NULL)
257  , lower_(NULL)
258  , upper_(NULL)
259  , status_(NULL)
260  , saveStatus_(NULL)
261  , savedKeyVariable_(NULL)
262  , backward_(NULL)
263  , backToPivotRow_(NULL)
264  , changeCost_(NULL)
265  , keyVariable_(NULL)
266  , next_(NULL)
267  , toIndex_(NULL)
268  , fromIndex_(NULL)
269  , model_(NULL)
270  , numberDualInfeasibilities_(0)
271  , numberPrimalInfeasibilities_(0)
272  , noCheck_(-1)
273  , numberSets_(0)
274  , saveNumber_(0)
275  , possiblePivotKey_(0)
276  , gubSlackIn_(-1)
277  , firstGub_(0)
278  , lastGub_(0)
279  , gubType_(0)
280{
281  setType(16);
282}
283
284//-------------------------------------------------------------------
285// Destructor
286//-------------------------------------------------------------------
287ClpGubMatrix::~ClpGubMatrix()
288{
289  delete[] start_;
290  delete[] end_;
291  delete[] lower_;
292  delete[] upper_;
293  delete[] status_;
294  delete[] saveStatus_;
295  delete[] savedKeyVariable_;
296  delete[] backward_;
297  delete[] backToPivotRow_;
298  delete[] changeCost_;
299  delete[] keyVariable_;
300  delete[] next_;
301  delete[] toIndex_;
302  delete[] fromIndex_;
303}
304
305//----------------------------------------------------------------
306// Assignment operator
307//-------------------------------------------------------------------
308ClpGubMatrix &
309ClpGubMatrix::operator=(const ClpGubMatrix &rhs)
310{
311  if (this != &rhs) {
312    ClpPackedMatrix::operator=(rhs);
313    delete[] start_;
314    delete[] end_;
315    delete[] lower_;
316    delete[] upper_;
317    delete[] status_;
318    delete[] saveStatus_;
319    delete[] savedKeyVariable_;
320    delete[] backward_;
321    delete[] backToPivotRow_;
322    delete[] changeCost_;
323    delete[] keyVariable_;
324    delete[] next_;
325    delete[] toIndex_;
326    delete[] fromIndex_;
327    numberSets_ = rhs.numberSets_;
328    saveNumber_ = rhs.saveNumber_;
329    possiblePivotKey_ = rhs.possiblePivotKey_;
330    gubSlackIn_ = rhs.gubSlackIn_;
331    start_ = ClpCopyOfArray(rhs.start_, numberSets_);
332    end_ = ClpCopyOfArray(rhs.end_, numberSets_);
333    lower_ = ClpCopyOfArray(rhs.lower_, numberSets_);
334    upper_ = ClpCopyOfArray(rhs.upper_, numberSets_);
335    status_ = ClpCopyOfArray(rhs.status_, numberSets_);
336    saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberSets_);
337    savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_, numberSets_);
338    int numberColumns = getNumCols();
339    backward_ = ClpCopyOfArray(rhs.backward_, numberColumns);
340    backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, numberColumns);
341    changeCost_ = ClpCopyOfArray(rhs.changeCost_, getNumRows() + numberSets_);
342    fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + numberSets_ + 1);
343    keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_);
344    // find longest set
345    int *longest = new int[numberSets_];
346    CoinZeroN(longest, numberSets_);
347    int j;
348    for (j = 0; j < numberColumns; j++) {
349      int iSet = backward_[j];
350      if (iSet >= 0)
351        longest[iSet]++;
352    }
353    int length = 0;
354    for (j = 0; j < numberSets_; j++)
355      length = CoinMax(length, longest[j]);
356    next_ = ClpCopyOfArray(rhs.next_, numberColumns + numberSets_ + 2 * length);
357    toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_);
358    sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
359    sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
360    sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
361    sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
362    infeasibilityWeight_ = rhs.infeasibilityWeight_;
363    numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
364    numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
365    noCheck_ = rhs.noCheck_;
366    firstGub_ = rhs.firstGub_;
367    lastGub_ = rhs.lastGub_;
368    gubType_ = rhs.gubType_;
369    model_ = rhs.model_;
370  }
371  return *this;
372}
373//-------------------------------------------------------------------
374// Clone
375//-------------------------------------------------------------------
376ClpMatrixBase *ClpGubMatrix::clone() const
377{
378  return new ClpGubMatrix(*this);
379}
380/* Subset clone (without gaps).  Duplicates are allowed
381   and order is as given */
382ClpMatrixBase *
383ClpGubMatrix::subsetClone(int numberRows, const int *whichRows,
384  int numberColumns,
385  const int *whichColumns) const
386{
387  return new ClpGubMatrix(*this, numberRows, whichRows,
388    numberColumns, whichColumns);
389}
390/* Returns a new matrix in reverse order without gaps
391   Is allowed to return NULL if doesn't want to have row copy */
392ClpMatrixBase *
393ClpGubMatrix::reverseOrderedCopy() const
394{
395  return NULL;
396}
397int ClpGubMatrix::hiddenRows() const
398{
399  return numberSets_;
400}
401/* Subset constructor (without gaps).  Duplicates are allowed
402   and order is as given */
403ClpGubMatrix::ClpGubMatrix(
404  const ClpGubMatrix &rhs,
405  int numberRows, const int *whichRows,
406  int numberColumns, const int *whichColumns)
407  : ClpPackedMatrix(rhs, numberRows, whichRows, numberColumns, whichColumns)
408{
409  // Assuming no gub rows deleted
410  // We also assume all sets in same order
411  // Get array with backward pointers
412  int numberColumnsOld = rhs.matrix_->getNumCols();
413  int *array = new int[numberColumnsOld];
414  int i;
415  for (i = 0; i < numberColumnsOld; i++)
416    array[i] = -1;
417  for (int iSet = 0; iSet < numberSets_; iSet++) {
418    for (int j = start_[iSet]; j < end_[iSet]; j++)
419      array[j] = iSet;
420  }
421  numberSets_ = -1;
422  int lastSet = -1;
423  bool inSet = false;
424  for (i = 0; i < numberColumns; i++) {
425    int iColumn = whichColumns[i];
426    int iSet = array[iColumn];
427    if (iSet < 0) {
428      inSet = false;
429    } else {
430      if (!inSet) {
431        // start of new set but check okay
432        if (iSet <= lastSet)
433          throw CoinError("overlapping or non-monotonic sets", "subset constructor", "ClpGubMatrix");
434        lastSet = iSet;
435        numberSets_++;
436        start_[numberSets_] = i;
437        end_[numberSets_] = i + 1;
438        lower_[numberSets_] = lower_[iSet];
439        upper_[numberSets_] = upper_[iSet];
440        inSet = true;
441      } else {
442        if (iSet < lastSet) {
443          throw CoinError("overlapping or non-monotonic sets", "subset constructor", "ClpGubMatrix");
444        } else if (iSet == lastSet) {
445          end_[numberSets_] = i + 1;
446        } else {
447          // new set
448          lastSet = iSet;
449          numberSets_++;
450          start_[numberSets_] = i;
451          end_[numberSets_] = i + 1;
452          lower_[numberSets_] = lower_[iSet];
453          upper_[numberSets_] = upper_[iSet];
454        }
455      }
456    }
457  }
458  delete[] array;
459  numberSets_++; // adjust
460  // Find type of gub
461  firstGub_ = numberColumns + 1;
462  lastGub_ = -1;
463  for (i = 0; i < numberColumns; i++) {
464    if (backward_[i] >= 0) {
465      firstGub_ = CoinMin(firstGub_, i);
466      lastGub_ = CoinMax(lastGub_, i);
467    }
468  }
469  if (lastGub_ > 0)
470    lastGub_++;
471  gubType_ = 0;
472  for (i = firstGub_; i < lastGub_; i++) {
473    if (backward_[i] < 0) {
474      gubType_ = 1;
475      break;
476    }
477  }
478
479  // Make sure key is feasible if only key in set
480}
481ClpGubMatrix::ClpGubMatrix(
482  const CoinPackedMatrix &rhs,
483  int numberRows, const int *whichRows,
484  int numberColumns, const int *whichColumns)
485  : ClpPackedMatrix(rhs, numberRows, whichRows, numberColumns, whichColumns)
486  , sumDualInfeasibilities_(0.0)
487  , sumPrimalInfeasibilities_(0.0)
488  , sumOfRelaxedDualInfeasibilities_(0.0)
489  , sumOfRelaxedPrimalInfeasibilities_(0.0)
490  , start_(NULL)
491  , end_(NULL)
492  , lower_(NULL)
493  , upper_(NULL)
494  , backward_(NULL)
495  , backToPivotRow_(NULL)
496  , changeCost_(NULL)
497  , keyVariable_(NULL)
498  , next_(NULL)
499  , toIndex_(NULL)
500  , fromIndex_(NULL)
501  , numberDualInfeasibilities_(0)
502  , numberPrimalInfeasibilities_(0)
503  , numberSets_(0)
504  , saveNumber_(0)
505  , possiblePivotKey_(0)
506  , gubSlackIn_(-1)
507  , firstGub_(0)
508  , lastGub_(0)
509  , gubType_(0)
510{
511  setType(16);
512}
513/* Return <code>x * A + y</code> in <code>z</code>.
514        Squashes small elements and knows about ClpSimplex */
515void ClpGubMatrix::transposeTimes(const ClpSimplex *model, double scalar,
516  const CoinIndexedVector *rowArray,
517  CoinIndexedVector *y,
518  CoinIndexedVector *columnArray) const
519{
520  columnArray->clear();
521  double *pi = rowArray->denseVector();
522  int numberNonZero = 0;
523  int *index = columnArray->getIndices();
524  double *array = columnArray->denseVector();
525  int numberInRowArray = rowArray->getNumElements();
526  // maybe I need one in OsiSimplex
527  double zeroTolerance = model->zeroTolerance();
528  int numberRows = model->numberRows();
529  ClpPackedMatrix *rowCopy = dynamic_cast< ClpPackedMatrix * >(model->rowCopy());
530  bool packed = rowArray->packedMode();
531  double factor = 0.3;
532  // We may not want to do by row if there may be cache problems
533  int numberColumns = model->numberColumns();
534  // It would be nice to find L2 cache size - for moment 512K
535  // Be slightly optimistic
536  if (numberColumns * sizeof(double) > 1000000) {
537    if (numberRows * 10 < numberColumns)
538      factor = 0.1;
539    else if (numberRows * 4 < numberColumns)
540      factor = 0.15;
541    else if (numberRows * 2 < numberColumns)
542      factor = 0.2;
543    //if (model->numberIterations()%50==0)
544    //printf("%d nonzero\n",numberInRowArray);
545  }
546  // reduce for gub
547  factor *= 0.5;
548  assert(!y->getNumElements());
549  if (numberInRowArray > factor * numberRows || !rowCopy) {
550    // do by column
551    int iColumn;
552    // get matrix data pointers
553    const int *row = matrix_->getIndices();
554    const CoinBigIndex *columnStart = matrix_->getVectorStarts();
555    const int *columnLength = matrix_->getVectorLengths();
556    const double *elementByColumn = matrix_->getElements();
557    const double *rowScale = model->rowScale();
558    int numberColumns = model->numberColumns();
559    int iSet = -1;
560    double djMod = 0.0;
561    if (packed) {
562      // need to expand pi into y
563      assert(y->capacity() >= numberRows);
564      double *piOld = pi;
565      pi = y->denseVector();
566      const int *whichRow = rowArray->getIndices();
567      int i;
568      if (!rowScale) {
569        // modify pi so can collapse to one loop
570        for (i = 0; i < numberInRowArray; i++) {
571          int iRow = whichRow[i];
572          pi[iRow] = scalar * piOld[i];
573        }
574        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
575          if (backward_[iColumn] != iSet) {
576            // get pi on gub row
577            iSet = backward_[iColumn];
578            if (iSet >= 0) {
579              int iBasic = keyVariable_[iSet];
580              if (iBasic < numberColumns) {
581                // get dj without
582                assert(model->getStatus(iBasic) == ClpSimplex::basic);
583                djMod = 0.0;
584                for (CoinBigIndex j = columnStart[iBasic];
585                     j < columnStart[iBasic] + columnLength[iBasic]; j++) {
586                  int jRow = row[j];
587                  djMod -= pi[jRow] * elementByColumn[j];
588                }
589              } else {
590                djMod = 0.0;
591              }
592            } else {
593              djMod = 0.0;
594            }
595          }
596          double value = -djMod;
597          CoinBigIndex j;
598          for (j = columnStart[iColumn];
599               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
600            int iRow = row[j];
601            value += pi[iRow] * elementByColumn[j];
602          }
603          if (fabs(value) > zeroTolerance) {
604            array[numberNonZero] = value;
605            index[numberNonZero++] = iColumn;
606          }
607        }
608      } else {
609        // scaled
610        // modify pi so can collapse to one loop
611        for (i = 0; i < numberInRowArray; i++) {
612          int iRow = whichRow[i];
613          pi[iRow] = scalar * piOld[i] * rowScale[iRow];
614        }
615        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
616          if (backward_[iColumn] != iSet) {
617            // get pi on gub row
618            iSet = backward_[iColumn];
619            if (iSet >= 0) {
620              int iBasic = keyVariable_[iSet];
621              if (iBasic < numberColumns) {
622                // get dj without
623                assert(model->getStatus(iBasic) == ClpSimplex::basic);
624                djMod = 0.0;
625                // scaled
626                for (CoinBigIndex j = columnStart[iBasic];
627                     j < columnStart[iBasic] + columnLength[iBasic]; j++) {
628                  int jRow = row[j];
629                  djMod -= pi[jRow] * elementByColumn[j] * rowScale[jRow];
630                }
631              } else {
632                djMod = 0.0;
633              }
634            } else {
635              djMod = 0.0;
636            }
637          }
638          double value = -djMod;
639          CoinBigIndex j;
640          const double *columnScale = model->columnScale();
641          for (j = columnStart[iColumn];
642               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
643            int iRow = row[j];
644            value += pi[iRow] * elementByColumn[j];
645          }
646          value *= columnScale[iColumn];
647          if (fabs(value) > zeroTolerance) {
648            array[numberNonZero] = value;
649            index[numberNonZero++] = iColumn;
650          }
651        }
652      }
653      // zero out
654      for (i = 0; i < numberInRowArray; i++) {
655        int iRow = whichRow[i];
656        pi[iRow] = 0.0;
657      }
658    } else {
659      // code later
660      assert(packed);
661      if (!rowScale) {
662        if (scalar == -1.0) {
663          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
664            double value = 0.0;
665            CoinBigIndex j;
666            for (j = columnStart[iColumn];
667                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
668              int iRow = row[j];
669              value += pi[iRow] * elementByColumn[j];
670            }
671            if (fabs(value) > zeroTolerance) {
672              index[numberNonZero++] = iColumn;
673              array[iColumn] = -value;
674            }
675          }
676        } else if (scalar == 1.0) {
677          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
678            double value = 0.0;
679            CoinBigIndex j;
680            for (j = columnStart[iColumn];
681                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
682              int iRow = row[j];
683              value += pi[iRow] * elementByColumn[j];
684            }
685            if (fabs(value) > zeroTolerance) {
686              index[numberNonZero++] = iColumn;
687              array[iColumn] = value;
688            }
689          }
690        } else {
691          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
692            double value = 0.0;
693            CoinBigIndex j;
694            for (j = columnStart[iColumn];
695                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
696              int iRow = row[j];
697              value += pi[iRow] * elementByColumn[j];
698            }
699            value *= scalar;
700            if (fabs(value) > zeroTolerance) {
701              index[numberNonZero++] = iColumn;
702              array[iColumn] = value;
703            }
704          }
705        }
706      } else {
707        // scaled
708        if (scalar == -1.0) {
709          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
710            double value = 0.0;
711            CoinBigIndex j;
712            const double *columnScale = model->columnScale();
713            for (j = columnStart[iColumn];
714                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
715              int iRow = row[j];
716              value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
717            }
718            value *= columnScale[iColumn];
719            if (fabs(value) > zeroTolerance) {
720              index[numberNonZero++] = iColumn;
721              array[iColumn] = -value;
722            }
723          }
724        } else if (scalar == 1.0) {
725          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
726            double value = 0.0;
727            CoinBigIndex j;
728            const double *columnScale = model->columnScale();
729            for (j = columnStart[iColumn];
730                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
731              int iRow = row[j];
732              value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
733            }
734            value *= columnScale[iColumn];
735            if (fabs(value) > zeroTolerance) {
736              index[numberNonZero++] = iColumn;
737              array[iColumn] = value;
738            }
739          }
740        } else {
741          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
742            double value = 0.0;
743            CoinBigIndex j;
744            const double *columnScale = model->columnScale();
745            for (j = columnStart[iColumn];
746                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
747              int iRow = row[j];
748              value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
749            }
750            value *= scalar * columnScale[iColumn];
751            if (fabs(value) > zeroTolerance) {
752              index[numberNonZero++] = iColumn;
753              array[iColumn] = value;
754            }
755          }
756        }
757      }
758    }
759    columnArray->setNumElements(numberNonZero);
760    y->setNumElements(0);
761  } else {
762    // do by row
763    transposeTimesByRow(model, scalar, rowArray, y, columnArray);
764  }
765  if (packed)
766    columnArray->setPackedMode(true);
767  if (0) {
768    columnArray->checkClean();
769    int numberNonZero = columnArray->getNumElements();
770    ;
771    int *index = columnArray->getIndices();
772    double *array = columnArray->denseVector();
773    int i;
774    for (i = 0; i < numberNonZero; i++) {
775      int j = index[i];
776      double value;
777      if (packed)
778        value = array[i];
779      else
780        value = array[j];
781      printf("Ti %d %d %g\n", i, j, value);
782    }
783  }
784}
785/* Return <code>x * A + y</code> in <code>z</code>.
786        Squashes small elements and knows about ClpSimplex */
787void ClpGubMatrix::transposeTimesByRow(const ClpSimplex *model, double scalar,
788  const CoinIndexedVector *rowArray,
789  CoinIndexedVector *y,
790  CoinIndexedVector *columnArray) const
791{
792  // Do packed part
793  ClpPackedMatrix::transposeTimesByRow(model, scalar, rowArray, y, columnArray);
794  if (numberSets_) {
795    /* what we need to do is do by row as normal but get list of sets touched
796             and then update those ones */
797    abort();
798  }
799}
800/* Return <code>x *A in <code>z</code> but
801   just for indices in y. */
802void ClpGubMatrix::subsetTransposeTimes(const ClpSimplex *model,
803  const CoinIndexedVector *rowArray,
804  const CoinIndexedVector *y,
805  CoinIndexedVector *columnArray) const
806{
807  columnArray->clear();
808  double *pi = rowArray->denseVector();
809  double *array = columnArray->denseVector();
810  int jColumn;
811  // get matrix data pointers
812  const int *row = matrix_->getIndices();
813  const CoinBigIndex *columnStart = matrix_->getVectorStarts();
814  const int *columnLength = matrix_->getVectorLengths();
815  const double *elementByColumn = matrix_->getElements();
816  const double *rowScale = model->rowScale();
817  int numberToDo = y->getNumElements();
818  const int *which = y->getIndices();
819  assert(!rowArray->packedMode());
820  columnArray->setPacked();
821  int numberTouched = 0;
822  if (!rowScale) {
823    for (jColumn = 0; jColumn < numberToDo; jColumn++) {
824      int iColumn = which[jColumn];
825      double value = 0.0;
826      CoinBigIndex j;
827      for (j = columnStart[iColumn];
828           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
829        int iRow = row[j];
830        value += pi[iRow] * elementByColumn[j];
831      }
832      array[jColumn] = value;
833      if (value) {
834        int iSet = backward_[iColumn];
835        if (iSet >= 0) {
836          int iBasic = keyVariable_[iSet];
837          if (iBasic == iColumn) {
838            toIndex_[iSet] = jColumn;
839            fromIndex_[numberTouched++] = iSet;
840          }
841        }
842      }
843    }
844  } else {
845    // scaled
846    for (jColumn = 0; jColumn < numberToDo; jColumn++) {
847      int iColumn = which[jColumn];
848      double value = 0.0;
849      CoinBigIndex j;
850      const double *columnScale = model->columnScale();
851      for (j = columnStart[iColumn];
852           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
853        int iRow = row[j];
854        value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
855      }
856      value *= columnScale[iColumn];
857      array[jColumn] = value;
858      if (value) {
859        int iSet = backward_[iColumn];
860        if (iSet >= 0) {
861          int iBasic = keyVariable_[iSet];
862          if (iBasic == iColumn) {
863            toIndex_[iSet] = jColumn;
864            fromIndex_[numberTouched++] = iSet;
865          }
866        }
867      }
868    }
869  }
870  // adjust djs
871  for (jColumn = 0; jColumn < numberToDo; jColumn++) {
872    int iColumn = which[jColumn];
873    int iSet = backward_[iColumn];
874    if (iSet >= 0) {
875      int kColumn = toIndex_[iSet];
876      if (kColumn >= 0)
877        array[jColumn] -= array[kColumn];
878    }
879  }
880  // and clear basic
881  for (int j = 0; j < numberTouched; j++) {
882    int iSet = fromIndex_[j];
883    int kColumn = toIndex_[iSet];
884    toIndex_[iSet] = -1;
885    array[kColumn] = 0.0;
886  }
887}
888/// returns number of elements in column part of basis,
889int ClpGubMatrix::countBasis(const int *whichColumn,
890  int &numberColumnBasic)
891{
892  int i;
893  int numberColumns = getNumCols();
894  const int *columnLength = matrix_->getVectorLengths();
895  int numberRows = getNumRows();
896  int numberBasic = 0;
897  CoinBigIndex numberElements = 0;
898  int lastSet = -1;
899  int key = -1;
900  int keyLength = -1;
901  double *work = new double[numberRows];
902  CoinZeroN(work, numberRows);
903  char *mark = new char[numberRows];
904  CoinZeroN(mark, numberRows);
905  const CoinBigIndex *columnStart = matrix_->getVectorStarts();
906  const int *row = matrix_->getIndices();
907  const double *elementByColumn = matrix_->getElements();
908  //ClpGubDynamicMatrix* gubx =
909  //dynamic_cast< ClpGubDynamicMatrix*>(this);
910  //int * id = gubx->id();
911  // just count
912  for (i = 0; i < numberColumnBasic; i++) {
913    int iColumn = whichColumn[i];
914    int iSet = backward_[iColumn];
915    int length = columnLength[iColumn];
916    if (iSet < 0 || keyVariable_[iSet] >= numberColumns) {
917      numberElements += length;
918      numberBasic++;
919      //printf("non gub - set %d id %d (column %d) nel %d\n",iSet,id[iColumn-20],iColumn,length);
920    } else {
921      // in gub set
922      if (iColumn != keyVariable_[iSet]) {
923        numberBasic++;
924        CoinBigIndex j;
925        // not key
926        if (lastSet < iSet) {
927          // erase work
928          if (key >= 0) {
929            for (j = columnStart[key]; j < columnStart[key] + keyLength; j++)
930              work[row[j]] = 0.0;
931          }
932          key = keyVariable_[iSet];
933          lastSet = iSet;
934          keyLength = columnLength[key];
935          for (j = columnStart[key]; j < columnStart[key] + keyLength; j++)
936            work[row[j]] = elementByColumn[j];
937        }
938        int extra = keyLength;
939        for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) {
940          int iRow = row[j];
941          double keyValue = work[iRow];
942          double value = elementByColumn[j];
943          if (!keyValue) {
944            if (fabs(value) > 1.0e-20)
945              extra++;
946          } else {
947            value -= keyValue;
948            if (fabs(value) <= 1.0e-20)
949              extra--;
950          }
951        }
952        numberElements += extra;
953        //printf("gub - set %d id %d (column %d) nel %d\n",iSet,id[iColumn-20],iColumn,extra);
954      }
955    }
956  }
957  delete[] work;
958  delete[] mark;
959  // update number of column basic
960  numberColumnBasic = numberBasic;
961#if COIN_BIG_INDEX
962  if (numberElements > COIN_INT_MAX) {
963    printf("Factorization too large\n");
964    abort();
965  }
966#endif
967  return static_cast< int >(numberElements);
968}
969void ClpGubMatrix::fillBasis(ClpSimplex *model,
970  const int *whichColumn,
971  int &numberColumnBasic,
972  int *indexRowU, int *start,
973  int *rowCount, int *columnCount,
974  CoinFactorizationDouble *elementU)
975{
976  int i;
977  int numberColumns = getNumCols();
978  const int *columnLength = matrix_->getVectorLengths();
979  int numberRows = getNumRows();
980  assert(next_ || !elementU);
981  CoinBigIndex numberElements = start[0];
982  int lastSet = -1;
983  int key = -1;
984  int keyLength = -1;
985  double *work = new double[numberRows];
986  CoinZeroN(work, numberRows);
987  char *mark = new char[numberRows];
988  CoinZeroN(mark, numberRows);
989  const CoinBigIndex *columnStart = matrix_->getVectorStarts();
990  const int *row = matrix_->getIndices();
991  const double *elementByColumn = matrix_->getElements();
992  const double *rowScale = model->rowScale();
993  int numberBasic = 0;
994  if (0) {
995    printf("%d basiccolumns\n", numberColumnBasic);
996    int i;
997    for (i = 0; i < numberSets_; i++) {
998      int k = keyVariable_[i];
999      if (k < numberColumns) {
1000        printf("key %d on set %d, %d elements\n", k, i, columnStart[k + 1] - columnStart[k]);
1001        for (int j = columnStart[k]; j < columnStart[k + 1]; j++)
1002          printf("row %d el %g\n", row[j], elementByColumn[j]);
1003      } else {
1004        printf("slack key on set %d\n", i);
1005      }
1006    }
1007  }
1008  // fill
1009  if (!rowScale) {
1010    // no scaling
1011    for (i = 0; i < numberColumnBasic; i++) {
1012      int iColumn = whichColumn[i];
1013      int iSet = backward_[iColumn];
1014      int length = columnLength[iColumn];
1015      if (0) {
1016        int k = iColumn;
1017        printf("column %d in set %d, %d elements\n", k, iSet, columnStart[k + 1] - columnStart[k]);
1018        for (int j = columnStart[k]; j < columnStart[k + 1]; j++)
1019          printf("row %d el %g\n", row[j], elementByColumn[j]);
1020      }
1021      CoinBigIndex j;
1022      if (iSet < 0 || keyVariable_[iSet] >= numberColumns) {
1023        for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1024          double value = elementByColumn[j];
1025          if (fabs(value) > 1.0e-20) {
1026            int iRow = row[j];
1027            indexRowU[numberElements] = iRow;
1028            rowCount[iRow]++;
1029            elementU[numberElements++] = value;
1030          }
1031        }
1032        // end of column
1033        columnCount[numberBasic] = numberElements - start[numberBasic];
1034        numberBasic++;
1035        start[numberBasic] = numberElements;
1036      } else {
1037        // in gub set
1038        if (iColumn != keyVariable_[iSet]) {
1039          // not key
1040          if (lastSet != iSet) {
1041            // erase work
1042            if (key >= 0) {
1043              for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
1044                int iRow = row[j];
1045                work[iRow] = 0.0;
1046                mark[iRow] = 0;
1047              }
1048            }
1049            key = keyVariable_[iSet];
1050            lastSet = iSet;
1051            keyLength = columnLength[key];
1052            for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
1053              int iRow = row[j];
1054              work[iRow] = elementByColumn[j];
1055              mark[iRow] = 1;
1056            }
1057          }
1058          for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) {
1059            int iRow = row[j];
1060            double value = elementByColumn[j];
1061            if (mark[iRow]) {
1062              mark[iRow] = 0;
1063              double keyValue = work[iRow];
1064              value -= keyValue;
1065            }
1066            if (fabs(value) > 1.0e-20) {
1067              indexRowU[numberElements] = iRow;
1068              rowCount[iRow]++;
1069              elementU[numberElements++] = value;
1070            }
1071          }
1072          for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
1073            int iRow = row[j];
1074            if (mark[iRow]) {
1075              double value = -work[iRow];
1076              if (fabs(value) > 1.0e-20) {
1077                indexRowU[numberElements] = iRow;
1078                rowCount[iRow]++;
1079                elementU[numberElements++] = value;
1080              }
1081            } else {
1082              // just put back mark
1083              mark[iRow] = 1;
1084            }
1085          }
1086          // end of column
1087          columnCount[numberBasic] = numberElements - start[numberBasic];
1088          numberBasic++;
1089          start[numberBasic] = numberElements;
1090        }
1091      }
1092    }
1093  } else {
1094    // scaling
1095    const double *columnScale = model->columnScale();
1096    for (i = 0; i < numberColumnBasic; i++) {
1097      int iColumn = whichColumn[i];
1098      int iSet = backward_[iColumn];
1099      int length = columnLength[iColumn];
1100      CoinBigIndex j;
1101      if (iSet < 0 || keyVariable_[iSet] >= numberColumns) {
1102        double scale = columnScale[iColumn];
1103        for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1104          int iRow = row[j];
1105          double value = elementByColumn[j] * scale * rowScale[iRow];
1106          if (fabs(value) > 1.0e-20) {
1107            indexRowU[numberElements] = iRow;
1108            rowCount[iRow]++;
1109            elementU[numberElements++] = value;
1110          }
1111        }
1112        // end of column
1113        columnCount[numberBasic] = numberElements - start[numberBasic];
1114        numberBasic++;
1115        start[numberBasic] = numberElements;
1116      } else {
1117        // in gub set
1118        if (iColumn != keyVariable_[iSet]) {
1119          double scale = columnScale[iColumn];
1120          // not key
1121          if (lastSet < iSet) {
1122            // erase work
1123            if (key >= 0) {
1124              for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
1125                int iRow = row[j];
1126                work[iRow] = 0.0;
1127                mark[iRow] = 0;
1128              }
1129            }
1130            key = keyVariable_[iSet];
1131            lastSet = iSet;
1132            keyLength = columnLength[key];
1133            double scale = columnScale[key];
1134            for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
1135              int iRow = row[j];
1136              work[iRow] = elementByColumn[j] * scale * rowScale[iRow];
1137              mark[iRow] = 1;
1138            }
1139          }
1140          for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) {
1141            int iRow = row[j];
1142            double value = elementByColumn[j] * scale * rowScale[iRow];
1143            if (mark[iRow]) {
1144              mark[iRow] = 0;
1145              double keyValue = work[iRow];
1146              value -= keyValue;
1147            }
1148            if (fabs(value) > 1.0e-20) {
1149              indexRowU[numberElements] = iRow;
1150              rowCount[iRow]++;
1151              elementU[numberElements++] = value;
1152            }
1153          }
1154          for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
1155            int iRow = row[j];
1156            if (mark[iRow]) {
1157              double value = -work[iRow];
1158              if (fabs(value) > 1.0e-20) {
1159                indexRowU[numberElements] = iRow;
1160                rowCount[iRow]++;
1161                elementU[numberElements++] = value;
1162              }
1163            } else {
1164              // just put back mark
1165              mark[iRow] = 1;
1166            }
1167          }
1168          // end of column
1169          columnCount[numberBasic] = numberElements - start[numberBasic];
1170          numberBasic++;
1171          start[numberBasic] = numberElements;
1172        }
1173      }
1174    }
1175  }
1176  delete[] work;
1177  delete[] mark;
1178  // update number of column basic
1179  numberColumnBasic = numberBasic;
1180}
1181/* Unpacks a column into an CoinIndexedvector
1182 */
1183void ClpGubMatrix::unpack(const ClpSimplex *model, CoinIndexedVector *rowArray,
1184  int iColumn) const
1185{
1186  assert(iColumn < model->numberColumns());
1187  // Do packed part
1188  ClpPackedMatrix::unpack(model, rowArray, iColumn);
1189  int iSet = backward_[iColumn];
1190  if (iSet >= 0) {
1191    int iBasic = keyVariable_[iSet];
1192    if (iBasic < model->numberColumns()) {
1193      add(model, rowArray, iBasic, -1.0);
1194    }
1195  }
1196}
1197/* Unpacks a column into a CoinIndexedVector
1198** in packed format
1199Note that model is NOT const.  Bounds and objective could
1200be modified if doing column generation (just for this variable) */
1201void ClpGubMatrix::unpackPacked(ClpSimplex *model,
1202  CoinIndexedVector *rowArray,
1203  int iColumn) const
1204{
1205  int numberColumns = model->numberColumns();
1206  if (iColumn < numberColumns) {
1207    // Do packed part
1208    ClpPackedMatrix::unpackPacked(model, rowArray, iColumn);
1209    int iSet = backward_[iColumn];
1210    if (iSet >= 0) {
1211      // columns are in order
1212      int iBasic = keyVariable_[iSet];
1213      if (iBasic < numberColumns) {
1214        int number = rowArray->getNumElements();
1215        const double *rowScale = model->rowScale();
1216        const int *row = matrix_->getIndices();
1217        const CoinBigIndex *columnStart = matrix_->getVectorStarts();
1218        const int *columnLength = matrix_->getVectorLengths();
1219        const double *elementByColumn = matrix_->getElements();
1220        double *array = rowArray->denseVector();
1221        int *index = rowArray->getIndices();
1222        CoinBigIndex i;
1223        int numberOld = number;
1224        int lastIndex = 0;
1225        int next = index[lastIndex];
1226        if (!rowScale) {
1227          for (i = columnStart[iBasic];
1228               i < columnStart[iBasic] + columnLength[iBasic]; i++) {
1229            int iRow = row[i];
1230            while (iRow > next) {
1231              lastIndex++;
1232              if (lastIndex == numberOld)
1233                next = matrix_->getNumRows();
1234              else
1235                next = index[lastIndex];
1236            }
1237            if (iRow < next) {
1238              array[number] = -elementByColumn[i];
1239              index[number++] = iRow;
1240            } else {
1241              assert(iRow == next);
1242              array[lastIndex] -= elementByColumn[i];
1243              if (!array[lastIndex])
1244                array[lastIndex] = 1.0e-100;
1245            }
1246          }
1247        } else {
1248          // apply scaling
1249          double scale = model->columnScale()[iBasic];
1250          for (i = columnStart[iBasic];
1251               i < columnStart[iBasic] + columnLength[iBasic]; i++) {
1252            int iRow = row[i];
1253            while (iRow > next) {
1254              lastIndex++;
1255              if (lastIndex == numberOld)
1256                next = matrix_->getNumRows();
1257              else
1258                next = index[lastIndex];
1259            }
1260            if (iRow < next) {
1261              array[number] = -elementByColumn[i] * scale * rowScale[iRow];
1262              index[number++] = iRow;
1263            } else {
1264              assert(iRow == next);
1265              array[lastIndex] -= elementByColumn[i] * scale * rowScale[iRow];
1266              if (!array[lastIndex])
1267                array[lastIndex] = 1.0e-100;
1268            }
1269          }
1270        }
1271        rowArray->setNumElements(number);
1272      }
1273    }
1274  } else {
1275    // key slack entering
1276    int iBasic = keyVariable_[gubSlackIn_];
1277    assert(iBasic < numberColumns);
1278    int number = 0;
1279    const double *rowScale = model->rowScale();
1280    const int *row = matrix_->getIndices();
1281    const CoinBigIndex *columnStart = matrix_->getVectorStarts();
1282    const int *columnLength = matrix_->getVectorLengths();
1283    const double *elementByColumn = matrix_->getElements();
1284    double *array = rowArray->denseVector();
1285    int *index = rowArray->getIndices();
1286    CoinBigIndex i;
1287    if (!rowScale) {
1288      for (i = columnStart[iBasic];
1289           i < columnStart[iBasic] + columnLength[iBasic]; i++) {
1290        int iRow = row[i];
1291        array[number] = elementByColumn[i];
1292        index[number++] = iRow;
1293      }
1294    } else {
1295      // apply scaling
1296      double scale = model->columnScale()[iBasic];
1297      for (i = columnStart[iBasic];
1298           i < columnStart[iBasic] + columnLength[iBasic]; i++) {
1299        int iRow = row[i];
1300        array[number] = elementByColumn[i] * scale * rowScale[iRow];
1301        index[number++] = iRow;
1302      }
1303    }
1304    rowArray->setNumElements(number);
1305    rowArray->setPacked();
1306  }
1307}
1308/* Adds multiple of a column into an CoinIndexedvector
1309      You can use quickAdd to add to vector */
1310void ClpGubMatrix::add(const ClpSimplex *model, CoinIndexedVector *rowArray,
1311  int iColumn, double multiplier) const
1312{
1313  assert(iColumn < model->numberColumns());
1314  // Do packed part
1315  ClpPackedMatrix::add(model, rowArray, iColumn, multiplier);
1316  int iSet = backward_[iColumn];
1317  if (iSet >= 0 && iColumn != keyVariable_[iSet]) {
1318    ClpPackedMatrix::add(model, rowArray, keyVariable_[iSet], -multiplier);
1319  }
1320}
1321/* Adds multiple of a column into an array */
1322void ClpGubMatrix::add(const ClpSimplex *model, double *array,
1323  int iColumn, double multiplier) const
1324{
1325  assert(iColumn < model->numberColumns());
1326  // Do packed part
1327  ClpPackedMatrix::add(model, array, iColumn, multiplier);
1328  if (iColumn < model->numberColumns()) {
1329    int iSet = backward_[iColumn];
1330    if (iSet >= 0 && iColumn != keyVariable_[iSet] && keyVariable_[iSet] < model->numberColumns()) {
1331      ClpPackedMatrix::add(model, array, keyVariable_[iSet], -multiplier);
1332    }
1333  }
1334}
1335// Partial pricing
1336void ClpGubMatrix::partialPricing(ClpSimplex *model, double startFraction, double endFraction,
1337  int &bestSequence, int &numberWanted)
1338{
1339  numberWanted = currentWanted_;
1340  if (numberSets_) {
1341    // Do packed part before gub
1342    int numberColumns = matrix_->getNumCols();
1343    double ratio = static_cast< double >(firstGub_) / static_cast< double >(numberColumns);
1344    ClpPackedMatrix::partialPricing(model, startFraction * ratio,
1345      endFraction * ratio, bestSequence, numberWanted);
1346    if (numberWanted || minimumGoodReducedCosts_ < -1) {
1347      // do gub
1348      const double *element = matrix_->getElements();
1349      const int *row = matrix_->getIndices();
1350      const CoinBigIndex *startColumn = matrix_->getVectorStarts();
1351      const int *length = matrix_->getVectorLengths();
1352      const double *rowScale = model->rowScale();
1353      const double *columnScale = model->columnScale();
1354      int iSequence;
1355      CoinBigIndex j;
1356      double tolerance = model->currentDualTolerance();
1357      double *reducedCost = model->djRegion();
1358      const double *duals = model->dualRowSolution();
1359      const double *cost = model->costRegion();
1360      double bestDj;
1361      int numberColumns = model->numberColumns();
1362      int numberRows = model->numberRows();
1363      if (bestSequence >= 0)
1364        bestDj = fabs(this->reducedCost(model, bestSequence));
1365      else
1366        bestDj = tolerance;
1367      int sequenceOut = model->sequenceOut();
1368      int saveSequence = bestSequence;
1369      int startG = firstGub_ + static_cast< int >(startFraction * (lastGub_ - firstGub_));
1370      int endG = firstGub_ + static_cast< int >(endFraction * (lastGub_ - firstGub_));
1371      endG = CoinMin(lastGub_, endG + 1);
1372      // If nothing found yet can go all the way to end
1373      int endAll = endG;
1374      if (bestSequence < 0 && !startG)
1375        endAll = lastGub_;
1376      int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_;
1377      int minNeg = minimumGoodReducedCosts_ == -1 ? 5 : minimumGoodReducedCosts_;
1378      int nSets = 0;
1379      int iSet = -1;
1380      double djMod = 0.0;
1381      double infeasibilityCost = model->infeasibilityCost();
1382      if (rowScale) {
1383        double bestDjMod = 0.0;
1384        // scaled
1385        for (iSequence = startG; iSequence < endAll; iSequence++) {
1386          if (numberWanted + minNeg < originalWanted_ && nSets > minSet) {
1387            // give up
1388            numberWanted = 0;
1389            break;
1390          } else if (iSequence == endG && bestSequence >= 0) {
1391            break;
1392          }
1393          if (backward_[iSequence] != iSet) {
1394            // get pi on gub row
1395            iSet = backward_[iSequence];
1396            if (iSet >= 0) {
1397              nSets++;
1398              int iBasic = keyVariable_[iSet];
1399              if (iBasic >= numberColumns) {
1400                djMod = -weight(iSet) * infeasibilityCost;
1401              } else {
1402                // get dj without
1403                assert(model->getStatus(iBasic) == ClpSimplex::basic);
1404                djMod = 0.0;
1405                // scaled
1406                for (j = startColumn[iBasic];
1407                     j < startColumn[iBasic] + length[iBasic]; j++) {
1408                  int jRow = row[j];
1409                  djMod -= duals[jRow] * element[j] * rowScale[jRow];
1410                }
1411                // allow for scaling
1412                djMod += cost[iBasic] / columnScale[iBasic];
1413                // See if gub slack possible - dj is djMod
1414                if (getStatus(iSet) == ClpSimplex::atLowerBound) {
1415                  double value = -djMod;
1416                  if (value > tolerance) {
1417                    numberWanted--;
1418                    if (value > bestDj) {
1419                      // check flagged variable and correct dj
1420                      if (!flagged(iSet)) {
1421                        bestDj = value;
1422                        bestSequence = numberRows + numberColumns + iSet;
1423                        bestDjMod = djMod;
1424                      } else {
1425                        // just to make sure we don't exit before got something
1426                        numberWanted++;
1427                        abort();
1428                      }
1429                    }
1430                  }
1431                } else if (getStatus(iSet) == ClpSimplex::atUpperBound) {
1432                  double value = djMod;
1433                  if (value > tolerance) {
1434                    numberWanted--;
1435                    if (value > bestDj) {
1436                      // check flagged variable and correct dj
1437                      if (!flagged(iSet)) {
1438                        bestDj = value;
1439                        bestSequence = numberRows + numberColumns + iSet;
1440                        bestDjMod = djMod;
1441                      } else {
1442                        // just to make sure we don't exit before got something
1443                        numberWanted++;
1444                        abort();
1445                      }
1446                    }
1447                  }
1448                }
1449              }
1450            } else {
1451              // not in set
1452              djMod = 0.0;
1453            }
1454          }
1455          if (iSequence != sequenceOut) {
1456            double value;
1457            ClpSimplex::Status status = model->getStatus(iSequence);
1458
1459            switch (status) {
1460
1461            case ClpSimplex::basic:
1462            case ClpSimplex::isFixed:
1463              break;
1464            case ClpSimplex::isFree:
1465            case ClpSimplex::superBasic:
1466              value = -djMod;
1467              // scaled
1468              for (j = startColumn[iSequence];
1469                   j < startColumn[iSequence] + length[iSequence]; j++) {
1470                int jRow = row[j];
1471                value -= duals[jRow] * element[j] * rowScale[jRow];
1472              }
1473              value = fabs(cost[iSequence] + value * columnScale[iSequence]);
1474              if (value > FREE_ACCEPT * tolerance) {
1475                numberWanted--;
1476                // we are going to bias towards free (but only if reasonable)
1477                value *= FREE_BIAS;
1478                if (value > bestDj) {
1479                  // check flagged variable and correct dj
1480                  if (!model->flagged(iSequence)) {
1481                    bestDj = value;
1482                    bestSequence = iSequence;
1483                    bestDjMod = djMod;
1484                  } else {
1485                    // just to make sure we don't exit before got something
1486                    numberWanted++;
1487                  }
1488                }
1489              }
1490              break;
1491            case ClpSimplex::atUpperBound:
1492              value = -djMod;
1493              // scaled
1494              for (j = startColumn[iSequence];
1495                   j < startColumn[iSequence] + length[iSequence]; j++) {
1496                int jRow = row[j];
1497                value -= duals[jRow] * element[j] * rowScale[jRow];
1498              }
1499              value = cost[iSequence] + value * columnScale[iSequence];
1500              if (value > tolerance) {
1501                numberWanted--;
1502                if (value > bestDj) {
1503                  // check flagged variable and correct dj
1504                  if (!model->flagged(iSequence)) {
1505                    bestDj = value;
1506                    bestSequence = iSequence;
1507                    bestDjMod = djMod;
1508                  } else {
1509                    // just to make sure we don't exit before got something
1510                    numberWanted++;
1511                  }
1512                }
1513              }
1514              break;
1515            case ClpSimplex::atLowerBound:
1516              value = -djMod;
1517              // scaled
1518              for (j = startColumn[iSequence];
1519                   j < startColumn[iSequence] + length[iSequence]; j++) {
1520                int jRow = row[j];
1521                value -= duals[jRow] * element[j] * rowScale[jRow];
1522              }
1523              value = -(cost[iSequence] + value * columnScale[iSequence]);
1524              if (value > tolerance) {
1525                numberWanted--;
1526                if (value > bestDj) {
1527                  // check flagged variable and correct dj
1528                  if (!model->flagged(iSequence)) {
1529                    bestDj = value;
1530                    bestSequence = iSequence;
1531                    bestDjMod = djMod;
1532                  } else {
1533                    // just to make sure we don't exit before got something
1534                    numberWanted++;
1535                  }
1536                }
1537              }
1538              break;
1539            }
1540          }
1541          if (!numberWanted)
1542            break;
1543        }
1544        if (bestSequence != saveSequence) {
1545          if (bestSequence < numberRows + numberColumns) {
1546            // recompute dj
1547            double value = bestDjMod;
1548            // scaled
1549            for (j = startColumn[bestSequence];
1550                 j < startColumn[bestSequence] + length[bestSequence]; j++) {
1551              int jRow = row[j];
1552              value -= duals[jRow] * element[j] * rowScale[jRow];
1553            }
1554            reducedCost[bestSequence] = cost[bestSequence] + value * columnScale[bestSequence];
1555            gubSlackIn_ = -1;
1556          } else {
1557            // slack - make last column
1558            gubSlackIn_ = bestSequence - numberRows - numberColumns;
1559            bestSequence = numberColumns + 2 * numberRows;
1560            reducedCost[bestSequence] = bestDjMod;
1561            model->setStatus(bestSequence, getStatus(gubSlackIn_));
1562            if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound)
1563              model->solutionRegion()[bestSequence] = upper_[gubSlackIn_];
1564            else
1565              model->solutionRegion()[bestSequence] = lower_[gubSlackIn_];
1566            model->lowerRegion()[bestSequence] = lower_[gubSlackIn_];
1567            model->upperRegion()[bestSequence] = upper_[gubSlackIn_];
1568            model->costRegion()[bestSequence] = 0.0;
1569          }
1570          savedBestSequence_ = bestSequence;
1571          savedBestDj_ = reducedCost[savedBestSequence_];
1572        }
1573      } else {
1574        double bestDjMod = 0.0;
1575        //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
1576        //     startG,endG,numberWanted);
1577        for (iSequence = startG; iSequence < endG; iSequence++) {
1578          if (numberWanted + minNeg < originalWanted_ && nSets > minSet) {
1579            // give up
1580            numberWanted = 0;
1581            break;
1582          } else if (iSequence == endG && bestSequence >= 0) {
1583            break;
1584          }
1585          if (backward_[iSequence] != iSet) {
1586            // get pi on gub row
1587            iSet = backward_[iSequence];
1588            if (iSet >= 0) {
1589              nSets++;
1590              int iBasic = keyVariable_[iSet];
1591              if (iBasic >= numberColumns) {
1592                djMod = -weight(iSet) * infeasibilityCost;
1593              } else {
1594                // get dj without
1595                assert(model->getStatus(iBasic) == ClpSimplex::basic);
1596                djMod = 0.0;
1597
1598                for (j = startColumn[iBasic];
1599                     j < startColumn[iBasic] + length[iBasic]; j++) {
1600                  int jRow = row[j];
1601                  djMod -= duals[jRow] * element[j];
1602                }
1603                djMod += cost[iBasic];
1604                // See if gub slack possible - dj is djMod
1605                if (getStatus(iSet) == ClpSimplex::atLowerBound) {
1606                  double value = -djMod;
1607                  if (value > tolerance) {
1608                    numberWanted--;
1609                    if (value > bestDj) {
1610                      // check flagged variable and correct dj
1611                      if (!flagged(iSet)) {
1612                        bestDj = value;
1613                        bestSequence = numberRows + numberColumns + iSet;
1614                        bestDjMod = djMod;
1615                      } else {
1616                        // just to make sure we don't exit before got something
1617                        numberWanted++;
1618                        abort();
1619                      }
1620                    }
1621                  }
1622                } else if (getStatus(iSet) == ClpSimplex::atUpperBound) {
1623                  double value = djMod;
1624                  if (value > tolerance) {
1625                    numberWanted--;
1626                    if (value > bestDj) {
1627                      // check flagged variable and correct dj
1628                      if (!flagged(iSet)) {
1629                        bestDj = value;
1630                        bestSequence = numberRows + numberColumns + iSet;
1631                        bestDjMod = djMod;
1632                      } else {
1633                        // just to make sure we don't exit before got something
1634                        numberWanted++;
1635                        abort();
1636                      }
1637                    }
1638                  }
1639                }
1640              }
1641            } else {
1642              // not in set
1643              djMod = 0.0;
1644            }
1645          }
1646          if (iSequence != sequenceOut) {
1647            double value;
1648            ClpSimplex::Status status = model->getStatus(iSequence);
1649
1650            switch (status) {
1651
1652            case ClpSimplex::basic:
1653            case ClpSimplex::isFixed:
1654              break;
1655            case ClpSimplex::isFree:
1656            case ClpSimplex::superBasic:
1657              value = cost[iSequence] - djMod;
1658              for (j = startColumn[iSequence];
1659                   j < startColumn[iSequence] + length[iSequence]; j++) {
1660                int jRow = row[j];
1661                value -= duals[jRow] * element[j];
1662              }
1663              value = fabs(value);
1664              if (value > FREE_ACCEPT * tolerance) {
1665                numberWanted--;
1666                // we are going to bias towards free (but only if reasonable)
1667                value *= FREE_BIAS;
1668                if (value > bestDj) {
1669                  // check flagged variable and correct dj
1670                  if (!model->flagged(iSequence)) {
1671                    bestDj = value;
1672                    bestSequence = iSequence;
1673                    bestDjMod = djMod;
1674                  } else {
1675                    // just to make sure we don't exit before got something
1676                    numberWanted++;
1677                  }
1678                }
1679              }
1680              break;
1681            case ClpSimplex::atUpperBound:
1682              value = cost[iSequence] - djMod;
1683              for (j = startColumn[iSequence];
1684                   j < startColumn[iSequence] + length[iSequence]; j++) {
1685                int jRow = row[j];
1686                value -= duals[jRow] * element[j];
1687              }
1688              if (value > tolerance) {
1689                numberWanted--;
1690                if (value > bestDj) {
1691                  // check flagged variable and correct dj
1692                  if (!model->flagged(iSequence)) {
1693                    bestDj = value;
1694                    bestSequence = iSequence;
1695                    bestDjMod = djMod;
1696                  } else {
1697                    // just to make sure we don't exit before got something
1698                    numberWanted++;
1699                  }
1700                }
1701              }
1702              break;
1703            case ClpSimplex::atLowerBound:
1704              value = cost[iSequence] - djMod;
1705              for (j = startColumn[iSequence];
1706                   j < startColumn[iSequence] + length[iSequence]; j++) {
1707                int jRow = row[j];
1708                value -= duals[jRow] * element[j];
1709              }
1710              value = -value;
1711              if (value > tolerance) {
1712                numberWanted--;
1713                if (value > bestDj) {
1714                  // check flagged variable and correct dj
1715                  if (!model->flagged(iSequence)) {
1716                    bestDj = value;
1717                    bestSequence = iSequence;
1718                    bestDjMod = djMod;
1719                  } else {
1720                    // just to make sure we don't exit before got something
1721                    numberWanted++;
1722                  }
1723                }
1724              }
1725              break;
1726            }
1727          }
1728          if (!numberWanted)
1729            break;
1730        }
1731        if (bestSequence != saveSequence) {
1732          if (bestSequence < numberRows + numberColumns) {
1733            // recompute dj
1734            double value = cost[bestSequence] - bestDjMod;
1735            for (j = startColumn[bestSequence];
1736                 j < startColumn[bestSequence] + length[bestSequence]; j++) {
1737              int jRow = row[j];
1738              value -= duals[jRow] * element[j];
1739            }
1740            //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,bestDjMod);
1741            reducedCost[bestSequence] = value;
1742            gubSlackIn_ = -1;
1743          } else {
1744            // slack - make last column
1745            gubSlackIn_ = bestSequence - numberRows - numberColumns;
1746            bestSequence = numberColumns + 2 * numberRows;
1747            reducedCost[bestSequence] = bestDjMod;
1748            //printf("price slack %d - gubpi %g\n",gubSlackIn_,bestDjMod);
1749            model->setStatus(bestSequence, getStatus(gubSlackIn_));
1750            if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound)
1751              model->solutionRegion()[bestSequence] = upper_[gubSlackIn_];
1752            else
1753              model->solutionRegion()[bestSequence] = lower_[gubSlackIn_];
1754            model->lowerRegion()[bestSequence] = lower_[gubSlackIn_];
1755            model->upperRegion()[bestSequence] = upper_[gubSlackIn_];
1756            model->costRegion()[bestSequence] = 0.0;
1757          }
1758        }
1759      }
1760      // See if may be finished
1761      if (startG == firstGub_ && bestSequence < 0)
1762        infeasibilityWeight_ = model_->infeasibilityCost();
1763      else if (bestSequence >= 0)
1764        infeasibilityWeight_ = -1.0;
1765    }
1766    if (numberWanted) {
1767      // Do packed part after gub
1768      double offset = static_cast< double >(lastGub_) / static_cast< double >(numberColumns);
1769      double ratio = static_cast< double >(numberColumns) / static_cast< double >(numberColumns) - offset;
1770      double start2 = offset + ratio * startFraction;
1771      double end2 = CoinMin(1.0, offset + ratio * endFraction + 1.0e-6);
1772      ClpPackedMatrix::partialPricing(model, start2, end2, bestSequence, numberWanted);
1773    }
1774  } else {
1775    // no gub
1776    ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
1777  }
1778  if (bestSequence >= 0)
1779    infeasibilityWeight_ = -1.0; // not optimal
1780  currentWanted_ = numberWanted;
1781}
1782/* expands an updated column to allow for extra rows which the main
1783   solver does not know about and returns number added.
1784*/
1785int ClpGubMatrix::extendUpdated(ClpSimplex *model, CoinIndexedVector *update, int mode)
1786{
1787  // I think we only need to bother about sets with two in basis or incoming set
1788  int number = update->getNumElements();
1789  double *array = update->denseVector();
1790  int *index = update->getIndices();
1791  int i;
1792  assert(!number || update->packedMode());
1793  int *pivotVariable = model->pivotVariable();
1794  int numberRows = model->numberRows();
1795  int numberColumns = model->numberColumns();
1796  int numberTotal = numberRows + numberColumns;
1797  int sequenceIn = model->sequenceIn();
1798  int returnCode = 0;
1799  int iSetIn;
1800  if (sequenceIn < numberColumns) {
1801    iSetIn = backward_[sequenceIn];
1802    gubSlackIn_ = -1; // in case set
1803  } else if (sequenceIn < numberRows + numberColumns) {
1804    iSetIn = -1;
1805    gubSlackIn_ = -1; // in case set
1806  } else {
1807    iSetIn = gubSlackIn_;
1808  }
1809  double *lower = model->lowerRegion();
1810  double *upper = model->upperRegion();
1811  double *cost = model->costRegion();
1812  double *solution = model->solutionRegion();
1813  int number2 = number;
1814  if (!mode) {
1815    double primalTolerance = model->primalTolerance();
1816    double infeasibilityCost = model->infeasibilityCost();
1817    // extend
1818    saveNumber_ = number;
1819    for (i = 0; i < number; i++) {
1820      int iRow = index[i];
1821      int iPivot = pivotVariable[iRow];
1822      if (iPivot < numberColumns) {
1823        int iSet = backward_[iPivot];
1824        if (iSet >= 0) {
1825          // two (or more) in set
1826          int iIndex = toIndex_[iSet];
1827          double otherValue = array[i];
1828          double value;
1829          if (iIndex < 0) {
1830            toIndex_[iSet] = number2;
1831            int iNew = number2 - number;
1832            fromIndex_[number2 - number] = iSet;
1833            iIndex = number2;
1834            index[number2] = numberRows + iNew;
1835            // do key stuff
1836            int iKey = keyVariable_[iSet];
1837            if (iKey < numberColumns) {
1838              // Save current cost of key
1839              changeCost_[number2 - number] = cost[iKey];
1840              if (iSet != iSetIn)
1841                value = 0.0;
1842              else if (iSetIn != gubSlackIn_)
1843                value = 1.0;
1844              else
1845                value = -1.0;
1846              pivotVariable[numberRows + iNew] = iKey;
1847              // Do I need to recompute?
1848              double sol;
1849              assert(getStatus(iSet) != ClpSimplex::basic);
1850              if (getStatus(iSet) == ClpSimplex::atLowerBound)
1851                sol = lower_[iSet];
1852              else
1853                sol = upper_[iSet];
1854              if ((gubType_ & 8) != 0) {
1855                int iColumn = next_[iKey];
1856                // sum all non-key variables
1857                while (iColumn >= 0) {
1858                  sol -= solution[iColumn];
1859                  iColumn = next_[iColumn];
1860                }
1861              } else {
1862                int stop = -(iKey + 1);
1863                int iColumn = next_[iKey];
1864                // sum all non-key variables
1865                while (iColumn != stop) {
1866                  if (iColumn < 0)
1867                    iColumn = -iColumn - 1;
1868                  sol -= solution[iColumn];
1869                  iColumn = next_[iColumn];
1870                }
1871              }
1872              solution[iKey] = sol;
1873              if (model->algorithm() > 0)
1874                model->nonLinearCost()->setOne(iKey, sol);
1875              //assert (fabs(sol-solution[iKey])<1.0e-3);
1876            } else {
1877              // gub slack is basic
1878              // Save current cost of key
1879              changeCost_[number2 - number] = -weight(iSet) * infeasibilityCost;
1880              otherValue = -otherValue; //allow for - sign on slack
1881              if (iSet != iSetIn)
1882                value = 0.0;
1883              else
1884                value = -1.0;
1885              pivotVariable[numberRows + iNew] = iNew + numberTotal;
1886              model->djRegion()[iNew + numberTotal] = 0.0;
1887              double sol = 0.0;
1888              if ((gubType_ & 8) != 0) {
1889                int iColumn = next_[iKey];
1890                // sum all non-key variables
1891                while (iColumn >= 0) {
1892                  sol += solution[iColumn];
1893                  iColumn = next_[iColumn];
1894                }
1895              } else {
1896                int stop = -(iKey + 1);
1897                int iColumn = next_[iKey];
1898                // sum all non-key variables
1899                while (iColumn != stop) {
1900                  if (iColumn < 0)
1901                    iColumn = -iColumn - 1;
1902                  sol += solution[iColumn];
1903                  iColumn = next_[iColumn];
1904                }
1905              }
1906              solution[iNew + numberTotal] = sol;
1907              // and do cost in nonLinearCost
1908              if (model->algorithm() > 0)
1909                model->nonLinearCost()->setOne(iNew + numberTotal, sol, lower_[iSet], upper_[iSet]);
1910              if (sol > upper_[iSet] + primalTolerance) {
1911                setAbove(iSet);
1912                lower[iNew + numberTotal] = upper_[iSet];
1913                upper[iNew + numberTotal] = COIN_DBL_MAX;
1914              } else if (sol < lower_[iSet] - primalTolerance) {
1915                setBelow(iSet);
1916                lower[iNew + numberTotal] = -COIN_DBL_MAX;
1917                upper[iNew + numberTotal] = lower_[iSet];
1918              } else {
1919                setFeasible(iSet);
1920                lower[iNew + numberTotal] = lower_[iSet];
1921                upper[iNew + numberTotal] = upper_[iSet];
1922              }
1923              cost[iNew + numberTotal] = weight(iSet) * infeasibilityCost;
1924            }
1925            number2++;
1926          } else {
1927            value = array[iIndex];
1928            int iKey = keyVariable_[iSet];
1929            if (iKey >= numberColumns)
1930              otherValue = -otherValue; //allow for - sign on slack
1931          }
1932          value -= otherValue;
1933          array[iIndex] = value;
1934        }
1935      }
1936    }
1937    if (iSetIn >= 0 && toIndex_[iSetIn] < 0) {
1938      // Do incoming
1939      update->setPacked(); // just in case no elements
1940      toIndex_[iSetIn] = number2;
1941      int iNew = number2 - number;
1942      fromIndex_[number2 - number] = iSetIn;
1943      // Save current cost of key
1944      double currentCost;
1945      int key = keyVariable_[iSetIn];
1946      if (key < numberColumns)
1947        currentCost = cost[key];
1948      else
1949        currentCost = -weight(iSetIn) * infeasibilityCost;
1950      changeCost_[number2 - number] = currentCost;
1951      index[number2] = numberRows + iNew;
1952      // do key stuff
1953      int iKey = keyVariable_[iSetIn];
1954      if (iKey < numberColumns) {
1955        if (gubSlackIn_ < 0)
1956          array[number2] = 1.0;
1957        else
1958          array[number2] = -1.0;
1959        pivotVariable[numberRows + iNew] = iKey;
1960        // Do I need to recompute?
1961        double sol;
1962        assert(getStatus(iSetIn) != ClpSimplex::basic);
1963        if (getStatus(iSetIn) == ClpSimplex::atLowerBound)
1964          sol = lower_[iSetIn];
1965        else
1966          sol = upper_[iSetIn];
1967        if ((gubType_ & 8) != 0) {
1968          int iColumn = next_[iKey];
1969          // sum all non-key variables
1970          while (iColumn >= 0) {
1971            sol -= solution[iColumn];
1972            iColumn = next_[iColumn];
1973          }
1974        } else {
1975          // bounds exist - sum over all except key
1976          int stop = -(iKey + 1);
1977          int iColumn = next_[iKey];
1978          // sum all non-key variables
1979          while (iColumn != stop) {
1980            if (iColumn < 0)
1981              iColumn = -iColumn - 1;
1982            sol -= solution[iColumn];
1983            iColumn = next_[iColumn];
1984          }
1985        }
1986        solution[iKey] = sol;
1987        if (model->algorithm() > 0)
1988          model->nonLinearCost()->setOne(iKey, sol);
1989        //assert (fabs(sol-solution[iKey])<1.0e-3);
1990      } else {
1991        // gub slack is basic
1992        array[number2] = -1.0;
1993        pivotVariable[numberRows + iNew] = iNew + numberTotal;
1994        model->djRegion()[iNew + numberTotal] = 0.0;
1995        double sol = 0.0;
1996        if ((gubType_ & 8) != 0) {
1997          int iColumn = next_[iKey];
1998          // sum all non-key variables
1999          while (iColumn >= 0) {
2000            sol += solution[iColumn];
2001            iColumn = next_[iColumn];
2002          }
2003        } else {
2004          // bounds exist - sum over all except key
2005          int stop = -(iKey + 1);
2006          int iColumn = next_[iKey];
2007          // sum all non-key variables
2008          while (iColumn != stop) {
2009            if (iColumn < 0)
2010              iColumn = -iColumn - 1;
2011            sol += solution[iColumn];
2012            iColumn = next_[iColumn];
2013          }
2014        }
2015        solution[iNew + numberTotal] = sol;
2016        // and do cost in nonLinearCost
2017        if (model->algorithm() > 0)
2018          model->nonLinearCost()->setOne(iNew + numberTotal, sol, lower_[iSetIn], upper_[iSetIn]);
2019        if (sol > upper_[iSetIn] + primalTolerance) {
2020          setAbove(iSetIn);
2021          lower[iNew + numberTotal] = upper_[iSetIn];
2022          upper[iNew + numberTotal] = COIN_DBL_MAX;
2023        } else if (sol < lower_[iSetIn] - primalTolerance) {
2024          setBelow(iSetIn);
2025          lower[iNew + numberTotal] = -COIN_DBL_MAX;
2026          upper[iNew + numberTotal] = lower_[iSetIn];
2027        } else {
2028          setFeasible(iSetIn);
2029          lower[iNew + numberTotal] = lower_[iSetIn];
2030          upper[iNew + numberTotal] = upper_[iSetIn];
2031        }
2032        cost[iNew + numberTotal] = weight(iSetIn) * infeasibilityCost;
2033      }
2034      number2++;
2035    }
2036    // mark end
2037    fromIndex_[number2 - number] = -1;
2038    returnCode = number2 - number;
2039    // make sure lower_ upper_ adjusted
2040    synchronize(model, 9);
2041  } else {
2042    // take off?
2043    if (number > saveNumber_) {
2044      // clear
2045      double theta = model->theta();
2046      double *solution = model->solutionRegion();
2047      for (i = saveNumber_; i < number; i++) {
2048        int iRow = index[i];
2049        int iColumn = pivotVariable[iRow];
2050#ifdef CLP_DEBUG_PRINT
2051        printf("Column %d (set %d) lower %g, upper %g - alpha %g - old value %g, new %g (theta %g)\n",
2052          iColumn, fromIndex_[i - saveNumber_], lower[iColumn], upper[iColumn], array[i],
2053          solution[iColumn], solution[iColumn] - model->theta() * array[i], model->theta());
2054#endif
2055        double value = array[i];
2056        array[i] = 0.0;
2057        int iSet = fromIndex_[i - saveNumber_];
2058        toIndex_[iSet] = -1;
2059        if (iSet == iSetIn && iColumn < numberColumns) {
2060          // update as may need value
2061          solution[iColumn] -= theta * value;
2062        }
2063      }
2064    }
2065#ifdef CLP_DEBUG
2066    for (i = 0; i < numberSets_; i++)
2067      assert(toIndex_[i] == -1);
2068#endif
2069    number2 = saveNumber_;
2070  }
2071  update->setNumElements(number2);
2072  return returnCode;
2073}
2074/*
2075     utility primal function for dealing with dynamic constraints
2076     mode=n see ClpGubMatrix.hpp for definition
2077     Remember to update here when settled down
2078*/
2079void ClpGubMatrix::primalExpanded(ClpSimplex *model, int mode)
2080{
2081  int numberColumns = model->numberColumns();
2082  switch (mode) {
2083    // If key variable then slot in gub rhs so will get correct contribution
2084  case 0: {
2085    int i;
2086    double *solution = model->solutionRegion();
2087    ClpSimplex::Status iStatus;
2088    for (i = 0; i < numberSets_; i++) {
2089      int iColumn = keyVariable_[i];
2090      if (iColumn < numberColumns) {
2091        // key is structural - where is slack
2092        iStatus = getStatus(i);
2093        assert(iStatus != ClpSimplex::basic);
2094        if (iStatus == ClpSimplex::atLowerBound)
2095          solution[iColumn] = lower_[i];
2096        else
2097          solution[iColumn] = upper_[i];
2098      }
2099    }
2100  } break;
2101  // Compute values of key variables
2102  case 1: {
2103    int i;
2104    double *solution = model->solutionRegion();
2105    //const int * columnLength = matrix_->getVectorLengths();
2106    //const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2107    //const int * row = matrix_->getIndices();
2108    //const double * elementByColumn = matrix_->getElements();
2109    //int * pivotVariable = model->pivotVariable();
2110    sumPrimalInfeasibilities_ = 0.0;
2111    numberPrimalInfeasibilities_ = 0;
2112    double primalTolerance = model->primalTolerance();
2113    double relaxedTolerance = primalTolerance;
2114    // we can't really trust infeasibilities if there is primal error
2115    double error = CoinMin(1.0e-2, model->largestPrimalError());
2116    // allow tolerance at least slightly bigger than standard
2117    relaxedTolerance = relaxedTolerance + error;
2118    // but we will be using difference
2119    relaxedTolerance -= primalTolerance;
2120    sumOfRelaxedPrimalInfeasibilities_ = 0.0;
2121    for (i = 0; i < numberSets_; i++) { // Could just be over basics (esp if no bounds)
2122      int kColumn = keyVariable_[i];
2123      double value = 0.0;
2124      if ((gubType_ & 8) != 0) {
2125        int iColumn = next_[kColumn];
2126        // sum all non-key variables
2127        while (iColumn >= 0) {
2128          value += solution[iColumn];
2129          iColumn = next_[iColumn];
2130        }
2131      } else {
2132        // bounds exist - sum over all except key
2133        int stop = -(kColumn + 1);
2134        int iColumn = next_[kColumn];
2135        // sum all non-key variables
2136        while (iColumn != stop) {
2137          if (iColumn < 0)
2138            iColumn = -iColumn - 1;
2139          value += solution[iColumn];
2140          iColumn = next_[iColumn];
2141        }
2142      }
2143      if (kColumn < numberColumns) {
2144        // make sure key is basic - so will be skipped in values pass
2145        model->setStatus(kColumn, ClpSimplex::basic);
2146        // feasibility will be done later
2147        assert(getStatus(i) != ClpSimplex::basic);
2148        if (getStatus(i) == ClpSimplex::atUpperBound)
2149          solution[kColumn] = upper_[i] - value;
2150        else
2151          solution[kColumn] = lower_[i] - value;
2152        //printf("Value of key structural %d for set %d is %g\n",kColumn,i,solution[kColumn]);
2153      } else {
2154        // slack is key
2155        assert(getStatus(i) == ClpSimplex::basic);
2156        double infeasibility = 0.0;
2157        if (value > upper_[i] + primalTolerance) {
2158          infeasibility = value - upper_[i] - primalTolerance;
2159          setAbove(i);
2160        } else if (value < lower_[i] - primalTolerance) {
2161          infeasibility = lower_[i] - value - primalTolerance;
2162          setBelow(i);
2163        } else {
2164          setFeasible(i);
2165        }
2166        //printf("Value of key slack for set %d is %g\n",i,value);
2167        if (infeasibility > 0.0) {
2168          sumPrimalInfeasibilities_ += infeasibility;
2169          if (infeasibility > relaxedTolerance)
2170            sumOfRelaxedPrimalInfeasibilities_ += infeasibility;
2171          numberPrimalInfeasibilities_++;
2172        }
2173      }
2174    }
2175  } break;
2176  // Report on infeasibilities of key variables
2177  case 2: {
2178    model->setSumPrimalInfeasibilities(model->sumPrimalInfeasibilities() + sumPrimalInfeasibilities_);
2179    model->setNumberPrimalInfeasibilities(model->numberPrimalInfeasibilities() + numberPrimalInfeasibilities_);
2180    model->setSumOfRelaxedPrimalInfeasibilities(model->sumOfRelaxedPrimalInfeasibilities() + sumOfRelaxedPrimalInfeasibilities_);
2181  } break;
2182  }
2183}
2184/*
2185     utility dual function for dealing with dynamic constraints
2186     mode=n see ClpGubMatrix.hpp for definition
2187     Remember to update here when settled down
2188*/
2189void ClpGubMatrix::dualExpanded(ClpSimplex *model,
2190  CoinIndexedVector *array,
2191  double * /*other*/, int mode)
2192{
2193  switch (mode) {
2194    // modify costs before transposeUpdate
2195  case 0: {
2196    int i;
2197    double *cost = model->costRegion();
2198    // not dual values yet
2199    //assert (!other);
2200    //double * work = array->denseVector();
2201    double infeasibilityCost = model->infeasibilityCost();
2202    int *pivotVariable = model->pivotVariable();
2203    int numberRows = model->numberRows();
2204    int numberColumns = model->numberColumns();
2205    for (i = 0; i < numberRows; i++) {
2206      int iPivot = pivotVariable[i];
2207      if (iPivot < numberColumns) {
2208        int iSet = backward_[iPivot];
2209        if (iSet >= 0) {
2210          int kColumn = keyVariable_[iSet];
2211          double costValue;
2212          if (kColumn < numberColumns) {
2213            // structural has cost
2214            costValue = cost[kColumn];
2215          } else {
2216            // slack is key
2217            assert(getStatus(iSet) == ClpSimplex::basic);
2218            // negative as -1.0 for slack
2219            costValue = -weight(iSet) * infeasibilityCost;
2220          }
2221          array->add(i, -costValue); // was work[i]-costValue
2222        }
2223      }
2224    }
2225  } break;
2226  // create duals for key variables (without check on dual infeasible)
2227  case 1: {
2228    // If key slack then dual 0.0 (if feasible)
2229    // dj for key is zero so that defines dual on set
2230    int i;
2231    double *dj = model->djRegion();
2232    int numberColumns = model->numberColumns();
2233    double infeasibilityCost = model->infeasibilityCost();
2234    for (i = 0; i < numberSets_; i++) {
2235      int kColumn = keyVariable_[i];
2236      if (kColumn < numberColumns) {
2237        // dj without set
2238        double value = dj[kColumn];
2239        // Now subtract out from all
2240        dj[kColumn] = 0.0;
2241        int iColumn = next_[kColumn];
2242        // modify all non-key variables
2243        while (iColumn >= 0) {
2244          dj[iColumn] -= value;
2245          iColumn = next_[iColumn];
2246        }
2247      } else {
2248        // slack key - may not be feasible
2249        assert(getStatus(i) == ClpSimplex::basic);
2250        // negative as -1.0 for slack
2251        double value = -weight(i) * infeasibilityCost;
2252        if (value) {
2253          int iColumn = next_[kColumn];
2254          // modify all non-key variables basic
2255          while (iColumn >= 0) {
2256            dj[iColumn] -= value;
2257            iColumn = next_[iColumn];
2258          }
2259        }
2260      }
2261    }
2262  } break;
2263  // as 1 but check slacks and compute djs
2264  case 2: {
2265    // If key slack then dual 0.0
2266    // If not then slack could be dual infeasible
2267    // dj for key is zero so that defines dual on set
2268    int i;
2269    // make sure fromIndex will not confuse pricing
2270    fromIndex_[0] = -1;
2271    possiblePivotKey_ = -1;
2272    // Create array
2273    int numberColumns = model->numberColumns();
2274    int *pivotVariable = model->pivotVariable();
2275    int numberRows = model->numberRows();
2276    for (i = 0; i < numberRows; i++) {
2277      int iPivot = pivotVariable[i];
2278      if (iPivot < numberColumns)
2279        backToPivotRow_[iPivot] = i;
2280    }
2281    if (noCheck_ >= 0) {
2282      if (infeasibilityWeight_ != model->infeasibilityCost()) {
2283        // don't bother checking
2284        sumDualInfeasibilities_ = 100.0;
2285        numberDualInfeasibilities_ = 1;
2286        sumOfRelaxedDualInfeasibilities_ = 100.0;
2287        return;
2288      }
2289    }
2290    double *dj = model->djRegion();
2291    double *dual = model->dualRowSolution();
2292    double *cost = model->costRegion();
2293    ClpSimplex::Status iStatus;
2294    const int *columnLength = matrix_->getVectorLengths();
2295    const CoinBigIndex *columnStart = matrix_->getVectorStarts();
2296    const int *row = matrix_->getIndices();
2297    const double *elementByColumn = matrix_->getElements();
2298    double infeasibilityCost = model->infeasibilityCost();
2299    sumDualInfeasibilities_ = 0.0;
2300    numberDualInfeasibilities_ = 0;
2301    double dualTolerance = model->dualTolerance();
2302    double relaxedTolerance = dualTolerance;
2303    // we can't really trust infeasibilities if there is dual error
2304    double error = CoinMin(1.0e-2, model->largestDualError());
2305    // allow tolerance at least slightly bigger than standard
2306    relaxedTolerance = relaxedTolerance + error;
2307    // but we will be using difference
2308    relaxedTolerance -= dualTolerance;
2309    sumOfRelaxedDualInfeasibilities_ = 0.0;
2310    for (i = 0; i < numberSets_; i++) {
2311      int kColumn = keyVariable_[i];
2312      if (kColumn < numberColumns) {
2313        // dj without set
2314        double value = cost[kColumn];
2315        for (CoinBigIndex j = columnStart[kColumn];
2316             j < columnStart[kColumn] + columnLength[kColumn]; j++) {
2317          int iRow = row[j];
2318          value -= dual[iRow] * elementByColumn[j];
2319        }
2320        // Now subtract out from all
2321        dj[kColumn] -= value;
2322        int stop = -(kColumn + 1);
2323        kColumn = next_[kColumn];
2324        while (kColumn != stop) {
2325          if (kColumn < 0)
2326            kColumn = -kColumn - 1;
2327          double djValue = dj[kColumn] - value;
2328          dj[kColumn] = djValue;
2329          ;
2330          double infeasibility = 0.0;
2331          iStatus = model->getStatus(kColumn);
2332          if (iStatus == ClpSimplex::atLowerBound) {
2333            if (djValue < -dualTolerance)
2334              infeasibility = -djValue - dualTolerance;
2335          } else if (iStatus == ClpSimplex::atUpperBound) {
2336            // at upper bound
2337            if (djValue > dualTolerance)
2338              infeasibility = djValue - dualTolerance;
2339          }
2340          if (infeasibility > 0.0) {
2341            sumDualInfeasibilities_ += infeasibility;
2342            if (infeasibility > relaxedTolerance)
2343              sumOfRelaxedDualInfeasibilities_ += infeasibility;
2344            numberDualInfeasibilities_++;
2345          }
2346          kColumn = next_[kColumn];
2347        }
2348        // check slack
2349        iStatus = getStatus(i);
2350        assert(iStatus != ClpSimplex::basic);
2351        double infeasibility = 0.0;
2352        // dj of slack is -(-1.0)value
2353        if (iStatus == ClpSimplex::atLowerBound) {
2354          if (value < -dualTolerance)
2355            infeasibility = -value - dualTolerance;
2356        } else if (iStatus == ClpSimplex::atUpperBound) {
2357          // at upper bound
2358          if (value > dualTolerance)
2359            infeasibility = value - dualTolerance;
2360        }
2361        if (infeasibility > 0.0) {
2362          sumDualInfeasibilities_ += infeasibility;
2363          if (infeasibility > relaxedTolerance)
2364            sumOfRelaxedDualInfeasibilities_ += infeasibility;
2365          numberDualInfeasibilities_++;
2366        }
2367      } else {
2368        // slack key - may not be feasible
2369        assert(getStatus(i) == ClpSimplex::basic);
2370        // negative as -1.0 for slack
2371        double value = -weight(i) * infeasibilityCost;
2372        if (value) {
2373          // Now subtract out from all
2374          int kColumn = i + numberColumns;
2375          int stop = -(kColumn + 1);
2376          kColumn = next_[kColumn];
2377          while (kColumn != stop) {
2378            if (kColumn < 0)
2379              kColumn = -kColumn - 1;
2380            double djValue = dj[kColumn] - value;
2381            dj[kColumn] = djValue;
2382            ;
2383            double infeasibility = 0.0;
2384            iStatus = model->getStatus(kColumn);
2385            if (iStatus == ClpSimplex::atLowerBound) {
2386              if (djValue < -dualTolerance)
2387                infeasibility = -djValue - dualTolerance;
2388            } else if (iStatus == ClpSimplex::atUpperBound) {
2389              // at upper bound
2390              if (djValue > dualTolerance)
2391                infeasibility = djValue - dualTolerance;
2392            }
2393            if (infeasibility > 0.0) {
2394              sumDualInfeasibilities_ += infeasibility;
2395              if (infeasibility > relaxedTolerance)
2396                sumOfRelaxedDualInfeasibilities_ += infeasibility;
2397              numberDualInfeasibilities_++;
2398            }
2399            kColumn = next_[kColumn];
2400          }
2401        }
2402      }
2403    }
2404    // and get statistics for column generation
2405    synchronize(model, 4);
2406    infeasibilityWeight_ = -1.0;
2407  } break;
2408  // Report on infeasibilities of key variables
2409  case 3: {
2410    model->setSumDualInfeasibilities(model->sumDualInfeasibilities() + sumDualInfeasibilities_);
2411    model->setNumberDualInfeasibilities(model->numberDualInfeasibilities() + numberDualInfeasibilities_);
2412    model->setSumOfRelaxedDualInfeasibilities(model->sumOfRelaxedDualInfeasibilities() + sumOfRelaxedDualInfeasibilities_);
2413  } break;
2414  // modify costs before transposeUpdate for partial pricing
2415  case 4: {
2416    // First compute new costs etc for interesting gubs
2417    int iLook = 0;
2418    int iSet = fromIndex_[0];
2419    double primalTolerance = model->primalTolerance();
2420    const double *cost = model->costRegion();
2421    double *solution = model->solutionRegion();
2422    double infeasibilityCost = model->infeasibilityCost();
2423    int numberColumns = model->numberColumns();
2424    int numberChanged = 0;
2425    int *pivotVariable = model->pivotVariable();
2426    while (iSet >= 0) {
2427      int key = keyVariable_[iSet];
2428      double value = 0.0;
2429      // sum over all except key
2430      if ((gubType_ & 8) != 0) {
2431        int iColumn = next_[key];
2432        // sum all non-key variables
2433        while (iColumn >= 0) {
2434          value += solution[iColumn];
2435          iColumn = next_[iColumn];
2436        }
2437      } else {
2438        // bounds exist - sum over all except key
2439        int stop = -(key + 1);
2440        int iColumn = next_[key];
2441        // sum all non-key variables
2442        while (iColumn != stop) {
2443          if (iColumn < 0)
2444            iColumn = -iColumn - 1;
2445          value += solution[iColumn];
2446          iColumn = next_[iColumn];
2447        }
2448      }
2449      double costChange;
2450      double oldCost = changeCost_[iLook];
2451      if (key < numberColumns) {
2452        assert(getStatus(iSet) != ClpSimplex::basic);
2453        double sol;
2454        if (getStatus(iSet) == ClpSimplex::atUpperBound)
2455          sol = upper_[iSet] - value;
2456        else
2457          sol = lower_[iSet] - value;
2458        solution[key] = sol;
2459        // fix up cost
2460        model->nonLinearCost()->setOne(key, sol);
2461#ifdef CLP_DEBUG_PRINT
2462        printf("yy Value of key structural %d for set %d is %g - cost %g old cost %g\n", key, iSet, sol,
2463          cost[key], oldCost);
2464#endif
2465        costChange = cost[key] - oldCost;
2466      } else {
2467        // slack is key
2468        if (value > upper_[iSet] + primalTolerance) {
2469          setAbove(iSet);
2470        } else if (value < lower_[iSet] - primalTolerance) {
2471          setBelow(iSet);
2472        } else {
2473          setFeasible(iSet);
2474        }
2475        // negative as -1.0 for slack
2476        costChange = -weight(iSet) * infeasibilityCost - oldCost;
2477#ifdef CLP_DEBUG_PRINT
2478        printf("yy Value of key slack for set %d is %g - cost %g old cost %g\n", iSet, value,
2479          weight(iSet) * infeasibilityCost, oldCost);
2480#endif
2481      }
2482      if (costChange) {
2483        fromIndex_[numberChanged] = iSet;
2484        toIndex_[iSet] = numberChanged;
2485        changeCost_[numberChanged++] = costChange;
2486      }
2487      iSet = fromIndex_[++iLook];
2488    }
2489    if (numberChanged || possiblePivotKey_ >= 0) {
2490      // first do those in list already
2491      int number = array->getNumElements();
2492      array->setPacked();
2493      int i;
2494      double *work = array->denseVector();
2495      int *which = array->getIndices();
2496      for (i = 0; i < number; i++) {
2497        int iRow = which[i];
2498        int iPivot = pivotVariable[iRow];
2499        if (iPivot < numberColumns) {
2500          int iSet = backward_[iPivot];
2501          if (iSet >= 0 && toIndex_[iSet] >= 0) {
2502            double newValue = work[i] + changeCost_[toIndex_[iSet]];
2503            if (!newValue)
2504              newValue = 1.0e-100;
2505            work[i] = newValue;
2506            // mark as done
2507            backward_[iPivot] = -1;
2508          }
2509        }
2510        if (possiblePivotKey_ == iRow) {
2511          double newValue = work[i] - model->dualIn();
2512          if (!newValue)
2513            newValue = 1.0e-100;
2514          work[i] = newValue;
2515          possiblePivotKey_ = -1;
2516        }
2517      }
2518      // now do rest and clean up
2519      for (i = 0; i < numberChanged; i++) {
2520        int iSet = fromIndex_[i];
2521        int key = keyVariable_[iSet];
2522        int iColumn = next_[key];
2523        double change = changeCost_[i];
2524        while (iColumn >= 0) {
2525          if (backward_[iColumn] >= 0) {
2526            int iRow = backToPivotRow_[iColumn];
2527            assert(iRow >= 0);
2528            work[number] = change;
2529            if (possiblePivotKey_ == iRow) {
2530              double newValue = work[number] - model->dualIn();
2531              if (!newValue)
2532                newValue = 1.0e-100;
2533              work[number] = newValue;
2534              possiblePivotKey_ = -1;
2535            }
2536            which[number++] = iRow;
2537          } else {
2538            // reset
2539            backward_[iColumn] = iSet;
2540          }
2541          iColumn = next_[iColumn];
2542        }
2543        toIndex_[iSet] = -1;
2544      }
2545      if (possiblePivotKey_ >= 0) {
2546        work[number] = -model->dualIn();
2547        which[number++] = possiblePivotKey_;
2548        possiblePivotKey_ = -1;
2549      }
2550      fromIndex_[0] = -1;
2551      array->setNumElements(number);
2552    }
2553  } break;
2554  }
2555}
2556// This is local to Gub to allow synchronization when status is good
2557int ClpGubMatrix::synchronize(ClpSimplex *, int)
2558{
2559  return 0;
2560}
2561/*
2562     general utility function for dealing with dynamic constraints
2563     mode=n see ClpGubMatrix.hpp for definition
2564     Remember to update here when settled down
2565*/
2566int ClpGubMatrix::generalExpanded(ClpSimplex *model, int mode, int &number)
2567{
2568  int returnCode = 0;
2569  int numberColumns = model->numberColumns();
2570  switch (mode) {
2571    // Fill in pivotVariable but not for key variables
2572  case 0: {
2573    if (!next_) {
2574      // do ordering
2575      assert(!rhsOffset_);
2576      // create and do gub crash
2577      useEffectiveRhs(model, false);
2578    }
2579    int i;
2580    int numberBasic = number;
2581    // Use different array so can build from true pivotVariable_
2582    //int * pivotVariable = model->pivotVariable();
2583    int *pivotVariable = model->rowArray(0)->getIndices();
2584    for (i = 0; i < numberColumns; i++) {
2585      if (model->getColumnStatus(i) == ClpSimplex::basic) {
2586        int iSet = backward_[i];
2587        if (iSet < 0 || i != keyVariable_[iSet])
2588          pivotVariable[numberBasic++] = i;
2589      }
2590    }
2591    number = numberBasic;
2592    if (model->numberIterations())
2593      assert(number == model->numberRows());
2594  } break;
2595  // Make all key variables basic
2596  case 1: {
2597    int i;
2598    for (i = 0; i < numberSets_; i++) {
2599      int iColumn = keyVariable_[i];
2600      if (iColumn < numberColumns)
2601        model->setColumnStatus(iColumn, ClpSimplex::basic);
2602    }
2603  } break;
2604  // Do initial extra rows + maximum basic
2605  case 2: {
2606    returnCode = getNumRows() + 1;
2607    number = model->numberRows() + numberSets_;
2608  } break;
2609  // Before normal replaceColumn
2610  case 3: {
2611    int sequenceIn = model->sequenceIn();
2612    int sequenceOut = model->sequenceOut();
2613    int numberColumns = model->numberColumns();
2614    int numberRows = model->numberRows();
2615    int pivotRow = model->pivotRow();
2616    if (gubSlackIn_ >= 0)
2617      assert(sequenceIn > numberRows + numberColumns);
2618    if (sequenceIn == sequenceOut)
2619      return -1;
2620    int iSetIn = -1;
2621    int iSetOut = -1;
2622    if (sequenceOut < numberColumns) {
2623      iSetOut = backward_[sequenceOut];
2624    } else if (sequenceOut >= numberRows + numberColumns) {
2625      assert(pivotRow >= numberRows);
2626      int iExtra = pivotRow - numberRows;
2627      assert(iExtra >= 0);
2628      if (iSetOut < 0)
2629        iSetOut = fromIndex_[iExtra];
2630      else
2631        assert(iSetOut == fromIndex_[iExtra]);
2632    }
2633    if (sequenceIn < numberColumns) {
2634      iSetIn = backward_[sequenceIn];
2635    } else if (gubSlackIn_ >= 0) {
2636      iSetIn = gubSlackIn_;
2637    }
2638    possiblePivotKey_ = -1;
2639    number = 0; // say do ordinary
2640    int *pivotVariable = model->pivotVariable();
2641    if (pivotRow >= numberRows) {
2642      int iExtra = pivotRow - numberRows;
2643      //const int * length = matrix_->getVectorLengths();
2644
2645      assert(sequenceOut >= numberRows + numberColumns || sequenceOut == keyVariable_[iSetOut]);
2646      int incomingColumn = sequenceIn; // to be used in updates
2647      if (iSetIn != iSetOut) {
2648        // We need to find a possible pivot for incoming
2649        // look through rowArray_[1]
2650        int n = model->rowArray(1)->getNumElements();
2651        int *which = model->rowArray(1)->getIndices();
2652        double *array = model->rowArray(1)->denseVector();
2653        double bestAlpha = 1.0e-5;
2654        //int shortest=numberRows+1;
2655        for (int i = 0; i < n; i++) {
2656          int iRow = which[i];
2657          int iPivot = pivotVariable[iRow];
2658          if (iPivot < numberColumns && backward_[iPivot] == iSetOut) {
2659            if (fabs(array[i]) > fabs(bestAlpha)) {
2660              bestAlpha = array[i];
2661              possiblePivotKey_ = iRow;
2662            }
2663          }
2664        }
2665        assert(possiblePivotKey_ >= 0); // could set returnCode=4
2666        number = 1;
2667        if (sequenceIn >= numberRows + numberColumns) {
2668          number = 3;
2669          // need swap as gub slack in and must become key
2670          // is this best way
2671          int key = keyVariable_[iSetIn];
2672          assert(key < numberColumns);
2673          // check other basic
2674          int iColumn = next_[key];
2675          // set new key to be used by unpack
2676          keyVariable_[iSetIn] = iSetIn + numberColumns;
2677          // change cost in changeCost
2678          {
2679            int iLook = 0;
2680            int iSet = fromIndex_[0];
2681            while (iSet >= 0) {
2682              if (iSet == iSetIn) {
2683                changeCost_[iLook] = 0.0;
2684                break;
2685              }
2686              iSet = fromIndex_[++iLook];
2687            }
2688          }
2689          while (iColumn >= 0) {
2690            if (iColumn != sequenceOut) {
2691              // need partial ftran and skip accuracy check in replaceColumn
2692#ifdef CLP_DEBUG_PRINT
2693              printf("TTTTTry 5\n");
2694#endif
2695              int iRow = backToPivotRow_[iColumn];
2696              assert(iRow >= 0);
2697              unpack(model, model->rowArray(3), iColumn);
2698              model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
2699              double alpha = model->rowArray(3)->denseVector()[iRow];
2700              //if (!alpha)
2701              //printf("zero alpha a\n");
2702              int updateStatus = model->factorization()->replaceColumn(model,
2703                model->rowArray(2),
2704                model->rowArray(3),
2705                iRow, alpha);
2706              returnCode = CoinMax(updateStatus, returnCode);
2707              model->rowArray(3)->clear();
2708              if (returnCode)
2709                break;
2710            }
2711            iColumn = next_[iColumn];
2712          }
2713          if (!returnCode) {
2714            // now factorization looks as if key is out
2715            // pivot back in
2716#ifdef CLP_DEBUG_PRINT
2717            printf("TTTTTry 6\n");
2718#endif
2719            unpack(model, model->rowArray(3), key);
2720            model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
2721            pivotRow = possiblePivotKey_;
2722            double alpha = model->rowArray(3)->denseVector()[pivotRow];
2723            //if (!alpha)
2724            //printf("zero alpha b\n");
2725            int updateStatus = model->factorization()->replaceColumn(model,
2726              model->rowArray(2),
2727              model->rowArray(3),
2728              pivotRow, alpha);
2729            returnCode = CoinMax(updateStatus, returnCode);
2730            model->rowArray(3)->clear();
2731          }
2732          // restore key
2733          keyVariable_[iSetIn] = key;
2734          // now alternate column can replace key on out
2735          incomingColumn = pivotVariable[possiblePivotKey_];
2736        } else {
2737#ifdef CLP_DEBUG_PRINT
2738          printf("TTTTTTry 4 %d\n", possiblePivotKey_);
2739#endif
2740          int updateStatus = model->factorization()->replaceColumn(model,
2741            model->rowArray(2),
2742            model->rowArray(1),
2743            possiblePivotKey_,
2744            bestAlpha);
2745          returnCode = CoinMax(updateStatus, returnCode);
2746          incomingColumn = pivotVariable[possiblePivotKey_];
2747        }
2748
2749        //returnCode=4; // need swap
2750      } else {
2751        // key swap
2752        number = -1;
2753      }
2754      int key = keyVariable_[iSetOut];
2755      if (key < numberColumns)
2756        assert(key == sequenceOut);
2757      // check if any other basic
2758      int iColumn = next_[key];
2759      if (returnCode)
2760        iColumn = -1; // skip if error on previous
2761      // set new key to be used by unpack
2762      if (incomingColumn < numberColumns)
2763        keyVariable_[iSetOut] = incomingColumn;
2764      else
2765        keyVariable_[iSetOut] = iSetIn + numberColumns;
2766      double *cost = model->costRegion();
2767      if (possiblePivotKey_ < 0) {
2768        double dj = model->djRegion()[sequenceIn] - cost[sequenceIn];
2769        changeCost_[iExtra] = -dj;
2770#ifdef CLP_DEBUG_PRINT
2771        printf("modifying changeCost %d by %g - cost %g\n", iExtra, dj, cost[sequenceIn]);
2772#endif
2773      }
2774      while (iColumn >= 0) {
2775        if (iColumn != incomingColumn) {
2776          number = -2;
2777          // need partial ftran and skip accuracy check in replaceColumn
2778#ifdef CLP_DEBUG_PRINT
2779          printf("TTTTTTry 1\n");
2780#endif
2781          int iRow = backToPivotRow_[iColumn];
2782          assert(iRow >= 0 && iRow < numberRows);
2783          unpack(model, model->rowArray(3), iColumn);
2784          model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
2785          double *array = model->rowArray(3)->denseVector();
2786          double alpha = array[iRow];
2787          //if (!alpha)
2788          //printf("zero alpha d\n");
2789          int updateStatus = model->factorization()->replaceColumn(model,
2790            model->rowArray(2),
2791            model->rowArray(3),
2792            iRow, alpha);
2793          returnCode = CoinMax(updateStatus, returnCode);
2794          model->rowArray(3)->clear();
2795          if (returnCode)
2796            break;
2797        }
2798        iColumn = next_[iColumn];
2799      }
2800      // restore key
2801      keyVariable_[iSetOut] = key;
2802    } else if (sequenceIn >= numberRows + numberColumns) {
2803      number = 2;
2804      //returnCode=4;
2805      // need swap as gub slack in and must become key
2806      // is this best way
2807      int key = keyVariable_[iSetIn];
2808      assert(key < numberColumns);
2809      // check other basic
2810      int iColumn = next_[key];
2811      // set new key to be used by unpack
2812      keyVariable_[iSetIn] = iSetIn + numberColumns;
2813      // change cost in changeCost
2814      {
2815        int iLook = 0;
2816        int iSet = fromIndex_[0];
2817        while (iSet >= 0) {
2818          if (iSet == iSetIn) {
2819            changeCost_[iLook] = 0.0;
2820            break;
2821          }
2822          iSet = fromIndex_[++iLook];
2823        }
2824      }
2825      while (iColumn >= 0) {
2826        if (iColumn != sequenceOut) {
2827          // need partial ftran and skip accuracy check in replaceColumn
2828#ifdef CLP_DEBUG_PRINT
2829          printf("TTTTTry 2\n");
2830#endif
2831          int iRow = backToPivotRow_[iColumn];
2832          assert(iRow >= 0);
2833          unpack(model, model->rowArray(3), iColumn);
2834          model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
2835          double alpha = model->rowArray(3)->denseVector()[iRow];
2836          //if (!alpha)
2837          //printf("zero alpha e\n");
2838          int updateStatus = model->factorization()->replaceColumn(model,
2839            model->rowArray(2),
2840            model->rowArray(3),
2841            iRow, alpha);
2842          returnCode = CoinMax(updateStatus, returnCode);
2843          model->rowArray(3)->clear();
2844          if (returnCode)
2845            break;
2846        }
2847        iColumn = next_[iColumn];
2848      }
2849      if (!returnCode) {
2850        // now factorization looks as if key is out
2851        // pivot back in
2852#ifdef CLP_DEBUG_PRINT
2853        printf("TTTTTry 3\n");
2854#endif
2855        unpack(model, model->rowArray(3), key);
2856        model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
2857        double alpha = model->rowArray(3)->denseVector()[pivotRow];
2858        //if (!alpha)
2859        //printf("zero alpha f\n");
2860        int updateStatus = model->factorization()->replaceColumn(model,
2861          model->rowArray(2),
2862          model->rowArray(3),
2863          pivotRow, alpha);
2864        returnCode = CoinMax(updateStatus, returnCode);
2865        model->rowArray(3)->clear();
2866      }
2867      // restore key
2868      keyVariable_[iSetIn] = key;
2869    } else {
2870      // normal - but might as well do here
2871      returnCode = model->factorization()->replaceColumn(model,
2872        model->rowArray(2),
2873        model->rowArray(1),
2874        model->pivotRow(),
2875        model->alpha());
2876    }
2877  }
2878#ifdef CLP_DEBUG_PRINT
2879    printf("Update type after %d - status %d - pivot row %d\n",
2880      number, returnCode, model->pivotRow());
2881#endif
2882    // see if column generation says time to re-factorize
2883    returnCode = CoinMax(returnCode, synchronize(model, 5));
2884    number = -1; // say no need for normal replaceColumn
2885    break;
2886  // To see if can dual or primal
2887  case 4: {
2888    returnCode = 1;
2889  } break;
2890  // save status
2891  case 5: {
2892    synchronize(model, 0);
2893    CoinMemcpyN(status_, numberSets_, saveStatus_);
2894    CoinMemcpyN(keyVariable_, numberSets_, savedKeyVariable_);
2895  } break;
2896  // restore status
2897  case 6: {
2898    CoinMemcpyN(saveStatus_, numberSets_, status_);
2899    CoinMemcpyN(savedKeyVariable_, numberSets_, keyVariable_);
2900    // restore firstAvailable_
2901    synchronize(model, 7);
2902    // redo next_
2903    int i;
2904    int *last = new int[numberSets_];
2905    for (i = 0; i < numberSets_; i++) {
2906      int iKey = keyVariable_[i];
2907      assert(iKey >= numberColumns || backward_[iKey] == i);
2908      last[i] = iKey;
2909      // make sure basic
2910      //if (iKey<numberColumns)
2911      //model->setStatus(iKey,ClpSimplex::basic);
2912    }
2913    for (i = 0; i < numberColumns; i++) {
2914      int iSet = backward_[i];
2915      if (iSet >= 0) {
2916        next_[last[iSet]] = i;
2917        last[iSet] = i;
2918      }
2919    }
2920    for (i = 0; i < numberSets_; i++) {
2921      next_[last[i]] = -(keyVariable_[i] + 1);
2922      redoSet(model, keyVariable_[i], keyVariable_[i], i);
2923    }
2924    delete[] last;
2925    // redo pivotVariable
2926    int *pivotVariable = model->pivotVariable();
2927    int iRow;
2928    int numberBasic = 0;
2929    int numberRows = model->numberRows();
2930    for (iRow = 0; iRow < numberRows; iRow++) {
2931      if (model->getRowStatus(iRow) == ClpSimplex::basic) {
2932        numberBasic++;
2933        pivotVariable[iRow] = iRow + numberColumns;
2934      } else {
2935        pivotVariable[iRow] = -1;
2936      }
2937    }
2938    i = 0;
2939    int iColumn;
2940    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2941      if (model->getStatus(iColumn) == ClpSimplex::basic) {
2942        int iSet = backward_[iColumn];
2943        if (iSet < 0 || keyVariable_[iSet] != iColumn) {
2944          while (pivotVariable[i] >= 0) {
2945            i++;
2946            assert(i < numberRows);
2947          }
2948          pivotVariable[i] = iColumn;
2949          backToPivotRow_[iColumn] = i;
2950          numberBasic++;
2951        }
2952      }
2953    }
2954    assert(numberBasic == numberRows);
2955    rhsOffset(model, true);
2956  } break;
2957  // flag a variable
2958  case 7: {
2959    assert(number == model->sequenceIn());
2960    synchronize(model, 1);
2961    synchronize(model, 8);
2962  } break;
2963  // unflag all variables
2964  case 8: {
2965    returnCode = synchronize(model, 2);
2966  } break;
2967  // redo costs in primal
2968  case 9: {
2969    returnCode = synchronize(model, 3);
2970  } break;
2971  // return 1 if there may be changing bounds on variable (column generation)
2972  case 10: {
2973    returnCode = synchronize(model, 6);
2974  } break;
2975  // make sure set is clean
2976  case 11: {
2977    assert(number == model->sequenceIn());
2978    returnCode = synchronize(model, 8);
2979  } break;
2980  default:
2981    break;
2982  }
2983  return returnCode;
2984}
2985// Sets up an effective RHS and does gub crash if needed
2986void ClpGubMatrix::useEffectiveRhs(ClpSimplex *model, bool cheapest)
2987{
2988  // Do basis - cheapest or slack if feasible (unless cheapest set)
2989  int longestSet = 0;
2990  int iSet;
2991  for (iSet = 0; iSet < numberSets_; iSet++)
2992    longestSet = CoinMax(longestSet, end_[iSet] - start_[iSet]);
2993
2994  double *upper = new double[longestSet + 1];
2995  double *cost = new double[longestSet + 1];
2996  double *lower = new double[longestSet + 1];
2997  double *solution = new double[longestSet + 1];
2998  assert(!next_);
2999  int numberColumns = getNumCols();
3000  const int *columnLength = matrix_->getVectorLengths();
3001  const double *columnLower = model->lowerRegion();
3002  const double *columnUpper = model->upperRegion();
3003  double *columnSolution = model->solutionRegion();
3004  const double *objective = model->costRegion();
3005  int numberRows = getNumRows();
3006  toIndex_ = new int[numberSets_];
3007  for (iSet = 0; iSet < numberSets_; iSet++)
3008    toIndex_[iSet] = -1;
3009  fromIndex_ = new int[getNumRows() + 1];
3010  double tolerance = model->primalTolerance();
3011  bool noNormalBounds = true;
3012  gubType_ &= ~8;
3013  bool gotBasis = false;
3014  for (iSet = 0; iSet < numberSets_; iSet++) {
3015    if (keyVariable_[iSet] < numberColumns)
3016      gotBasis = true;
3017    CoinBigIndex j;
3018    CoinBigIndex iStart = start_[iSet];
3019    CoinBigIndex iEnd = end_[iSet];
3020    for (j = iStart; j < iEnd; j++) {
3021      if (columnLower[j] && columnLower[j] > -1.0e20)
3022        noNormalBounds = false;
3023      if (columnUpper[j] && columnUpper[j] < 1.0e20)
3024        noNormalBounds = false;
3025    }
3026  }
3027  if (noNormalBounds)
3028    gubType_ |= 8;
3029  if (!gotBasis) {
3030    for (iSet = 0; iSet < numberSets_; iSet++) {
3031      CoinBigIndex j;
3032      int numberBasic = 0;
3033      int iBasic = -1;
3034      CoinBigIndex iStart = start_[iSet];
3035      CoinBigIndex iEnd = end_[iSet];
3036      // find one with smallest length
3037      int smallest = numberRows + 1;
3038      double value = 0.0;
3039      for (j = iStart; j < iEnd; j++) {
3040        if (model->getStatus(j) == ClpSimplex::basic) {
3041          if (columnLength[j] < smallest) {
3042            smallest = columnLength[j];
3043            iBasic = j;
3044          }
3045          numberBasic++;
3046        }
3047        value += columnSolution[j];
3048      }
3049      bool done = false;
3050      if (numberBasic > 1 || (numberBasic == 1 && getStatus(iSet) == ClpSimplex::basic)) {
3051        if (getStatus(iSet) == ClpSimplex::basic)
3052          iBasic = iSet + numberColumns; // slack key - use
3053        done = true;
3054      } else if (numberBasic == 1) {
3055        // see if can be key
3056        double thisSolution = columnSolution[iBasic];
3057        if (thisSolution > columnUpper[iBasic]) {
3058          value -= thisSolution - columnUpper[iBasic];
3059          thisSolution = columnUpper[iBasic];
3060          columnSolution[iBasic] = thisSolution;
3061        }
3062        if (thisSolution < columnLower[iBasic]) {
3063          value -= thisSolution - columnLower[iBasic];
3064          thisSolution = columnLower[iBasic];
3065          columnSolution[iBasic] = thisSolution;
3066        }
3067        // try setting slack to a bound
3068        assert(upper_[iSet] < 1.0e20 || lower_[iSet] > -1.0e20);
3069        double cost1 = COIN_DBL_MAX;
3070        int whichBound = -1;
3071        if (upper_[iSet] < 1.0e20) {
3072          // try slack at ub
3073          double newBasic = thisSolution + upper_[iSet] - value;
3074          if (newBasic >= columnLower[iBasic] - tolerance && newBasic <= columnUpper[iBasic] + tolerance) {
3075            // can go
3076            whichBound = 1;
3077            cost1 = newBasic * objective[iBasic];
3078            // But if exact then may be good solution
3079            if (fabs(upper_[iSet] - value) < tolerance)
3080              cost1 = -COIN_DBL_MAX;
3081          }
3082        }
3083        if (lower_[iSet] > -1.0e20) {
3084          // try slack at lb
3085          double newBasic = thisSolution + lower_[iSet] - value;
3086          if (newBasic >= columnLower[iBasic] - tolerance && newBasic <= columnUpper[iBasic] + tolerance) {
3087            // can go but is it cheaper
3088            double cost2 = newBasic * objective[iBasic];
3089            // But if exact then may be good solution
3090            if (fabs(lower_[iSet] - value) < tolerance)
3091              cost2 = -COIN_DBL_MAX;
3092            if (cost2 < cost1)
3093              whichBound = 0;
3094          }
3095        }
3096        if (whichBound != -1) {
3097          // key
3098          done = true;
3099          if (whichBound) {
3100            // slack to upper
3101            columnSolution[iBasic] = thisSolution + upper_[iSet] - value;
3102            setStatus(iSet, ClpSimplex::atUpperBound);
3103          } else {
3104            // slack to lower
3105            columnSolution[iBasic] = thisSolution + lower_[iSet] - value;
3106            setStatus(iSet, ClpSimplex::atLowerBound);
3107          }
3108        }
3109      }
3110      if (!done) {
3111        if (!cheapest) {
3112          // see if slack can be key
3113          if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
3114            done = true;
3115            setStatus(iSet, ClpSimplex::basic);
3116            iBasic = iSet + numberColumns;
3117          }
3118        }
3119        if (!done) {
3120          // set non basic if there was one
3121          if (iBasic >= 0)
3122            model->setStatus(iBasic, ClpSimplex::atLowerBound);
3123          // find cheapest
3124          int numberInSet = iEnd - iStart;
3125          CoinMemcpyN(columnLower + iStart, numberInSet, lower);
3126          CoinMemcpyN(columnUpper + iStart, numberInSet, upper);
3127          CoinMemcpyN(columnSolution + iStart, numberInSet, solution);
3128          // and slack
3129          iBasic = numberInSet;
3130          solution[iBasic] = -value;
3131          lower[iBasic] = -upper_[iSet];
3132          upper[iBasic] = -lower_[iSet];
3133          int kphase;
3134          if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
3135            // feasible
3136            kphase = 1;
3137            cost[iBasic] = 0.0;
3138            CoinMemcpyN(objective + iStart, numberInSet, cost);
3139          } else {
3140            // infeasible
3141            kphase = 0;
3142            // remember bounds are flipped so opposite to natural
3143            if (value < lower_[iSet] - tolerance)
3144              cost[iBasic] = 1.0;
3145            else
3146              cost[iBasic] = -1.0;
3147            CoinZeroN(cost, numberInSet);
3148          }
3149          double dualTolerance = model->dualTolerance();
3150          for (int iphase = kphase; iphase < 2; iphase++) {
3151            if (iphase) {
3152              cost[numberInSet] = 0.0;
3153              CoinMemcpyN(objective + iStart, numberInSet, cost);
3154            }
3155            // now do one row lp
3156            bool improve = true;
3157            while (improve) {
3158              improve = false;
3159              double dual = cost[iBasic];
3160              int chosen = -1;
3161              double best = dualTolerance;
3162              int way = 0;
3163              for (int i = 0; i <= numberInSet; i++) {
3164                double dj = cost[i] - dual;
3165                double improvement = 0.0;
3166                if (iphase || i < numberInSet)
3167                  assert(solution[i] >= lower[i] && solution[i] <= upper[i]);
3168                if (dj > dualTolerance)
3169                  improvement = dj * (solution[i] - lower[i]);
3170                else if (dj < -dualTolerance)
3171                  improvement = dj * (solution[i] - upper[i]);
3172                if (improvement > best) {
3173                  best = improvement;
3174                  chosen = i;
3175                  if (dj < 0.0) {
3176                    way = 1;
3177                  } else {
3178                    way = -1;
3179                  }
3180                }
3181              }
3182              if (chosen >= 0) {
3183                improve = true;
3184                // now see how far
3185                if (way > 0) {
3186                  // incoming increasing so basic decreasing
3187                  // if phase 0 then go to nearest bound
3188                  double distance = upper[chosen] - solution[chosen];
3189                  double basicDistance;
3190                  if (!iphase) {
3191                    assert(iBasic == numberInSet);
3192                    assert(solution[iBasic] > upper[iBasic]);
3193                    basicDistance = solution[iBasic] - upper[iBasic];
3194                  } else {
3195                    basicDistance = solution[iBasic] - lower[iBasic];
3196                  }
3197                  // need extra coding for unbounded
3198                  assert(CoinMin(distance, basicDistance) < 1.0e20);
3199                  if (distance > basicDistance) {
3200                    // incoming becomes basic
3201                    solution[chosen] += basicDistance;
3202                    if (!iphase)
3203                      solution[iBasic] = upper[iBasic];
3204                    else
3205                      solution[iBasic] = lower[iBasic];
3206                    iBasic = chosen;
3207                  } else {
3208                    // flip
3209                    solution[chosen] = upper[chosen];
3210                    solution[iBasic] -= distance;
3211                  }
3212                } else {
3213                  // incoming decreasing so basic increasing
3214                  // if phase 0 then go to nearest bound
3215                  double distance = solution[chosen] - lower[chosen];
3216                  double basicDistance;
3217                  if (!iphase) {
3218                    assert(iBasic == numberInSet);
3219                    assert(solution[iBasic] < lower[iBasic]);
3220                    basicDistance = lower[iBasic] - solution[iBasic];
3221                  } else {
3222                    basicDistance = upper[iBasic] - solution[iBasic];
3223                  }
3224                  // need extra coding for unbounded - for now just exit
3225                  if (CoinMin(distance, basicDistance) > 1.0e20) {
3226                    printf("unbounded on set %d\n", iSet);
3227                    iphase = 1;
3228                    iBasic = numberInSet;
3229                    break;
3230                  }
3231                  if (distance > basicDistance) {
3232                    // incoming becomes basic
3233                    solution[chosen] -= basicDistance;
3234                    if (!iphase)
3235                      solution[iBasic] = lower[iBasic];
3236                    else
3237                      solution[iBasic] = upper[iBasic];
3238                    iBasic = chosen;
3239                  } else {
3240                    // flip
3241                    solution[chosen] = lower[chosen];
3242                    solution[iBasic] += distance;
3243                  }
3244                }
3245                if (!iphase) {
3246                  if (iBasic < numberInSet)
3247                    break; // feasible
3248                  else if (solution[iBasic] >= lower[iBasic] && solution[iBasic] <= upper[iBasic])
3249                    break; // feasible (on flip)
3250                }
3251              }
3252            }
3253          }
3254          // convert iBasic back and do bounds
3255          if (iBasic == numberInSet) {
3256            // slack basic
3257            setStatus(iSet, ClpSimplex::basic);
3258            iBasic = iSet + numberColumns;
3259          } else {
3260            iBasic += start_[iSet];
3261            model->setStatus(iBasic, ClpSimplex::basic);
3262            // remember bounds flipped
3263            if (upper[numberInSet] == lower[numberInSet])
3264              setStatus(iSet, ClpSimplex::isFixed);
3265            else if (solution[numberInSet] == upper[numberInSet])
3266              setStatus(iSet, ClpSimplex::atLowerBound);
3267            else if (solution[numberInSet] == lower[numberInSet])
3268              setStatus(iSet, ClpSimplex::atUpperBound);
3269            else
3270              abort();
3271          }
3272          for (j = iStart; j < iEnd; j++) {
3273            if (model->getStatus(j) != ClpSimplex::basic) {
3274              int inSet = j - iStart;
3275              columnSolution[j] = solution[inSet];
3276              if (upper[inSet] == lower[inSet])
3277                model->setStatus(j, ClpSimplex::isFixed);
3278              else if (solution[inSet] == upper[inSet])
3279                model->setStatus(j, ClpSimplex::atUpperBound);
3280              else if (solution[inSet] == lower[inSet])
3281                model->setStatus(j, ClpSimplex::atLowerBound);
3282            }
3283          }
3284        }
3285      }
3286      keyVariable_[iSet] = iBasic;
3287    }
3288  }
3289  delete[] lower;
3290  delete[] solution;
3291  delete[] upper;
3292  delete[] cost;
3293  // make sure matrix is in good shape
3294  matrix_->orderMatrix();
3295  // create effective rhs
3296  delete[] rhsOffset_;
3297  rhsOffset_ = new double[numberRows];
3298  delete[] next_;
3299  next_ = new int[numberColumns + numberSets_ + 2 * longestSet];
3300  char *mark = new char[numberColumns];
3301  memset(mark, 0, numberColumns);
3302  for (int iColumn = 0; iColumn < numberColumns; iColumn++)
3303    next_[iColumn] = COIN_INT_MAX;
3304  int i;
3305  int *keys = new int[numberSets_];
3306  for (i = 0; i < numberSets_; i++)
3307    keys[i] = COIN_INT_MAX;
3308  // set up chains
3309  for (i = 0; i < numberColumns; i++) {
3310    if (model->getStatus(i) == ClpSimplex::basic)
3311      mark[i] = 1;
3312    int iSet = backward_[i];
3313    if (iSet >= 0) {
3314      int iNext = keys[iSet];
3315      next_[i] = iNext;
3316      keys[iSet] = i;
3317    }
3318  }
3319  for (i = 0; i < numberSets_; i++) {
3320    int j;
3321    if (getStatus(i) != ClpSimplex::basic) {
3322      // make sure fixed if it is
3323      if (upper_[i] == lower_[i])
3324        setStatus(i, ClpSimplex::isFixed);
3325      // slack not key - choose one with smallest length
3326      int smallest = numberRows + 1;
3327      int key = -1;
3328      j = keys[i];
3329      if (j != COIN_INT_MAX) {
3330        while (1) {
3331          if (mark[j] && columnLength[j] < smallest && !gotBasis) {
3332            key = j;
3333            smallest = columnLength[j];
3334          }
3335          if (next_[j] != COIN_INT_MAX) {
3336            j = next_[j];
3337          } else {
3338            // correct end
3339            next_[j] = -(keys[i] + 1);
3340            break;
3341          }
3342        }
3343      } else {
3344        next_[i + numberColumns] = -(numberColumns + i + 1);
3345      }
3346      if (gotBasis)
3347        key = keyVariable_[i];
3348      if (key >= 0) {
3349        keyVariable_[i] = key;
3350      } else {
3351        // nothing basic - make slack key
3352        //((ClpGubMatrix *)this)->setStatus(i,ClpSimplex::basic);
3353        // fudge to avoid const problem
3354        status_[i] = 1;
3355      }
3356    } else {
3357      // slack key
3358      keyVariable_[i] = numberColumns + i;
3359      int j;
3360      double sol = 0.0;
3361      j = keys[i];
3362      if (j != COIN_INT_MAX) {
3363        while (1) {
3364          sol += columnSolution[j];
3365          if (next_[j] != COIN_INT_MAX) {
3366            j = next_[j];
3367          } else {
3368            // correct end
3369            next_[j] = -(keys[i] + 1);
3370            break;
3371          }
3372        }
3373      } else {
3374        next_[i + numberColumns] = -(numberColumns + i + 1);
3375      }
3376      if (sol > upper_[i] + tolerance) {
3377        setAbove(i);
3378      } else if (sol < lower_[i] - tolerance) {
3379        setBelow(i);
3380      } else {
3381        setFeasible(i);
3382      }
3383    }
3384    // Create next_
3385    int key = keyVariable_[i];
3386    redoSet(model, key, keys[i], i);
3387  }
3388  delete[] keys;
3389  delete[] mark;
3390  rhsOffset(model, true);
3391}
3392// redoes next_ for a set.
3393void ClpGubMatrix::redoSet(ClpSimplex *model, int newKey, int oldKey, int iSet)
3394{
3395  int numberColumns = model->numberColumns();
3396  int *save = next_ + numberColumns + numberSets_;
3397  int number = 0;
3398  int stop = -(oldKey + 1);
3399  int j = next_[oldKey];
3400  while (j != stop) {
3401    if (j < 0)
3402      j = -j - 1;
3403    if (j != newKey)
3404      save[number++] = j;
3405    j = next_[j];
3406  }
3407  // and add oldkey
3408  if (newKey != oldKey)
3409    save[number++] = oldKey;
3410  // now do basic
3411  int lastMarker = -(newKey + 1);
3412  keyVariable_[iSet] = newKey;
3413  next_[newKey] = lastMarker;
3414  int last = newKey;
3415  for (j = 0; j < number; j++) {
3416    int iColumn = save[j];
3417    if (iColumn < numberColumns) {
3418      if (model->getStatus(iColumn) == ClpSimplex::basic) {
3419        next_[last] = iColumn;
3420        next_[iColumn] = lastMarker;
3421        last = iColumn;
3422      }
3423    }
3424  }
3425  // now add in non-basic
3426  for (j = 0; j < number; j++) {
3427    int iColumn = save[j];
3428    if (iColumn < numberColumns) {
3429      if (model->getStatus(iColumn) != ClpSimplex::basic) {
3430        next_[last] = -(iColumn + 1);
3431        next_[iColumn] = lastMarker;
3432        last = iColumn;
3433      }
3434    }
3435  }
3436}
3437/* Returns effective RHS if it is being used.  This is used for long problems
3438   or big gub or anywhere where going through full columns is
3439   expensive.  This may re-compute */
3440double *
3441ClpGubMatrix::rhsOffset(ClpSimplex *model, bool forceRefresh, bool
3442#ifdef CLP_DEBUG
3443                                                                check
3444#endif
3445)
3446{
3447  //forceRefresh=true;
3448  if (rhsOffset_) {
3449#ifdef CLP_DEBUG
3450    if (check) {
3451      // no need - but check anyway
3452      // zero out basic
3453      int numberRows = model->numberRows();
3454      int numberColumns = model->numberColumns();
3455      double *solution = new double[numberColumns];
3456      double *rhs = new double[numberRows];
3457      CoinMemcpyN(model->solutionRegion(), numberColumns, solution);
3458      CoinZeroN(rhs, numberRows);
3459      int iRow;
3460      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3461        if (model->getColumnStatus(iColumn) == ClpSimplex::basic)
3462          solution[iColumn] = 0.0;
3463      }
3464      for (int iSet = 0; iSet < numberSets_; iSet++) {
3465        int iColumn = keyVariable_[iSet];
3466        if (iColumn < numberColumns)
3467          solution[iColumn] = 0.0;
3468      }
3469      times(-1.0, solution, rhs);
3470      delete[] solution;
3471      const double *columnSolution = model->solutionRegion();
3472      // and now subtract out non basic
3473      ClpSimplex::Status iStatus;
3474      for (int iSet = 0; iSet < numberSets_; iSet++) {
3475        int iColumn = keyVariable_[iSet];
3476        if (iColumn < numberColumns) {
3477          double b = 0.0;
3478          // key is structural - where is slack
3479          iStatus = getStatus(iSet);
3480          assert(iStatus != ClpSimplex::basic);
3481          if (iStatus == ClpSimplex::atLowerBound)
3482            b = lower_[iSet];
3483          else
3484            b = upper_[iSet];
3485          // subtract out others at bounds
3486          if ((gubType_ & 8) == 0) {
3487            int stop = -(iColumn + 1);
3488            int jColumn = next_[iColumn];
3489            // sum all non-basic variables - first skip basic
3490            while (jColumn >= 0)
3491              jColumn = next_[jColumn];
3492            while (jColumn != stop) {
3493              assert(jColumn < 0);
3494              jColumn = -jColumn - 1;
3495              b -= columnSolution[jColumn];
3496              jColumn = next_[jColumn];
3497            }
3498          }
3499          // subtract out
3500          ClpPackedMatrix::add(model, rhs, iColumn, -b);
3501        }
3502      }
3503      for (iRow = 0; iRow < numberRows; iRow++) {
3504        if (fabs(rhs[iRow] - rhsOffset_[iRow]) > 1.0e-3)
3505          printf("** bad effective %d - true %g old %g\n", iRow, rhs[iRow], rhsOffset_[iRow]);
3506      }
3507      delete[] rhs;
3508    }
3509#endif
3510    if (forceRefresh || (refreshFrequency_ && model->numberIterations() >= lastRefresh_ + refreshFrequency_)) {
3511      // zero out basic
3512      int numberRows = model->numberRows();
3513      int numberColumns = model->numberColumns();
3514      double *solution = new double[numberColumns];
3515      CoinMemcpyN(model->solutionRegion(), numberColumns, solution);
3516      CoinZeroN(rhsOffset_, numberRows);
3517      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3518        if (model->getColumnStatus(iColumn) == ClpSimplex::basic)
3519          solution[iColumn] = 0.0;
3520      }
3521      int iSet;
3522      for (iSet = 0; iSet < numberSets_; iSet++) {
3523        int iColumn = keyVariable_[iSet];
3524        if (iColumn < numberColumns)
3525          solution[iColumn] = 0.0;
3526      }
3527      times(-1.0, solution, rhsOffset_);
3528      delete[] solution;
3529      lastRefresh_ = model->numberIterations();
3530      const double *columnSolution = model->solutionRegion();
3531      // and now subtract out non basic
3532      ClpSimplex::Status iStatus;
3533      for (iSet = 0; iSet < numberSets_; iSet++) {
3534        int iColumn = keyVariable_[iSet];
3535        if (iColumn < numberColumns) {
3536          double b = 0.0;
3537          // key is structural - where is slack
3538          iStatus = getStatus(iSet);
3539          assert(iStatus != ClpSimplex::basic);
3540          if (iStatus == ClpSimplex::atLowerBound)
3541            b = lower_[iSet];
3542          else
3543            b = upper_[iSet];
3544          // subtract out others at bounds
3545          if ((gubType_ & 8) == 0) {
3546            int stop = -(iColumn + 1);
3547            int jColumn = next_[iColumn];
3548            // sum all non-basic variables - first skip basic
3549            while (jColumn >= 0)
3550              jColumn = next_[jColumn];
3551            while (jColumn != stop) {
3552              assert(jColumn < 0);
3553              jColumn = -jColumn - 1;
3554              b -= columnSolution[jColumn];
3555              jColumn = next_[jColumn];
3556            }
3557          }
3558          // subtract out
3559          if (b)
3560            ClpPackedMatrix::add(model, rhsOffset_, iColumn, -b);
3561        }
3562      }
3563    }
3564  }
3565  return rhsOffset_;
3566}
3567/*
3568   update information for a pivot (and effective rhs)
3569*/
3570int ClpGubMatrix::updatePivot(ClpSimplex *model, double oldInValue, double /*oldOutValue*/)
3571{
3572  int sequenceIn = model->sequenceIn();
3573  int sequenceOut = model->sequenceOut();
3574  double *solution = model->solutionRegion();
3575  int numberColumns = model->numberColumns();
3576  int numberRows = model->numberRows();
3577  int pivotRow = model->pivotRow();
3578  int iSetIn;
3579  // Correct sequence in
3580  trueSequenceIn_ = sequenceIn;
3581  if (sequenceIn < numberColumns) {
3582    iSetIn = backward_[sequenceIn];
3583  } else if (sequenceIn < numberColumns + numberRows) {
3584    iSetIn = -1;
3585  } else {
3586    iSetIn = gubSlackIn_;
3587    trueSequenceIn_ = numberColumns + numberRows + iSetIn;
3588  }
3589  int iSetOut = -1;
3590  trueSequenceOut_ = sequenceOut;
3591  if (sequenceOut < numberColumns) {
3592    iSetOut = backward_[sequenceOut];
3593  } else if (sequenceOut >= numberRows + numberColumns) {
3594    assert(pivotRow >= numberRows);
3595    int iExtra = pivotRow - numberRows;
3596    assert(iExtra >= 0);
3597    if (iSetOut < 0)
3598      iSetOut = fromIndex_[iExtra];
3599    else
3600      assert(iSetOut == fromIndex_[iExtra]);
3601    trueSequenceOut_ = numberColumns + numberRows + iSetOut;
3602  }
3603  if (rhsOffset_) {
3604    // update effective rhs
3605    if (sequenceIn == sequenceOut) {
3606      assert(sequenceIn < numberRows + numberColumns); // should be easy to deal with
3607      if (sequenceIn < numberColumns)
3608        add(model, rhsOffset_, sequenceIn, oldInValue - solution[sequenceIn]);
3609    } else {
3610      if (sequenceIn < numberColumns) {
3611        // we need to test if WILL be key
3612        ClpPackedMatrix::add(model, rhsOffset_, sequenceIn, oldInValue);
3613        if (iSetIn >= 0) {
3614          // old contribution to rhsOffset_
3615          int key = keyVariable_[iSetIn];
3616          if (key < numberColumns) {
3617            double oldB = 0.0;
3618            ClpSimplex::Status iStatus = getStatus(iSetIn);
3619            if (iStatus == ClpSimplex::atLowerBound)
3620              oldB = lower_[iSetIn];
3621            else
3622              oldB = upper_[iSetIn];
3623            // subtract out others at bounds
3624            if ((gubType_ & 8) == 0) {
3625              int stop = -(key + 1);
3626              int iColumn = next_[key];
3627              // skip basic
3628              while (iColumn >= 0)
3629                iColumn = next_[iColumn];
3630              // sum all non-key variables
3631              while (iColumn != stop) {
3632                assert(iColumn < 0);
3633                iColumn = -iColumn - 1;
3634                if (iColumn == sequenceIn)
3635                  oldB -= oldInValue;
3636                else if (iColumn != sequenceOut)
3637                  oldB -= solution[iColumn];
3638                iColumn = next_[iColumn];
3639              }
3640            }
3641            if (oldB)
3642              ClpPackedMatrix::add(model, rhsOffset_, key, oldB);
3643          }
3644        }
3645      } else if (sequenceIn < numberRows + numberColumns) {
3646        //rhsOffset_[sequenceIn-numberColumns] -= oldInValue;
3647      } else {
3648#ifdef CLP_DEBUG_PRINT
3649        printf("** in is key slack %d\n", sequenceIn);
3650#endif
3651        // old contribution to rhsOffset_
3652        int key = keyVariable_[iSetIn];
3653        if (key < numberColumns) {
3654          double oldB = 0.0;
3655          ClpSimplex::Status iStatus = getStatus(iSetIn);
3656          if (iStatus == ClpSimplex::atLowerBound)
3657            oldB = lower_[iSetIn];
3658          else
3659            oldB = upper_[iSetIn];
3660          // subtract out others at bounds
3661          if ((gubType_ & 8) == 0) {
3662            int stop = -(key + 1);
3663            int iColumn = next_[key];
3664            // skip basic
3665            while (iColumn >= 0)
3666              iColumn = next_[iColumn];
3667            // sum all non-key variables
3668            while (iColumn != stop) {
3669              assert(iColumn < 0);
3670              iColumn = -iColumn - 1;
3671              if (iColumn != sequenceOut)
3672                oldB -= solution[iColumn];
3673              iColumn = next_[iColumn];
3674            }
3675          }
3676          if (oldB)
3677            ClpPackedMatrix::add(model, rhsOffset_, key, oldB);
3678        }
3679      }
3680      if (sequenceOut < numberColumns) {
3681        ClpPackedMatrix::add(model, rhsOffset_, sequenceOut, -solution[sequenceOut]);
3682        if (iSetOut >= 0) {
3683          // old contribution to rhsOffset_
3684          int key = keyVariable_[iSetOut];
3685          if (key < numberColumns && iSetIn != iSetOut) {
3686            double oldB = 0.0;
3687            ClpSimplex::Status iStatus = getStatus(iSetOut);
3688            if (iStatus == ClpSimplex::atLowerBound)
3689              oldB = lower_[iSetOut];
3690            else
3691              oldB = upper_[iSetOut];
3692            // subtract out others at bounds
3693            if ((gubType_ & 8) == 0) {
3694              int stop = -(key + 1);
3695              int iColumn = next_[key];
3696              // skip basic
3697              while (iColumn >= 0)
3698                iColumn = next_[iColumn];
3699              // sum all non-key variables
3700              while (iColumn != stop) {
3701                assert(iColumn < 0);
3702                iColumn = -iColumn - 1;
3703                if (iColumn == sequenceIn)
3704                  oldB -= oldInValue;
3705                else if (iColumn != sequenceOut)
3706                  oldB -= solution[iColumn];
3707                iColumn = next_[iColumn];
3708              }
3709            }
3710            if (oldB)
3711              ClpPackedMatrix::add(model, rhsOffset_, key, oldB);
3712          }
3713        }
3714      } else if (sequenceOut < numberRows + numberColumns) {
3715        //rhsOffset_[sequenceOut-numberColumns] -= -solution[sequenceOut];
3716      } else {
3717#ifdef CLP_DEBUG_PRINT
3718        printf("** out is key slack %d\n", sequenceOut);
3719#endif
3720        assert(pivotRow >= numberRows);
3721      }
3722    }
3723  }
3724  int *pivotVariable = model->pivotVariable();
3725  // may need to deal with key
3726  // Also need coding to mark/allow key slack entering
3727  if (pivotRow >= numberRows) {
3728    assert(sequenceOut >= numberRows + numberColumns || sequenceOut == keyVariable_[iSetOut]);
3729#ifdef CLP_DEBUG_PRINT
3730    if (sequenceIn >= numberRows + numberColumns)
3731      printf("key slack %d in, set out %d\n", gubSlackIn_, iSetOut);
3732    printf("** danger - key out for set %d in %d (set %d)\n", iSetOut, sequenceIn,
3733      iSetIn);
3734#endif
3735    // if slack out mark correctly
3736    if (sequenceOut >= numberRows + numberColumns) {
3737      double value = model->valueOut();
3738      if (value == upper_[iSetOut]) {
3739        setStatus(iSetOut, ClpSimplex::atUpperBound);
3740      } else if (value == lower_[iSetOut]) {
3741        setStatus(iSetOut, ClpSimplex::atLowerBound);
3742      } else {
3743        if (fabs(value - upper_[iSetOut]) < fabs(value - lower_[iSetOut])) {
3744          setStatus(iSetOut, ClpSimplex::atUpperBound);
3745        } else {
3746          setStatus(iSetOut, ClpSimplex::atLowerBound);
3747        }
3748      }
3749      if (upper_[iSetOut] == lower_[iSetOut])
3750        setStatus(iSetOut, ClpSimplex::isFixed);
3751      setFeasible(iSetOut);
3752    }
3753    if (iSetOut == iSetIn) {
3754      // key swap
3755      int key;
3756      if (sequenceIn >= numberRows + numberColumns) {
3757        key = numberColumns + iSetIn;
3758        setStatus(iSetIn, ClpSimplex::basic);
3759      } else {
3760        key = sequenceIn;
3761      }
3762      redoSet(model, key, keyVariable_[iSetIn], iSetIn);
3763    } else {
3764      // key was chosen
3765      assert(possiblePivotKey_ >= 0 && possiblePivotKey_ < numberRows);
3766      int key = pivotVariable[possiblePivotKey_];
3767      // and set incoming here
3768      if (sequenceIn >= numberRows + numberColumns) {
3769        // slack in - so use old key
3770        sequenceIn = keyVariable_[iSetIn];
3771        model->setStatus(sequenceIn, ClpSimplex::basic);
3772        setStatus(iSetIn, ClpSimplex::basic);
3773        redoSet(model, iSetIn + numberColumns, keyVariable_[iSetIn], iSetIn);
3774      }
3775      //? do not do if iSetIn<0 ? as will be done later
3776      pivotVariable[possiblePivotKey_] = sequenceIn;
3777      if (sequenceIn < numberColumns)
3778        backToPivotRow_[sequenceIn] = possiblePivotKey_;
3779      redoSet(model, key, keyVariable_[iSetOut], iSetOut);
3780    }
3781  } else {
3782    if (sequenceOut < numberColumns) {
3783      if (iSetIn >= 0 && iSetOut == iSetIn) {
3784        // key not out - only problem is if slack in
3785        int key;
3786        if (sequenceIn >= numberRows + numberColumns) {
3787          key = numberColumns + iSetIn;
3788          setStatus(iSetIn, ClpSimplex::basic);
3789          assert(pivotRow < numberRows);
3790          // must swap with current key
3791          int key = keyVariable_[iSetIn];
3792          model->setStatus(key, ClpSimplex::basic);
3793          pivotVariable[pivotRow] = key;
3794          backToPivotRow_[key] = pivotRow;
3795        } else {
3796          key = keyVariable_[iSetIn];
3797        }
3798        redoSet(model, key, keyVariable_[iSetIn], iSetIn);
3799      } else if (iSetOut >= 0) {
3800        // just redo set
3801        int key = keyVariable_[iSetOut];
3802        ;
3803        redoSet(model, key, keyVariable_[iSetOut], iSetOut);
3804      }
3805    }
3806  }
3807  if (iSetIn >= 0 && iSetIn != iSetOut) {
3808    int key = keyVariable_[iSetIn];
3809    if (sequenceIn == numberColumns + 2 * numberRows) {
3810      // key slack in
3811      assert(pivotRow < numberRows);
3812      // must swap with current key
3813      model->setStatus(key, ClpSimplex::basic);
3814      pivotVariable[pivotRow] = key;
3815      backToPivotRow_[key] = pivotRow;
3816      setStatus(iSetIn, ClpSimplex::basic);
3817      key = iSetIn + numberColumns;
3818    }
3819    // redo set to allow for new one
3820    redoSet(model, key, keyVariable_[iSetIn], iSetIn);
3821  }
3822  // update pivot
3823  if (sequenceIn < numberColumns) {
3824    if (pivotRow < numberRows) {
3825      backToPivotRow_[sequenceIn] = pivotRow;
3826    } else {
3827      if (possiblePivotKey_ >= 0) {
3828        assert(possiblePivotKey_ < numberRows);
3829        backToPivotRow_[sequenceIn] = possiblePivotKey_;
3830        pivotVariable[possiblePivotKey_] = sequenceIn;
3831      }
3832    }
3833  } else if (sequenceIn >= numberRows + numberColumns) {
3834    // key in - something should have been done before
3835    int key = keyVariable_[iSetIn];
3836    assert(key == numberColumns + iSetIn);
3837    //pivotVariable[pivotRow]=key;
3838    //backToPivotRow_[key]=pivotRow;
3839    //model->setStatus(key,ClpSimplex::basic);
3840    //key=numberColumns+iSetIn;
3841    setStatus(iSetIn, ClpSimplex::basic);
3842    redoSet(model, key, keyVariable_[iSetIn], iSetIn);
3843  }
3844#ifdef CLP_DEBUG
3845  {
3846    char *xx = new char[numberColumns + numberRows];
3847    memset(xx, 0, numberRows + numberColumns);
3848    for (int i = 0; i < numberRows; i++) {
3849      int iPivot = pivotVariable[i];
3850      assert(iPivot < numberRows + numberColumns);
3851      assert(!xx[iPivot]);
3852      xx[iPivot] = 1;
3853      if (iPivot < numberColumns) {
3854        int iBack = backToPivotRow_[iPivot];
3855        assert(i == iBack);
3856      }
3857    }
3858    delete[] xx;
3859  }
3860#endif
3861  if (rhsOffset_) {
3862    // update effective rhs
3863    if (sequenceIn != sequenceOut) {
3864      if (sequenceIn < numberColumns) {
3865        if (iSetIn >= 0) {
3866          // new contribution to rhsOffset_
3867          int key = keyVariable_[iSetIn];
3868          if (key < numberColumns) {
3869            double newB = 0.0;
3870            ClpSimplex::Status iStatus = getStatus(iSetIn);
3871            if (iStatus == ClpSimplex::atLowerBound)
3872              newB = lower_[iSetIn];
3873            else
3874              newB = upper_[iSetIn];
3875            // subtract out others at bounds
3876            if ((gubType_ & 8) == 0) {
3877              int stop = -(key + 1);
3878              int iColumn = next_[key];
3879              // skip basic
3880              while (iColumn >= 0)
3881                iColumn = next_[iColumn];
3882              // sum all non-key variables
3883              while (iColumn != stop) {
3884                assert(iColumn < 0);
3885                iColumn = -iColumn - 1;
3886                newB -= solution[iColumn];
3887                iColumn = next_[iColumn];
3888              }
3889            }
3890            if (newB)
3891              ClpPackedMatrix::add(model, rhsOffset_, key, -newB);
3892          }
3893        }
3894      }
3895      if (iSetOut >= 0) {
3896        // new contribution to rhsOffset_
3897        int key = keyVariable_[iSetOut];
3898        if (key < numberColumns && iSetIn != iSetOut) {
3899          double newB = 0.0;
3900          ClpSimplex::Status iStatus = getStatus(iSetOut);
3901          if (iStatus == ClpSimplex::atLowerBound)
3902            newB = lower_[iSetOut];
3903          else
3904            newB = upper_[iSetOut];
3905          // subtract out others at bounds
3906          if ((gubType_ & 8) == 0) {
3907            int stop = -(key + 1);
3908            int iColumn = next_[key];
3909            // skip basic
3910            while (iColumn >= 0)
3911              iColumn = next_[iColumn];
3912            // sum all non-key variables
3913            while (iColumn != stop) {
3914              assert(iColumn < 0);
3915              iColumn = -iColumn - 1;
3916              newB -= solution[iColumn];
3917              iColumn = next_[iColumn];
3918            }
3919          }
3920          if (newB)
3921            ClpPackedMatrix::add(model, rhsOffset_, key, -newB);
3922        }
3923      }
3924    }
3925  }
3926#ifdef CLP_DEBUG
3927  // debug
3928  {
3929    int i;
3930    char *xxxx = new char[numberColumns];
3931    memset(xxxx, 0, numberColumns);
3932    for (i = 0; i < numberRows; i++) {
3933      int iPivot = pivotVariable[i];
3934      assert(model->getStatus(iPivot) == ClpSimplex::basic);
3935      if (iPivot < numberColumns && backward_[iPivot] >= 0)
3936        xxxx[iPivot] = 1;
3937    }
3938    double primalTolerance = model->primalTolerance();
3939    for (i = 0; i < numberSets_; i++) {
3940      int key = keyVariable_[i];
3941      double value = 0.0;
3942      // sum over all except key
3943      int iColumn = next_[key];
3944      // sum all non-key variables
3945      int k = 0;
3946      int stop = -(key + 1);
3947      while (iColumn != stop) {
3948        if (iColumn < 0)
3949          iColumn = -iColumn - 1;
3950        value += solution[iColumn];
3951        k++;
3952        assert(k < 100);
3953        assert(backward_[iColumn] == i);
3954        iColumn = next_[iColumn];
3955      }
3956      iColumn = next_[key];
3957      if (key < numberColumns) {
3958        // feasibility will be done later
3959        assert(getStatus(i) != ClpSimplex::basic);
3960        double sol;
3961        if (getStatus(i) == ClpSimplex::atUpperBound)
3962          sol = upper_[i] - value;
3963        else
3964          sol = lower_[i] - value;
3965        //printf("xx Value of key structural %d for set %d is %g - cost %g\n",key,i,sol,
3966        //     cost[key]);
3967        //if (fabs(sol-solution[key])>1.0e-3)
3968        //printf("** stored value was %g\n",solution[key]);
3969      } else {
3970        // slack is key
3971        double infeasibility = 0.0;
3972        if (value > upper_[i] + primalTolerance) {
3973          infeasibility = value - upper_[i] - primalTolerance;
3974          //setAbove(i);
3975        } else if (value < lower_[i] - primalTolerance) {
3976          infeasibility = lower_[i] - value - primalTolerance;
3977          //setBelow(i);
3978        } else {
3979          //setFeasible(i);
3980        }
3981        //printf("xx Value of key slack for set %d is %g\n",i,value);
3982      }
3983      while (iColumn >= 0) {
3984        assert(xxxx[iColumn]);
3985        xxxx[iColumn] = 0;
3986        iColumn = next_[iColumn];
3987      }
3988    }
3989    for (i = 0; i < numberColumns; i++) {
3990      if (i < numberColumns && backward_[i] >= 0) {
3991        assert(!xxxx[i] || i == keyVariable_[backward_[i]]);
3992      }
3993    }
3994    delete[] xxxx;
3995  }
3996#endif
3997  return 0;
3998}
3999// Switches off dj checking each factorization (for BIG models)
4000void ClpGubMatrix::switchOffCheck()
4001{
4002  noCheck_ = 0;
4003  infeasibilityWeight_ = 0.0;
4004}
4005// Correct sequence in and out to give true value
4006void ClpGubMatrix::correctSequence(const ClpSimplex * /*model*/, int &sequenceIn, int &sequenceOut)
4007{
4008  if (sequenceIn != -999) {
4009    sequenceIn = trueSequenceIn_;
4010    sequenceOut = trueSequenceOut_;
4011  }
4012}
4013
4014/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
4015*/
Note: See TracBrowser for help on using the repository browser.