source: trunk/ClpGubMatrix.cpp @ 451

Last change on this file since 451 was 451, checked in by forrest, 16 years ago

Changes to make a difficult problem faster (dual) + quadratic

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 118.6 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/* If element NULL returns number of elements in column part of basis,
892   If not NULL fills in as well */
893CoinBigIndex
894ClpGubMatrix::fillBasis(ClpSimplex * model,
895                           const int * whichColumn, 
896                           int numberBasic,
897                           int & numberColumnBasic,
898                           int * indexRowU, int * indexColumnU,
899                           double * elementU) 
900{
901  int i;
902  int numberColumns = getNumCols();
903  const int * columnLength = matrix_->getVectorLengths(); 
904  int numberRows = getNumRows();
905  int saveNumberBasic=numberBasic;
906  assert (next_ ||!elementU) ;
907  CoinBigIndex numberElements=0;
908  int lastSet=-1;
909  int key=-1;
910  int keyLength=-1;
911  double * work = new double[numberRows];
912  CoinZeroN(work,numberRows);
913  char * mark = new char[numberRows];
914  CoinZeroN(mark,numberRows);
915  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
916  const int * row = matrix_->getIndices();
917  const double * elementByColumn = matrix_->getElements();
918  const double * rowScale = model->rowScale();
919  if (elementU!=NULL) {
920    if (0) {
921      printf("%d basic, %d columns\n",numberBasic,numberColumnBasic);
922      int i;
923      for (i=0;i<numberSets_;i++) {
924        int k=keyVariable_[i];
925        if (k<numberColumns) {
926          printf("key %d on set %d, %d elements\n",k,i,columnStart[k+1]-columnStart[k]);
927          for (int j=columnStart[k];j<columnStart[k+1];j++)
928            printf("row %d el %g\n",row[j],elementByColumn[j]);
929        } else {
930          printf("slack key on set %d\n",i);
931        }
932      }
933    }
934    // fill
935    if (!rowScale) {
936      // no scaling
937      for (i=0;i<numberColumnBasic;i++) {
938        int iColumn = whichColumn[i];
939        int iSet = backward_[iColumn];
940        int length = columnLength[iColumn];
941        if (0) {
942          int k=iColumn;
943          printf("column %d in set %d, %d elements\n",k,iSet,columnStart[k+1]-columnStart[k]);
944          for (int j=columnStart[k];j<columnStart[k+1];j++)
945            printf("row %d el %g\n",row[j],elementByColumn[j]);
946        }
947        CoinBigIndex j;
948        if (iSet<0||keyVariable_[iSet]>=numberColumns) {
949          for (j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
950            double value = elementByColumn[j];
951            if (fabs(value)>1.0e-20) {
952              int iRow = row[j];
953              indexRowU[numberElements]=iRow;
954              indexColumnU[numberElements]=numberBasic;
955              elementU[numberElements++]=value;
956            }
957          }
958          numberBasic++;
959        } else {
960          // in gub set
961          if (iColumn!=keyVariable_[iSet]) {
962            // not key
963            if (lastSet!=iSet) {
964              // erase work
965              if (key>=0) {
966                for (j=columnStart[key];j<columnStart[key]+keyLength;j++) {
967                  int iRow=row[j];
968                  work[iRow]=0.0;
969                  mark[iRow]=0;
970                }
971              }
972              key=keyVariable_[iSet];
973              lastSet=iSet;
974              keyLength = columnLength[key];
975              for (j=columnStart[key];j<columnStart[key]+keyLength;j++) {
976                int iRow=row[j];
977                work[iRow]=elementByColumn[j];
978                mark[iRow]=1;
979              }
980            }
981            for (j=columnStart[iColumn];j<columnStart[iColumn]+length;j++) {
982              int iRow = row[j];
983              double value=elementByColumn[j];
984              if (mark[iRow]) {
985                mark[iRow]=0;
986                double keyValue = work[iRow];
987                value -= keyValue;
988              }
989              if (fabs(value)>1.0e-20) {
990                indexRowU[numberElements]=iRow;
991                indexColumnU[numberElements]=numberBasic;
992                elementU[numberElements++]=value;
993              }
994            }
995            for (j=columnStart[key];j<columnStart[key]+keyLength;j++) {
996              int iRow = row[j];
997              if (mark[iRow]) {
998                double value = -work[iRow];
999                if (fabs(value)>1.0e-20) {
1000                  indexRowU[numberElements]=iRow;
1001                  indexColumnU[numberElements]=numberBasic;
1002                  elementU[numberElements++]=value;
1003                }
1004              } else {
1005                // just put back mark
1006                mark[iRow]=1;
1007              }
1008            }
1009            numberBasic++;
1010          }
1011        }
1012      }
1013    } else {
1014      // scaling
1015      const double * columnScale = model->columnScale();
1016      for (i=0;i<numberColumnBasic;i++) {
1017        int iColumn = whichColumn[i];
1018        int iSet = backward_[iColumn];
1019        int length = columnLength[iColumn];
1020        CoinBigIndex j;
1021        if (iSet<0||keyVariable_[iSet]>=numberColumns) {
1022          double scale = columnScale[iColumn];
1023          for (j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
1024            int iRow = row[j];
1025            double value = elementByColumn[j]*scale*rowScale[iRow];
1026            if (fabs(value)>1.0e-20) {
1027              indexRowU[numberElements]=iRow;
1028              indexColumnU[numberElements]=numberBasic;
1029              elementU[numberElements++]=value;
1030            }
1031          }
1032          numberBasic++;
1033        } else {
1034          // in gub set
1035          if (iColumn!=keyVariable_[iSet]) {
1036            double scale = columnScale[iColumn];
1037            // not key
1038            if (lastSet<iSet) {
1039              // erase work
1040              if (key>=0) {
1041                for (j=columnStart[key];j<columnStart[key]+keyLength;j++) {
1042                  int iRow=row[j];
1043                  work[iRow]=0.0;
1044                  mark[iRow]=0;
1045                }
1046              }
1047              key=keyVariable_[iSet];
1048              lastSet=iSet;
1049              keyLength = columnLength[key];
1050              double scale = columnScale[key];
1051              for (j=columnStart[key];j<columnStart[key]+keyLength;j++) {
1052                int iRow=row[j];
1053                work[iRow]=elementByColumn[j]*scale*rowScale[iRow];
1054                mark[iRow]=1;
1055              }
1056            }
1057            for (j=columnStart[iColumn];j<columnStart[iColumn]+length;j++) {
1058              int iRow = row[j];
1059              double value=elementByColumn[j]*scale*rowScale[iRow];
1060              if (mark[iRow]) {
1061                mark[iRow]=0;
1062                double keyValue = work[iRow];
1063                value -= keyValue;
1064              }
1065              if (fabs(value)>1.0e-20) {
1066                indexRowU[numberElements]=iRow;
1067                indexColumnU[numberElements]=numberBasic;
1068                elementU[numberElements++]=value;
1069              }
1070            }
1071            for (j=columnStart[key];j<columnStart[key]+keyLength;j++) {
1072              int iRow = row[j];
1073              if (mark[iRow]) {
1074                double value = -work[iRow];
1075                if (fabs(value)>1.0e-20) {
1076                  indexRowU[numberElements]=iRow;
1077                  indexColumnU[numberElements]=numberBasic;
1078                  elementU[numberElements++]=value;
1079                }
1080              } else {
1081                // just put back mark
1082                mark[iRow]=1;
1083              }
1084            }
1085            numberBasic++;
1086          }
1087        }
1088      }
1089    }
1090  } else {
1091    // just count
1092    if (!rowScale) {
1093      for (i=0;i<numberColumnBasic;i++) {
1094        int iColumn = whichColumn[i];
1095        int iSet = backward_[iColumn];
1096        int length = columnLength[iColumn];
1097        if (iSet<0||keyVariable_[iSet]>=numberColumns) {
1098          numberElements += length;
1099          numberBasic++;
1100        } else {
1101          // in gub set
1102          if (iColumn!=keyVariable_[iSet]) {
1103            numberBasic++;
1104            CoinBigIndex j;
1105            // not key
1106            if (lastSet<iSet) {
1107              // erase work
1108              if (key>=0) {
1109                for (j=columnStart[key];j<columnStart[key]+keyLength;j++)
1110                  work[row[j]]=0.0;
1111              }
1112              key=keyVariable_[iSet];
1113              lastSet=iSet;
1114              keyLength = columnLength[key];
1115              for (j=columnStart[key];j<columnStart[key]+keyLength;j++)
1116                work[row[j]]=elementByColumn[j];
1117            }
1118            int extra=keyLength;
1119            for (j=columnStart[iColumn];j<columnStart[iColumn]+length;j++) {
1120              int iRow = row[j];
1121              double keyValue = work[iRow];
1122              double value=elementByColumn[j];
1123              if (!keyValue) {
1124                if (fabs(value)>1.0e-20)
1125                  extra++;
1126              } else {
1127                value -= keyValue;
1128                if (fabs(value)<=1.0e-20)
1129                  extra--;
1130              }
1131            }
1132            numberElements+=extra;
1133          }
1134        }
1135      }
1136    } else {
1137      // scaled
1138      const double * columnScale = model->columnScale();
1139      for (i=0;i<numberColumnBasic;i++) {
1140        int iColumn = whichColumn[i];
1141        int iSet = backward_[iColumn];
1142        int length = columnLength[iColumn];
1143        if (iSet<0||keyVariable_[iSet]>=numberColumns) {
1144          numberElements += length;
1145          numberBasic++;
1146        } else {
1147          // in gub set
1148          if (iColumn!=keyVariable_[iSet]) {
1149            numberBasic++;
1150            CoinBigIndex j;
1151            double scale = columnScale[iColumn];
1152            // not key
1153            if (lastSet<iSet) {
1154              // erase work
1155              if (key>=0) {
1156                for (j=columnStart[key];j<columnStart[key]+keyLength;j++)
1157                  work[row[j]]=0.0;
1158              }
1159              key=keyVariable_[iSet];
1160              lastSet=iSet;
1161              keyLength = columnLength[key];
1162              double scale = columnScale[key];
1163              for (j=columnStart[key];j<columnStart[key]+keyLength;j++) {
1164                int iRow = row[j];
1165                work[iRow]=elementByColumn[j]*scale*rowScale[iRow];
1166              }
1167            }
1168            int extra=keyLength;
1169            for (j=columnStart[iColumn];j<columnStart[iColumn]+length;j++) {
1170              int iRow = row[j];
1171              double keyValue = work[iRow];
1172              double value=elementByColumn[j]*scale*rowScale[iRow];
1173              if (!keyValue) {
1174                if (fabs(value)>1.0e-20)
1175                  extra++;
1176              } else {
1177                value -= keyValue;
1178                if (fabs(value)<=1.0e-20)
1179                  extra--;
1180              }
1181            }
1182            numberElements+=extra;
1183          }
1184        }
1185      }
1186    }
1187  }
1188  delete [] work;
1189  delete [] mark;
1190  // update number of column basic
1191  numberColumnBasic = numberBasic-saveNumberBasic;
1192  return numberElements;
1193}
1194/* Unpacks a column into an CoinIndexedvector
1195 */
1196void 
1197ClpGubMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
1198                   int iColumn) const 
1199{
1200  assert (iColumn<model->numberColumns());
1201  // Do packed part
1202  ClpPackedMatrix::unpack(model,rowArray,iColumn);
1203  int iSet = backward_[iColumn];
1204  if (iSet>=0) {
1205    int iBasic = keyVariable_[iSet];
1206    if (iBasic <model->numberColumns()) {
1207      add(model,rowArray,iBasic,-1.0);
1208    }
1209  }
1210}
1211/* Unpacks a column into a CoinIndexedVector
1212** in packed format
1213Note that model is NOT const.  Bounds and objective could
1214be modified if doing column generation (just for this variable) */
1215void 
1216ClpGubMatrix::unpackPacked(ClpSimplex * model,
1217                            CoinIndexedVector * rowArray,
1218                            int iColumn) const
1219{
1220  int numberColumns = model->numberColumns();
1221  if (iColumn<numberColumns) {
1222    // Do packed part
1223    ClpPackedMatrix::unpackPacked(model,rowArray,iColumn);
1224    int iSet = backward_[iColumn];
1225    if (iSet>=0) {
1226      // columns are in order
1227      int iBasic = keyVariable_[iSet];
1228      if (iBasic <numberColumns) {
1229        int number = rowArray->getNumElements();
1230        const double * rowScale = model->rowScale();
1231        const int * row = matrix_->getIndices();
1232        const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1233        const int * columnLength = matrix_->getVectorLengths(); 
1234        const double * elementByColumn = matrix_->getElements();
1235        double * array = rowArray->denseVector();
1236        int * index = rowArray->getIndices();
1237        CoinBigIndex i;
1238        int numberOld=number;
1239        int lastIndex=0;
1240        int next=index[lastIndex];
1241        if (!rowScale) {
1242          for (i=columnStart[iBasic];
1243               i<columnStart[iBasic]+columnLength[iBasic];i++) {
1244            int iRow = row[i];
1245            while (iRow>next) {
1246              lastIndex++;
1247              if (lastIndex==numberOld)
1248                next=matrix_->getNumRows();
1249              else
1250                next=index[lastIndex];
1251            }
1252            if (iRow<next) {
1253              array[number]=-elementByColumn[i];
1254              index[number++]=iRow;
1255            } else {
1256              assert (iRow==next);
1257              array[lastIndex] -= elementByColumn[i];
1258              if (!array[lastIndex])
1259                array[lastIndex]=1.0e-100;
1260            }
1261          }
1262        } else {
1263          // apply scaling
1264          double scale = model->columnScale()[iBasic];
1265          for (i=columnStart[iBasic];
1266               i<columnStart[iBasic]+columnLength[iBasic];i++) {
1267            int iRow = row[i];
1268            while (iRow>next) {
1269              lastIndex++;
1270              if (lastIndex==numberOld)
1271                next=matrix_->getNumRows();
1272              else
1273                next=index[lastIndex];
1274            }
1275            if (iRow<next) {
1276              array[number]=-elementByColumn[i]*scale*rowScale[iRow];
1277              index[number++]=iRow;
1278            } else {
1279              assert (iRow==next);
1280              array[lastIndex] -=elementByColumn[i]*scale*rowScale[iRow];
1281              if (!array[lastIndex])
1282                array[lastIndex]=1.0e-100;
1283            }
1284          }
1285        }
1286        rowArray->setNumElements(number);
1287      }
1288    }
1289  } else {
1290    // key slack entering
1291    int iBasic = keyVariable_[gubSlackIn_];
1292    assert (iBasic <numberColumns);
1293    int number = 0;
1294    const double * rowScale = model->rowScale();
1295    const int * row = matrix_->getIndices();
1296    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1297    const int * columnLength = matrix_->getVectorLengths(); 
1298    const double * elementByColumn = matrix_->getElements();
1299    double * array = rowArray->denseVector();
1300    int * index = rowArray->getIndices();
1301    CoinBigIndex i;
1302    if (!rowScale) {
1303      for (i=columnStart[iBasic];
1304           i<columnStart[iBasic]+columnLength[iBasic];i++) {
1305        int iRow = row[i];
1306        array[number]=elementByColumn[i];
1307        index[number++]=iRow;
1308      }
1309    } else {
1310      // apply scaling
1311      double scale = model->columnScale()[iBasic];
1312      for (i=columnStart[iBasic];
1313           i<columnStart[iBasic]+columnLength[iBasic];i++) {
1314        int iRow = row[i];
1315        array[number]=elementByColumn[i]*scale*rowScale[iRow];
1316        index[number++]=iRow;
1317      }
1318    }
1319    rowArray->setNumElements(number);
1320    rowArray->setPacked();
1321  }
1322}
1323/* Adds multiple of a column into an CoinIndexedvector
1324      You can use quickAdd to add to vector */
1325void 
1326ClpGubMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
1327                   int iColumn, double multiplier) const 
1328{
1329  assert (iColumn<model->numberColumns());
1330  // Do packed part
1331  ClpPackedMatrix::add(model,rowArray,iColumn,multiplier);
1332  int iSet = backward_[iColumn];
1333  if (iSet>=0&&iColumn!=keyVariable_[iSet]) {
1334    ClpPackedMatrix::add(model,rowArray,keyVariable_[iSet],-multiplier);
1335  }
1336}
1337/* Adds multiple of a column into an array */
1338void 
1339ClpGubMatrix::add(const ClpSimplex * model,double * array,
1340                   int iColumn, double multiplier) const 
1341{
1342  assert (iColumn<model->numberColumns());
1343  // Do packed part
1344  ClpPackedMatrix::add(model,array,iColumn,multiplier);
1345  if (iColumn<model->numberColumns()) {
1346    int iSet = backward_[iColumn];
1347    if (iSet>=0&&iColumn!=keyVariable_[iSet]&&keyVariable_[iSet]<model->numberColumns()) {
1348      ClpPackedMatrix::add(model,array,keyVariable_[iSet],-multiplier);
1349    }
1350  }
1351}
1352// Partial pricing
1353void 
1354ClpGubMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
1355                              int & bestSequence, int & numberWanted)
1356{
1357  numberWanted=currentWanted_;
1358  if (numberSets_) {
1359    // Do packed part before gub
1360    int numberColumns = matrix_->getNumCols();
1361    double ratio = ((double) firstGub_)/((double) numberColumns);
1362    ClpPackedMatrix::partialPricing(model,startFraction*ratio,
1363                                    endFraction*ratio,bestSequence,numberWanted);
1364    if (numberWanted||minimumGoodReducedCosts_<-1) {
1365      // do gub
1366      const double * element =matrix_->getElements();
1367      const int * row = matrix_->getIndices();
1368      const CoinBigIndex * startColumn = matrix_->getVectorStarts();
1369      const int * length = matrix_->getVectorLengths();
1370      const double * rowScale = model->rowScale();
1371      const double * columnScale = model->columnScale();
1372      int iSequence;
1373      CoinBigIndex j;
1374      double tolerance=model->currentDualTolerance();
1375      double * reducedCost = model->djRegion();
1376      const double * duals = model->dualRowSolution();
1377      const double * cost = model->costRegion();
1378      double bestDj;
1379      int numberColumns = model->numberColumns();
1380      int numberRows = model->numberRows();
1381      if (bestSequence>=0)
1382        bestDj = fabs(this->reducedCost(model,bestSequence));
1383      else
1384        bestDj=tolerance;
1385      int sequenceOut = model->sequenceOut();
1386      int saveSequence = bestSequence;
1387      int startG = firstGub_+ (int) (startFraction* (lastGub_-firstGub_));
1388      int endG = firstGub_+ (int) (endFraction* (lastGub_-firstGub_));
1389      endG = CoinMin(lastGub_,endG+1);
1390      // If nothing found yet can go all the way to end
1391      int endAll = endG;
1392      if (bestSequence<0&&!startG)
1393        endAll = lastGub_;
1394      int minSet = minimumObjectsScan_<0 ? 5 : minimumObjectsScan_; 
1395      int minNeg = minimumGoodReducedCosts_==-1 ? 5 : minimumGoodReducedCosts_;
1396      int nSets=0;
1397      int iSet = -1;
1398      double djMod=0.0;
1399      double infeasibilityCost = model->infeasibilityCost();
1400      if (rowScale) {
1401        double bestDjMod=0.0;
1402        // scaled
1403        for (iSequence=startG;iSequence<endAll;iSequence++) {
1404          if (numberWanted+minNeg<originalWanted_&&nSets>minSet) {
1405            // give up
1406            numberWanted=0;
1407            break;
1408          } else if (iSequence==endG&&bestSequence>=0) {
1409            break;
1410          }
1411          if (backward_[iSequence]!=iSet) {
1412            // get pi on gub row
1413            iSet = backward_[iSequence];
1414            if (iSet>=0) {
1415              nSets++;
1416              int iBasic = keyVariable_[iSet];
1417              if (iBasic>=numberColumns) {
1418                djMod = - weight(iSet)*infeasibilityCost;
1419              } else {
1420                // get dj without
1421                assert (model->getStatus(iBasic)==ClpSimplex::basic);
1422                djMod=0.0;
1423                // scaled
1424                for (j=startColumn[iBasic];
1425                     j<startColumn[iBasic]+length[iBasic];j++) {
1426                  int jRow=row[j];
1427                  djMod -= duals[jRow]*element[j]*rowScale[jRow];
1428                }
1429                // allow for scaling
1430                djMod +=  cost[iBasic]/columnScale[iBasic];
1431                // See if gub slack possible - dj is djMod
1432                if (getStatus(iSet)==ClpSimplex::atLowerBound) {
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                } else if (getStatus(iSet)==ClpSimplex::atUpperBound) {
1450                  double value = djMod;
1451                  if (value>tolerance) {
1452                    numberWanted--;
1453                    if (value>bestDj) {
1454                      // check flagged variable and correct dj
1455                      if (!flagged(iSet)) {
1456                        bestDj=value;
1457                        bestSequence = numberRows+numberColumns+iSet;
1458                        bestDjMod = djMod;
1459                      } else {
1460                        // just to make sure we don't exit before got something
1461                        numberWanted++;
1462                        abort();
1463                      }
1464                    }
1465                  }
1466                }
1467              }
1468            } else {
1469              // not in set
1470              djMod=0.0;
1471            }
1472          }
1473          if (iSequence!=sequenceOut) {
1474            double value;
1475            ClpSimplex::Status status = model->getStatus(iSequence);
1476           
1477            switch(status) {
1478             
1479            case ClpSimplex::basic:
1480            case ClpSimplex::isFixed:
1481              break;
1482            case ClpSimplex::isFree:
1483            case ClpSimplex::superBasic:
1484              value=-djMod;
1485              // scaled
1486              for (j=startColumn[iSequence];
1487                   j<startColumn[iSequence]+length[iSequence];j++) {
1488                int jRow=row[j];
1489                value -= duals[jRow]*element[j]*rowScale[jRow];
1490              }
1491              value = fabs(cost[iSequence] +value*columnScale[iSequence]);
1492              if (value>FREE_ACCEPT*tolerance) {
1493                numberWanted--;
1494                // we are going to bias towards free (but only if reasonable)
1495                value *= FREE_BIAS;
1496                if (value>bestDj) {
1497                  // check flagged variable and correct dj
1498                  if (!model->flagged(iSequence)) {
1499                    bestDj=value;
1500                    bestSequence = iSequence;
1501                    bestDjMod = djMod;
1502                  } else {
1503                    // just to make sure we don't exit before got something
1504                    numberWanted++;
1505                  }
1506                }
1507              }
1508              break;
1509            case ClpSimplex::atUpperBound:
1510              value=-djMod;
1511              // scaled
1512              for (j=startColumn[iSequence];
1513                   j<startColumn[iSequence]+length[iSequence];j++) {
1514                int jRow=row[j];
1515                value -= duals[jRow]*element[j]*rowScale[jRow];
1516              }
1517              value = cost[iSequence] +value*columnScale[iSequence];
1518              if (value>tolerance) {
1519                numberWanted--;
1520                if (value>bestDj) {
1521                  // check flagged variable and correct dj
1522                  if (!model->flagged(iSequence)) {
1523                    bestDj=value;
1524                    bestSequence = iSequence;
1525                    bestDjMod = djMod;
1526                  } else {
1527                    // just to make sure we don't exit before got something
1528                    numberWanted++;
1529                  }
1530                }
1531              }
1532              break;
1533            case ClpSimplex::atLowerBound:
1534              value=-djMod;
1535              // scaled
1536              for (j=startColumn[iSequence];
1537                   j<startColumn[iSequence]+length[iSequence];j++) {
1538                int jRow=row[j];
1539                value -= duals[jRow]*element[j]*rowScale[jRow];
1540              }
1541              value = -(cost[iSequence] +value*columnScale[iSequence]);
1542              if (value>tolerance) {
1543                numberWanted--;
1544                if (value>bestDj) {
1545                  // check flagged variable and correct dj
1546                  if (!model->flagged(iSequence)) {
1547                    bestDj=value;
1548                    bestSequence = iSequence;
1549                    bestDjMod = djMod;
1550                  } else {
1551                    // just to make sure we don't exit before got something
1552                    numberWanted++;
1553                  }
1554                }
1555              }
1556              break;
1557            }
1558          }
1559          if (!numberWanted)
1560            break;
1561        }
1562        if (bestSequence!=saveSequence) {
1563          if (bestSequence<numberRows+numberColumns) {
1564            // recompute dj
1565            double value=bestDjMod;
1566            // scaled
1567            for (j=startColumn[bestSequence];
1568                 j<startColumn[bestSequence]+length[bestSequence];j++) {
1569              int jRow=row[j];
1570              value -= duals[jRow]*element[j]*rowScale[jRow];
1571            }
1572            reducedCost[bestSequence] = cost[bestSequence] +value*columnScale[bestSequence];
1573            gubSlackIn_=-1;
1574          } else {
1575            // slack - make last column
1576            gubSlackIn_= bestSequence - numberRows - numberColumns;
1577            bestSequence = numberColumns + 2*numberRows;
1578            reducedCost[bestSequence] = bestDjMod;
1579            model->setStatus(bestSequence,getStatus(gubSlackIn_));
1580            if (getStatus(gubSlackIn_)==ClpSimplex::atUpperBound)
1581              model->solutionRegion()[bestSequence] = upper_[gubSlackIn_];
1582            else
1583              model->solutionRegion()[bestSequence] = lower_[gubSlackIn_];
1584            model->lowerRegion()[bestSequence] = lower_[gubSlackIn_];
1585            model->upperRegion()[bestSequence] = upper_[gubSlackIn_];
1586            model->costRegion()[bestSequence] = 0.0;
1587          }
1588          savedBestSequence_ = bestSequence;
1589          savedBestDj_ = reducedCost[savedBestSequence_];
1590        }
1591      } else {
1592        double bestDjMod=0.0;
1593        //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
1594        //     startG,endG,numberWanted);
1595        for (iSequence=startG;iSequence<endG;iSequence++) {
1596          if (numberWanted+minNeg<originalWanted_&&nSets>minSet) {
1597            // give up
1598            numberWanted=0;
1599            break;
1600          } else if (iSequence==endG&&bestSequence>=0) {
1601            break;
1602          }
1603          if (backward_[iSequence]!=iSet) {
1604            // get pi on gub row
1605            iSet = backward_[iSequence];
1606            if (iSet>=0) {
1607              nSets++;
1608              int iBasic = keyVariable_[iSet];
1609              if (iBasic>=numberColumns) {
1610                djMod = - weight(iSet)*infeasibilityCost;
1611              } else {
1612                // get dj without
1613                assert (model->getStatus(iBasic)==ClpSimplex::basic);
1614                djMod=0.0;
1615               
1616                for (j=startColumn[iBasic];
1617                     j<startColumn[iBasic]+length[iBasic];j++) {
1618                  int jRow=row[j];
1619                  djMod -= duals[jRow]*element[j];
1620                }
1621                djMod += cost[iBasic];
1622                // See if gub slack possible - dj is djMod
1623                if (getStatus(iSet)==ClpSimplex::atLowerBound) {
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                } else if (getStatus(iSet)==ClpSimplex::atUpperBound) {
1641                  double value = djMod;
1642                  if (value>tolerance) {
1643                    numberWanted--;
1644                    if (value>bestDj) {
1645                      // check flagged variable and correct dj
1646                      if (!flagged(iSet)) {
1647                        bestDj=value;
1648                        bestSequence = numberRows+numberColumns+iSet;
1649                        bestDjMod = djMod;
1650                      } else {
1651                        // just to make sure we don't exit before got something
1652                        numberWanted++;
1653                        abort();
1654                      }
1655                    }
1656                  }
1657                }
1658              }
1659            } else {
1660              // not in set
1661              djMod=0.0;
1662            }
1663          }
1664          if (iSequence!=sequenceOut) {
1665            double value;
1666            ClpSimplex::Status status = model->getStatus(iSequence);
1667           
1668            switch(status) {
1669             
1670            case ClpSimplex::basic:
1671            case ClpSimplex::isFixed:
1672              break;
1673            case ClpSimplex::isFree:
1674            case ClpSimplex::superBasic:
1675              value=cost[iSequence]-djMod;
1676              for (j=startColumn[iSequence];
1677                   j<startColumn[iSequence]+length[iSequence];j++) {
1678                int jRow=row[j];
1679                value -= duals[jRow]*element[j];
1680              }
1681              value = fabs(value);
1682              if (value>FREE_ACCEPT*tolerance) {
1683                numberWanted--;
1684                // we are going to bias towards free (but only if reasonable)
1685                value *= FREE_BIAS;
1686                if (value>bestDj) {
1687                  // check flagged variable and correct dj
1688                  if (!model->flagged(iSequence)) {
1689                    bestDj=value;
1690                    bestSequence = iSequence;
1691                    bestDjMod = djMod;
1692                  } else {
1693                    // just to make sure we don't exit before got something
1694                    numberWanted++;
1695                  }
1696                }
1697              }
1698              break;
1699            case ClpSimplex::atUpperBound:
1700              value=cost[iSequence]-djMod;
1701              for (j=startColumn[iSequence];
1702                   j<startColumn[iSequence]+length[iSequence];j++) {
1703                int jRow=row[j];
1704                value -= duals[jRow]*element[j];
1705              }
1706              if (value>tolerance) {
1707                numberWanted--;
1708                if (value>bestDj) {
1709                  // check flagged variable and correct dj
1710                  if (!model->flagged(iSequence)) {
1711                    bestDj=value;
1712                    bestSequence = iSequence;
1713                    bestDjMod = djMod;
1714                  } else {
1715                    // just to make sure we don't exit before got something
1716                    numberWanted++;
1717                  }
1718                }
1719              }
1720              break;
1721            case ClpSimplex::atLowerBound:
1722              value=cost[iSequence]-djMod;
1723              for (j=startColumn[iSequence];
1724                   j<startColumn[iSequence]+length[iSequence];j++) {
1725                int jRow=row[j];
1726                value -= duals[jRow]*element[j];
1727              }
1728              value = -value;
1729              if (value>tolerance) {
1730                numberWanted--;
1731                if (value>bestDj) {
1732                  // check flagged variable and correct dj
1733                  if (!model->flagged(iSequence)) {
1734                    bestDj=value;
1735                    bestSequence = iSequence;
1736                    bestDjMod = djMod;
1737                  } else {
1738                    // just to make sure we don't exit before got something
1739                    numberWanted++;
1740                  }
1741                }
1742              }
1743              break;
1744            }
1745          }
1746          if (!numberWanted)
1747            break;
1748        }
1749        if (bestSequence!=saveSequence) {
1750          if (bestSequence<numberRows+numberColumns) {
1751            // recompute dj
1752            double value=cost[bestSequence]-bestDjMod;
1753            for (j=startColumn[bestSequence];
1754                 j<startColumn[bestSequence]+length[bestSequence];j++) {
1755              int jRow=row[j];
1756              value -= duals[jRow]*element[j];
1757            }
1758            //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,bestDjMod);
1759            reducedCost[bestSequence] = value;
1760            gubSlackIn_=-1;
1761          } else {
1762            // slack - make last column
1763            gubSlackIn_= bestSequence - numberRows - numberColumns;
1764            bestSequence = numberColumns + 2*numberRows;
1765            reducedCost[bestSequence] = bestDjMod;
1766            //printf("price slack %d - gubpi %g\n",gubSlackIn_,bestDjMod);
1767            model->setStatus(bestSequence,getStatus(gubSlackIn_));
1768            if (getStatus(gubSlackIn_)==ClpSimplex::atUpperBound)
1769              model->solutionRegion()[bestSequence] = upper_[gubSlackIn_];
1770            else
1771              model->solutionRegion()[bestSequence] = lower_[gubSlackIn_];
1772            model->lowerRegion()[bestSequence] = lower_[gubSlackIn_];
1773            model->upperRegion()[bestSequence] = upper_[gubSlackIn_];
1774            model->costRegion()[bestSequence] = 0.0;
1775          }
1776        }
1777      }
1778      // See if may be finished
1779      if (startG==firstGub_&&bestSequence<0)
1780        infeasibilityWeight_=model_->infeasibilityCost();
1781      else if (bestSequence>=0) 
1782        infeasibilityWeight_=-1.0;
1783    }
1784    if (numberWanted) {
1785      // Do packed part after gub
1786      double offset = ((double) lastGub_)/((double) numberColumns); 
1787      double ratio = ((double) numberColumns)/((double) numberColumns)-offset;
1788      double start2 = offset + ratio*startFraction;
1789      double end2 = CoinMin(1.0,offset + ratio*endFraction+1.0e-6);
1790      ClpPackedMatrix::partialPricing(model,start2,end2,bestSequence,numberWanted);
1791    }
1792  } else {
1793    // no gub
1794    ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
1795  }
1796  if (bestSequence>=0)
1797    infeasibilityWeight_=-1.0; // not optimal
1798  currentWanted_=numberWanted;
1799}
1800/* expands an updated column to allow for extra rows which the main
1801   solver does not know about and returns number added.
1802*/
1803int 
1804ClpGubMatrix::extendUpdated(ClpSimplex * model,CoinIndexedVector * update, int mode)
1805{
1806  // I think we only need to bother about sets with two in basis or incoming set
1807  int number = update->getNumElements();
1808  double * array = update->denseVector();
1809  int * index = update->getIndices();
1810  int i;
1811  assert (!number||update->packedMode());
1812  int * pivotVariable = model->pivotVariable();
1813  int numberRows = model->numberRows();
1814  int numberColumns = model->numberColumns();
1815  int numberTotal = numberRows+numberColumns;
1816  int sequenceIn = model->sequenceIn();
1817  int returnCode=0;
1818  int iSetIn;
1819  if (sequenceIn<numberColumns) {
1820    iSetIn = backward_[sequenceIn];
1821    gubSlackIn_ = -1; // in case set
1822  } else if (sequenceIn<numberRows+numberColumns) {
1823    iSetIn = -1;
1824    gubSlackIn_ = -1; // in case set
1825  } else {
1826    iSetIn = gubSlackIn_;
1827  }
1828  double * lower = model->lowerRegion();
1829  double * upper = model->upperRegion();
1830  double * cost = model->costRegion();
1831  double * solution = model->solutionRegion();
1832  int number2=number;
1833  if (!mode) {
1834    double primalTolerance = model->primalTolerance();
1835    double infeasibilityCost = model->infeasibilityCost();
1836    // extend
1837    saveNumber_ = number;
1838    for (i=0;i<number;i++) {
1839      int iRow = index[i];
1840      int iPivot = pivotVariable[iRow];
1841      if (iPivot<numberColumns) {
1842        int iSet = backward_[iPivot];
1843        if (iSet>=0) {
1844          // two (or more) in set
1845          int iIndex =toIndex_[iSet];
1846          double otherValue=array[i];
1847          double value;
1848          if (iIndex<0) {
1849            toIndex_[iSet]=number2;
1850            int iNew = number2-number;
1851            fromIndex_[number2-number]=iSet;
1852            iIndex=number2;
1853            index[number2]=numberRows+iNew;
1854            // do key stuff
1855            int iKey = keyVariable_[iSet];
1856            if (iKey<numberColumns) {
1857              // Save current cost of key
1858              changeCost_[number2-number] = cost[iKey];
1859              if (iSet!=iSetIn)
1860                value = 0.0;
1861              else if (iSetIn!=gubSlackIn_)
1862                value = 1.0;
1863              else
1864                value =-1.0;
1865              pivotVariable[numberRows+iNew]=iKey;
1866              // Do I need to recompute?
1867              double sol;
1868              assert (getStatus(iSet)!=ClpSimplex::basic);
1869              if (getStatus(iSet)==ClpSimplex::atLowerBound)
1870                sol = lower_[iSet];
1871              else
1872                sol = upper_[iSet];
1873              if ((gubType_&8)!=0) {
1874                int iColumn =next_[iKey];
1875                // sum all non-key variables
1876                while(iColumn>=0) {
1877                  sol -= solution[iColumn];
1878                  iColumn=next_[iColumn];
1879                }
1880              } else {
1881                int stop = -(iKey+1);
1882                int iColumn =next_[iKey];
1883                // sum all non-key variables
1884                while(iColumn!=stop) {
1885                  if (iColumn<0)
1886                    iColumn = -iColumn-1;
1887                  sol -= solution[iColumn];
1888                  iColumn=next_[iColumn];
1889                }
1890              }
1891              solution[iKey]=sol;
1892              if (model->algorithm()>0)
1893                model->nonLinearCost()->setOne(iKey,sol);
1894              //assert (fabs(sol-solution[iKey])<1.0e-3);
1895            } else {
1896              // gub slack is basic
1897              // Save current cost of key
1898              changeCost_[number2-number]= -weight(iSet)*infeasibilityCost;
1899              otherValue = - otherValue; //allow for - sign on slack
1900              if (iSet!=iSetIn)
1901                value = 0.0;
1902              else
1903                value = -1.0;
1904              pivotVariable[numberRows+iNew]=iNew+numberTotal;
1905              model->djRegion()[iNew+numberTotal]=0.0;
1906              double sol=0.0;
1907              if ((gubType_&8)!=0) {
1908                int iColumn =next_[iKey];
1909                // sum all non-key variables
1910                while(iColumn>=0) {
1911                  sol += solution[iColumn];
1912                  iColumn=next_[iColumn];
1913                }
1914              } else {
1915                int stop = -(iKey+1);
1916                int iColumn =next_[iKey];
1917                // sum all non-key variables
1918                while(iColumn!=stop) {
1919                  if (iColumn<0)
1920                    iColumn = -iColumn-1;
1921                  sol += solution[iColumn];
1922                  iColumn=next_[iColumn];
1923                }
1924              }
1925              solution[iNew+numberTotal]=sol;
1926              // and do cost in nonLinearCost
1927              if (model->algorithm()>0)
1928                model->nonLinearCost()->setOne(iNew+numberTotal,sol,lower_[iSet],upper_[iSet]);
1929              if (sol>upper_[iSet]+primalTolerance) {
1930                setAbove(iSet);
1931                lower[iNew+numberTotal]=upper_[iSet];
1932                upper[iNew+numberTotal]=COIN_DBL_MAX;
1933              } else if (sol<lower_[iSet]-primalTolerance) {
1934                setBelow(iSet);
1935                lower[iNew+numberTotal]=-COIN_DBL_MAX;
1936                upper[iNew+numberTotal]=lower_[iSet];
1937              } else {
1938                setFeasible(iSet);
1939                lower[iNew+numberTotal]=lower_[iSet];
1940                upper[iNew+numberTotal]=upper_[iSet];
1941              }
1942              cost[iNew+numberTotal]=weight(iSet)*infeasibilityCost;
1943            }
1944            number2++;
1945          } else {
1946            value = array[iIndex];
1947            int iKey = keyVariable_[iSet];
1948            if (iKey>=numberColumns) 
1949              otherValue = - otherValue; //allow for - sign on slack
1950          }
1951          value -= otherValue;
1952          array[iIndex]=value;
1953        }
1954      }
1955    }
1956    if (iSetIn>=0&&toIndex_[iSetIn]<0) {
1957      // Do incoming
1958      update->setPacked(); // just in case no elements
1959      toIndex_[iSetIn]=number2;
1960      int iNew = number2-number;
1961      fromIndex_[number2-number]=iSetIn;
1962      // Save current cost of key
1963      double currentCost;
1964      int key=keyVariable_[iSetIn];
1965      if (key<numberColumns) 
1966        currentCost = cost[key];
1967      else
1968        currentCost = -weight(iSetIn)*infeasibilityCost;
1969      changeCost_[number2-number]=currentCost;
1970      index[number2]=numberRows+iNew;
1971      // do key stuff
1972      int iKey = keyVariable_[iSetIn];
1973      if (iKey<numberColumns) {
1974        if (gubSlackIn_<0)
1975          array[number2]=1.0;
1976        else
1977          array[number2]=-1.0;
1978        pivotVariable[numberRows+iNew]=iKey;
1979        // Do I need to recompute?
1980        double sol;
1981        assert (getStatus(iSetIn)!=ClpSimplex::basic);
1982        if (getStatus(iSetIn)==ClpSimplex::atLowerBound)
1983          sol = lower_[iSetIn];
1984        else
1985          sol = upper_[iSetIn];
1986        if ((gubType_&8)!=0) {
1987          int iColumn =next_[iKey];
1988          // sum all non-key variables
1989          while(iColumn>=0) {
1990            sol -= solution[iColumn];
1991            iColumn=next_[iColumn];
1992          }
1993        } else {
1994          // bounds exist - sum over all except key
1995          int stop = -(iKey+1);
1996          int iColumn =next_[iKey];
1997          // sum all non-key variables
1998          while(iColumn!=stop) {
1999            if (iColumn<0)
2000              iColumn = -iColumn-1;
2001            sol -= solution[iColumn];
2002            iColumn=next_[iColumn];
2003          }
2004        }
2005        solution[iKey]=sol;
2006        if (model->algorithm()>0)
2007          model->nonLinearCost()->setOne(iKey,sol);
2008        //assert (fabs(sol-solution[iKey])<1.0e-3);
2009      } else {
2010        // gub slack is basic
2011        array[number2]=-1.0;
2012        pivotVariable[numberRows+iNew]=iNew+numberTotal;
2013        model->djRegion()[iNew+numberTotal]=0.0;
2014        double sol=0.0;
2015        if ((gubType_&8)!=0) {
2016          int iColumn =next_[iKey];
2017          // sum all non-key variables
2018          while(iColumn>=0) {
2019            sol += solution[iColumn];
2020            iColumn=next_[iColumn];
2021          }
2022        } else {
2023          // bounds exist - sum over all except key
2024          int stop = -(iKey+1);
2025          int iColumn =next_[iKey];
2026          // sum all non-key variables
2027          while(iColumn!=stop) {
2028            if (iColumn<0)
2029              iColumn = -iColumn-1;
2030            sol += solution[iColumn];
2031            iColumn=next_[iColumn];
2032          }
2033        }
2034        solution[iNew+numberTotal]=sol;
2035        // and do cost in nonLinearCost
2036        if (model->algorithm()>0)
2037          model->nonLinearCost()->setOne(iNew+numberTotal,sol,lower_[iSetIn],upper_[iSetIn]);
2038        if (sol>upper_[iSetIn]+primalTolerance) {
2039          setAbove(iSetIn);
2040          lower[iNew+numberTotal]=upper_[iSetIn];
2041          upper[iNew+numberTotal]=COIN_DBL_MAX;
2042        } else if (sol<lower_[iSetIn]-primalTolerance) {
2043          setBelow(iSetIn);
2044          lower[iNew+numberTotal]=-COIN_DBL_MAX;
2045          upper[iNew+numberTotal]=lower_[iSetIn];
2046        } else {
2047          setFeasible(iSetIn);
2048          lower[iNew+numberTotal]=lower_[iSetIn];
2049          upper[iNew+numberTotal]=upper_[iSetIn];
2050        }
2051        cost[iNew+numberTotal]=weight(iSetIn)*infeasibilityCost;
2052      }
2053      number2++;
2054    }
2055    // mark end
2056    fromIndex_[number2-number]=-1;
2057    returnCode = number2-number;
2058    // make sure lower_ upper_ adjusted
2059    synchronize(model,9);
2060  } else {
2061    // take off?
2062    if (number>saveNumber_) {
2063      // clear
2064      double theta = model->theta();
2065      double * solution = model->solutionRegion();
2066      for (i=saveNumber_;i<number;i++) {
2067        int iRow = index[i];
2068        int iColumn = pivotVariable[iRow];
2069#ifdef CLP_DEBUG_PRINT
2070        printf("Column %d (set %d) lower %g, upper %g - alpha %g - old value %g, new %g (theta %g)\n",
2071               iColumn,fromIndex_[i-saveNumber_],lower[iColumn],upper[iColumn],array[i],
2072               solution[iColumn],solution[iColumn]-model->theta()*array[i],model->theta());
2073#endif
2074        double value = array[i];
2075        array[i]=0.0;
2076        int iSet=fromIndex_[i-saveNumber_];
2077        toIndex_[iSet]=-1;
2078        if (iSet==iSetIn&&iColumn<numberColumns) {
2079          // update as may need value
2080          solution[iColumn] -= theta*value;
2081        }
2082      }
2083    }
2084#ifdef CLP_DEBUG
2085    for (i=0;i<numberSets_;i++)
2086      assert(toIndex_[i]==-1);
2087#endif
2088    number2= saveNumber_;
2089  }
2090  update->setNumElements(number2);
2091  return returnCode;
2092}
2093/*
2094     utility primal function for dealing with dynamic constraints
2095     mode=n see ClpGubMatrix.hpp for definition
2096     Remember to update here when settled down
2097*/
2098void 
2099ClpGubMatrix::primalExpanded(ClpSimplex * model,int mode)
2100{
2101  int numberColumns = model->numberColumns();
2102  switch (mode) {
2103    // If key variable then slot in gub rhs so will get correct contribution
2104  case 0:
2105    {
2106      int i;
2107      double * solution = model->solutionRegion();
2108      ClpSimplex::Status iStatus;
2109      for (i=0;i<numberSets_;i++) {
2110        int iColumn = keyVariable_[i];
2111        if (iColumn<numberColumns) {
2112          // key is structural - where is slack
2113          iStatus = getStatus(i);
2114          assert (iStatus!=ClpSimplex::basic);
2115          if (iStatus==ClpSimplex::atLowerBound)
2116            solution[iColumn]=lower_[i];
2117          else
2118            solution[iColumn]=upper_[i];
2119        }
2120      }
2121    }
2122    break;
2123    // Compute values of key variables
2124  case 1:
2125    {
2126      int i;
2127      double * solution = model->solutionRegion();
2128      ClpSimplex::Status iStatus;
2129      //const int * columnLength = matrix_->getVectorLengths();
2130      //const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2131      //const int * row = matrix_->getIndices();
2132      //const double * elementByColumn = matrix_->getElements();
2133      //int * pivotVariable = model->pivotVariable();
2134      sumPrimalInfeasibilities_=0.0;
2135      numberPrimalInfeasibilities_=0;
2136      double primalTolerance = model->primalTolerance();
2137      double relaxedTolerance=primalTolerance;
2138      // we can't really trust infeasibilities if there is primal error
2139      double error = CoinMin(1.0e-3,model->largestPrimalError());
2140      // allow tolerance at least slightly bigger than standard
2141      relaxedTolerance = relaxedTolerance +  error;
2142      // but we will be using difference
2143      relaxedTolerance -= primalTolerance;
2144      sumOfRelaxedPrimalInfeasibilities_ = 0.0;
2145      for (i=0;i<numberSets_;i++) { // Could just be over basics (esp if no bounds)
2146        int kColumn = keyVariable_[i];
2147        double value=0.0;
2148        if ((gubType_&8)!=0) {
2149          int iColumn =next_[kColumn];
2150          // sum all non-key variables
2151          while(iColumn>=0) {
2152            value+=solution[iColumn];
2153            iColumn=next_[iColumn];
2154          }
2155        } else {
2156          // bounds exist - sum over all except key
2157          int stop = -(kColumn+1);
2158          int iColumn =next_[kColumn];
2159          // sum all non-key variables
2160          while(iColumn!=stop) {
2161            if (iColumn<0)
2162              iColumn = -iColumn-1;
2163            value += solution[iColumn];
2164            iColumn=next_[iColumn];
2165          }
2166        }
2167        if (kColumn<numberColumns) {
2168          // make sure key is basic - so will be skipped in values pass
2169          model->setStatus(kColumn,ClpSimplex::basic);
2170          // feasibility will be done later
2171          assert (getStatus(i)!=ClpSimplex::basic);
2172          if (getStatus(i)==ClpSimplex::atUpperBound)
2173            solution[kColumn] = upper_[i]-value;
2174          else
2175            solution[kColumn] = lower_[i]-value;
2176          //printf("Value of key structural %d for set %d is %g\n",kColumn,i,solution[kColumn]);
2177        } else {
2178          // slack is key
2179          iStatus = getStatus(i);
2180          assert (iStatus==ClpSimplex::basic);
2181          double infeasibility=0.0;
2182          if (value>upper_[i]+primalTolerance) {
2183            infeasibility=value-upper_[i]-primalTolerance;
2184            setAbove(i);
2185          } else if (value<lower_[i]-primalTolerance) {
2186            infeasibility=lower_[i]-value-primalTolerance ;
2187            setBelow(i);
2188          } else {
2189            setFeasible(i);
2190          }
2191          //printf("Value of key slack for set %d is %g\n",i,value);
2192          if (infeasibility>0.0) {
2193            sumPrimalInfeasibilities_ += infeasibility;
2194            if (infeasibility>relaxedTolerance) 
2195              sumOfRelaxedPrimalInfeasibilities_ += infeasibility;
2196            numberPrimalInfeasibilities_ ++;
2197          }
2198        }
2199      }
2200    }
2201    break;
2202    // Report on infeasibilities of key variables
2203  case 2:
2204    {
2205      model->setSumPrimalInfeasibilities(model->sumPrimalInfeasibilities()+
2206                                         sumPrimalInfeasibilities_);
2207      model->setNumberPrimalInfeasibilities(model->numberPrimalInfeasibilities()+
2208                                         numberPrimalInfeasibilities_);
2209      model->setSumOfRelaxedPrimalInfeasibilities(model->sumOfRelaxedPrimalInfeasibilities()+
2210                                         sumOfRelaxedPrimalInfeasibilities_);
2211    }
2212    break;
2213  }
2214}
2215/*
2216     utility dual function for dealing with dynamic constraints
2217     mode=n see ClpGubMatrix.hpp for definition
2218     Remember to update here when settled down
2219*/
2220void 
2221ClpGubMatrix::dualExpanded(ClpSimplex * model,
2222                            CoinIndexedVector * array,
2223                            double * other,int mode)
2224{
2225  switch (mode) {
2226    // modify costs before transposeUpdate
2227  case 0:
2228    {
2229      int i;
2230      double * cost = model->costRegion();
2231      ClpSimplex::Status iStatus;
2232      // not dual values yet
2233      assert (!other);
2234      //double * work = array->denseVector();
2235      double infeasibilityCost = model->infeasibilityCost();
2236      int * pivotVariable = model->pivotVariable();
2237      int numberRows = model->numberRows();
2238      int numberColumns = model->numberColumns();
2239      for (i=0;i<numberRows;i++) {
2240        int iPivot = pivotVariable[i];
2241        if (iPivot<numberColumns) {
2242          int iSet = backward_[iPivot];
2243          if (iSet>=0) {
2244            int kColumn = keyVariable_[iSet];
2245            double costValue;
2246            if (kColumn<numberColumns) {
2247              // structural has cost
2248              costValue = cost[kColumn];
2249            } else {
2250              // slack is key
2251              iStatus = getStatus(iSet);
2252              assert (iStatus==ClpSimplex::basic);
2253              // negative as -1.0 for slack
2254              costValue=-weight(iSet)*infeasibilityCost;
2255            }
2256            array->add(i,-costValue); // was work[i]-costValue
2257          }
2258        }
2259      }
2260    }
2261    break;
2262    // create duals for key variables (without check on dual infeasible)
2263  case 1:
2264    {
2265      // If key slack then dual 0.0 (if feasible)
2266      // dj for key is zero so that defines dual on set
2267      int i;
2268      double * dj = model->djRegion();
2269      int numberColumns = model->numberColumns();
2270      double infeasibilityCost = model->infeasibilityCost();
2271      for (i=0;i<numberSets_;i++) {
2272        int kColumn = keyVariable_[i];
2273        if (kColumn<numberColumns) {
2274          // dj without set
2275          double value = dj[kColumn];
2276          // Now subtract out from all
2277          dj[kColumn] =0.0;
2278          int iColumn =next_[kColumn];
2279          // modify all non-key variables
2280          while(iColumn>=0) {
2281            dj[iColumn]-=value;
2282            iColumn=next_[iColumn];
2283          }
2284        } else {
2285          // slack key - may not be feasible
2286          assert (getStatus(i)==ClpSimplex::basic);
2287          // negative as -1.0 for slack
2288          double value=-weight(i)*infeasibilityCost;
2289          if (value) {
2290            int iColumn =next_[kColumn];
2291            // modify all non-key variables basic
2292            while(iColumn>=0) {
2293              dj[iColumn]-=value;
2294              iColumn=next_[iColumn];
2295            }
2296          }
2297        }
2298      }
2299    }
2300    break;
2301    // as 1 but check slacks and compute djs
2302  case 2:
2303    {
2304      // If key slack then dual 0.0
2305      // If not then slack could be dual infeasible
2306      // dj for key is zero so that defines dual on set
2307      int i;
2308      // make sure fromIndex will not confuse pricing
2309      fromIndex_[0]=-1;
2310      possiblePivotKey_=-1;
2311      // Create array
2312      int numberColumns = model->numberColumns();
2313      int * pivotVariable = model->pivotVariable();
2314      int numberRows = model->numberRows();
2315      for (i=0;i<numberRows;i++) {
2316        int iPivot = pivotVariable[i];
2317        if (iPivot<numberColumns)
2318          backToPivotRow_[iPivot]=i;
2319      }
2320      if (noCheck_>=0) {
2321        if (infeasibilityWeight_!=model->infeasibilityCost()) {
2322          // don't bother checking
2323          sumDualInfeasibilities_=100.0;
2324          numberDualInfeasibilities_=1;
2325          sumOfRelaxedDualInfeasibilities_ =100.0;
2326          return;
2327        }
2328      }
2329      double * dj = model->djRegion();
2330      double * dual = model->dualRowSolution();
2331      double * cost = model->costRegion();
2332      ClpSimplex::Status iStatus;
2333      const int * columnLength = matrix_->getVectorLengths(); 
2334      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2335      const int * row = matrix_->getIndices();
2336      const double * elementByColumn = matrix_->getElements();
2337      double infeasibilityCost = model->infeasibilityCost();
2338      sumDualInfeasibilities_=0.0;
2339      numberDualInfeasibilities_=0;
2340      double dualTolerance = model->dualTolerance();
2341      double relaxedTolerance=dualTolerance;
2342      // we can't really trust infeasibilities if there is dual error
2343      double error = CoinMin(1.0e-3,model->largestDualError());
2344      // allow tolerance at least slightly bigger than standard
2345      relaxedTolerance = relaxedTolerance +  error;
2346      // but we will be using difference
2347      relaxedTolerance -= dualTolerance;
2348      sumOfRelaxedDualInfeasibilities_ = 0.0;
2349      for (i=0;i<numberSets_;i++) {
2350        int kColumn = keyVariable_[i];
2351        if (kColumn<numberColumns) {
2352          // dj without set
2353          double value = cost[kColumn];
2354          for (CoinBigIndex j=columnStart[kColumn];
2355               j<columnStart[kColumn]+columnLength[kColumn];j++) {
2356            int iRow = row[j];
2357            value -= dual[iRow]*elementByColumn[j];
2358          }
2359          // Now subtract out from all
2360          dj[kColumn] -= value;
2361          int stop = -(kColumn+1);
2362          kColumn = next_[kColumn];
2363          while (kColumn!=stop) {
2364            if (kColumn<0)
2365              kColumn = -kColumn-1;
2366            double djValue = dj[kColumn]-value;
2367            dj[kColumn] = djValue;;
2368            double infeasibility=0.0;
2369            iStatus = model->getStatus(kColumn);
2370            if (iStatus==ClpSimplex::atLowerBound) {
2371              if (djValue<-dualTolerance) 
2372                infeasibility=-djValue-dualTolerance;
2373            } else if (iStatus==ClpSimplex::atUpperBound) {
2374              // at upper bound
2375              if (djValue>dualTolerance) 
2376                infeasibility=djValue-dualTolerance;
2377            }
2378            if (infeasibility>0.0) {
2379              sumDualInfeasibilities_ += infeasibility;
2380              if (infeasibility>relaxedTolerance) 
2381                sumOfRelaxedDualInfeasibilities_ += infeasibility;
2382              numberDualInfeasibilities_ ++;
2383            }
2384            kColumn = next_[kColumn];
2385          }
2386          // check slack
2387          iStatus = getStatus(i);
2388          assert (iStatus!=ClpSimplex::basic);
2389          double infeasibility=0.0;
2390          // dj of slack is -(-1.0)value
2391          if (iStatus==ClpSimplex::atLowerBound) {
2392            if (value<-dualTolerance) 
2393              infeasibility=-value-dualTolerance;
2394          } else if (iStatus==ClpSimplex::atUpperBound) {
2395            // at upper bound
2396            if (value>dualTolerance) 
2397              infeasibility=value-dualTolerance;
2398          }
2399          if (infeasibility>0.0) {
2400            sumDualInfeasibilities_ += infeasibility;
2401            if (infeasibility>relaxedTolerance) 
2402              sumOfRelaxedDualInfeasibilities_ += infeasibility;
2403            numberDualInfeasibilities_ ++;
2404          }
2405        } else {
2406          // slack key - may not be feasible
2407          assert (getStatus(i)==ClpSimplex::basic);
2408          // negative as -1.0 for slack
2409          double value=-weight(i)*infeasibilityCost;
2410          if (value) {
2411          // Now subtract out from all
2412            int kColumn = i+numberColumns;
2413            int stop = -(kColumn+1);
2414            kColumn = next_[kColumn];
2415            while (kColumn!=stop) {
2416              if (kColumn<0)
2417                kColumn = -kColumn-1;
2418              double djValue = dj[kColumn]-value;
2419              dj[kColumn] = djValue;;
2420              double infeasibility=0.0;
2421              iStatus = model->getStatus(kColumn);
2422              if (iStatus==ClpSimplex::atLowerBound) {
2423                if (djValue<-dualTolerance) 
2424                  infeasibility=-djValue-dualTolerance;
2425              } else if (iStatus==ClpSimplex::atUpperBound) {
2426                // at upper bound
2427                if (djValue>dualTolerance) 
2428                  infeasibility=djValue-dualTolerance;
2429              }
2430              if (infeasibility>0.0) {
2431                sumDualInfeasibilities_ += infeasibility;
2432                if (infeasibility>relaxedTolerance) 
2433                  sumOfRelaxedDualInfeasibilities_ += infeasibility;
2434                numberDualInfeasibilities_ ++;
2435              }
2436              kColumn = next_[kColumn];
2437            }
2438          }
2439        }
2440      }
2441      // and get statistics for column generation
2442      synchronize(model,4);
2443      infeasibilityWeight_=-1.0;
2444    }
2445    break;
2446    // Report on infeasibilities of key variables
2447  case 3:
2448    {
2449      model->setSumDualInfeasibilities(model->sumDualInfeasibilities()+
2450                                         sumDualInfeasibilities_);
2451      model->setNumberDualInfeasibilities(model->numberDualInfeasibilities()+
2452                                         numberDualInfeasibilities_);
2453      model->setSumOfRelaxedDualInfeasibilities(model->sumOfRelaxedDualInfeasibilities()+
2454                                         sumOfRelaxedDualInfeasibilities_);
2455    }
2456    break;
2457    // modify costs before transposeUpdate for partial pricing
2458  case 4:
2459    {
2460      // First compute new costs etc for interesting gubs
2461      int iLook=0;
2462      int iSet=fromIndex_[0];
2463      double primalTolerance = model->primalTolerance();
2464      const double * cost = model->costRegion();
2465      double * solution = model->solutionRegion();
2466      double infeasibilityCost = model->infeasibilityCost();
2467      int numberColumns = model->numberColumns();
2468      int numberChanged=0;
2469      int * pivotVariable = model->pivotVariable();
2470      while (iSet>=0) {
2471        int key=keyVariable_[iSet];
2472        double value=0.0;
2473        // sum over all except key
2474        if ((gubType_&8)!=0) {
2475          int iColumn =next_[key];
2476          // sum all non-key variables
2477          while(iColumn>=0) {
2478            value += solution[iColumn];
2479            iColumn=next_[iColumn];
2480          }
2481        } else {
2482          // bounds exist - sum over all except key
2483          int stop = -(key+1);
2484          int iColumn =next_[key];
2485          // sum all non-key variables
2486          while(iColumn!=stop) {
2487            if (iColumn<0)
2488              iColumn = -iColumn-1;
2489            value += solution[iColumn];
2490            iColumn=next_[iColumn];
2491          }
2492        }
2493        double costChange;
2494        double oldCost = changeCost_[iLook];
2495        if (key<numberColumns) {
2496          assert (getStatus(iSet)!=ClpSimplex::basic);
2497          double sol;
2498          if (getStatus(iSet)==ClpSimplex::atUpperBound)
2499            sol = upper_[iSet]-value;
2500          else
2501            sol = lower_[iSet]-value;
2502          solution[key]=sol;
2503          // fix up cost
2504          model->nonLinearCost()->setOne(key,sol);
2505#ifdef CLP_DEBUG_PRINT
2506          printf("yy Value of key structural %d for set %d is %g - cost %g old cost %g\n",key,iSet,sol,
2507                 cost[key],oldCost);
2508#endif
2509          costChange = cost[key]-oldCost;
2510        } else {
2511          // slack is key
2512          if (value>upper_[iSet]+primalTolerance) {
2513            setAbove(iSet);
2514          } else if (value<lower_[iSet]-primalTolerance) {
2515            setBelow(iSet);
2516          } else {
2517            setFeasible(iSet);
2518          }
2519          // negative as -1.0 for slack
2520          costChange=-weight(iSet)*infeasibilityCost-oldCost;
2521#ifdef CLP_DEBUG_PRINT
2522          printf("yy Value of key slack for set %d is %g - cost %g old cost %g\n",iSet,value,
2523                 weight(iSet)*infeasibilityCost,oldCost);
2524#endif
2525        }
2526        if (costChange) {
2527          fromIndex_[numberChanged]=iSet;
2528          toIndex_[iSet]=numberChanged;
2529          changeCost_[numberChanged++]=costChange;
2530        }
2531        iSet = fromIndex_[++iLook];
2532      }
2533      if (numberChanged||possiblePivotKey_>=0) {
2534        // first do those in list already
2535        int number = array->getNumElements();
2536        array->setPacked();
2537        int i;
2538        double * work = array->denseVector();
2539        int * which = array->getIndices();
2540        for (i=0;i<number;i++) {
2541          int iRow = which[i];
2542          int iPivot = pivotVariable[iRow];
2543          if (iPivot<numberColumns) {
2544            int iSet = backward_[iPivot];
2545            if (iSet>=0&&toIndex_[iSet]>=0) {
2546              double newValue = work[i]+changeCost_[toIndex_[iSet]];
2547              if (!newValue)
2548                newValue =1.0e-100;
2549              work[i]=newValue;
2550              // mark as done
2551              backward_[iPivot]=-1;
2552            }
2553          }
2554          if (possiblePivotKey_==iRow) {
2555            double newValue = work[i]-model->dualIn();
2556            if (!newValue)
2557              newValue =1.0e-100;
2558            work[i]=newValue;
2559            possiblePivotKey_=-1;
2560          }
2561        }
2562        // now do rest and clean up
2563        for (i=0;i<numberChanged;i++) {
2564          int iSet = fromIndex_[i];
2565          int key=keyVariable_[iSet];
2566          int iColumn = next_[key];
2567          double change = changeCost_[i];
2568          while (iColumn>=0) {
2569            if (backward_[iColumn]>=0) {
2570              int iRow = backToPivotRow_[iColumn];
2571              assert (iRow>=0);
2572              work[number] = change;
2573              if (possiblePivotKey_==iRow) {
2574                double newValue = work[number]-model->dualIn();
2575                if (!newValue)
2576                  newValue =1.0e-100;
2577                work[number]=newValue;
2578                possiblePivotKey_=-1;
2579              }
2580              which[number++]=iRow;
2581            } else {
2582              // reset
2583              backward_[iColumn]=iSet;
2584            }
2585            iColumn=next_[iColumn];
2586          }
2587          toIndex_[iSet]=-1;
2588        }
2589        if (possiblePivotKey_>=0) {
2590          work[number] = -model->dualIn();
2591          which[number++]=possiblePivotKey_;
2592          possiblePivotKey_=-1;
2593        }
2594        fromIndex_[0]=-1;
2595        array->setNumElements(number);
2596      }
2597    }
2598    break;
2599  }
2600}
2601// This is local to Gub to allow synchronization when status is good
2602int 
2603ClpGubMatrix::synchronize(ClpSimplex * model, int mode)
2604{
2605  return 0;
2606}
2607/*
2608     general utility function for dealing with dynamic constraints
2609     mode=n see ClpGubMatrix.hpp for definition
2610     Remember to update here when settled down
2611*/
2612int
2613ClpGubMatrix::generalExpanded(ClpSimplex * model,int mode,int &number)
2614{
2615  int returnCode=0;
2616  int numberColumns = model->numberColumns();
2617  switch (mode) {
2618    // Fill in pivotVariable but not for key variables
2619  case 0:
2620    {
2621      if (!next_ ) {
2622        // do ordering
2623        assert (!rhsOffset_);
2624        // create and do gub crash
2625        useEffectiveRhs(model,false);
2626      }
2627      int i;
2628      int numberBasic=number;
2629      // Use different array so can build from true pivotVariable_
2630      //int * pivotVariable = model->pivotVariable();
2631      int * pivotVariable = model->rowArray(0)->getIndices();
2632      for (i=0;i<numberColumns;i++) {
2633        if (model->getColumnStatus(i) == ClpSimplex::basic) {
2634          int iSet = backward_[i];
2635          if (iSet<0||i!=keyVariable_[iSet])
2636            pivotVariable[numberBasic++]=i;
2637        }
2638      }
2639      number = numberBasic;
2640      if (model->numberIterations())
2641        assert (number==model->numberRows());
2642    }
2643    break;
2644    // Make all key variables basic
2645  case 1:
2646    {
2647      int i;
2648      for (i=0;i<numberSets_;i++) {
2649        int iColumn = keyVariable_[i];
2650        if (iColumn<numberColumns)
2651          model->setColumnStatus(iColumn,ClpSimplex::basic);
2652      }
2653    }
2654    break;
2655    // Do initial extra rows + maximum basic
2656  case 2:
2657    {
2658      returnCode= getNumRows()+1;
2659      number = model->numberRows()+numberSets_;
2660    }
2661    break;
2662    // Before normal replaceColumn
2663  case 3:
2664    {
2665      int sequenceIn = model->sequenceIn();
2666      int sequenceOut = model->sequenceOut();
2667      int numberColumns = model->numberColumns();
2668      int numberRows = model->numberRows();
2669      int pivotRow = model->pivotRow();
2670      if (gubSlackIn_>=0)
2671        assert (sequenceIn>numberRows+numberColumns);
2672      if (sequenceIn==sequenceOut) 
2673        return -1;
2674      int iSetIn=-1;
2675      int iSetOut=-1;
2676      if (sequenceOut<numberColumns) {
2677        iSetOut = backward_[sequenceOut];
2678      } else if (sequenceOut>=numberRows+numberColumns) {
2679        assert (pivotRow>=numberRows);
2680        int iExtra = pivotRow-numberRows;
2681        assert (iExtra>=0);
2682        if (iSetOut<0)
2683          iSetOut = fromIndex_[iExtra];
2684        else
2685          assert(iSetOut == fromIndex_[iExtra]);
2686      }
2687      if (sequenceIn<numberColumns) {
2688        iSetIn = backward_[sequenceIn];
2689      } else if (gubSlackIn_>=0) {
2690        iSetIn = gubSlackIn_;
2691      }
2692      possiblePivotKey_=-1;
2693      number=0; // say do ordinary
2694      int * pivotVariable = model->pivotVariable();
2695      if (pivotRow>=numberRows) {
2696        int iExtra = pivotRow-numberRows;
2697        //const int * length = matrix_->getVectorLengths();
2698       
2699        assert (sequenceOut>=numberRows+numberColumns||
2700                sequenceOut==keyVariable_[iSetOut]);
2701        int incomingColumn = sequenceIn; // to be used in updates
2702        if (iSetIn!=iSetOut) {
2703          // We need to find a possible pivot for incoming
2704          // look through rowArray_[1]
2705          int n = model->rowArray(1)->getNumElements();
2706          int * which = model->rowArray(1)->getIndices();
2707          double * array = model->rowArray(1)->denseVector();
2708          double bestAlpha = 1.0e-5;
2709          //int shortest=numberRows+1;
2710          for (int i=0;i<n;i++) {
2711            int iRow = which[i];
2712            int iPivot = pivotVariable[iRow];
2713            if (iPivot<numberColumns&&backward_[iPivot]==iSetOut) {
2714              if (fabs(array[i])>fabs(bestAlpha)) {
2715                bestAlpha = array[i];
2716                possiblePivotKey_=iRow;
2717              }
2718            }
2719          }
2720          assert (possiblePivotKey_>=0); // could set returnCode=4
2721          number=1;
2722          if (sequenceIn>=numberRows+numberColumns) {
2723            number=3;
2724            // need swap as gub slack in and must become key
2725            // is this best way
2726            int key=keyVariable_[iSetIn];
2727            assert (key<numberColumns);
2728            // check other basic
2729            int iColumn = next_[key];
2730            // set new key to be used by unpack
2731            keyVariable_[iSetIn]=iSetIn+numberColumns;
2732            // change cost in changeCost
2733            {
2734              int iLook=0;
2735              int iSet=fromIndex_[0];
2736              while (iSet>=0) {
2737                if (iSet==iSetIn) {
2738                  changeCost_[iLook]=0.0;
2739                  break;
2740                }
2741                iSet = fromIndex_[++iLook];
2742              }
2743            }
2744            while (iColumn>=0) {
2745              if (iColumn!=sequenceOut) {
2746                // need partial ftran and skip accuracy check in replaceColumn
2747#ifdef CLP_DEBUG_PRINT
2748                printf("TTTTTry 5\n");
2749#endif
2750                int iRow = backToPivotRow_[iColumn];
2751                assert (iRow>=0);
2752                unpack(model,model->rowArray(3),iColumn);
2753                model->factorization()->updateColumnFT(model->rowArray(2),model->rowArray(3));
2754                double alpha = model->rowArray(3)->denseVector()[iRow];
2755                //if (!alpha)
2756                //printf("zero alpha a\n");
2757                int updateStatus = model->factorization()->replaceColumn(model,
2758                                                                         model->rowArray(2),
2759                                                                         model->rowArray(3),
2760                                                                         iRow, alpha);
2761                returnCode = CoinMax(updateStatus, returnCode);
2762                model->rowArray(3)->clear();
2763                if (returnCode)
2764                  break;
2765              }
2766              iColumn=next_[iColumn];
2767            }
2768            if (!returnCode) {
2769              // now factorization looks as if key is out
2770              // pivot back in
2771#ifdef CLP_DEBUG_PRINT
2772              printf("TTTTTry 6\n");
2773#endif
2774              unpack(model,model->rowArray(3),key);
2775              model->factorization()->updateColumnFT(model->rowArray(2),model->rowArray(3));
2776              pivotRow = possiblePivotKey_;
2777              double alpha = model->rowArray(3)->denseVector()[pivotRow];
2778              //if (!alpha)
2779              //printf("zero alpha b\n");
2780              int updateStatus = model->factorization()->replaceColumn(model,
2781                                                                       model->rowArray(2),
2782                                                                       model->rowArray(3),
2783                                                                       pivotRow, alpha);
2784              returnCode = CoinMax(updateStatus, returnCode);
2785              model->rowArray(3)->clear();
2786            }
2787            // restore key
2788            keyVariable_[iSetIn]=key;
2789            // now alternate column can replace key on out
2790            incomingColumn = pivotVariable[possiblePivotKey_];
2791          } else {
2792#ifdef CLP_DEBUG_PRINT
2793            printf("TTTTTTry 4 %d\n",possiblePivotKey_);
2794#endif
2795            int updateStatus = model->factorization()->replaceColumn(model,
2796                                                                     model->rowArray(2),
2797                                                                     model->rowArray(1),
2798                                                                     possiblePivotKey_, 
2799                                                                     bestAlpha);
2800            returnCode = CoinMax(updateStatus, returnCode);
2801            incomingColumn = pivotVariable[possiblePivotKey_];
2802          }
2803         
2804          //returnCode=4; // need swap
2805        } else {
2806          // key swap
2807          number=-1;
2808        }
2809        int key=keyVariable_[iSetOut];
2810        if (key<numberColumns)
2811          assert(key==sequenceOut);
2812        // check if any other basic
2813        int iColumn = next_[key];
2814        if (returnCode)
2815          iColumn = -1; // skip if error on previous
2816        // set new key to be used by unpack
2817        if (incomingColumn<numberColumns)
2818          keyVariable_[iSetOut]=incomingColumn;
2819        else
2820          keyVariable_[iSetOut]=iSetIn+numberColumns;
2821        double * cost = model->costRegion();
2822        if (possiblePivotKey_<0) {
2823          double dj = model->djRegion()[sequenceIn]-cost[sequenceIn];
2824          changeCost_[iExtra] = -dj;
2825#ifdef CLP_DEBUG_PRINT
2826          printf("modifying changeCost %d by %g - cost %g\n",iExtra, dj,cost[sequenceIn]);
2827#endif
2828        }
2829        while (iColumn>=0) {
2830          if (iColumn!=incomingColumn) {
2831            number=-2;
2832            // need partial ftran and skip accuracy check in replaceColumn
2833#ifdef CLP_DEBUG_PRINT
2834            printf("TTTTTTry 1\n");
2835#endif
2836            int iRow = backToPivotRow_[iColumn];
2837            assert (iRow>=0&&iRow<numberRows);
2838            unpack(model,model->rowArray(3),iColumn);
2839            model->factorization()->updateColumnFT(model->rowArray(2),model->rowArray(3));
2840            double * array = model->rowArray(3)->denseVector();
2841            double alpha = array[iRow];
2842            //if (!alpha)
2843            //printf("zero alpha d\n");
2844            int updateStatus = model->factorization()->replaceColumn(model,
2845                                                                     model->rowArray(2),
2846                                                                     model->rowArray(3),
2847                                                                     iRow, alpha);
2848            returnCode = CoinMax(updateStatus, returnCode);
2849            model->rowArray(3)->clear();
2850            if (returnCode)
2851              break;
2852          }
2853          iColumn=next_[iColumn];
2854        }
2855        // restore key
2856        keyVariable_[iSetOut]=key;
2857      } else if (sequenceIn>=numberRows+numberColumns) {
2858        number=2;
2859        //returnCode=4;
2860        // need swap as gub slack in and must become key
2861        // is this best way
2862        int key=keyVariable_[iSetIn];
2863        assert (key<numberColumns);
2864        // check other basic
2865        int iColumn = next_[key];
2866        // set new key to be used by unpack
2867        keyVariable_[iSetIn]=iSetIn+numberColumns;
2868        // change cost in changeCost
2869        {
2870          int iLook=0;
2871          int iSet=fromIndex_[0];
2872          while (iSet>=0) {
2873            if (iSet==iSetIn) {
2874              changeCost_[iLook]=0.0;
2875              break;
2876            }
2877            iSet = fromIndex_[++iLook];
2878          }
2879        }
2880        while (iColumn>=0) {
2881          if (iColumn!=sequenceOut) {
2882            // need partial ftran and skip accuracy check in replaceColumn
2883#ifdef CLP_DEBUG_PRINT
2884            printf("TTTTTry 2\n");
2885#endif
2886            int iRow = backToPivotRow_[iColumn];
2887            assert (iRow>=0);
2888            unpack(model,model->rowArray(3),iColumn);
2889            model->factorization()->updateColumnFT(model->rowArray(2),model->rowArray(3));
2890            double alpha = model->rowArray(3)->denseVector()[iRow];
2891            //if (!alpha)
2892            //printf("zero alpha e\n");
2893            int updateStatus = model->factorization()->replaceColumn(model,
2894                                                                     model->rowArray(2),
2895                                                                     model->rowArray(3),
2896                                                                     iRow, alpha);
2897            returnCode = CoinMax(updateStatus, returnCode);
2898            model->rowArray(3)->clear();
2899            if (returnCode)
2900              break;
2901          }
2902          iColumn=next_[iColumn];
2903        }
2904        if (!returnCode) {
2905          // now factorization looks as if key is out
2906          // pivot back in
2907#ifdef CLP_DEBUG_PRINT
2908          printf("TTTTTry 3\n");
2909#endif
2910          unpack(model,model->rowArray(3),key);
2911          model->factorization()->updateColumnFT(model->rowArray(2),model->rowArray(3));
2912          double alpha = model->rowArray(3)->denseVector()[pivotRow];
2913          //if (!alpha)
2914          //printf("zero alpha f\n");
2915          int updateStatus = model->factorization()->replaceColumn(model,
2916                                                                   model->rowArray(2),
2917                                                                   model->rowArray(3),
2918                                                                   pivotRow, alpha);
2919          returnCode = CoinMax(updateStatus, returnCode);
2920          model->rowArray(3)->clear();
2921        }
2922        // restore key
2923        keyVariable_[iSetIn]=key;
2924      } else {
2925        // normal - but might as well do here
2926        returnCode = model->factorization()->replaceColumn(model,
2927                                                           model->rowArray(2),
2928                                                           model->rowArray(1),
2929                                                           model->pivotRow(),
2930                                                           model->alpha());
2931      }
2932    }
2933#ifdef CLP_DEBUG_PRINT
2934    printf("Update type after %d - status %d - pivot row %d\n",
2935           number,returnCode,model->pivotRow());
2936#endif
2937    // see if column generation says time to re-factorize
2938    returnCode = CoinMax(returnCode,synchronize(model,5));
2939    number=-1; // say no need for normal replaceColumn
2940    break;
2941    // To see if can dual or primal
2942  case 4:
2943    {
2944      returnCode= 1;
2945    }
2946    break;
2947    // save status
2948  case 5:
2949    {
2950      synchronize(model,0);
2951      memcpy(saveStatus_,status_,numberSets_);
2952      memcpy(savedKeyVariable_,keyVariable_,numberSets_*sizeof(int));
2953    }
2954    break;
2955    // restore status
2956  case 6:
2957    {
2958      memcpy(status_,saveStatus_,numberSets_);
2959      memcpy(keyVariable_,savedKeyVariable_,numberSets_*sizeof(int));
2960      // restore firstAvailable_
2961      synchronize(model,7);
2962      // redo next_
2963      int i;
2964      int * last = new int[numberSets_];
2965      for (i=0;i<numberSets_;i++) {
2966        int iKey=keyVariable_[i];
2967        assert(iKey>=numberColumns||backward_[iKey]==i);
2968        last[i]=iKey;
2969        // make sure basic
2970        //if (iKey<numberColumns)
2971        //model->setStatus(iKey,ClpSimplex::basic);
2972      }
2973      for (i=0;i<numberColumns;i++){
2974        int iSet = backward_[i];
2975        if (iSet>=0) {
2976          next_[last[iSet]]=i;
2977          last[iSet]=i;
2978        }
2979      }
2980      for (i=0;i<numberSets_;i++) {
2981        next_[last[i]]=-(keyVariable_[i]+1);
2982        redoSet(model,keyVariable_[i],keyVariable_[i],i);
2983      }
2984      delete [] last;
2985      // redo pivotVariable
2986      int * pivotVariable = model->pivotVariable();
2987      int iRow;
2988      int numberBasic=0;
2989      int numberRows = model->numberRows();
2990      for (iRow=0;iRow<numberRows;iRow++) {
2991        if (model->getRowStatus(iRow)==ClpSimplex::basic) {
2992          numberBasic++;
2993          pivotVariable[iRow]=iRow+numberColumns;
2994        } else {
2995          pivotVariable[iRow]=-1;
2996        }
2997      }
2998      i=0;
2999      int iColumn;
3000      for (iColumn=0;iColumn<numberColumns;iColumn++) {
3001        if (model->getStatus(iColumn)==ClpSimplex::basic) {
3002          int iSet = backward_[iColumn];
3003          if (iSet<0||keyVariable_[iSet]!=iColumn) {
3004            while (pivotVariable[i]>=0) {
3005              i++;
3006              assert (i<numberRows);
3007            }
3008            pivotVariable[i]=iColumn;
3009            backToPivotRow_[iColumn]=i;
3010            numberBasic++;
3011          }
3012        }
3013      }
3014      assert (numberBasic==numberRows);
3015      rhsOffset(model,true);
3016    }
3017    break;
3018    // flag a variable
3019  case 7:
3020    {
3021      assert (number==model->sequenceIn());
3022      synchronize(model,1);
3023      synchronize(model,8);
3024    }
3025    break;
3026    // unflag all variables
3027  case 8:
3028    {
3029      returnCode=synchronize(model,2);
3030    }
3031    break;
3032    // redo costs in primal
3033  case 9:
3034    {
3035      returnCode=synchronize(model,3);
3036    }
3037    break;
3038    // return 1 if there may be changing bounds on variable (column generation)
3039  case 10:
3040    {
3041      returnCode=synchronize(model,6);
3042    }
3043    break;
3044    // make sure set is clean
3045  case 11:
3046    {
3047      assert (number==model->sequenceIn());
3048      returnCode=synchronize(model,8);
3049    }
3050    break;
3051  default:
3052    break;
3053  }
3054  return returnCode;
3055}
3056// Sets up an effective RHS and does gub crash if needed
3057void 
3058ClpGubMatrix::useEffectiveRhs(ClpSimplex * model, bool cheapest)
3059{
3060  // Do basis - cheapest or slack if feasible (unless cheapest set)
3061  int longestSet=0;
3062  int iSet;
3063  for (iSet=0;iSet<numberSets_;iSet++) 
3064    longestSet = CoinMax(longestSet,end_[iSet]-start_[iSet]);
3065   
3066  double * upper = new double[longestSet+1];
3067  double * cost = new double[longestSet+1];
3068  double * lower = new double[longestSet+1];
3069  double * solution = new double[longestSet+1];
3070  assert (!next_);
3071  int numberColumns = getNumCols();
3072  const int * columnLength = matrix_->getVectorLengths(); 
3073  const double * columnLower = model->lowerRegion();
3074  const double * columnUpper = model->upperRegion();
3075  double * columnSolution = model->solutionRegion();
3076  const double * objective = model->costRegion();
3077  int numberRows = getNumRows();
3078  toIndex_ = new int[numberSets_];
3079  for (iSet=0;iSet<numberSets_;iSet++) 
3080    toIndex_[iSet]=-1;
3081  fromIndex_ = new int [getNumRows()+1];
3082  double tolerance = model->primalTolerance();
3083  bool noNormalBounds=true;
3084  gubType_ &= ~8;
3085  bool gotBasis=false;
3086  for (iSet=0;iSet<numberSets_;iSet++) {
3087    if (keyVariable_[iSet]<numberColumns)
3088      gotBasis=true;
3089    CoinBigIndex j;
3090    CoinBigIndex iStart = start_[iSet];
3091    CoinBigIndex iEnd=end_[iSet];
3092    for (j=iStart;j<iEnd;j++) {
3093      if (columnLower[j]&&columnLower[j]>-1.0e20)
3094        noNormalBounds=false;
3095      if (columnUpper[j]&&columnUpper[j]<1.0e20)
3096        noNormalBounds=false;
3097    }
3098  }
3099  if (noNormalBounds)
3100    gubType_ |= 8;
3101  if (!gotBasis) {
3102    for (iSet=0;iSet<numberSets_;iSet++) {
3103      CoinBigIndex j;
3104      int numberBasic=0;
3105      int iBasic=-1;
3106      CoinBigIndex iStart = start_[iSet];
3107      CoinBigIndex iEnd=end_[iSet];
3108      // find one with smallest length
3109      int smallest=numberRows+1;
3110      double value=0.0;
3111      for (j=iStart;j<iEnd;j++) {
3112        if (model->getStatus(j)== ClpSimplex::basic) {
3113          if (columnLength[j]<smallest) {
3114            smallest=columnLength[j];
3115            iBasic=j;
3116          }
3117          numberBasic++;
3118        }
3119        value += columnSolution[j];
3120      }
3121      bool done=false;
3122      if (numberBasic>1||(numberBasic==1&&getStatus(iSet)==ClpSimplex::basic)) {
3123        if (getStatus(iSet)==ClpSimplex::basic) 
3124          iBasic = iSet+numberColumns;// slack key - use
3125        done=true;
3126      } else if (numberBasic==1) {
3127        // see if can be key
3128        double thisSolution = columnSolution[iBasic];
3129        if (thisSolution>columnUpper[iBasic]) {
3130          value -= thisSolution-columnUpper[iBasic];
3131          thisSolution = columnUpper[iBasic];
3132          columnSolution[iBasic]=thisSolution;
3133        }
3134        if (thisSolution<columnLower[iBasic]) {
3135          value -= thisSolution-columnLower[iBasic];
3136          thisSolution = columnLower[iBasic];
3137          columnSolution[iBasic]=thisSolution;
3138        }
3139        // try setting slack to a bound
3140        assert (upper_[iSet]<1.0e20||lower_[iSet]>-1.0e20);
3141        double cost1 = COIN_DBL_MAX;
3142        int whichBound=-1;
3143        if (upper_[iSet]<1.0e20) {
3144          // try slack at ub
3145          double newBasic = thisSolution +upper_[iSet]-value;
3146          if (newBasic>=columnLower[iBasic]-tolerance&&
3147              newBasic<=columnUpper[iBasic]+tolerance) {
3148            // can go
3149            whichBound=1;
3150            cost1 = newBasic*objective[iBasic];
3151            // But if exact then may be good solution
3152            if (fabs(upper_[iSet]-value)<tolerance)
3153              cost1=-COIN_DBL_MAX;
3154          }
3155        }
3156        if (lower_[iSet]>-1.0e20) {
3157          // try slack at lb
3158          double newBasic = thisSolution +lower_[iSet]-value;
3159          if (newBasic>=columnLower[iBasic]-tolerance&&
3160              newBasic<=columnUpper[iBasic]+tolerance) {
3161            // can go but is it cheaper
3162            double cost2 = newBasic*objective[iBasic];
3163            // But if exact then may be good solution
3164            if (fabs(lower_[iSet]-value)<tolerance)
3165              cost2=-COIN_DBL_MAX;
3166            if (cost2<cost1)
3167              whichBound=0;
3168          }
3169        }
3170        if (whichBound!=-1) {
3171          // key
3172          done=true;
3173          if (whichBound) {
3174            // slack to upper
3175            columnSolution[iBasic]=thisSolution + upper_[iSet]-value;
3176            setStatus(iSet,ClpSimplex::atUpperBound);
3177          } else {
3178            // slack to lower
3179            columnSolution[iBasic]=thisSolution + lower_[iSet]-value;
3180            setStatus(iSet,ClpSimplex::atLowerBound);
3181          }
3182        }
3183      }
3184      if (!done) {
3185        if (!cheapest) {
3186          // see if slack can be key
3187          if (value>=lower_[iSet]-tolerance&&value<=upper_[iSet]+tolerance) {
3188            done=true;
3189            setStatus(iSet,ClpSimplex::basic);
3190            iBasic=iSet+numberColumns;
3191          }
3192        }
3193        if (!done) {
3194          // set non basic if there was one
3195          if (iBasic>=0)
3196            model->setStatus(iBasic,ClpSimplex::atLowerBound);
3197          // find cheapest
3198          int numberInSet = iEnd-iStart;
3199          CoinMemcpyN(columnLower+iStart,numberInSet,lower);
3200          CoinMemcpyN(columnUpper+iStart,numberInSet,upper);
3201          CoinMemcpyN(columnSolution+iStart,numberInSet,solution);
3202          // and slack
3203          iBasic=numberInSet;
3204          solution[iBasic]=-value;
3205          lower[iBasic]=-upper_[iSet];
3206          upper[iBasic]=-lower_[iSet];
3207          int kphase;
3208          if (value>=lower_[iSet]-tolerance&&value<=upper_[iSet]+tolerance) {
3209            // feasible
3210            kphase=1;
3211            cost[iBasic]=0.0;
3212            CoinMemcpyN(objective+iStart,numberInSet,cost);
3213          } else {
3214            // infeasible
3215            kphase=0;
3216            // remember bounds are flipped so opposite to natural
3217            if (value<lower_[iSet]-tolerance)
3218              cost[iBasic]=1.0;
3219            else
3220              cost[iBasic]=-1.0;
3221            CoinZeroN(cost,numberInSet);
3222          }
3223          double dualTolerance =model->dualTolerance();
3224          for (int iphase =kphase;iphase<2;iphase++) {
3225            if (iphase) {
3226              cost[numberInSet]=0.0;
3227              CoinMemcpyN(objective+iStart,numberInSet,cost);
3228            }
3229            // now do one row lp
3230            bool improve=true;
3231            while (improve) {
3232              improve=false;
3233              double dual = cost[iBasic];
3234              int chosen =-1;
3235              double best=dualTolerance;
3236              int way=0;
3237              for (int i=0;i<=numberInSet;i++) {
3238                double dj = cost[i]-dual;
3239                double improvement =0.0;
3240                double distance=0.0;
3241                if (iphase||i<numberInSet)
3242                  assert (solution[i]>=lower[i]&&solution[i]<=upper[i]);
3243                if (dj>dualTolerance)
3244                  improvement = dj*(solution[i]-lower[i]);
3245                else if (dj<-dualTolerance)
3246                  improvement = dj*(solution[i]-upper[i]);
3247                if (improvement>best) {
3248                  best=improvement;
3249                  chosen=i;
3250                  if (dj<0.0) {
3251                    way = 1;
3252                    distance = upper[i]-solution[i];
3253                  } else {
3254                    way = -1;
3255                    distance = solution[i]-lower[i];
3256                  }
3257                }
3258              }
3259              if (chosen>=0) {
3260                improve=true;
3261                // now see how far
3262                if (way>0) {
3263                  // incoming increasing so basic decreasing
3264                  // if phase 0 then go to nearest bound
3265                  double distance=upper[chosen]-solution[chosen];
3266                  double basicDistance;
3267                  if (!iphase) {
3268                    assert (iBasic==numberInSet);
3269                    assert (solution[iBasic]>upper[iBasic]);
3270                    basicDistance = solution[iBasic]-upper[iBasic];
3271                  } else {
3272                    basicDistance = solution[iBasic]-lower[iBasic];
3273                  }
3274                  // need extra coding for unbounded
3275                  assert (CoinMin(distance,basicDistance)<1.0e20);
3276                  if (distance>basicDistance) {
3277                    // incoming becomes basic
3278                    solution[chosen] += basicDistance;
3279                    if (!iphase) 
3280                      solution[iBasic]=upper[iBasic];
3281                    else 
3282                      solution[iBasic]=lower[iBasic];
3283                    iBasic = chosen;
3284                  } else {
3285                    // flip
3286                    solution[chosen]=upper[chosen];
3287                    solution[iBasic] -= distance;
3288                  }
3289                } else {
3290                  // incoming decreasing so basic increasing
3291                  // if phase 0 then go to nearest bound
3292                  double distance=solution[chosen]-lower[chosen];
3293                  double basicDistance;
3294                  if (!iphase) {
3295                    assert (iBasic==numberInSet);
3296                    assert (solution[iBasic]<lower[iBasic]);
3297                    basicDistance = lower[iBasic]-solution[iBasic];
3298                  } else {
3299                    basicDistance = upper[iBasic]-solution[iBasic];
3300                  }
3301                  // need extra coding for unbounded - for now just exit
3302                  if (CoinMin(distance,basicDistance)>1.0e20) {
3303                    printf("unbounded on set %d\n",iSet);
3304                    iphase=1;
3305                    iBasic=numberInSet;
3306                    break;
3307                  }
3308                  if (distance>basicDistance) {
3309                    // incoming becomes basic
3310                    solution[chosen] -= basicDistance;
3311                    if (!iphase) 
3312                      solution[iBasic]=lower[iBasic];
3313                    else 
3314                      solution[iBasic]=upper[iBasic];
3315                    iBasic = chosen;
3316                  } else {
3317                    // flip
3318                    solution[chosen]=lower[chosen];
3319                    solution[iBasic] += distance;
3320                  }
3321                }
3322                if (!iphase) {
3323                  if(iBasic<numberInSet)
3324                    break; // feasible
3325                  else if (solution[iBasic]>=lower[iBasic]&&
3326                           solution[iBasic]<=upper[iBasic])
3327                    break; // feasible (on flip)
3328                }
3329              }
3330            }
3331          }
3332          // convert iBasic back and do bounds
3333          if (iBasic==numberInSet) {
3334            // slack basic
3335            setStatus(iSet,ClpSimplex::basic);
3336            iBasic=iSet+numberColumns;
3337          } else {
3338            iBasic += start_[iSet];
3339            model->setStatus(iBasic,ClpSimplex::basic);
3340            // remember bounds flipped
3341            if (upper[numberInSet]==lower[numberInSet]) 
3342              setStatus(iSet,ClpSimplex::isFixed);
3343            else if (solution[numberInSet]==upper[numberInSet])
3344              setStatus(iSet,ClpSimplex::atLowerBound);
3345            else if (solution[numberInSet]==lower[numberInSet])
3346              setStatus(iSet,ClpSimplex::atUpperBound);
3347            else 
3348              abort();
3349          }
3350          for (j=iStart;j<iEnd;j++) {
3351            if (model->getStatus(j)!=ClpSimplex::basic) {
3352              int inSet=j-iStart;
3353              columnSolution[j]=solution[inSet];
3354              if (upper[inSet]==lower[inSet]) 
3355                model->setStatus(j,ClpSimplex::isFixed);
3356              else if (solution[inSet]==upper[inSet])
3357                model->setStatus(j,ClpSimplex::atUpperBound);
3358              else if (solution[inSet]==lower[inSet])
3359                model->setStatus(j,ClpSimplex::atLowerBound);
3360            }
3361          }
3362        }
3363      } 
3364      keyVariable_[iSet]=iBasic;
3365    }
3366  }
3367  delete [] lower;
3368  delete [] solution;
3369  delete [] upper;
3370  delete [] cost;
3371  // make sure matrix is in good shape
3372  matrix_->orderMatrix();
3373  // create effective rhs
3374  delete [] rhsOffset_;
3375  rhsOffset_ = new double[numberRows];
3376  delete [] next_;
3377  next_ = new int[numberColumns+numberSets_+2*longestSet];
3378  char * mark = new char[numberColumns];
3379  memset(mark,0,numberColumns);
3380  for (int iColumn=0;iColumn<numberColumns;iColumn++) 
3381    next_[iColumn]=INT_MAX;
3382  int i;
3383  int * keys = new int[numberSets_];
3384  for (i=0;i<numberSets_;i++) 
3385    keys[i]=INT_MAX;
3386  // set up chains
3387  for (i=0;i<numberColumns;i++){
3388    if (model->getStatus(i)==ClpSimplex::basic) 
3389      mark[i]=1;
3390    int iSet = backward_[i];
3391    if (iSet>=0) {
3392      int iNext = keys[iSet];
3393      next_[i]=iNext;
3394      keys[iSet]=i;
3395    }
3396  }
3397  for (i=0;i<numberSets_;i++) {
3398    int j;
3399    if (getStatus(i)!=ClpSimplex::basic) {
3400      // make sure fixed if it is
3401      if (upper_[i]==lower_[i])
3402        setStatus(i,ClpSimplex::isFixed);
3403      // slack not key - choose one with smallest length
3404      int smallest=numberRows+1;
3405      int key=-1;
3406      j = keys[i];
3407      if (j!=INT_MAX) {
3408        while (1) {
3409          if (mark[j]&&columnLength[j]<smallest&&!gotBasis) {
3410            key=j;
3411            smallest=columnLength[j];
3412          }
3413          if (next_[j]!=INT_MAX) {
3414            j = next_[j];
3415          } else {
3416            // correct end
3417            next_[j]=-(keys[i]+1);
3418            break;
3419          }
3420        }
3421      } else {
3422        next_[i+numberColumns] = -(numberColumns+i+1);
3423      } 
3424      if (gotBasis)
3425        key =keyVariable_[i];
3426      if (key>=0) {
3427        keyVariable_[i]=key;
3428      } else {
3429        // nothing basic - make slack key
3430        //((ClpGubMatrix *)this)->setStatus(i,ClpSimplex::basic);
3431        // fudge to avoid const problem
3432        status_[i]=1;
3433      }
3434    } else {
3435      // slack key
3436      keyVariable_[i]=numberColumns+i;
3437      int j;
3438      double sol=0.0;
3439      j = keys[i];
3440      if (j!=INT_MAX) {
3441        while (1) {
3442          sol += columnSolution[j];
3443          if (next_[j]!=INT_MAX) {
3444            j = next_[j];
3445          } else {
3446            // correct end
3447            next_[j]=-(keys[i]+1);
3448            break;
3449          }
3450        }
3451      } else {
3452        next_[i+numberColumns] = -(numberColumns+i+1);
3453      }
3454      if (sol>upper_[i]+tolerance) {
3455        setAbove(i);
3456      } else if (sol<lower_[i]-tolerance) {
3457        setBelow(i);
3458      } else {
3459        setFeasible(i);
3460      }
3461    }
3462    // Create next_
3463    int key = keyVariable_[i];
3464    redoSet(model,key,keys[i],i);
3465  }
3466  delete [] keys;
3467  delete [] mark;
3468  rhsOffset(model,true);
3469}
3470// redoes next_ for a set. 
3471void 
3472ClpGubMatrix::redoSet(ClpSimplex * model, int newKey, int oldKey, int iSet)
3473{
3474  int numberColumns = model->numberColumns();
3475  int * save = next_+numberColumns+numberSets_;
3476  int number=0;
3477  int stop = -(oldKey+1);
3478  int j=next_[oldKey];
3479  while (j!=stop) {
3480    if (j<0)
3481      j = -j-1;
3482    if (j!=newKey)
3483      save[number++]=j;
3484    j=next_[j];
3485  }
3486  // and add oldkey
3487  if (newKey!=oldKey)
3488    save[number++]=oldKey;
3489  // now do basic
3490  int lastMarker = -(newKey+1);
3491  keyVariable_[iSet]=newKey;
3492  next_[newKey]=lastMarker;
3493  int last = newKey;
3494  for ( j=0;j<number;j++) {
3495    int iColumn=save[j];
3496    if (iColumn<numberColumns) {
3497      if (model->getStatus(iColumn)==ClpSimplex::basic) {
3498        next_[last]=iColumn;
3499        next_[iColumn]=lastMarker;
3500        last = iColumn;
3501      }
3502    }
3503  }
3504  // now add in non-basic
3505  for ( j=0;j<number;j++) {
3506    int iColumn=save[j];
3507    if (iColumn<numberColumns) {
3508      if (model->getStatus(iColumn)!=ClpSimplex::basic) {
3509        next_[last]=-(iColumn+1);
3510        next_[iColumn]=lastMarker;
3511        last = iColumn;
3512      }
3513    }
3514  }
3515
3516}
3517/* Returns effective RHS if it is being used.  This is used for long problems
3518   or big gub or anywhere where going through full columns is
3519   expensive.  This may re-compute */
3520double * 
3521ClpGubMatrix::rhsOffset(ClpSimplex * model,bool forceRefresh,bool check)
3522{
3523  //forceRefresh=true;
3524  if (rhsOffset_) {
3525#ifdef CLP_DEBUG
3526    if (check) {
3527      // no need - but check anyway
3528      // zero out basic
3529      int numberRows = model->numberRows();
3530      int numberColumns = model->numberColumns();
3531      double * solution = new double [numberColumns];
3532      double * rhs = new double[numberRows];
3533      CoinMemcpyN(model->solutionRegion(),numberColumns,solution);
3534      CoinZeroN(rhs,numberRows);
3535      int iRow;
3536      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
3537        if (model->getColumnStatus(iColumn)==ClpSimplex::basic)
3538          solution[iColumn]=0.0;
3539      }
3540      for (int iSet=0;iSet<numberSets_;iSet++) {
3541        int iColumn = keyVariable_[iSet];
3542        if (iColumn<numberColumns) 
3543          solution[iColumn]=0.0;
3544      }
3545      times(-1.0,solution,rhs);
3546      delete [] solution;
3547      const double * columnSolution = model->solutionRegion();
3548      // and now subtract out non basic
3549      ClpSimplex::Status iStatus;
3550      for (int iSet=0;iSet<numberSets_;iSet++) {
3551        int iColumn = keyVariable_[iSet];
3552        if (iColumn<numberColumns) {
3553          double b=0.0;
3554          // key is structural - where is slack
3555          iStatus = getStatus(iSet);
3556          assert (iStatus!=ClpSimplex::basic);
3557          if (iStatus==ClpSimplex::atLowerBound)
3558            b=lower_[iSet];
3559          else
3560            b=upper_[iSet];
3561          // subtract out others at bounds
3562          if ((gubType_&8)==0) {
3563            int stop = -(iColumn+1);
3564            int jColumn =next_[iColumn];
3565            // sum all non-basic variables - first skip basic
3566            while(jColumn>=0) 
3567              jColumn=next_[jColumn];
3568            while(jColumn!=stop) {
3569              assert (jColumn<0);
3570              jColumn = -jColumn-1;
3571              b -= columnSolution[jColumn];
3572              jColumn=next_[jColumn];
3573            }
3574          }
3575          // subtract out
3576          ClpPackedMatrix::add(model,rhs,iColumn,-b);
3577        }
3578      }
3579      for (iRow=0;iRow<numberRows;iRow++) {
3580        if (fabs(rhs[iRow]-rhsOffset_[iRow])>1.0e-3)
3581          printf("** bad effective %d - true %g old %g\n",iRow,rhs[iRow],rhsOffset_[iRow]);
3582      }
3583      delete [] rhs;
3584    }
3585#endif
3586    if (forceRefresh||(refreshFrequency_&&model->numberIterations()>=
3587                       lastRefresh_+refreshFrequency_)) {
3588      // zero out basic
3589      int numberRows = model->numberRows();
3590      int numberColumns = model->numberColumns();
3591      double * solution = new double [numberColumns];
3592      CoinMemcpyN(model->solutionRegion(),numberColumns,solution);
3593      CoinZeroN(rhsOffset_,numberRows);
3594      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
3595        if (model->getColumnStatus(iColumn)==ClpSimplex::basic)
3596          solution[iColumn]=0.0;
3597      }
3598      for (int iSet=0;iSet<numberSets_;iSet++) {
3599        int iColumn = keyVariable_[iSet];
3600        if (iColumn<numberColumns) 
3601          solution[iColumn]=0.0;
3602      }
3603      times(-1.0,solution,rhsOffset_);
3604      delete [] solution;
3605      lastRefresh_ = model->numberIterations();
3606      const double * columnSolution = model->solutionRegion();
3607      // and now subtract out non basic
3608      ClpSimplex::Status iStatus;
3609      for (int iSet=0;iSet<numberSets_;iSet++) {
3610        int iColumn = keyVariable_[iSet];
3611        if (iColumn<numberColumns) {
3612          double b=0.0;
3613          // key is structural - where is slack
3614          iStatus = getStatus(iSet);
3615          assert (iStatus!=ClpSimplex::basic);
3616          if (iStatus==ClpSimplex::atLowerBound)
3617            b=lower_[iSet];
3618          else
3619            b=upper_[iSet];
3620          // subtract out others at bounds
3621          if ((gubType_&8)==0) {
3622            int stop = -(iColumn+1);
3623            int jColumn =next_[iColumn];
3624            // sum all non-basic variables - first skip basic
3625            while(jColumn>=0) 
3626              jColumn=next_[jColumn];
3627            while(jColumn!=stop) {
3628              assert (jColumn<0);
3629              jColumn = -jColumn-1;
3630              b -= columnSolution[jColumn];
3631              jColumn=next_[jColumn];
3632            }
3633          }
3634          // subtract out
3635          if (b)
3636            ClpPackedMatrix::add(model,rhsOffset_,iColumn,-b);
3637        }
3638      }
3639    }
3640  }
3641  return rhsOffset_;
3642}
3643/*
3644   update information for a pivot (and effective rhs)
3645*/
3646int 
3647ClpGubMatrix::updatePivot(ClpSimplex * model,double oldInValue, double oldOutValue)
3648{
3649  int sequenceIn = model->sequenceIn();
3650  int sequenceOut = model->sequenceOut();
3651  double * solution = model->solutionRegion();
3652  int numberColumns = model->numberColumns();
3653  int numberRows = model->numberRows();
3654  int pivotRow = model->pivotRow();
3655  int iSetIn;
3656  // Correct sequence in
3657  trueSequenceIn_=sequenceIn;
3658  if (sequenceIn<numberColumns) {
3659    iSetIn = backward_[sequenceIn];
3660  } else if (sequenceIn<numberColumns+numberRows) {
3661    iSetIn = -1;
3662  } else {
3663    iSetIn = gubSlackIn_;
3664    trueSequenceIn_ = numberColumns+numberRows+iSetIn;
3665  }
3666  int iSetOut=-1;
3667  trueSequenceOut_=sequenceOut;
3668  if (sequenceOut<numberColumns) {
3669    iSetOut = backward_[sequenceOut];
3670  } else if (sequenceOut>=numberRows+numberColumns) {
3671    assert (pivotRow>=numberRows);
3672    int iExtra = pivotRow-numberRows;
3673    assert (iExtra>=0);
3674    if (iSetOut<0)
3675      iSetOut = fromIndex_[iExtra];
3676    else
3677      assert(iSetOut == fromIndex_[iExtra]);
3678    trueSequenceOut_ = numberColumns+numberRows+iSetOut;
3679  }
3680  if (rhsOffset_) {
3681    // update effective rhs
3682    if (sequenceIn==sequenceOut) {
3683      assert (sequenceIn<numberRows+numberColumns); // should be easy to deal with
3684      if (sequenceIn<numberColumns)
3685        add(model,rhsOffset_,sequenceIn,oldInValue-solution[sequenceIn]);
3686    } else {
3687      if (sequenceIn<numberColumns) {
3688        // we need to test if WILL be key
3689        ClpPackedMatrix::add(model,rhsOffset_,sequenceIn,oldInValue);
3690        if (iSetIn>=0)  {
3691          // old contribution to rhsOffset_
3692          int key = keyVariable_[iSetIn];
3693          if (key<numberColumns) {
3694            double oldB=0.0;
3695            ClpSimplex::Status iStatus = getStatus(iSetIn);
3696            if (iStatus==ClpSimplex::atLowerBound)
3697              oldB=lower_[iSetIn];
3698            else
3699              oldB=upper_[iSetIn];
3700            // subtract out others at bounds
3701            if ((gubType_&8)==0) {
3702              int stop = -(key+1);
3703              int iColumn =next_[key];
3704              // skip basic
3705              while (iColumn>=0)
3706                iColumn = next_[iColumn];
3707              // sum all non-key variables
3708              while(iColumn!=stop) {
3709                assert (iColumn<0);
3710                iColumn = -iColumn-1;
3711                if (iColumn == sequenceIn) 
3712                  oldB -= oldInValue;
3713                else if ( iColumn != sequenceOut )
3714                  oldB -= solution[iColumn];
3715                iColumn=next_[iColumn];
3716              }
3717            }
3718            if (oldB)
3719              ClpPackedMatrix::add(model,rhsOffset_,key,oldB);
3720          }
3721        }
3722      } else if (sequenceIn<numberRows+numberColumns) {
3723        //rhsOffset_[sequenceIn-numberColumns] -= oldInValue;
3724      } else {
3725#ifdef CLP_DEBUG_PRINT
3726        printf("** in is key slack %d\n",sequenceIn);
3727#endif
3728        // old contribution to rhsOffset_
3729        int key = keyVariable_[iSetIn];
3730        if (key<numberColumns) {
3731          double oldB=0.0;
3732          ClpSimplex::Status iStatus = getStatus(iSetIn);
3733          if (iStatus==ClpSimplex::atLowerBound)
3734            oldB=lower_[iSetIn];
3735          else
3736            oldB=upper_[iSetIn];
3737          // subtract out others at bounds
3738          if ((gubType_&8)==0) {
3739            int stop = -(key+1);
3740            int iColumn =next_[key];
3741            // skip basic
3742            while (iColumn>=0)
3743              iColumn = next_[iColumn];
3744            // sum all non-key variables
3745            while(iColumn!=stop) {
3746              assert (iColumn<0);
3747              iColumn = -iColumn-1;
3748              if ( iColumn != sequenceOut )
3749                oldB -= solution[iColumn];
3750              iColumn=next_[iColumn];
3751            }
3752          }
3753          if (oldB)
3754            ClpPackedMatrix::add(model,rhsOffset_,key,oldB);
3755        }
3756      }
3757      if (sequenceOut<numberColumns) {
3758        ClpPackedMatrix::add(model,rhsOffset_,sequenceOut,-solution[sequenceOut]);
3759        if (iSetOut>=0) {
3760          // old contribution to rhsOffset_
3761          int key = keyVariable_[iSetOut];
3762          if (key<numberColumns&&iSetIn!=iSetOut) {
3763            double oldB=0.0;
3764            ClpSimplex::Status iStatus = getStatus(iSetOut);
3765            if (iStatus==ClpSimplex::atLowerBound)
3766              oldB=lower_[iSetOut];
3767            else
3768              oldB=upper_[iSetOut];
3769            // subtract out others at bounds
3770            if ((gubType_&8)==0) {
3771              int stop = -(key+1);
3772              int iColumn =next_[key];
3773              // skip basic
3774              while (iColumn>=0)
3775                iColumn = next_[iColumn];
3776              // sum all non-key variables
3777              while(iColumn!=stop) {
3778                assert (iColumn<0);
3779                iColumn = -iColumn-1;
3780                if (iColumn == sequenceIn) 
3781                  oldB -= oldInValue;
3782                else if ( iColumn != sequenceOut )
3783                  oldB -= solution[iColumn];
3784                iColumn=next_[iColumn];
3785              }
3786            }
3787            if (oldB)
3788              ClpPackedMatrix::add(model,rhsOffset_,key,oldB);
3789          }
3790        }
3791      } else if (sequenceOut<numberRows+numberColumns) {
3792        //rhsOffset_[sequenceOut-numberColumns] -= -solution[sequenceOut];
3793      } else {
3794#ifdef CLP_DEBUG_PRINT
3795        printf("** out is key slack %d\n",sequenceOut);
3796#endif
3797        assert (pivotRow>=numberRows);
3798      }
3799    }
3800  }
3801  int * pivotVariable = model->pivotVariable();
3802  // may need to deal with key
3803  // Also need coding to mark/allow key slack entering
3804  if (pivotRow>=numberRows) {
3805    assert (sequenceOut>=numberRows+numberColumns||sequenceOut==keyVariable_[iSetOut]);
3806#ifdef CLP_DEBUG_PRINT
3807    if (sequenceIn>=numberRows+numberColumns)
3808      printf("key slack %d in, set out %d\n",gubSlackIn_,iSetOut);
3809    printf("** danger - key out for set %d in %d (set %d)\n",iSetOut,sequenceIn,
3810           iSetIn);
3811#endif   
3812    // if slack out mark correctly
3813    if (sequenceOut>=numberRows+numberColumns) {
3814      double value=model->valueOut();
3815      if (value==upper_[iSetOut]) {
3816        setStatus(iSetOut,ClpSimplex::atUpperBound);
3817      } else if (value==lower_[iSetOut]) {
3818        setStatus(iSetOut,ClpSimplex::atLowerBound);
3819      } else {
3820        if (fabs(value-upper_[iSetOut])<
3821            fabs(value-lower_[iSetOut])) {
3822          setStatus(iSetOut,ClpSimplex::atUpperBound);
3823        } else {
3824          setStatus(iSetOut,ClpSimplex::atLowerBound);
3825        }
3826      }
3827      if (upper_[iSetOut]==lower_[iSetOut])
3828        setStatus(iSetOut,ClpSimplex::isFixed);
3829      setFeasible(iSetOut);
3830    }
3831    if (iSetOut==iSetIn) {
3832      // key swap
3833      int key;
3834      if (sequenceIn>=numberRows+numberColumns) {
3835        key = numberColumns+iSetIn;
3836        setStatus(iSetIn,ClpSimplex::basic);
3837      } else {
3838        key = sequenceIn;
3839      }
3840      redoSet(model,key,keyVariable_[iSetIn],iSetIn);
3841    } else {
3842      // key was chosen
3843      assert (possiblePivotKey_>=0&&possiblePivotKey_<numberRows);
3844      int key=pivotVariable[possiblePivotKey_];
3845      // and set incoming here
3846      if (sequenceIn>=numberRows+numberColumns) {
3847        // slack in - so use old key
3848        sequenceIn = keyVariable_[iSetIn];
3849        model->setStatus(sequenceIn,ClpSimplex::basic);
3850        setStatus(iSetIn,ClpSimplex::basic);
3851        redoSet(model,iSetIn+numberColumns,keyVariable_[iSetIn],iSetIn);
3852      }
3853      //? do not do if iSetIn<0 ? as will be done later
3854      pivotVariable[possiblePivotKey_]=sequenceIn;
3855      if (sequenceIn<numberColumns)
3856        backToPivotRow_[sequenceIn]=possiblePivotKey_;
3857      redoSet(model,key,keyVariable_[iSetOut],iSetOut);
3858    }
3859  } else {
3860    if (sequenceOut<numberColumns) {
3861      if (iSetIn>=0&&iSetOut==iSetIn) {
3862        // key not out - only problem is if slack in
3863        int key;
3864        if (sequenceIn>=numberRows+numberColumns) {
3865          key = numberColumns+iSetIn;
3866          setStatus(iSetIn,ClpSimplex::basic);
3867          assert (pivotRow<numberRows);
3868          // must swap with current key
3869          int key=keyVariable_[iSetIn];
3870          model->setStatus(key,ClpSimplex::basic);
3871          pivotVariable[pivotRow]=key;
3872          backToPivotRow_[key]=pivotRow;
3873        } else {
3874          key = keyVariable_[iSetIn];
3875        }
3876        redoSet(model,key,keyVariable_[iSetIn],iSetIn);
3877      } else if (iSetOut>=0) {
3878        // just redo set
3879        int key=keyVariable_[iSetOut];;
3880        redoSet(model,key,keyVariable_[iSetOut],iSetOut);
3881      }
3882    }
3883  }
3884  if (iSetIn>=0&&iSetIn!=iSetOut) {
3885    int key=keyVariable_[iSetIn];
3886    if (sequenceIn == numberColumns+2*numberRows) {
3887      // key slack in
3888      assert (pivotRow<numberRows);
3889      // must swap with current key
3890      model->setStatus(key,ClpSimplex::basic);
3891      pivotVariable[pivotRow]=key;
3892      backToPivotRow_[key]=pivotRow;
3893      setStatus(iSetIn,ClpSimplex::basic);
3894      key = iSetIn+numberColumns;
3895    }
3896    // redo set to allow for new one
3897    redoSet(model,key,keyVariable_[iSetIn],iSetIn);
3898  }
3899  // update pivot
3900  if (sequenceIn<numberColumns) {
3901    if (pivotRow<numberRows) {
3902      backToPivotRow_[sequenceIn]=pivotRow;
3903    } else {
3904      if (possiblePivotKey_>=0) {
3905        assert (possiblePivotKey_<numberRows);
3906        backToPivotRow_[sequenceIn]=possiblePivotKey_;
3907        pivotVariable[possiblePivotKey_]=sequenceIn;
3908      }
3909    }
3910  } else if (sequenceIn>=numberRows+numberColumns) {
3911    // key in - something should have been done before
3912    int key = keyVariable_[iSetIn];
3913    assert (key==numberColumns+iSetIn);
3914    //pivotVariable[pivotRow]=key;
3915    //backToPivotRow_[key]=pivotRow;
3916    //model->setStatus(key,ClpSimplex::basic);
3917    //key=numberColumns+iSetIn;
3918    setStatus(iSetIn,ClpSimplex::basic);
3919    redoSet(model,key,keyVariable_[iSetIn],iSetIn);
3920  }
3921#ifdef CLP_DEBUG
3922  {
3923    char * xx = new char[numberColumns+numberRows];
3924    memset(xx,0,numberRows+numberColumns);
3925    for (int i=0;i<numberRows;i++) {
3926      int iPivot = pivotVariable[i];
3927      assert (iPivot<numberRows+numberColumns);
3928      assert (!xx[iPivot]);
3929      xx[iPivot]=1;
3930      if (iPivot<numberColumns) {
3931        int iBack = backToPivotRow_[iPivot];
3932        assert (i==iBack);
3933      }
3934    }
3935    delete [] xx;
3936  }
3937#endif
3938  if (rhsOffset_) {
3939    // update effective rhs
3940    if (sequenceIn!=sequenceOut) {
3941      if (sequenceIn<numberColumns) {
3942        if (iSetIn>=0) {
3943          // new contribution to rhsOffset_
3944          int key = keyVariable_[iSetIn];
3945          if (key<numberColumns) {
3946            double newB=0.0;
3947            ClpSimplex::Status iStatus = getStatus(iSetIn);
3948            if (iStatus==ClpSimplex::atLowerBound)
3949              newB=lower_[iSetIn];
3950            else
3951              newB=upper_[iSetIn];
3952            // subtract out others at bounds
3953            if ((gubType_&8)==0) {
3954              int stop = -(key+1);
3955              int iColumn =next_[key];
3956              // skip basic
3957              while (iColumn>=0)
3958                iColumn = next_[iColumn];
3959              // sum all non-key variables
3960              while(iColumn!=stop) {
3961                assert (iColumn<0);
3962                iColumn = -iColumn-1;
3963                newB -= solution[iColumn];
3964                iColumn=next_[iColumn];
3965              }
3966            }
3967            if (newB)
3968              ClpPackedMatrix::add(model,rhsOffset_,key,-newB);
3969          }
3970        }
3971      }
3972      if (iSetOut>=0) {
3973        // new contribution to rhsOffset_
3974        int key = keyVariable_[iSetOut];
3975        if (key<numberColumns&&iSetIn!=iSetOut) {
3976          double newB=0.0;
3977          ClpSimplex::Status iStatus = getStatus(iSetOut);
3978          if (iStatus==ClpSimplex::atLowerBound)
3979            newB=lower_[iSetOut];
3980          else
3981            newB=upper_[iSetOut];
3982          // subtract out others at bounds
3983          if ((gubType_&8)==0) {
3984            int stop = -(key+1);
3985            int iColumn =next_[key];
3986            // skip basic
3987            while (iColumn>=0)
3988              iColumn = next_[iColumn];
3989            // sum all non-key variables
3990            while(iColumn!=stop) {
3991              assert (iColumn<0);
3992              iColumn = -iColumn-1;
3993              newB -= solution[iColumn];
3994              iColumn=next_[iColumn];
3995            }
3996          }
3997          if (newB)
3998            ClpPackedMatrix::add(model,rhsOffset_,key,-newB);
3999        }
4000      }
4001    }
4002  }
4003#ifdef CLP_DEBUG
4004  // debug
4005  {
4006    int i;
4007    char * xxxx = new char[numberColumns];
4008    memset(xxxx,0,numberColumns);
4009    for (i=0;i<numberRows;i++) {
4010      int iPivot = pivotVariable[i];
4011      assert (model->getStatus(iPivot)==ClpSimplex::basic);
4012      if (iPivot<numberColumns&&backward_[iPivot]>=0)
4013        xxxx[iPivot]=1;
4014    }
4015    double primalTolerance = model->primalTolerance();
4016    for (i=0;i<numberSets_;i++) {
4017      int key=keyVariable_[i];
4018      double value=0.0;
4019      // sum over all except key
4020      int iColumn =next_[key];
4021      // sum all non-key variables
4022      int k=0;
4023      int stop = -(key+1);
4024      while (iColumn!=stop) {
4025        if (iColumn<0)
4026          iColumn = -iColumn-1;
4027        value += solution[iColumn];
4028        k++;
4029        assert (k<100);
4030        assert (backward_[iColumn]==i);
4031        iColumn=next_[iColumn];
4032      }
4033      iColumn = next_[key];
4034      if (key<numberColumns) {
4035        // feasibility will be done later
4036        assert (getStatus(i)!=ClpSimplex::basic);
4037        double sol;
4038        if (getStatus(i)==ClpSimplex::atUpperBound)
4039          sol = upper_[i]-value;
4040        else
4041          sol = lower_[i]-value;
4042        //printf("xx Value of key structural %d for set %d is %g - cost %g\n",key,i,sol,
4043        //     cost[key]);
4044        //if (fabs(sol-solution[key])>1.0e-3)
4045        //printf("** stored value was %g\n",solution[key]);
4046      } else {
4047        // slack is key
4048        double infeasibility=0.0;
4049        if (value>upper_[i]+primalTolerance) {
4050          infeasibility=value-upper_[i]-primalTolerance;
4051          //setAbove(i);
4052        } else if (value<lower_[i]-primalTolerance) {
4053          infeasibility=lower_[i]-value-primalTolerance ;
4054          //setBelow(i);
4055        } else {
4056          //setFeasible(i);
4057        }
4058        //printf("xx Value of key slack for set %d is %g\n",i,value);
4059      }
4060      while (iColumn>=0) {
4061        assert (xxxx[iColumn]);
4062        xxxx[iColumn]=0;
4063        iColumn=next_[iColumn];
4064      }
4065    }
4066    for (i=0;i<numberColumns;i++) {
4067      if (i<numberColumns&&backward_[i]>=0) {
4068        assert (!xxxx[i]||i==keyVariable_[backward_[i]]);
4069      }
4070    }
4071    delete [] xxxx;
4072  }
4073#endif
4074  return 0;
4075}
4076// Switches off dj checking each factorization (for BIG models)
4077void 
4078ClpGubMatrix::switchOffCheck()
4079{
4080  noCheck_=0;
4081  infeasibilityWeight_=0.0;
4082}
4083// Correct sequence in and out to give true value
4084void 
4085ClpGubMatrix::correctSequence(int & sequenceIn, int & sequenceOut) const
4086{
4087  sequenceIn = trueSequenceIn_;
4088  sequenceOut = trueSequenceOut_;
4089}
Note: See TracBrowser for help on using the repository browser.