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

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

changes to try and make more stable

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