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

Last change on this file since 754 was 754, checked in by andreasw, 14 years ago

first version

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