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

Last change on this file since 1402 was 1402, checked in by forrest, 12 years ago

get rid of compiler warnings

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