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

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

out some printf statements

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