source: stable/1.15/Clp/src/ClpPrimalColumnSteepest.cpp @ 1955

Last change on this file since 1955 was 1955, checked in by forrest, 6 years ago

fix assert if variable wrongly marked

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