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

Last change on this file since 1732 was 1732, checked in by forrest, 8 years ago

various fixes plus slightly weighted pricing plus lagomory switches

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