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

Last change on this file since 1054 was 1054, checked in by forrest, 13 years ago

changes to compile under gcc 4.3

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