source: trunk/ClpGubMatrix.cpp @ 437

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

to CoinMax? and CoinMin?

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