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

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

various fixes plus slightly weighted pricing plus lagomory switches

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