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

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

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