source: stable/1.15/Clp/src/ClpGubMatrix.cpp @ 1931

Last change on this file since 1931 was 1931, checked in by stefan, 7 years ago

sync with trunk rev 1930

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 186.3 KB
Line 
1/* $Id: ClpGubMatrix.cpp 1931 2013-04-06 20:44:29Z stefan $ */
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          //const int * columnLength = matrix_->getVectorLengths();
2119          //const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2120          //const int * row = matrix_->getIndices();
2121          //const double * elementByColumn = matrix_->getElements();
2122          //int * pivotVariable = model->pivotVariable();
2123          sumPrimalInfeasibilities_ = 0.0;
2124          numberPrimalInfeasibilities_ = 0;
2125          double primalTolerance = model->primalTolerance();
2126          double relaxedTolerance = primalTolerance;
2127          // we can't really trust infeasibilities if there is primal error
2128          double error = CoinMin(1.0e-2, model->largestPrimalError());
2129          // allow tolerance at least slightly bigger than standard
2130          relaxedTolerance = relaxedTolerance +  error;
2131          // but we will be using difference
2132          relaxedTolerance -= primalTolerance;
2133          sumOfRelaxedPrimalInfeasibilities_ = 0.0;
2134          for (i = 0; i < numberSets_; i++) { // Could just be over basics (esp if no bounds)
2135               int kColumn = keyVariable_[i];
2136               double value = 0.0;
2137               if ((gubType_ & 8) != 0) {
2138                    int iColumn = next_[kColumn];
2139                    // sum all non-key variables
2140                    while(iColumn >= 0) {
2141                         value += solution[iColumn];
2142                         iColumn = next_[iColumn];
2143                    }
2144               } else {
2145                    // bounds exist - sum over all except key
2146                    int stop = -(kColumn + 1);
2147                    int iColumn = next_[kColumn];
2148                    // sum all non-key variables
2149                    while(iColumn != stop) {
2150                         if (iColumn < 0)
2151                              iColumn = -iColumn - 1;
2152                         value += solution[iColumn];
2153                         iColumn = next_[iColumn];
2154                    }
2155               }
2156               if (kColumn < numberColumns) {
2157                    // make sure key is basic - so will be skipped in values pass
2158                    model->setStatus(kColumn, ClpSimplex::basic);
2159                    // feasibility will be done later
2160                    assert (getStatus(i) != ClpSimplex::basic);
2161                    if (getStatus(i) == ClpSimplex::atUpperBound)
2162                         solution[kColumn] = upper_[i] - value;
2163                    else
2164                         solution[kColumn] = lower_[i] - value;
2165                    //printf("Value of key structural %d for set %d is %g\n",kColumn,i,solution[kColumn]);
2166               } else {
2167                    // slack is key
2168                    assert (getStatus(i) == ClpSimplex::basic);
2169                    double infeasibility = 0.0;
2170                    if (value > upper_[i] + primalTolerance) {
2171                         infeasibility = value - upper_[i] - primalTolerance;
2172                         setAbove(i);
2173                    } else if (value < lower_[i] - primalTolerance) {
2174                         infeasibility = lower_[i] - value - primalTolerance ;
2175                         setBelow(i);
2176                    } else {
2177                         setFeasible(i);
2178                    }
2179                    //printf("Value of key slack for set %d is %g\n",i,value);
2180                    if (infeasibility > 0.0) {
2181                         sumPrimalInfeasibilities_ += infeasibility;
2182                         if (infeasibility > relaxedTolerance)
2183                              sumOfRelaxedPrimalInfeasibilities_ += infeasibility;
2184                         numberPrimalInfeasibilities_ ++;
2185                    }
2186               }
2187          }
2188     }
2189     break;
2190     // Report on infeasibilities of key variables
2191     case 2: {
2192          model->setSumPrimalInfeasibilities(model->sumPrimalInfeasibilities() +
2193                                             sumPrimalInfeasibilities_);
2194          model->setNumberPrimalInfeasibilities(model->numberPrimalInfeasibilities() +
2195                                                numberPrimalInfeasibilities_);
2196          model->setSumOfRelaxedPrimalInfeasibilities(model->sumOfRelaxedPrimalInfeasibilities() +
2197                    sumOfRelaxedPrimalInfeasibilities_);
2198     }
2199     break;
2200     }
2201}
2202/*
2203     utility dual function for dealing with dynamic constraints
2204     mode=n see ClpGubMatrix.hpp for definition
2205     Remember to update here when settled down
2206*/
2207void
2208ClpGubMatrix::dualExpanded(ClpSimplex * model,
2209                           CoinIndexedVector * array,
2210                           double * /*other*/, int mode)
2211{
2212     switch (mode) {
2213          // modify costs before transposeUpdate
2214     case 0: {
2215          int i;
2216          double * cost = model->costRegion();
2217          // not dual values yet
2218          //assert (!other);
2219          //double * work = array->denseVector();
2220          double infeasibilityCost = model->infeasibilityCost();
2221          int * pivotVariable = model->pivotVariable();
2222          int numberRows = model->numberRows();
2223          int numberColumns = model->numberColumns();
2224          for (i = 0; i < numberRows; i++) {
2225               int iPivot = pivotVariable[i];
2226               if (iPivot < numberColumns) {
2227                    int iSet = backward_[iPivot];
2228                    if (iSet >= 0) {
2229                         int kColumn = keyVariable_[iSet];
2230                         double costValue;
2231                         if (kColumn < numberColumns) {
2232                              // structural has cost
2233                              costValue = cost[kColumn];
2234                         } else {
2235                              // slack is key
2236                              assert (getStatus(iSet) == ClpSimplex::basic);
2237                              // negative as -1.0 for slack
2238                              costValue = -weight(iSet) * infeasibilityCost;
2239                         }
2240                         array->add(i, -costValue); // was work[i]-costValue
2241                    }
2242               }
2243          }
2244     }
2245     break;
2246     // create duals for key variables (without check on dual infeasible)
2247     case 1: {
2248          // If key slack then dual 0.0 (if feasible)
2249          // dj for key is zero so that defines dual on set
2250          int i;
2251          double * dj = model->djRegion();
2252          int numberColumns = model->numberColumns();
2253          double infeasibilityCost = model->infeasibilityCost();
2254          for (i = 0; i < numberSets_; i++) {
2255               int kColumn = keyVariable_[i];
2256               if (kColumn < numberColumns) {
2257                    // dj without set
2258                    double value = dj[kColumn];
2259                    // Now subtract out from all
2260                    dj[kColumn] = 0.0;
2261                    int iColumn = next_[kColumn];
2262                    // modify all non-key variables
2263                    while(iColumn >= 0) {
2264                         dj[iColumn] -= value;
2265                         iColumn = next_[iColumn];
2266                    }
2267               } else {
2268                    // slack key - may not be feasible
2269                    assert (getStatus(i) == ClpSimplex::basic);
2270                    // negative as -1.0 for slack
2271                    double value = -weight(i) * infeasibilityCost;
2272                    if (value) {
2273                         int iColumn = next_[kColumn];
2274                         // modify all non-key variables basic
2275                         while(iColumn >= 0) {
2276                              dj[iColumn] -= value;
2277                              iColumn = next_[iColumn];
2278                         }
2279                    }
2280               }
2281          }
2282     }
2283     break;
2284     // as 1 but check slacks and compute djs
2285     case 2: {
2286          // If key slack then dual 0.0
2287          // If not then slack could be dual infeasible
2288          // dj for key is zero so that defines dual on set
2289          int i;
2290          // make sure fromIndex will not confuse pricing
2291          fromIndex_[0] = -1;
2292          possiblePivotKey_ = -1;
2293          // Create array
2294          int numberColumns = model->numberColumns();
2295          int * pivotVariable = model->pivotVariable();
2296          int numberRows = model->numberRows();
2297          for (i = 0; i < numberRows; i++) {
2298               int iPivot = pivotVariable[i];
2299               if (iPivot < numberColumns)
2300                    backToPivotRow_[iPivot] = i;
2301          }
2302          if (noCheck_ >= 0) {
2303               if (infeasibilityWeight_ != model->infeasibilityCost()) {
2304                    // don't bother checking
2305                    sumDualInfeasibilities_ = 100.0;
2306                    numberDualInfeasibilities_ = 1;
2307                    sumOfRelaxedDualInfeasibilities_ = 100.0;
2308                    return;
2309               }
2310          }
2311          double * dj = model->djRegion();
2312          double * dual = model->dualRowSolution();
2313          double * cost = model->costRegion();
2314          ClpSimplex::Status iStatus;
2315          const int * columnLength = matrix_->getVectorLengths();
2316          const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2317          const int * row = matrix_->getIndices();
2318          const double * elementByColumn = matrix_->getElements();
2319          double infeasibilityCost = model->infeasibilityCost();
2320          sumDualInfeasibilities_ = 0.0;
2321          numberDualInfeasibilities_ = 0;
2322          double dualTolerance = model->dualTolerance();
2323          double relaxedTolerance = dualTolerance;
2324          // we can't really trust infeasibilities if there is dual error
2325          double error = CoinMin(1.0e-2, model->largestDualError());
2326          // allow tolerance at least slightly bigger than standard
2327          relaxedTolerance = relaxedTolerance +  error;
2328          // but we will be using difference
2329          relaxedTolerance -= dualTolerance;
2330          sumOfRelaxedDualInfeasibilities_ = 0.0;
2331          for (i = 0; i < numberSets_; i++) {
2332               int kColumn = keyVariable_[i];
2333               if (kColumn < numberColumns) {
2334                    // dj without set
2335                    double value = cost[kColumn];
2336                    for (CoinBigIndex j = columnStart[kColumn];
2337                              j < columnStart[kColumn] + columnLength[kColumn]; j++) {
2338                         int iRow = row[j];
2339                         value -= dual[iRow] * elementByColumn[j];
2340                    }
2341                    // Now subtract out from all
2342                    dj[kColumn] -= value;
2343                    int stop = -(kColumn + 1);
2344                    kColumn = next_[kColumn];
2345                    while (kColumn != stop) {
2346                         if (kColumn < 0)
2347                              kColumn = -kColumn - 1;
2348                         double djValue = dj[kColumn] - value;
2349                         dj[kColumn] = djValue;;
2350                         double infeasibility = 0.0;
2351                         iStatus = model->getStatus(kColumn);
2352                         if (iStatus == ClpSimplex::atLowerBound) {
2353                              if (djValue < -dualTolerance)
2354                                   infeasibility = -djValue - dualTolerance;
2355                         } else if (iStatus == ClpSimplex::atUpperBound) {
2356                              // at upper bound
2357                              if (djValue > dualTolerance)
2358                                   infeasibility = djValue - dualTolerance;
2359                         }
2360                         if (infeasibility > 0.0) {
2361                              sumDualInfeasibilities_ += infeasibility;
2362                              if (infeasibility > relaxedTolerance)
2363                                   sumOfRelaxedDualInfeasibilities_ += infeasibility;
2364                              numberDualInfeasibilities_ ++;
2365                         }
2366                         kColumn = next_[kColumn];
2367                    }
2368                    // check slack
2369                    iStatus = getStatus(i);
2370                    assert (iStatus != ClpSimplex::basic);
2371                    double infeasibility = 0.0;
2372                    // dj of slack is -(-1.0)value
2373                    if (iStatus == ClpSimplex::atLowerBound) {
2374                         if (value < -dualTolerance)
2375                              infeasibility = -value - dualTolerance;
2376                    } else if (iStatus == ClpSimplex::atUpperBound) {
2377                         // at upper bound
2378                         if (value > dualTolerance)
2379                              infeasibility = value - dualTolerance;
2380                    }
2381                    if (infeasibility > 0.0) {
2382                         sumDualInfeasibilities_ += infeasibility;
2383                         if (infeasibility > relaxedTolerance)
2384                              sumOfRelaxedDualInfeasibilities_ += infeasibility;
2385                         numberDualInfeasibilities_ ++;
2386                    }
2387               } else {
2388                    // slack key - may not be feasible
2389                    assert (getStatus(i) == ClpSimplex::basic);
2390                    // negative as -1.0 for slack
2391                    double value = -weight(i) * infeasibilityCost;
2392                    if (value) {
2393                         // Now subtract out from all
2394                         int kColumn = i + numberColumns;
2395                         int stop = -(kColumn + 1);
2396                         kColumn = next_[kColumn];
2397                         while (kColumn != stop) {
2398                              if (kColumn < 0)
2399                                   kColumn = -kColumn - 1;
2400                              double djValue = dj[kColumn] - value;
2401                              dj[kColumn] = djValue;;
2402                              double infeasibility = 0.0;
2403                              iStatus = model->getStatus(kColumn);
2404                              if (iStatus == ClpSimplex::atLowerBound) {
2405                                   if (djValue < -dualTolerance)
2406                                        infeasibility = -djValue - dualTolerance;
2407                              } else if (iStatus == ClpSimplex::atUpperBound) {
2408                                   // at upper bound
2409                                   if (djValue > dualTolerance)
2410                                        infeasibility = djValue - dualTolerance;
2411                              }
2412                              if (infeasibility > 0.0) {
2413                                   sumDualInfeasibilities_ += infeasibility;
2414                                   if (infeasibility > relaxedTolerance)
2415                                        sumOfRelaxedDualInfeasibilities_ += infeasibility;
2416                                   numberDualInfeasibilities_ ++;
2417                              }
2418                              kColumn = next_[kColumn];
2419                         }
2420                    }
2421               }
2422          }
2423          // and get statistics for column generation
2424          synchronize(model, 4);
2425          infeasibilityWeight_ = -1.0;
2426     }
2427     break;
2428     // Report on infeasibilities of key variables
2429     case 3: {
2430          model->setSumDualInfeasibilities(model->sumDualInfeasibilities() +
2431                                           sumDualInfeasibilities_);
2432          model->setNumberDualInfeasibilities(model->numberDualInfeasibilities() +
2433                                              numberDualInfeasibilities_);
2434          model->setSumOfRelaxedDualInfeasibilities(model->sumOfRelaxedDualInfeasibilities() +
2435                    sumOfRelaxedDualInfeasibilities_);
2436     }
2437     break;
2438     // modify costs before transposeUpdate for partial pricing
2439     case 4: {
2440          // First compute new costs etc for interesting gubs
2441          int iLook = 0;
2442          int iSet = fromIndex_[0];
2443          double primalTolerance = model->primalTolerance();
2444          const double * cost = model->costRegion();
2445          double * solution = model->solutionRegion();
2446          double infeasibilityCost = model->infeasibilityCost();
2447          int numberColumns = model->numberColumns();
2448          int numberChanged = 0;
2449          int * pivotVariable = model->pivotVariable();
2450          while (iSet >= 0) {
2451               int key = keyVariable_[iSet];
2452               double value = 0.0;
2453               // sum over all except key
2454               if ((gubType_ & 8) != 0) {
2455                    int iColumn = next_[key];
2456                    // sum all non-key variables
2457                    while(iColumn >= 0) {
2458                         value += solution[iColumn];
2459                         iColumn = next_[iColumn];
2460                    }
2461               } else {
2462                    // bounds exist - sum over all except key
2463                    int stop = -(key + 1);
2464                    int iColumn = next_[key];
2465                    // sum all non-key variables
2466                    while(iColumn != stop) {
2467                         if (iColumn < 0)
2468                              iColumn = -iColumn - 1;
2469                         value += solution[iColumn];
2470                         iColumn = next_[iColumn];
2471                    }
2472               }
2473               double costChange;
2474               double oldCost = changeCost_[iLook];
2475               if (key < numberColumns) {
2476                    assert (getStatus(iSet) != ClpSimplex::basic);
2477                    double sol;
2478                    if (getStatus(iSet) == ClpSimplex::atUpperBound)
2479                         sol = upper_[iSet] - value;
2480                    else
2481                         sol = lower_[iSet] - value;
2482                    solution[key] = sol;
2483                    // fix up cost
2484                    model->nonLinearCost()->setOne(key, sol);
2485#ifdef CLP_DEBUG_PRINT
2486                    printf("yy Value of key structural %d for set %d is %g - cost %g old cost %g\n", key, iSet, sol,
2487                           cost[key], oldCost);
2488#endif
2489                    costChange = cost[key] - oldCost;
2490               } else {
2491                    // slack is key
2492                    if (value > upper_[iSet] + primalTolerance) {
2493                         setAbove(iSet);
2494                    } else if (value < lower_[iSet] - primalTolerance) {
2495                         setBelow(iSet);
2496                    } else {
2497                         setFeasible(iSet);
2498                    }
2499                    // negative as -1.0 for slack
2500                    costChange = -weight(iSet) * infeasibilityCost - oldCost;
2501#ifdef CLP_DEBUG_PRINT
2502                    printf("yy Value of key slack for set %d is %g - cost %g old cost %g\n", iSet, value,
2503                           weight(iSet)*infeasibilityCost, oldCost);
2504#endif
2505               }
2506               if (costChange) {
2507                    fromIndex_[numberChanged] = iSet;
2508                    toIndex_[iSet] = numberChanged;
2509                    changeCost_[numberChanged++] = costChange;
2510               }
2511               iSet = fromIndex_[++iLook];
2512          }
2513          if (numberChanged || possiblePivotKey_ >= 0) {
2514               // first do those in list already
2515               int number = array->getNumElements();
2516               array->setPacked();
2517               int i;
2518               double * work = array->denseVector();
2519               int * which = array->getIndices();
2520               for (i = 0; i < number; i++) {
2521                    int iRow = which[i];
2522                    int iPivot = pivotVariable[iRow];
2523                    if (iPivot < numberColumns) {
2524                         int iSet = backward_[iPivot];
2525                         if (iSet >= 0 && toIndex_[iSet] >= 0) {
2526                              double newValue = work[i] + changeCost_[toIndex_[iSet]];
2527                              if (!newValue)
2528                                   newValue = 1.0e-100;
2529                              work[i] = newValue;
2530                              // mark as done
2531                              backward_[iPivot] = -1;
2532                         }
2533                    }
2534                    if (possiblePivotKey_ == iRow) {
2535                         double newValue = work[i] - model->dualIn();
2536                         if (!newValue)
2537                              newValue = 1.0e-100;
2538                         work[i] = newValue;
2539                         possiblePivotKey_ = -1;
2540                    }
2541               }
2542               // now do rest and clean up
2543               for (i = 0; i < numberChanged; i++) {
2544                    int iSet = fromIndex_[i];
2545                    int key = keyVariable_[iSet];
2546                    int iColumn = next_[key];
2547                    double change = changeCost_[i];
2548                    while (iColumn >= 0) {
2549                         if (backward_[iColumn] >= 0) {
2550                              int iRow = backToPivotRow_[iColumn];
2551                              assert (iRow >= 0);
2552                              work[number] = change;
2553                              if (possiblePivotKey_ == iRow) {
2554                                   double newValue = work[number] - model->dualIn();
2555                                   if (!newValue)
2556                                        newValue = 1.0e-100;
2557                                   work[number] = newValue;
2558                                   possiblePivotKey_ = -1;
2559                              }
2560                              which[number++] = iRow;
2561                         } else {
2562                              // reset
2563                              backward_[iColumn] = iSet;
2564                         }
2565                         iColumn = next_[iColumn];
2566                    }
2567                    toIndex_[iSet] = -1;
2568               }
2569               if (possiblePivotKey_ >= 0) {
2570                    work[number] = -model->dualIn();
2571                    which[number++] = possiblePivotKey_;
2572                    possiblePivotKey_ = -1;
2573               }
2574               fromIndex_[0] = -1;
2575               array->setNumElements(number);
2576          }
2577     }
2578     break;
2579     }
2580}
2581// This is local to Gub to allow synchronization when status is good
2582int
2583ClpGubMatrix::synchronize(ClpSimplex *, int)
2584{
2585     return 0;
2586}
2587/*
2588     general utility function for dealing with dynamic constraints
2589     mode=n see ClpGubMatrix.hpp for definition
2590     Remember to update here when settled down
2591*/
2592int
2593ClpGubMatrix::generalExpanded(ClpSimplex * model, int mode, int &number)
2594{
2595     int returnCode = 0;
2596     int numberColumns = model->numberColumns();
2597     switch (mode) {
2598          // Fill in pivotVariable but not for key variables
2599     case 0: {
2600          if (!next_ ) {
2601               // do ordering
2602               assert (!rhsOffset_);
2603               // create and do gub crash
2604               useEffectiveRhs(model, false);
2605          }
2606          int i;
2607          int numberBasic = number;
2608          // Use different array so can build from true pivotVariable_
2609          //int * pivotVariable = model->pivotVariable();
2610          int * pivotVariable = model->rowArray(0)->getIndices();
2611          for (i = 0; i < numberColumns; i++) {
2612               if (model->getColumnStatus(i) == ClpSimplex::basic) {
2613                    int iSet = backward_[i];
2614                    if (iSet < 0 || i != keyVariable_[iSet])
2615                         pivotVariable[numberBasic++] = i;
2616               }
2617          }
2618          number = numberBasic;
2619          if (model->numberIterations())
2620               assert (number == model->numberRows());
2621     }
2622     break;
2623     // Make all key variables basic
2624     case 1: {
2625          int i;
2626          for (i = 0; i < numberSets_; i++) {
2627               int iColumn = keyVariable_[i];
2628               if (iColumn < numberColumns)
2629                    model->setColumnStatus(iColumn, ClpSimplex::basic);
2630          }
2631     }
2632     break;
2633     // Do initial extra rows + maximum basic
2634     case 2: {
2635          returnCode = getNumRows() + 1;
2636          number = model->numberRows() + numberSets_;
2637     }
2638     break;
2639     // Before normal replaceColumn
2640     case 3: {
2641          int sequenceIn = model->sequenceIn();
2642          int sequenceOut = model->sequenceOut();
2643          int numberColumns = model->numberColumns();
2644          int numberRows = model->numberRows();
2645          int pivotRow = model->pivotRow();
2646          if (gubSlackIn_ >= 0)
2647               assert (sequenceIn > numberRows + numberColumns);
2648          if (sequenceIn == sequenceOut)
2649               return -1;
2650          int iSetIn = -1;
2651          int iSetOut = -1;
2652          if (sequenceOut < numberColumns) {
2653               iSetOut = backward_[sequenceOut];
2654          } else if (sequenceOut >= numberRows + numberColumns) {
2655               assert (pivotRow >= numberRows);
2656               int iExtra = pivotRow - numberRows;
2657               assert (iExtra >= 0);
2658               if (iSetOut < 0)
2659                    iSetOut = fromIndex_[iExtra];
2660               else
2661                    assert(iSetOut == fromIndex_[iExtra]);
2662          }
2663          if (sequenceIn < numberColumns) {
2664               iSetIn = backward_[sequenceIn];
2665          } else if (gubSlackIn_ >= 0) {
2666               iSetIn = gubSlackIn_;
2667          }
2668          possiblePivotKey_ = -1;
2669          number = 0; // say do ordinary
2670          int * pivotVariable = model->pivotVariable();
2671          if (pivotRow >= numberRows) {
2672               int iExtra = pivotRow - numberRows;
2673               //const int * length = matrix_->getVectorLengths();
2674
2675               assert (sequenceOut >= numberRows + numberColumns ||
2676                       sequenceOut == keyVariable_[iSetOut]);
2677               int incomingColumn = sequenceIn; // to be used in updates
2678               if (iSetIn != iSetOut) {
2679                    // We need to find a possible pivot for incoming
2680                    // look through rowArray_[1]
2681                    int n = model->rowArray(1)->getNumElements();
2682                    int * which = model->rowArray(1)->getIndices();
2683                    double * array = model->rowArray(1)->denseVector();
2684                    double bestAlpha = 1.0e-5;
2685                    //int shortest=numberRows+1;
2686                    for (int i = 0; i < n; i++) {
2687                         int iRow = which[i];
2688                         int iPivot = pivotVariable[iRow];
2689                         if (iPivot < numberColumns && backward_[iPivot] == iSetOut) {
2690                              if (fabs(array[i]) > fabs(bestAlpha)) {
2691                                   bestAlpha = array[i];
2692                                   possiblePivotKey_ = iRow;
2693                              }
2694                         }
2695                    }
2696                    assert (possiblePivotKey_ >= 0); // could set returnCode=4
2697                    number = 1;
2698                    if (sequenceIn >= numberRows + numberColumns) {
2699                         number = 3;
2700                         // need swap as gub slack in and must become key
2701                         // is this best way
2702                         int key = keyVariable_[iSetIn];
2703                         assert (key < numberColumns);
2704                         // check other basic
2705                         int iColumn = next_[key];
2706                         // set new key to be used by unpack
2707                         keyVariable_[iSetIn] = iSetIn + numberColumns;
2708                         // change cost in changeCost
2709                         {
2710                              int iLook = 0;
2711                              int iSet = fromIndex_[0];
2712                              while (iSet >= 0) {
2713                                   if (iSet == iSetIn) {
2714                                        changeCost_[iLook] = 0.0;
2715                                        break;
2716                                   }
2717                                   iSet = fromIndex_[++iLook];
2718                              }
2719                         }
2720                         while (iColumn >= 0) {
2721                              if (iColumn != sequenceOut) {
2722                                   // need partial ftran and skip accuracy check in replaceColumn
2723#ifdef CLP_DEBUG_PRINT
2724                                   printf("TTTTTry 5\n");
2725#endif
2726                                   int iRow = backToPivotRow_[iColumn];
2727                                   assert (iRow >= 0);
2728                                   unpack(model, model->rowArray(3), iColumn);
2729                                   model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
2730                                   double alpha = model->rowArray(3)->denseVector()[iRow];
2731                                   //if (!alpha)
2732                                   //printf("zero alpha a\n");
2733                                   int updateStatus = model->factorization()->replaceColumn(model,
2734                                                      model->rowArray(2),
2735                                                      model->rowArray(3),
2736                                                      iRow, alpha);
2737                                   returnCode = CoinMax(updateStatus, returnCode);
2738                                   model->rowArray(3)->clear();
2739                                   if (returnCode)
2740                                        break;
2741                              }
2742                              iColumn = next_[iColumn];
2743                         }
2744                         if (!returnCode) {
2745                              // now factorization looks as if key is out
2746                              // pivot back in
2747#ifdef CLP_DEBUG_PRINT
2748                              printf("TTTTTry 6\n");
2749#endif
2750                              unpack(model, model->rowArray(3), key);
2751                              model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
2752                              pivotRow = possiblePivotKey_;
2753                              double alpha = model->rowArray(3)->denseVector()[pivotRow];
2754                              //if (!alpha)
2755                              //printf("zero alpha b\n");
2756                              int updateStatus = model->factorization()->replaceColumn(model,
2757                                                 model->rowArray(2),
2758                                                 model->rowArray(3),
2759                                                 pivotRow, alpha);
2760                              returnCode = CoinMax(updateStatus, returnCode);
2761                              model->rowArray(3)->clear();
2762                         }
2763                         // restore key
2764                         keyVariable_[iSetIn] = key;
2765                         // now alternate column can replace key on out
2766                         incomingColumn = pivotVariable[possiblePivotKey_];
2767                    } else {
2768#ifdef CLP_DEBUG_PRINT
2769                         printf("TTTTTTry 4 %d\n", possiblePivotKey_);
2770#endif
2771                         int updateStatus = model->factorization()->replaceColumn(model,
2772                                            model->rowArray(2),
2773                                            model->rowArray(1),
2774                                            possiblePivotKey_,
2775                                            bestAlpha);
2776                         returnCode = CoinMax(updateStatus, returnCode);
2777                         incomingColumn = pivotVariable[possiblePivotKey_];
2778                    }
2779
2780                    //returnCode=4; // need swap
2781               } else {
2782                    // key swap
2783                    number = -1;
2784               }
2785               int key = keyVariable_[iSetOut];
2786               if (key < numberColumns)
2787                    assert(key == sequenceOut);
2788               // check if any other basic
2789               int iColumn = next_[key];
2790               if (returnCode)
2791                    iColumn = -1; // skip if error on previous
2792               // set new key to be used by unpack
2793               if (incomingColumn < numberColumns)
2794                    keyVariable_[iSetOut] = incomingColumn;
2795               else
2796                    keyVariable_[iSetOut] = iSetIn + numberColumns;
2797               double * cost = model->costRegion();
2798               if (possiblePivotKey_ < 0) {
2799                    double dj = model->djRegion()[sequenceIn] - cost[sequenceIn];
2800                    changeCost_[iExtra] = -dj;
2801#ifdef CLP_DEBUG_PRINT
2802                    printf("modifying changeCost %d by %g - cost %g\n", iExtra, dj, cost[sequenceIn]);
2803#endif
2804               }
2805               while (iColumn >= 0) {
2806                    if (iColumn != incomingColumn) {
2807                         number = -2;
2808                         // need partial ftran and skip accuracy check in replaceColumn
2809#ifdef CLP_DEBUG_PRINT
2810                         printf("TTTTTTry 1\n");
2811#endif
2812                         int iRow = backToPivotRow_[iColumn];
2813                         assert (iRow >= 0 && iRow < numberRows);
2814                         unpack(model, model->rowArray(3), iColumn);
2815                         model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
2816                         double * array = model->rowArray(3)->denseVector();
2817                         double alpha = array[iRow];
2818                         //if (!alpha)
2819                         //printf("zero alpha d\n");
2820                         int updateStatus = model->factorization()->replaceColumn(model,
2821                                            model->rowArray(2),
2822                                            model->rowArray(3),
2823                                            iRow, alpha);
2824                         returnCode = CoinMax(updateStatus, returnCode);
2825                         model->rowArray(3)->clear();
2826                         if (returnCode)
2827                              break;
2828                    }
2829                    iColumn = next_[iColumn];
2830               }
2831               // restore key
2832               keyVariable_[iSetOut] = key;
2833          } else if (sequenceIn >= numberRows + numberColumns) {
2834               number = 2;
2835               //returnCode=4;
2836               // need swap as gub slack in and must become key
2837               // is this best way
2838               int key = keyVariable_[iSetIn];
2839               assert (key < numberColumns);
2840               // check other basic
2841               int iColumn = next_[key];
2842               // set new key to be used by unpack
2843               keyVariable_[iSetIn] = iSetIn + numberColumns;
2844               // change cost in changeCost
2845               {
2846                    int iLook = 0;
2847                    int iSet = fromIndex_[0];
2848                    while (iSet >= 0) {
2849                         if (iSet == iSetIn) {
2850                              changeCost_[iLook] = 0.0;
2851                              break;
2852                         }
2853                         iSet = fromIndex_[++iLook];
2854                    }
2855               }
2856               while (iColumn >= 0) {
2857                    if (iColumn != sequenceOut) {
2858                         // need partial ftran and skip accuracy check in replaceColumn
2859#ifdef CLP_DEBUG_PRINT
2860                         printf("TTTTTry 2\n");
2861#endif
2862                         int iRow = backToPivotRow_[iColumn];
2863                         assert (iRow >= 0);
2864                         unpack(model, model->rowArray(3), iColumn);
2865                         model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
2866                         double alpha = model->rowArray(3)->denseVector()[iRow];
2867                         //if (!alpha)
2868                         //printf("zero alpha e\n");
2869                         int updateStatus = model->factorization()->replaceColumn(model,
2870                                            model->rowArray(2),
2871                                            model->rowArray(3),
2872                                            iRow, alpha);
2873                         returnCode = CoinMax(updateStatus, returnCode);
2874                         model->rowArray(3)->clear();
2875                         if (returnCode)
2876                              break;
2877                    }
2878                    iColumn = next_[iColumn];
2879               }
2880               if (!returnCode) {
2881                    // now factorization looks as if key is out
2882                    // pivot back in
2883#ifdef CLP_DEBUG_PRINT
2884                    printf("TTTTTry 3\n");
2885#endif
2886                    unpack(model, model->rowArray(3), key);
2887                    model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
2888                    double alpha = model->rowArray(3)->denseVector()[pivotRow];
2889                    //if (!alpha)
2890                    //printf("zero alpha f\n");
2891                    int updateStatus = model->factorization()->replaceColumn(model,
2892                                       model->rowArray(2),
2893                                       model->rowArray(3),
2894                                       pivotRow, alpha);
2895                    returnCode = CoinMax(updateStatus, returnCode);
2896                    model->rowArray(3)->clear();
2897               }
2898               // restore key
2899               keyVariable_[iSetIn] = key;
2900          } else {
2901               // normal - but might as well do here
2902               returnCode = model->factorization()->replaceColumn(model,
2903                            model->rowArray(2),
2904                            model->rowArray(1),
2905                            model->pivotRow(),
2906                            model->alpha());
2907          }
2908     }
2909#ifdef CLP_DEBUG_PRINT
2910     printf("Update type after %d - status %d - pivot row %d\n",
2911            number, returnCode, model->pivotRow());
2912#endif
2913     // see if column generation says time to re-factorize
2914     returnCode = CoinMax(returnCode, synchronize(model, 5));
2915     number = -1; // say no need for normal replaceColumn
2916     break;
2917     // To see if can dual or primal
2918     case 4: {
2919          returnCode = 1;
2920     }
2921     break;
2922     // save status
2923     case 5: {
2924          synchronize(model, 0);
2925          CoinMemcpyN(status_, numberSets_, saveStatus_);
2926          CoinMemcpyN(keyVariable_, numberSets_, savedKeyVariable_);
2927     }
2928     break;
2929     // restore status
2930     case 6: {
2931          CoinMemcpyN(saveStatus_, numberSets_, status_);
2932          CoinMemcpyN(savedKeyVariable_, numberSets_, keyVariable_);
2933          // restore firstAvailable_
2934          synchronize(model, 7);
2935          // redo next_
2936          int i;
2937          int * last = new int[numberSets_];
2938          for (i = 0; i < numberSets_; i++) {
2939               int iKey = keyVariable_[i];
2940               assert(iKey >= numberColumns || backward_[iKey] == i);
2941               last[i] = iKey;
2942               // make sure basic
2943               //if (iKey<numberColumns)
2944               //model->setStatus(iKey,ClpSimplex::basic);
2945          }
2946          for (i = 0; i < numberColumns; i++) {
2947               int iSet = backward_[i];
2948               if (iSet >= 0) {
2949                    next_[last[iSet]] = i;
2950                    last[iSet] = i;
2951               }
2952          }
2953          for (i = 0; i < numberSets_; i++) {
2954               next_[last[i]] = -(keyVariable_[i] + 1);
2955               redoSet(model, keyVariable_[i], keyVariable_[i], i);
2956          }
2957          delete [] last;
2958          // redo pivotVariable
2959          int * pivotVariable = model->pivotVariable();
2960          int iRow;
2961          int numberBasic = 0;
2962          int numberRows = model->numberRows();
2963          for (iRow = 0; iRow < numberRows; iRow++) {
2964               if (model->getRowStatus(iRow) == ClpSimplex::basic) {
2965                    numberBasic++;
2966                    pivotVariable[iRow] = iRow + numberColumns;
2967               } else {
2968                    pivotVariable[iRow] = -1;
2969               }
2970          }
2971          i = 0;
2972          int iColumn;
2973          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2974               if (model->getStatus(iColumn) == ClpSimplex::basic) {
2975                    int iSet = backward_[iColumn];
2976                    if (iSet < 0 || keyVariable_[iSet] != iColumn) {
2977                         while (pivotVariable[i] >= 0) {
2978                              i++;
2979                              assert (i < numberRows);
2980                         }
2981                         pivotVariable[i] = iColumn;
2982                         backToPivotRow_[iColumn] = i;
2983                         numberBasic++;
2984                    }
2985               }
2986          }
2987          assert (numberBasic == numberRows);
2988          rhsOffset(model, true);
2989     }
2990     break;
2991     // flag a variable
2992     case 7: {
2993          assert (number == model->sequenceIn());
2994          synchronize(model, 1);
2995          synchronize(model, 8);
2996     }
2997     break;
2998     // unflag all variables
2999     case 8: {
3000          returnCode = synchronize(model, 2);
3001     }
3002     break;
3003     // redo costs in primal
3004     case 9: {
3005          returnCode = synchronize(model, 3);
3006     }
3007     break;
3008     // return 1 if there may be changing bounds on variable (column generation)
3009     case 10: {
3010          returnCode = synchronize(model, 6);
3011     }
3012     break;
3013     // make sure set is clean
3014     case 11: {
3015          assert (number == model->sequenceIn());
3016          returnCode = synchronize(model, 8);
3017     }
3018     break;
3019     default:
3020          break;
3021     }
3022     return returnCode;
3023}
3024// Sets up an effective RHS and does gub crash if needed
3025void
3026ClpGubMatrix::useEffectiveRhs(ClpSimplex * model, bool cheapest)
3027{
3028     // Do basis - cheapest or slack if feasible (unless cheapest set)
3029     int longestSet = 0;
3030     int iSet;
3031     for (iSet = 0; iSet < numberSets_; iSet++)
3032          longestSet = CoinMax(longestSet, end_[iSet] - start_[iSet]);
3033
3034     double * upper = new double[longestSet+1];
3035     double * cost = new double[longestSet+1];
3036     double * lower = new double[longestSet+1];
3037     double * solution = new double[longestSet+1];
3038     assert (!next_);
3039     int numberColumns = getNumCols();
3040     const int * columnLength = matrix_->getVectorLengths();
3041     const double * columnLower = model->lowerRegion();
3042     const double * columnUpper = model->upperRegion();
3043     double * columnSolution = model->solutionRegion();
3044     const double * objective = model->costRegion();
3045     int numberRows = getNumRows();
3046     toIndex_ = new int[numberSets_];
3047     for (iSet = 0; iSet < numberSets_; iSet++)
3048          toIndex_[iSet] = -1;
3049     fromIndex_ = new int [getNumRows()+1];
3050     double tolerance = model->primalTolerance();
3051     bool noNormalBounds = true;
3052     gubType_ &= ~8;
3053     bool gotBasis = false;
3054     for (iSet = 0; iSet < numberSets_; iSet++) {
3055          if (keyVariable_[iSet] < numberColumns)
3056               gotBasis = true;
3057          CoinBigIndex j;
3058          CoinBigIndex iStart = start_[iSet];
3059          CoinBigIndex iEnd = end_[iSet];
3060          for (j = iStart; j < iEnd; j++) {
3061               if (columnLower[j] && columnLower[j] > -1.0e20)
3062                    noNormalBounds = false;
3063               if (columnUpper[j] && columnUpper[j] < 1.0e20)
3064                    noNormalBounds = false;
3065          }
3066     }
3067     if (noNormalBounds)
3068          gubType_ |= 8;
3069     if (!gotBasis) {
3070          for (iSet = 0; iSet < numberSets_; iSet++) {
3071               CoinBigIndex j;
3072               int numberBasic = 0;
3073               int iBasic = -1;
3074               CoinBigIndex iStart = start_[iSet];
3075               CoinBigIndex iEnd = end_[iSet];
3076               // find one with smallest length
3077               int smallest = numberRows + 1;
3078               double value = 0.0;
3079               for (j = iStart; j < iEnd; j++) {
3080                    if (model->getStatus(j) == ClpSimplex::basic) {
3081                         if (columnLength[j] < smallest) {
3082                              smallest = columnLength[j];
3083                              iBasic = j;
3084                         }
3085                         numberBasic++;
3086                    }
3087                    value += columnSolution[j];
3088               }
3089               bool done = false;
3090               if (numberBasic > 1 || (numberBasic == 1 && getStatus(iSet) == ClpSimplex::basic)) {
3091                    if (getStatus(iSet) == ClpSimplex::basic)
3092                         iBasic = iSet + numberColumns; // slack key - use
3093                    done = true;
3094               } else if (numberBasic == 1) {
3095                    // see if can be key
3096                    double thisSolution = columnSolution[iBasic];
3097                    if (thisSolution > columnUpper[iBasic]) {
3098                         value -= thisSolution - columnUpper[iBasic];
3099                         thisSolution = columnUpper[iBasic];
3100                         columnSolution[iBasic] = thisSolution;
3101                    }
3102                    if (thisSolution < columnLower[iBasic]) {
3103                         value -= thisSolution - columnLower[iBasic];
3104                         thisSolution = columnLower[iBasic];
3105                         columnSolution[iBasic] = thisSolution;
3106                    }
3107                    // try setting slack to a bound
3108                    assert (upper_[iSet] < 1.0e20 || lower_[iSet] > -1.0e20);
3109                    double cost1 = COIN_DBL_MAX;
3110                    int whichBound = -1;
3111                    if (upper_[iSet] < 1.0e20) {
3112                         // try slack at ub
3113                         double newBasic = thisSolution + upper_[iSet] - value;
3114                         if (newBasic >= columnLower[iBasic] - tolerance &&
3115                                   newBasic <= columnUpper[iBasic] + tolerance) {
3116                              // can go
3117                              whichBound = 1;
3118                              cost1 = newBasic * objective[iBasic];
3119                              // But if exact then may be good solution
3120                              if (fabs(upper_[iSet] - value) < tolerance)
3121                                   cost1 = -COIN_DBL_MAX;
3122                         }
3123                    }
3124                    if (lower_[iSet] > -1.0e20) {
3125                         // try slack at lb
3126                         double newBasic = thisSolution + lower_[iSet] - value;
3127                         if (newBasic >= columnLower[iBasic] - tolerance &&
3128                                   newBasic <= columnUpper[iBasic] + tolerance) {
3129                              // can go but is it cheaper
3130                              double cost2 = newBasic * objective[iBasic];
3131                              // But if exact then may be good solution
3132                              if (fabs(lower_[iSet] - value) < tolerance)
3133                                   cost2 = -COIN_DBL_MAX;
3134                              if (cost2 < cost1)
3135                                   whichBound = 0;
3136                         }
3137                    }
3138                    if (whichBound != -1) {
3139                         // key
3140                         done = true;
3141                         if (whichBound) {
3142                              // slack to upper
3143                              columnSolution[iBasic] = thisSolution + upper_[iSet] - value;
3144                              setStatus(iSet, ClpSimplex::atUpperBound);
3145                         } else {
3146                              // slack to lower
3147                              columnSolution[iBasic] = thisSolution + lower_[iSet] - value;
3148                              setStatus(iSet, ClpSimplex::atLowerBound);
3149                         }
3150                    }
3151               }
3152               if (!done) {
3153                    if (!cheapest) {
3154                         // see if slack can be key
3155                         if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
3156                              done = true;
3157                              setStatus(iSet, ClpSimplex::basic);
3158                              iBasic = iSet + numberColumns;
3159                         }
3160                    }
3161                    if (!done) {
3162                         // set non basic if there was one
3163                         if (iBasic >= 0)
3164                              model->setStatus(iBasic, ClpSimplex::atLowerBound);
3165                         // find cheapest
3166                         int numberInSet = iEnd - iStart;
3167                         CoinMemcpyN(columnLower + iStart, numberInSet, lower);
3168                         CoinMemcpyN(columnUpper + iStart, numberInSet, upper);
3169                         CoinMemcpyN(columnSolution + iStart, numberInSet, solution);
3170                         // and slack
3171                         iBasic = numberInSet;
3172                         solution[iBasic] = -value;
3173                         lower[iBasic] = -upper_[iSet];
3174                         upper[iBasic] = -lower_[iSet];
3175                         int kphase;
3176                         if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
3177                              // feasible
3178                              kphase = 1;
3179                              cost[iBasic] = 0.0;
3180                              CoinMemcpyN(objective + iStart, numberInSet, cost);
3181                         } else {
3182                              // infeasible
3183                              kphase = 0;
3184                              // remember bounds are flipped so opposite to natural
3185                              if (value < lower_[iSet] - tolerance)
3186                                   cost[iBasic] = 1.0;
3187                              else
3188                                   cost[iBasic] = -1.0;
3189                              CoinZeroN(cost, numberInSet);
3190                         }
3191                         double dualTolerance = model->dualTolerance();
3192                         for (int iphase = kphase; iphase < 2; iphase++) {
3193                              if (iphase) {
3194                                   cost[numberInSet] = 0.0;
3195                                   CoinMemcpyN(objective + iStart, numberInSet, cost);
3196                              }
3197                              // now do one row lp
3198                              bool improve = true;
3199                              while (improve) {
3200                                   improve = false;
3201                                   double dual = cost[iBasic];
3202                                   int chosen = -1;
3203                                   double best = dualTolerance;
3204                                   int way = 0;
3205                                   for (int i = 0; i <= numberInSet; i++) {
3206                                        double dj = cost[i] - dual;
3207                                        double improvement = 0.0;
3208                                        if (iphase || i < numberInSet)
3209                                             assert (solution[i] >= lower[i] && solution[i] <= upper[i]);
3210                                        if (dj > dualTolerance)
3211                                             improvement = dj * (solution[i] - lower[i]);
3212                                        else if (dj < -dualTolerance)
3213                                             improvement = dj * (solution[i] - upper[i]);
3214                                        if (improvement > best) {
3215                                             best = improvement;
3216                                             chosen = i;
3217                                             if (dj < 0.0) {
3218                                                  way = 1;
3219                                             } else {
3220                                                  way = -1;
3221                                             }
3222                                        }
3223                                   }
3224                                   if (chosen >= 0) {
3225                                        improve = true;
3226                                        // now see how far
3227                                        if (way > 0) {
3228                                             // incoming increasing so basic decreasing
3229                                             // if phase 0 then go to nearest bound
3230                                             double distance = upper[chosen] - solution[chosen];
3231                                             double basicDistance;
3232                                             if (!iphase) {
3233                                                  assert (iBasic == numberInSet);
3234                                                  assert (solution[iBasic] > upper[iBasic]);
3235                                                  basicDistance = solution[iBasic] - upper[iBasic];
3236                                             } else {
3237                                                  basicDistance = solution[iBasic] - lower[iBasic];
3238                                             }
3239                                             // need extra coding for unbounded
3240                                             assert (CoinMin(distance, basicDistance) < 1.0e20);
3241                                             if (distance > basicDistance) {
3242                                                  // incoming becomes basic
3243                                                  solution[chosen] += basicDistance;
3244                                                  if (!iphase)
3245                                                       solution[iBasic] = upper[iBasic];
3246                                                  else
3247                                                       solution[iBasic] = lower[iBasic];
3248                                                  iBasic = chosen;
3249                                             } else {
3250                                                  // flip
3251                                                  solution[chosen] = upper[chosen];
3252                                                  solution[iBasic] -= distance;
3253                                             }
3254                                        } else {
3255                                             // incoming decreasing so basic increasing
3256                                             // if phase 0 then go to nearest bound
3257                                             double distance = solution[chosen] - lower[chosen];
3258                                             double basicDistance;
3259                                             if (!iphase) {
3260                                                  assert (iBasic == numberInSet);
3261                                                  assert (solution[iBasic] < lower[iBasic]);
3262                                                  basicDistance = lower[iBasic] - solution[iBasic];
3263                                             } else {
3264                                                  basicDistance = upper[iBasic] - solution[iBasic];
3265                                             }
3266                                             // need extra coding for unbounded - for now just exit
3267                                             if (CoinMin(distance, basicDistance) > 1.0e20) {
3268                                                  printf("unbounded on set %d\n", iSet);
3269                                                  iphase = 1;
3270                                                  iBasic = numberInSet;
3271                                                  break;
3272                                             }
3273                                             if (distance > basicDistance) {
3274                                                  // incoming becomes basic
3275                                                  solution[chosen] -= basicDistance;
3276                                                  if (!iphase)
3277                                                       solution[iBasic] = lower[iBasic];
3278                                                  else
3279                                                       solution[iBasic] = upper[iBasic];
3280                                                  iBasic = chosen;
3281                                             } else {
3282                                                  // flip
3283                                                  solution[chosen] = lower[chosen];
3284                                                  solution[iBasic] += distance;
3285                                             }
3286                                        }
3287                                        if (!iphase) {
3288                                             if(iBasic < numberInSet)
3289                                                  break; // feasible
3290                                             else if (solution[iBasic] >= lower[iBasic] &&
3291                                                       solution[iBasic] <= upper[iBasic])
3292                                                  break; // feasible (on flip)
3293                                        }
3294                                   }
3295                              }
3296                         }
3297                         // convert iBasic back and do bounds
3298                         if (iBasic == numberInSet) {
3299                              // slack basic
3300                              setStatus(iSet, ClpSimplex::basic);
3301                              iBasic = iSet + numberColumns;
3302                         } else {
3303                              iBasic += start_[iSet];
3304                              model->setStatus(iBasic, ClpSimplex::basic);
3305                              // remember bounds flipped
3306                              if (upper[numberInSet] == lower[numberInSet])
3307                                   setStatus(iSet, ClpSimplex::isFixed);
3308                              else if (solution[numberInSet] == upper[numberInSet])
3309                                   setStatus(iSet, ClpSimplex::atLowerBound);
3310                              else if (solution[numberInSet] == lower[numberInSet])
3311                                   setStatus(iSet, ClpSimplex::atUpperBound);
3312                              else
3313                                   abort();
3314                         }
3315                         for (j = iStart; j < iEnd; j++) {
3316                              if (model->getStatus(j) != ClpSimplex::basic) {
3317                                   int inSet = j - iStart;
3318                                   columnSolution[j] = solution[inSet];
3319                                   if (upper[inSet] == lower[inSet])
3320                                        model->setStatus(j, ClpSimplex::isFixed);
3321                                   else if (solution[inSet] == upper[inSet])
3322                                        model->setStatus(j, ClpSimplex::atUpperBound);
3323                                   else if (solution[inSet] == lower[inSet])
3324                                        model->setStatus(j, ClpSimplex::atLowerBound);
3325                              }
3326                         }
3327                    }
3328               }
3329               keyVariable_[iSet] = iBasic;
3330          }
3331     }
3332     delete [] lower;
3333     delete [] solution;
3334     delete [] upper;
3335     delete [] cost;
3336     // make sure matrix is in good shape
3337     matrix_->orderMatrix();
3338     // create effective rhs
3339     delete [] rhsOffset_;
3340     rhsOffset_ = new double[numberRows];
3341     delete [] next_;
3342     next_ = new int[numberColumns+numberSets_+2*longestSet];
3343     char * mark = new char[numberColumns];
3344     memset(mark, 0, numberColumns);
3345     for (int iColumn = 0; iColumn < numberColumns; iColumn++)
3346          next_[iColumn] = COIN_INT_MAX;
3347     int i;
3348     int * keys = new int[numberSets_];
3349     for (i = 0; i < numberSets_; i++)
3350          keys[i] = COIN_INT_MAX;
3351     // set up chains
3352     for (i = 0; i < numberColumns; i++) {
3353          if (model->getStatus(i) == ClpSimplex::basic)
3354               mark[i] = 1;
3355          int iSet = backward_[i];
3356          if (iSet >= 0) {
3357               int iNext = keys[iSet];
3358               next_[i] = iNext;
3359               keys[iSet] = i;
3360          }
3361     }
3362     for (i = 0; i < numberSets_; i++) {
3363          int j;
3364          if (getStatus(i) != ClpSimplex::basic) {
3365               // make sure fixed if it is
3366               if (upper_[i] == lower_[i])
3367                    setStatus(i, ClpSimplex::isFixed);
3368               // slack not key - choose one with smallest length
3369               int smallest = numberRows + 1;
3370               int key = -1;
3371               j = keys[i];
3372               if (j != COIN_INT_MAX) {
3373                    while (1) {
3374                         if (mark[j] && columnLength[j] < smallest && !gotBasis) {
3375                              key = j;
3376                              smallest = columnLength[j];
3377                         }
3378                         if (next_[j] != COIN_INT_MAX) {
3379                              j = next_[j];
3380                         } else {
3381                              // correct end
3382                              next_[j] = -(keys[i] + 1);
3383                              break;
3384                         }
3385                    }
3386               } else {
3387                    next_[i+numberColumns] = -(numberColumns + i + 1);
3388               }
3389               if (gotBasis)
3390                    key = keyVariable_[i];
3391               if (key >= 0) {
3392                    keyVariable_[i] = key;
3393               } else {
3394                    // nothing basic - make slack key
3395                    //((ClpGubMatrix *)this)->setStatus(i,ClpSimplex::basic);
3396                    // fudge to avoid const problem
3397                    status_[i] = 1;
3398               }
3399          } else {
3400               // slack key
3401               keyVariable_[i] = numberColumns + i;
3402               int j;
3403               double sol = 0.0;
3404               j = keys[i];
3405               if (j != COIN_INT_MAX) {
3406                    while (1) {
3407                         sol += columnSolution[j];
3408                         if (next_[j] != COIN_INT_MAX) {
3409                              j = next_[j];
3410                         } else {
3411                              // correct end
3412                              next_[j] = -(keys[i] + 1);
3413                              break;
3414                         }
3415                    }
3416               } else {
3417                    next_[i+numberColumns] = -(numberColumns + i + 1);
3418               }
3419               if (sol > upper_[i] + tolerance) {
3420                    setAbove(i);
3421               } else if (sol < lower_[i] - tolerance) {
3422                    setBelow(i);
3423               } else {
3424                    setFeasible(i);
3425               }
3426          }
3427          // Create next_
3428          int key = keyVariable_[i];
3429          redoSet(model, key, keys[i], i);
3430     }
3431     delete [] keys;
3432     delete [] mark;
3433     rhsOffset(model, true);
3434}
3435// redoes next_ for a set.
3436void
3437ClpGubMatrix::redoSet(ClpSimplex * model, int newKey, int oldKey, int iSet)
3438{
3439     int numberColumns = model->numberColumns();
3440     int * save = next_ + numberColumns + numberSets_;
3441     int number = 0;
3442     int stop = -(oldKey + 1);
3443     int j = next_[oldKey];
3444     while (j != stop) {
3445          if (j < 0)
3446               j = -j - 1;
3447          if (j != newKey)
3448               save[number++] = j;
3449          j = next_[j];
3450     }
3451     // and add oldkey
3452     if (newKey != oldKey)
3453          save[number++] = oldKey;
3454     // now do basic
3455     int lastMarker = -(newKey + 1);
3456     keyVariable_[iSet] = newKey;
3457     next_[newKey] = lastMarker;
3458     int last = newKey;
3459     for ( j = 0; j < number; j++) {
3460          int iColumn = save[j];
3461          if (iColumn < numberColumns) {
3462               if (model->getStatus(iColumn) == ClpSimplex::basic) {
3463                    next_[last] = iColumn;
3464                    next_[iColumn] = lastMarker;
3465                    last = iColumn;
3466               }
3467          }
3468     }
3469     // now add in non-basic
3470     for ( j = 0; j < number; j++) {
3471          int iColumn = save[j];
3472          if (iColumn < numberColumns) {
3473               if (model->getStatus(iColumn) != ClpSimplex::basic) {
3474                    next_[last] = -(iColumn + 1);
3475                    next_[iColumn] = lastMarker;
3476                    last = iColumn;
3477               }
3478          }
3479     }
3480
3481}
3482/* Returns effective RHS if it is being used.  This is used for long problems
3483   or big gub or anywhere where going through full columns is
3484   expensive.  This may re-compute */
3485double *
3486ClpGubMatrix::rhsOffset(ClpSimplex * model, bool forceRefresh, bool
3487#ifdef CLP_DEBUG
3488                        check
3489#endif
3490                       )
3491{
3492     //forceRefresh=true;
3493     if (rhsOffset_) {
3494#ifdef CLP_DEBUG
3495          if (check) {
3496               // no need - but check anyway
3497               // zero out basic
3498               int numberRows = model->numberRows();
3499               int numberColumns = model->numberColumns();
3500               double * solution = new double [numberColumns];
3501               double * rhs = new double[numberRows];
3502               CoinMemcpyN(model->solutionRegion(), numberColumns, solution);
3503               CoinZeroN(rhs, numberRows);
3504               int iRow;
3505               for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3506                    if (model->getColumnStatus(iColumn) == ClpSimplex::basic)
3507                         solution[iColumn] = 0.0;
3508               }
3509               for (int iSet = 0; iSet < numberSets_; iSet++) {
3510                    int iColumn = keyVariable_[iSet];
3511                    if (iColumn < numberColumns)
3512                         solution[iColumn] = 0.0;
3513               }
3514               times(-1.0, solution, rhs);
3515               delete [] solution;
3516               const double * columnSolution = model->solutionRegion();
3517               // and now subtract out non basic
3518               ClpSimplex::Status iStatus;
3519               for (int iSet = 0; iSet < numberSets_; iSet++) {
3520                    int iColumn = keyVariable_[iSet];
3521                    if (iColumn < numberColumns) {
3522                         double b = 0.0;
3523                         // key is structural - where is slack
3524                         iStatus = getStatus(iSet);
3525                         assert (iStatus != ClpSimplex::basic);
3526                         if (iStatus == ClpSimplex::atLowerBound)
3527                              b = lower_[iSet];
3528                         else
3529                              b = upper_[iSet];
3530                         // subtract out others at bounds
3531                         if ((gubType_ & 8) == 0) {
3532                              int stop = -(iColumn + 1);
3533                              int jColumn = next_[iColumn];
3534                              // sum all non-basic variables - first skip basic
3535                              while(jColumn >= 0)
3536                                   jColumn = next_[jColumn];
3537                              while(jColumn != stop) {
3538                                   assert (jColumn < 0);
3539                                   jColumn = -jColumn - 1;
3540                                   b -= columnSolution[jColumn];
3541                                   jColumn = next_[jColumn];
3542                              }
3543                         }
3544                         // subtract out
3545                         ClpPackedMatrix::add(model, rhs, iColumn, -b);
3546                    }
3547               }
3548               for (iRow = 0; iRow < numberRows; iRow++) {
3549                    if (fabs(rhs[iRow] - rhsOffset_[iRow]) > 1.0e-3)
3550                         printf("** bad effective %d - true %g old %g\n", iRow, rhs[iRow], rhsOffset_[iRow]);
3551               }
3552               delete [] rhs;
3553          }
3554#endif
3555          if (forceRefresh || (refreshFrequency_ && model->numberIterations() >=
3556                               lastRefresh_ + refreshFrequency_)) {
3557               // zero out basic
3558               int numberRows = model->numberRows();
3559               int numberColumns = model->numberColumns();
3560               double * solution = new double [numberColumns];
3561               CoinMemcpyN(model->solutionRegion(), numberColumns, solution);
3562               CoinZeroN(rhsOffset_, numberRows);
3563               for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3564                    if (model->getColumnStatus(iColumn) == ClpSimplex::basic)
3565                         solution[iColumn] = 0.0;
3566               }
3567               int iSet;
3568               for ( iSet = 0; iSet < numberSets_; iSet++) {
3569                    int iColumn = keyVariable_[iSet];
3570                    if (iColumn < numberColumns)
3571                         solution[iColumn] = 0.0;
3572               }
3573               times(-1.0, solution, rhsOffset_);
3574               delete [] solution;
3575               lastRefresh_ = model->numberIterations();
3576               const double * columnSolution = model->solutionRegion();
3577               // and now subtract out non basic
3578               ClpSimplex::Status iStatus;
3579               for ( iSet = 0; iSet < numberSets_; iSet++) {
3580                    int iColumn = keyVariable_[iSet];
3581                    if (iColumn < numberColumns) {
3582                         double b = 0.0;
3583                         // key is structural - where is slack
3584                         iStatus = getStatus(iSet);
3585                         assert (iStatus != ClpSimplex::basic);
3586                         if (iStatus == ClpSimplex::atLowerBound)
3587                              b = lower_[iSet];
3588                         else
3589                              b = upper_[iSet];
3590                         // subtract out others at bounds
3591                         if ((gubType_ & 8) == 0) {
3592                              int stop = -(iColumn + 1);
3593                              int jColumn = next_[iColumn];
3594                              // sum all non-basic variables - first skip basic
3595                              while(jColumn >= 0)
3596                                   jColumn = next_[jColumn];
3597                              while(jColumn != stop) {
3598                                   assert (jColumn < 0);
3599                                   jColumn = -jColumn - 1;
3600                                   b -= columnSolution[jColumn];
3601                                   jColumn = next_[jColumn];
3602                              }
3603                         }
3604                         // subtract out
3605                         if (b)
3606                              ClpPackedMatrix::add(model, rhsOffset_, iColumn, -b);
3607                    }
3608               }
3609          }
3610     }
3611     return rhsOffset_;
3612}
3613/*
3614   update information for a pivot (and effective rhs)
3615*/
3616int
3617ClpGubMatrix::updatePivot(ClpSimplex * model, double oldInValue, double /*oldOutValue*/)
3618{
3619     int sequenceIn = model->sequenceIn();
3620     int sequenceOut = model->sequenceOut();
3621     double * solution = model->solutionRegion();
3622     int numberColumns = model->numberColumns();
3623     int numberRows = model->numberRows();
3624     int pivotRow = model->pivotRow();
3625     int iSetIn;
3626     // Correct sequence in
3627     trueSequenceIn_ = sequenceIn;
3628     if (sequenceIn < numberColumns) {
3629          iSetIn = backward_[sequenceIn];
3630     } else if (sequenceIn < numberColumns + numberRows) {
3631          iSetIn = -1;
3632     } else {
3633          iSetIn = gubSlackIn_;
3634          trueSequenceIn_ = numberColumns + numberRows + iSetIn;
3635     }
3636     int iSetOut = -1;
3637     trueSequenceOut_ = sequenceOut;
3638     if (sequenceOut < numberColumns) {
3639          iSetOut = backward_[sequenceOut];
3640     } else if (sequenceOut >= numberRows + numberColumns) {
3641          assert (pivotRow >= numberRows);
3642          int iExtra = pivotRow - numberRows;
3643          assert (iExtra >= 0);
3644          if (iSetOut < 0)
3645               iSetOut = fromIndex_[iExtra];
3646          else
3647               assert(iSetOut == fromIndex_[iExtra]);
3648          trueSequenceOut_ = numberColumns + numberRows + iSetOut;
3649     }
3650     if (rhsOffset_) {
3651          // update effective rhs
3652          if (sequenceIn == sequenceOut) {
3653               assert (sequenceIn < numberRows + numberColumns); // should be easy to deal with
3654               if (sequenceIn < numberColumns)
3655                    add(model, rhsOffset_, sequenceIn, oldInValue - solution[sequenceIn]);
3656          } else {
3657               if (sequenceIn < numberColumns) {
3658                    // we need to test if WILL be key
3659                    ClpPackedMatrix::add(model, rhsOffset_, sequenceIn, oldInValue);
3660                    if (iSetIn >= 0)  {
3661                         // old contribution to rhsOffset_
3662                         int key = keyVariable_[iSetIn];
3663                         if (key < numberColumns) {
3664                              double oldB = 0.0;
3665                              ClpSimplex::Status iStatus = getStatus(iSetIn);
3666                              if (iStatus == ClpSimplex::atLowerBound)
3667                                   oldB = lower_[iSetIn];
3668                              else
3669                                   oldB = upper_[iSetIn];
3670                              // subtract out others at bounds
3671                              if ((gubType_ & 8) == 0) {
3672                                   int stop = -(key + 1);
3673                                   int iColumn = next_[key];
3674                                   // skip basic
3675                                   while (iColumn >= 0)
3676                                        iColumn = next_[iColumn];
3677                                   // sum all non-key variables
3678                                   while(iColumn != stop) {
3679                                        assert (iColumn < 0);
3680                                        iColumn = -iColumn - 1;
3681                                        if (iColumn == sequenceIn)
3682                                             oldB -= oldInValue;
3683                                        else if ( iColumn != sequenceOut )
3684                                             oldB -= solution[iColumn];
3685                                        iColumn = next_[iColumn];
3686                                   }
3687                              }
3688                              if (oldB)
3689                                   ClpPackedMatrix::add(model, rhsOffset_, key, oldB);
3690                         }
3691                    }
3692               } else if (sequenceIn < numberRows + numberColumns) {
3693                    //rhsOffset_[sequenceIn-numberColumns] -= oldInValue;
3694               } else {
3695#ifdef CLP_DEBUG_PRINT
3696                    printf("** in is key slack %d\n", sequenceIn);
3697#endif
3698                    // old contribution to rhsOffset_
3699                    int key = keyVariable_[iSetIn];
3700                    if (key < numberColumns) {
3701                         double oldB = 0.0;
3702                         ClpSimplex::Status iStatus = getStatus(iSetIn);
3703                         if (iStatus == ClpSimplex::atLowerBound)
3704                              oldB = lower_[iSetIn];
3705                         else
3706                              oldB = upper_[iSetIn];
3707                         // subtract out others at bounds
3708                         if ((gubType_ & 8) == 0) {
3709                              int stop = -(key + 1);
3710                              int iColumn = next_[key];
3711                              // skip basic
3712                              while (iColumn >= 0)
3713                                   iColumn = next_[iColumn];
3714                              // sum all non-key variables
3715                              while(iColumn != stop) {
3716                                   assert (iColumn < 0);
3717                                   iColumn = -iColumn - 1;
3718                                   if ( iColumn != sequenceOut )
3719                                        oldB -= solution[iColumn];
3720                                   iColumn = next_[iColumn];
3721                              }
3722                         }
3723                         if (oldB)
3724                              ClpPackedMatrix::add(model, rhsOffset_, key, oldB);
3725                    }
3726               }
3727               if (sequenceOut < numberColumns) {
3728                    ClpPackedMatrix::add(model, rhsOffset_, sequenceOut, -solution[sequenceOut]);
3729                    if (iSetOut >= 0) {
3730                         // old contribution to rhsOffset_
3731                         int key = keyVariable_[iSetOut];
3732                         if (key < numberColumns && iSetIn != iSetOut) {
3733                              double oldB = 0.0;
3734                              ClpSimplex::Status iStatus = getStatus(iSetOut);
3735                              if (iStatus == ClpSimplex::atLowerBound)
3736                                   oldB = lower_[iSetOut];
3737                              else
3738                                   oldB = upper_[iSetOut];
3739                              // subtract out others at bounds
3740                              if ((gubType_ & 8) == 0) {
3741                                   int stop = -(key + 1);
3742                                   int iColumn = next_[key];
3743                                   // skip basic
3744                                   while (iColumn >= 0)
3745                                        iColumn = next_[iColumn];
3746                                   // sum all non-key variables
3747                                   while(iColumn != stop) {
3748                                        assert (iColumn < 0);
3749                                        iColumn = -iColumn - 1;
3750                                        if (iColumn == sequenceIn)
3751                                             oldB -= oldInValue;
3752                                        else if ( iColumn != sequenceOut )
3753                                             oldB -= solution[iColumn];
3754                                        iColumn = next_[iColumn];
3755                                   }
3756                              }
3757                              if (oldB)
3758                                   ClpPackedMatrix::add(model, rhsOffset_, key, oldB);
3759                         }
3760                    }
3761               } else if (sequenceOut < numberRows + numberColumns) {
3762                    //rhsOffset_[sequenceOut-numberColumns] -= -solution[sequenceOut];
3763               } else {
3764#ifdef CLP_DEBUG_PRINT
3765                    printf("** out is key slack %d\n", sequenceOut);
3766#endif
3767                    assert (pivotRow >= numberRows);
3768               }
3769          }
3770     }
3771     int * pivotVariable = model->pivotVariable();
3772     // may need to deal with key
3773     // Also need coding to mark/allow key slack entering
3774     if (pivotRow >= numberRows) {
3775          assert (sequenceOut >= numberRows + numberColumns || sequenceOut == keyVariable_[iSetOut]);
3776#ifdef CLP_DEBUG_PRINT
3777          if (sequenceIn >= numberRows + numberColumns)
3778               printf("key slack %d in, set out %d\n", gubSlackIn_, iSetOut);
3779          printf("** danger - key out for set %d in %d (set %d)\n", iSetOut, sequenceIn,
3780                 iSetIn);
3781#endif
3782          // if slack out mark correctly
3783          if (sequenceOut >= numberRows + numberColumns) {
3784               double value = model->valueOut();
3785               if (value == upper_[iSetOut]) {
3786                    setStatus(iSetOut, ClpSimplex::atUpperBound);
3787               } else if (value == lower_[iSetOut]) {
3788                    setStatus(iSetOut, ClpSimplex::atLowerBound);
3789               } else {
3790                    if (fabs(value - upper_[iSetOut]) <
3791                              fabs(value - lower_[iSetOut])) {
3792                         setStatus(iSetOut, ClpSimplex::atUpperBound);
3793                    } else {
3794                         setStatus(iSetOut, ClpSimplex::atLowerBound);
3795                    }
3796               }
3797               if (upper_[iSetOut] == lower_[iSetOut])
3798                    setStatus(iSetOut, ClpSimplex::isFixed);
3799               setFeasible(iSetOut);
3800          }
3801          if (iSetOut == iSetIn) {
3802               // key swap
3803               int key;
3804               if (sequenceIn >= numberRows + numberColumns) {
3805                    key = numberColumns + iSetIn;
3806                    setStatus(iSetIn, ClpSimplex::basic);
3807               } else {
3808                    key = sequenceIn;
3809               }
3810               redoSet(model, key, keyVariable_[iSetIn], iSetIn);
3811          } else {
3812               // key was chosen
3813               assert (possiblePivotKey_ >= 0 && possiblePivotKey_ < numberRows);
3814               int key = pivotVariable[possiblePivotKey_];
3815               // and set incoming here
3816               if (sequenceIn >= numberRows + numberColumns) {
3817                    // slack in - so use old key
3818                    sequenceIn = keyVariable_[iSetIn];
3819                    model->setStatus(sequenceIn, ClpSimplex::basic);
3820                    setStatus(iSetIn, ClpSimplex::basic);
3821                    redoSet(model, iSetIn + numberColumns, keyVariable_[iSetIn], iSetIn);
3822               }
3823               //? do not do if iSetIn<0 ? as will be done later
3824               pivotVariable[possiblePivotKey_] = sequenceIn;
3825               if (sequenceIn < numberColumns)
3826                    backToPivotRow_[sequenceIn] = possiblePivotKey_;
3827               redoSet(model, key, keyVariable_[iSetOut], iSetOut);
3828          }
3829     } else {
3830          if (sequenceOut < numberColumns) {
3831               if (iSetIn >= 0 && iSetOut == iSetIn) {
3832                    // key not out - only problem is if slack in
3833                    int key;
3834                    if (sequenceIn >= numberRows + numberColumns) {
3835                         key = numberColumns + iSetIn;
3836                         setStatus(iSetIn, ClpSimplex::basic);
3837                         assert (pivotRow < numberRows);
3838                         // must swap with current key
3839                         int key = keyVariable_[iSetIn];
3840                         model->setStatus(key, ClpSimplex::basic);
3841                         pivotVariable[pivotRow] = key;
3842                         backToPivotRow_[key] = pivotRow;
3843                    } else {
3844                         key = keyVariable_[iSetIn];
3845                    }
3846                    redoSet(model, key, keyVariable_[iSetIn], iSetIn);
3847               } else if (iSetOut >= 0) {
3848                    // just redo set
3849                    int key = keyVariable_[iSetOut];;
3850                    redoSet(model, key, keyVariable_[iSetOut], iSetOut);
3851               }
3852          }
3853     }
3854     if (iSetIn >= 0 && iSetIn != iSetOut) {
3855          int key = keyVariable_[iSetIn];
3856          if (sequenceIn == numberColumns + 2 * numberRows) {
3857               // key slack in
3858               assert (pivotRow < numberRows);
3859               // must swap with current key
3860               model->setStatus(key, ClpSimplex::basic);
3861               pivotVariable[pivotRow] = key;
3862               backToPivotRow_[key] = pivotRow;
3863               setStatus(iSetIn, ClpSimplex::basic);
3864               key = iSetIn + numberColumns;
3865          }
3866          // redo set to allow for new one
3867          redoSet(model, key, keyVariable_[iSetIn], iSetIn);
3868     }
3869     // update pivot
3870     if (sequenceIn < numberColumns) {
3871          if (pivotRow < numberRows) {
3872               backToPivotRow_[sequenceIn] = pivotRow;
3873          } else {
3874               if (possiblePivotKey_ >= 0) {
3875                    assert (possiblePivotKey_ < numberRows);
3876                    backToPivotRow_[sequenceIn] = possiblePivotKey_;
3877                    pivotVariable[possiblePivotKey_] = sequenceIn;
3878               }
3879          }
3880     } else if (sequenceIn >= numberRows + numberColumns) {
3881          // key in - something should have been done before
3882          int key = keyVariable_[iSetIn];
3883          assert (key == numberColumns + iSetIn);
3884          //pivotVariable[pivotRow]=key;
3885          //backToPivotRow_[key]=pivotRow;
3886          //model->setStatus(key,ClpSimplex::basic);
3887          //key=numberColumns+iSetIn;
3888          setStatus(iSetIn, ClpSimplex::basic);
3889          redoSet(model, key, keyVariable_[iSetIn], iSetIn);
3890     }
3891#ifdef CLP_DEBUG
3892     {
3893          char * xx = new char[numberColumns+numberRows];
3894          memset(xx, 0, numberRows + numberColumns);
3895          for (int i = 0; i < numberRows; i++) {
3896               int iPivot = pivotVariable[i];
3897               assert (iPivot < numberRows + numberColumns);
3898               assert (!xx[iPivot]);
3899               xx[iPivot] = 1;
3900               if (iPivot < numberColumns) {
3901                    int iBack = backToPivotRow_[iPivot];
3902                    assert (i == iBack);
3903               }
3904          }
3905          delete [] xx;
3906     }
3907#endif
3908     if (rhsOffset_) {
3909          // update effective rhs
3910          if (sequenceIn != sequenceOut) {
3911               if (sequenceIn < numberColumns) {
3912                    if (iSetIn >= 0) {
3913                         // new contribution to rhsOffset_
3914                         int key = keyVariable_[iSetIn];
3915                         if (key < numberColumns) {
3916                              double newB = 0.0;
3917                              ClpSimplex::Status iStatus = getStatus(iSetIn);
3918                              if (iStatus == ClpSimplex::atLowerBound)
3919                                   newB = lower_[iSetIn];
3920                              else
3921                                   newB = upper_[iSetIn];
3922                              // subtract out others at bounds
3923                              if ((gubType_ & 8) == 0) {
3924                                   int stop = -(key + 1);
3925                                   int iColumn = next_[key];
3926                                   // skip basic
3927                                   while (iColumn >= 0)
3928                                        iColumn = next_[iColumn];
3929                                   // sum all non-key variables
3930                                   while(iColumn != stop) {
3931                                        assert (iColumn < 0);
3932                                        iColumn = -iColumn - 1;
3933                                        newB -= solution[iColumn];
3934                                        iColumn = next_[iColumn];
3935                                   }
3936                              }
3937                              if (newB)
3938                                   ClpPackedMatrix::add(model, rhsOffset_, key, -newB);
3939                         }
3940                    }
3941               }
3942               if (iSetOut >= 0) {
3943                    // new contribution to rhsOffset_
3944                    int key = keyVariable_[iSetOut];
3945                    if (key < numberColumns && iSetIn != iSetOut) {
3946                         double newB = 0.0;
3947                         ClpSimplex::Status iStatus = getStatus(iSetOut);
3948                         if (iStatus == ClpSimplex::atLowerBound)
3949                              newB = lower_[iSetOut];
3950                         else
3951                              newB = upper_[iSetOut];
3952                         // subtract out others at bounds
3953                         if ((gubType_ & 8) == 0) {
3954                              int stop = -(key + 1);
3955                              int iColumn = next_[key];
3956                              // skip basic
3957                              while (iColumn >= 0)
3958                                   iColumn = next_[iColumn];
3959                              // sum all non-key variables
3960                              while(iColumn != stop) {
3961                                   assert (iColumn < 0);
3962                                   iColumn = -iColumn - 1;
3963                                   newB -= solution[iColumn];
3964                                   iColumn = next_[iColumn];
3965                              }
3966                         }
3967                         if (newB)
3968                              ClpPackedMatrix::add(model, rhsOffset_, key, -newB);
3969                    }
3970               }
3971          }
3972     }
3973#ifdef CLP_DEBUG
3974     // debug
3975     {
3976          int i;
3977          char * xxxx = new char[numberColumns];
3978          memset(xxxx, 0, numberColumns);
3979          for (i = 0; i < numberRows; i++) {
3980               int iPivot = pivotVariable[i];
3981               assert (model->getStatus(iPivot) == ClpSimplex::basic);
3982               if (iPivot < numberColumns && backward_[iPivot] >= 0)
3983                    xxxx[iPivot] = 1;
3984          }
3985          double primalTolerance = model->primalTolerance();
3986          for (i = 0; i < numberSets_; i++) {
3987               int key = keyVariable_[i];
3988               double value = 0.0;
3989               // sum over all except key
3990               int iColumn = next_[key];
3991               // sum all non-key variables
3992               int k = 0;
3993               int stop = -(key + 1);
3994               while (iColumn != stop) {
3995                    if (iColumn < 0)
3996                         iColumn = -iColumn - 1;
3997                    value += solution[iColumn];
3998                    k++;
3999                    assert (k < 100);
4000                    assert (backward_[iColumn] == i);
4001                    iColumn = next_[iColumn];
4002               }
4003               iColumn = next_[key];
4004               if (key < numberColumns) {
4005                    // feasibility will be done later
4006                    assert (getStatus(i) != ClpSimplex::basic);
4007                    double sol;
4008                    if (getStatus(i) == ClpSimplex::atUpperBound)
4009                         sol = upper_[i] - value;
4010                    else
4011                         sol = lower_[i] - value;
4012                    //printf("xx Value of key structural %d for set %d is %g - cost %g\n",key,i,sol,
4013                    //     cost[key]);
4014                    //if (fabs(sol-solution[key])>1.0e-3)
4015                    //printf("** stored value was %g\n",solution[key]);
4016               } else {
4017                    // slack is key
4018                    double infeasibility = 0.0;
4019                    if (value > upper_[i] + primalTolerance) {
4020                         infeasibility = value - upper_[i] - primalTolerance;
4021                         //setAbove(i);
4022                    } else if (value < lower_[i] - primalTolerance) {
4023                         infeasibility = lower_[i] - value - primalTolerance ;
4024                         //setBelow(i);
4025                    } else {
4026                         //setFeasible(i);
4027                    }
4028                    //printf("xx Value of key slack for set %d is %g\n",i,value);
4029               }
4030               while (iColumn >= 0) {
4031                    assert (xxxx[iColumn]);
4032                    xxxx[iColumn] = 0;
4033                    iColumn = next_[iColumn];
4034               }
4035          }
4036          for (i = 0; i < numberColumns; i++) {
4037               if (i < numberColumns && backward_[i] >= 0) {
4038                    assert (!xxxx[i] || i == keyVariable_[backward_[i]]);
4039               }
4040          }
4041          delete [] xxxx;
4042     }
4043#endif
4044     return 0;
4045}
4046// Switches off dj checking each factorization (for BIG models)
4047void
4048ClpGubMatrix::switchOffCheck()
4049{
4050     noCheck_ = 0;
4051     infeasibilityWeight_ = 0.0;
4052}
4053// Correct sequence in and out to give true value
4054void
4055ClpGubMatrix::correctSequence(const ClpSimplex * /*model*/, int & sequenceIn, int & sequenceOut)
4056{
4057     if (sequenceIn != -999) {
4058          sequenceIn = trueSequenceIn_;
4059          sequenceOut = trueSequenceOut_;
4060     }
4061}
Note: See TracBrowser for help on using the repository browser.