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

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

fix for ray and warning messages

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