source: trunk/Clp/src/AbcDualRowSteepest.cpp @ 2385

Last change on this file since 2385 was 2385, checked in by unxusr, 4 months ago

formatting

  • Property svn:keywords set to Id
File size: 48.4 KB
Line 
1/* $Id: AbcDualRowSteepest.cpp 2385 2019-01-06 19:43:06Z unxusr $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others, Copyright (C) 2012, FasterCoin.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinPragma.hpp"
7#include "AbcSimplex.hpp"
8#include "ClpMessage.hpp"
9#include "AbcDualRowSteepest.hpp"
10#include "CoinIndexedVector.hpp"
11#include "AbcSimplexFactorization.hpp"
12#include "CoinHelperFunctions.hpp"
13#include "CoinAbcHelperFunctions.hpp"
14#include <cstdio>
15//#############################################################################
16// Constructors / Destructor / Assignment
17//#############################################################################
18//#define CLP_DEBUG 2
19//-------------------------------------------------------------------
20// Default Constructor
21//-------------------------------------------------------------------
22AbcDualRowSteepest::AbcDualRowSteepest(int mode)
23  : AbcDualRowPivot()
24  , norm_(0.0)
25  , factorizationRatio_(10.0)
26  , state_(-1)
27  , mode_(mode)
28  , persistence_(normal)
29  , weights_(NULL)
30  , infeasible_(NULL)
31  , savedWeights_(NULL)
32{
33  type_ = 2 + 64 * mode;
34  //printf("XX %x AbcDualRowSteepest constructor\n",this);
35}
36
37//-------------------------------------------------------------------
38// Copy constructor
39//-------------------------------------------------------------------
40AbcDualRowSteepest::AbcDualRowSteepest(const AbcDualRowSteepest &rhs)
41  : AbcDualRowPivot(rhs)
42{
43  //printf("XX %x AbcDualRowSteepest constructor from %x\n",this,&rhs);
44  norm_ = rhs.norm_;
45  factorizationRatio_ = rhs.factorizationRatio_;
46  state_ = rhs.state_;
47  mode_ = rhs.mode_;
48  persistence_ = rhs.persistence_;
49  model_ = rhs.model_;
50  if ((model_ && model_->whatsChanged() & 1) != 0) {
51    int number = model_->numberRows();
52    if (rhs.savedWeights_)
53      number = CoinMin(number, rhs.savedWeights_->capacity());
54    if (rhs.infeasible_) {
55      infeasible_ = new CoinIndexedVector(rhs.infeasible_);
56    } else {
57      infeasible_ = NULL;
58    }
59    if (rhs.weights_) {
60      weights_ = new CoinIndexedVector(rhs.weights_);
61    } else {
62      weights_ = NULL;
63    }
64    if (rhs.savedWeights_) {
65      savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
66    } else {
67      savedWeights_ = NULL;
68    }
69  } else {
70    infeasible_ = NULL;
71    weights_ = NULL;
72    savedWeights_ = NULL;
73  }
74}
75
76//-------------------------------------------------------------------
77// Destructor
78//-------------------------------------------------------------------
79AbcDualRowSteepest::~AbcDualRowSteepest()
80{
81  //printf("XX %x AbcDualRowSteepest destructor\n",this);
82  delete weights_;
83  delete infeasible_;
84  delete savedWeights_;
85}
86
87//----------------------------------------------------------------
88// Assignment operator
89//-------------------------------------------------------------------
90AbcDualRowSteepest &
91AbcDualRowSteepest::operator=(const AbcDualRowSteepest &rhs)
92{
93  if (this != &rhs) {
94    AbcDualRowPivot::operator=(rhs);
95    norm_ = rhs.norm_;
96    factorizationRatio_ = rhs.factorizationRatio_;
97    state_ = rhs.state_;
98    mode_ = rhs.mode_;
99    persistence_ = rhs.persistence_;
100    model_ = rhs.model_;
101    delete weights_;
102    delete infeasible_;
103    delete savedWeights_;
104    assert(model_);
105    int number = model_->numberRows();
106    if (rhs.savedWeights_)
107      number = CoinMin(number, rhs.savedWeights_->capacity());
108    if (rhs.infeasible_ != NULL) {
109      infeasible_ = new CoinIndexedVector(rhs.infeasible_);
110    } else {
111      infeasible_ = NULL;
112    }
113    if (rhs.weights_ != NULL) {
114      weights_ = new CoinIndexedVector(rhs.weights_);
115    } else {
116      weights_ = NULL;
117    }
118    if (rhs.savedWeights_ != NULL) {
119      savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
120    } else {
121      savedWeights_ = NULL;
122    }
123  }
124  return *this;
125}
126// Fill most values
127void AbcDualRowSteepest::fill(const AbcDualRowSteepest &rhs)
128{
129  norm_ = rhs.norm_;
130  factorizationRatio_ = rhs.factorizationRatio_;
131  state_ = rhs.state_;
132  mode_ = rhs.mode_;
133  persistence_ = rhs.persistence_;
134  assert(model_->numberRows() == rhs.model_->numberRows());
135  model_ = rhs.model_;
136  assert(model_);
137  int number = model_->numberRows();
138  if (rhs.savedWeights_)
139    number = CoinMin(number, rhs.savedWeights_->capacity());
140  if (rhs.infeasible_ != NULL) {
141    if (!infeasible_)
142      infeasible_ = new CoinIndexedVector(rhs.infeasible_);
143    else
144      *infeasible_ = *rhs.infeasible_;
145  } else {
146    delete infeasible_;
147    infeasible_ = NULL;
148  }
149  if (rhs.weights_ != NULL) {
150    if (!weights_)
151      weights_ = new CoinIndexedVector(rhs.weights_);
152  } else {
153    delete weights_;
154    weights_ = NULL;
155  }
156  if (rhs.savedWeights_ != NULL) {
157    if (!savedWeights_)
158      savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
159    else
160      *savedWeights_ = *rhs.savedWeights_;
161  } else {
162    delete savedWeights_;
163    savedWeights_ = NULL;
164  }
165}
166#ifdef CHECK_NUMBER_WANTED
167static int xTimes = 0;
168static int xWanted = 0;
169#endif
170#if ABC_PARALLEL == 2
171#define DO_REDUCE 2
172#ifdef DO_REDUCE
173#if DO_REDUCE == 1
174#include <cilk/reducer_max.h>
175static void choose(int &chosenRow, double &largest, int n,
176  const int *index, const double *infeas,
177  const double *weights, double tolerance)
178{
179  cilk::reducer_max_index< int, double > maximumIndex(chosenRow, largest);
180#pragma cilk_grainsize = 128
181  cilk_for(int i = 0; i < n; i++)
182  {
183    int iRow = index[i];
184    double value = infeas[iRow];
185    if (value > tolerance) {
186      double thisWeight = CoinMin(weights[iRow], 1.0e50);
187      maximumIndex.calc_max(iRow, value / thisWeight);
188    }
189  }
190  chosenRow = maximumIndex.get_index();
191  largest = maximumIndex.get_value();
192}
193#else
194static void choose(AbcDualRowSteepest *steepest,
195  int &chosenRowSave, double &largestSave, int first, int last,
196  double tolerance)
197{
198  if (last - first > 256) {
199    int mid = (last + first) >> 1;
200    int chosenRow2 = chosenRowSave;
201    double largest2 = largestSave;
202    cilk_spawn choose(steepest, chosenRow2, largest2, first, mid,
203      tolerance);
204    choose(steepest, chosenRowSave, largestSave, mid, last,
205      tolerance);
206    cilk_sync;
207    if (largest2 > largestSave) {
208      largestSave = largest2;
209      chosenRowSave = chosenRow2;
210    }
211  } else {
212    const int *index = steepest->infeasible()->getIndices();
213    const double *infeas = steepest->infeasible()->denseVector();
214    const double *weights = steepest->weights()->denseVector();
215    double largest = largestSave;
216    int chosenRow = chosenRowSave;
217    if ((steepest->model()->stateOfProblem() & VALUES_PASS2) == 0) {
218      for (int i = first; i < last; i++) {
219        int iRow = index[i];
220        double value = infeas[iRow];
221        if (value > tolerance) {
222          double thisWeight = CoinMin(weights[iRow], 1.0e50);
223          if (value > largest * thisWeight) {
224            largest = value / thisWeight;
225            chosenRow = iRow;
226          }
227        }
228      }
229    } else {
230      const double *fakeDjs = steepest->model()->fakeDjs();
231      const int *pivotVariable = steepest->model()->pivotVariable();
232      for (int i = first; i < last; i++) {
233        int iRow = index[i];
234        double value = infeas[iRow];
235        if (value > tolerance) {
236          int iSequence = pivotVariable[iRow];
237          value *= (fabs(fakeDjs[iSequence]) + 1.0e-6);
238          double thisWeight = CoinMin(weights[iRow], 1.0e50);
239          /*
240            Ideas
241            always use fake
242            use fake until chosen
243            if successful would move djs to basic region?
244            how to switch off if cilking - easiest at saveWeights
245           */
246          if (value > largest * thisWeight) {
247            largest = value / thisWeight;
248            chosenRow = iRow;
249          }
250        }
251      }
252    }
253    chosenRowSave = chosenRow;
254    largestSave = largest;
255  }
256}
257static void choose2(AbcDualRowSteepest *steepest,
258  int &chosenRowSave, double &largestSave, int first, int last,
259  double tolerance)
260{
261  if (last - first > 256) {
262    int mid = (last + first) >> 1;
263    int chosenRow2 = chosenRowSave;
264    double largest2 = largestSave;
265    cilk_spawn choose2(steepest, chosenRow2, largest2, first, mid,
266      tolerance);
267    choose2(steepest, chosenRowSave, largestSave, mid, last,
268      tolerance);
269    cilk_sync;
270    if (largest2 > largestSave) {
271      largestSave = largest2;
272      chosenRowSave = chosenRow2;
273    }
274  } else {
275    const int *index = steepest->infeasible()->getIndices();
276    const double *infeas = steepest->infeasible()->denseVector();
277    double largest = largestSave;
278    int chosenRow = chosenRowSave;
279    for (int i = first; i < last; i++) {
280      int iRow = index[i];
281      double value = infeas[iRow];
282      if (value > largest) {
283        largest = value;
284        chosenRow = iRow;
285      }
286    }
287    chosenRowSave = chosenRow;
288    largestSave = largest;
289  }
290}
291#endif
292#endif
293#endif
294// Returns pivot row, -1 if none
295int AbcDualRowSteepest::pivotRow()
296{
297  assert(model_);
298  double *COIN_RESTRICT infeas = infeasible_->denseVector();
299  int *COIN_RESTRICT index = infeasible_->getIndices();
300  int number = infeasible_->getNumElements();
301  int lastPivotRow = model_->lastPivotRow();
302  assert(lastPivotRow < model_->numberRows());
303  double tolerance = model_->currentPrimalTolerance();
304  // we can't really trust infeasibilities if there is primal error
305  // this coding has to mimic coding in checkPrimalSolution
306  double error = CoinMin(1.0e-2, model_->largestPrimalError());
307  // allow tolerance at least slightly bigger than standard
308  tolerance = tolerance + error;
309  // But cap
310  tolerance = CoinMin(1000.0, tolerance);
311  tolerance *= tolerance; // as we are using squares
312  bool toleranceChanged = false;
313  const double *COIN_RESTRICT solutionBasic = model_->solutionBasic();
314  const double *COIN_RESTRICT lowerBasic = model_->lowerBasic();
315  const double *COIN_RESTRICT upperBasic = model_->upperBasic();
316  // do last pivot row one here
317  if (lastPivotRow >= 0 && lastPivotRow < model_->numberRows()) {
318    double value = solutionBasic[lastPivotRow];
319    double lower = lowerBasic[lastPivotRow];
320    double upper = upperBasic[lastPivotRow];
321    if (value > upper + tolerance) {
322      value -= upper;
323      value *= value;
324      // store square in list
325      if (infeas[lastPivotRow])
326        infeas[lastPivotRow] = value; // already there
327      else
328        infeasible_->quickAdd(lastPivotRow, value);
329    } else if (value < lower - tolerance) {
330      value -= lower;
331      value *= value;
332      // store square in list
333      if (infeas[lastPivotRow])
334        infeas[lastPivotRow] = value; // already there
335      else
336        infeasible_->add(lastPivotRow, value);
337    } else {
338      // feasible - was it infeasible - if so set tiny
339      if (infeas[lastPivotRow])
340        infeas[lastPivotRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
341    }
342    number = infeasible_->getNumElements();
343  }
344  if (model_->numberIterations() < model_->lastBadIteration() + 200) {
345    // we can't really trust infeasibilities if there is dual error
346    if (model_->largestDualError() > model_->largestPrimalError()) {
347      tolerance *= CoinMin(model_->largestDualError() / model_->largestPrimalError(), 1000.0);
348      toleranceChanged = true;
349    }
350  }
351  int numberWanted;
352  if (mode_ < 2) {
353    numberWanted = number + 1;
354  } else {
355    if (model_->factorization()->pivots() == 0) {
356      int numberElements = model_->factorization()->numberElements();
357      int numberRows = model_->numberRows();
358#if 0
359      int numberSlacks = model_->factorization()->numberSlacks();
360      factorizationRatio_ = static_cast<double> (numberElements+1) /
361        static_cast<double> (numberRows+1-numberSlacks);
362#else
363      factorizationRatio_ = static_cast< double >(numberElements) / static_cast< double >(numberRows);
364#endif
365      //printf("%d iterations, numberElements %d ratio %g\n",
366      //     model_->numberIterations(),numberElements,factorizationRatio_);
367      // Bias so less likely to go to steepest if large
368      double weights[] = { 1.0, 2.0, 5.0 };
369      double modifier;
370      if (numberRows < 100000)
371        modifier = weights[0];
372      else if (numberRows < 200000)
373        modifier = weights[1];
374      else
375        modifier = weights[2];
376      if (mode_ == 2 && factorizationRatio_ > modifier) {
377        // switch from dantzig to steepest
378        mode_ = 3;
379        saveWeights(model_, 5);
380      }
381    }
382#if 0
383    // look at this and modify
384    numberWanted = number+1;
385    if (factorizationRatio_ < 1.0) {
386      numberWanted = CoinMax(2000, numberRows / 20);
387#if 0
388    } else if (factorizationRatio_ > 10.0) {
389      double ratio = number * (factorizationRatio_ / 80.0);
390      if (ratio > number)
391        numberWanted = number + 1;
392      else
393        numberWanted = CoinMax(2000, static_cast<int> (ratio));
394#else
395    } else {
396      //double multiplier=CoinMin(factorizationRatio_*0.1,1.0)
397#endif
398    }
399#else
400    numberWanted = CoinMax(2000, number / 8);
401    if (factorizationRatio_ < 1.5) {
402      numberWanted = CoinMax(2000, number / 20);
403    } else if (factorizationRatio_ > 5.0) {
404      numberWanted = number + 1;
405    }
406#endif
407    numberWanted = CoinMin(numberWanted, number + 1);
408  }
409  if (model_->largestPrimalError() > 1.0e-3)
410    numberWanted = number + 1; // be safe
411#ifdef CHECK_NUMBER_WANTED
412  xWanted += numberWanted;
413  xTimes++;
414  if ((xTimes % 1000) == 0)
415    printf("time %d average wanted %g\n", xTimes, static_cast< double >(xWanted) / xTimes);
416#endif
417  int iPass;
418  // Setup two passes
419  int start[4];
420  start[1] = number;
421  start[2] = 0;
422  double dstart = static_cast< double >(number) * model_->randomNumberGenerator()->randomDouble();
423  start[0] = static_cast< int >(dstart);
424  start[3] = start[0];
425  //double largestWeight=0.0;
426  //double smallestWeight=1.0e100;
427  const double *weights = weights_->denseVector();
428  const int *pivotVariable = model_->pivotVariable();
429  double savePivotWeight = weights_->denseVector()[lastPivotRow];
430  weights_->denseVector()[lastPivotRow] *= 1.0e10;
431  double largest = 0.0;
432  int chosenRow = -1;
433  int saveNumberWanted = numberWanted;
434#ifdef DO_REDUCE
435  bool doReduce = true;
436  int lastChosen = -1;
437  double lastLargest = 0.0;
438#endif
439  for (iPass = 0; iPass < 2; iPass++) {
440    int endThis = start[2 * iPass + 1];
441    int startThis = start[2 * iPass];
442    while (startThis < endThis) {
443      int end = CoinMin(endThis, startThis + numberWanted);
444#ifdef DO_REDUCE
445      if (doReduce) {
446#if DO_REDUCE == 1
447        choose(chosenRow, largest, end - startThis,
448          index + startThis, infeas, weights, tolerance);
449#else
450        if (mode_ != 2)
451          choose(this, chosenRow, largest, startThis, end, tolerance);
452        else
453          choose2(this, chosenRow, largest, startThis, end, tolerance);
454#endif
455        if (chosenRow != lastChosen) {
456          assert(chosenRow >= 0);
457          if (model_->flagged(pivotVariable[chosenRow]) || (solutionBasic[chosenRow] <= upperBasic[chosenRow] + tolerance && solutionBasic[chosenRow] >= lowerBasic[chosenRow] - tolerance)) {
458            doReduce = false;
459            chosenRow = lastChosen;
460            largest = lastLargest;
461          } else {
462            lastChosen = chosenRow;
463            lastLargest = largest;
464          }
465        }
466      }
467      if (!doReduce) {
468#endif
469        for (int i = startThis; i < end; i++) {
470          int iRow = index[i];
471          double value = infeas[iRow];
472          if (value > tolerance) {
473            //#define OUT_EQ
474#ifdef OUT_EQ
475            if (upper[iRow] == lower[iRow])
476              value *= 2.0;
477#endif
478            double thisWeight = CoinMin(weights[iRow], 1.0e50);
479            //largestWeight = CoinMax(largestWeight,weight);
480            //smallestWeight = CoinMin(smallestWeight,weight);
481            //double dubious = dubiousWeights_[iRow];
482            //weight *= dubious;
483            //if (value>2.0*largest*weight||(value>0.5*largest*weight&&value*largestWeight>dubious*largest*weight)) {
484            if (value > largest * thisWeight) {
485              // make last pivot row last resort choice
486              //if (iRow == lastPivotRow) {
487              //if (value * 1.0e-10 < largest * thisWeight)
488              //  continue;
489              //else
490              //  value *= 1.0e-10;
491              //}
492              if (!model_->flagged(pivotVariable[iRow])) {
493                //#define CLP_DEBUG 3
494#ifdef CLP_DEBUG
495                double value2 = 0.0;
496                if (solutionBasic[iRow] > upperBasic[iRow] + tolerance)
497                  value2 = solutionBasic[iRow] - upperBasic[iRow];
498                else if (solutionBasic[iRow] < lowerBasic[iRow] - tolerance)
499                  value2 = solutionBasic[iRow] - lowerBasic[iRow];
500                assert(fabs(value2 * value2 - infeas[iRow]) < 1.0e-8 * CoinMin(value2 * value2, infeas[iRow]));
501#endif
502                if (solutionBasic[iRow] > upperBasic[iRow] + tolerance || solutionBasic[iRow] < lowerBasic[iRow] - tolerance) {
503                  chosenRow = iRow;
504                  largest = value / thisWeight;
505                  //largestWeight = dubious;
506                }
507              }
508            }
509          }
510        }
511#ifdef DO_REDUCE
512      }
513#endif
514      numberWanted -= (end - startThis);
515      if (!numberWanted) {
516        if (chosenRow >= 0)
517          break;
518        else
519          numberWanted = (saveNumberWanted + 1) >> 1;
520      }
521      startThis = end;
522    }
523    if (!numberWanted) {
524      if (chosenRow >= 0)
525        break;
526      else
527        numberWanted = (saveNumberWanted + 1) >> 1;
528    }
529  }
530  weights_->denseVector()[lastPivotRow] = savePivotWeight;
531  //printf("smallest %g largest %g\n",smallestWeight,largestWeight);
532  if (chosenRow < 0 && toleranceChanged) {
533    // won't line up with checkPrimalSolution - do again
534    double saveError = model_->largestDualError();
535    model_->setLargestDualError(0.0);
536    // can't loop
537    chosenRow = pivotRow();
538    model_->setLargestDualError(saveError);
539  }
540  if (chosenRow < 0 && lastPivotRow < 0) {
541    int nLeft = 0;
542    for (int i = 0; i < number; i++) {
543      int iRow = index[i];
544      if (fabs(infeas[iRow]) > 1.0e-50) {
545        index[nLeft++] = iRow;
546      } else {
547        infeas[iRow] = 0.0;
548      }
549    }
550    infeasible_->setNumElements(nLeft);
551    model_->setNumberPrimalInfeasibilities(nLeft);
552  }
553  if (true && (model_->stateOfProblem() & VALUES_PASS2) != 0) {
554    if (chosenRow >= 0) {
555      double *fakeDjs = model_->fakeDjs();
556      //const int * pivotVariable = model_->pivotVariable();
557      fakeDjs[pivotVariable[chosenRow]] = 0.0;
558    }
559  }
560  return chosenRow;
561}
562//#define PRINT_WEIGHTS_FREQUENCY
563#ifdef PRINT_WEIGHTS_FREQUENCY
564int xweights1 = 0;
565int xweights2 = 0;
566int xweights3 = 0;
567#endif
568/* Updates weights and returns pivot alpha.
569   Also does FT update */
570double
571AbcDualRowSteepest::updateWeights(CoinIndexedVector &input,
572  CoinIndexedVector &updateColumn)
573{
574  double *COIN_RESTRICT weights = weights_->denseVector();
575#if CLP_DEBUG > 2
576  // Very expensive debug
577  {
578    int numberRows = model_->numberRows();
579    CoinIndexedVector temp;
580    temp.reserve(numberRows);
581    double *array = temp.denseVector();
582    int *which = temp.getIndices();
583    for (int i = 0; i < numberRows; i++) {
584      double value = 0.0;
585      array[i] = 1.0;
586      which[0] = i;
587      temp.setNumElements(1);
588      model_->factorization()->updateColumnTranspose(temp);
589      int number = temp.getNumElements();
590      for (int j = 0; j < number; j++) {
591        int iRow = which[j];
592        value += array[iRow] * array[iRow];
593        array[iRow] = 0.0;
594      }
595      temp.setNumElements(0);
596      double w = CoinMax(weights[i], value) * .1;
597      if (fabs(weights[i] - value) > w) {
598        printf("%d old %g, true %g\n", i, weights[i], value);
599        weights[i] = value; // to reduce printout
600      }
601    }
602  }
603#endif
604  double alpha = 0.0;
605  double norm = 0.0;
606  int i;
607  double *COIN_RESTRICT work2 = input.denseVector();
608  int numberNonZero = input.getNumElements();
609  int *COIN_RESTRICT which = input.getIndices();
610
611  //compute norm
612  for (int i = 0; i < numberNonZero; i++) {
613    int iRow = which[i];
614    double value = work2[iRow];
615    norm += value * value;
616  }
617  // Do FT update
618  if (true) {
619    model_->factorization()->updateTwoColumnsFT(updateColumn, input);
620#if CLP_DEBUG > 3
621    printf("REGION after %d els\n", input.getNumElements());
622    input.print();
623#endif
624  } else {
625    // Leave as old way
626    model_->factorization()->updateColumnFT(updateColumn);
627    model_->factorization()->updateColumn(input);
628  }
629  //#undef CLP_DEBUG
630  numberNonZero = input.getNumElements();
631  int pivotRow = model_->lastPivotRow();
632  assert(model_->pivotRow() == model_->lastPivotRow());
633#ifdef CLP_DEBUG
634  if (model_->logLevel() > 4 && fabs(norm - weights[pivotRow]) > 1.0e-2 * (1.0 + norm))
635    printf("on row %d, true weight %g, old %g\n",
636      pivotRow, sqrt(norm), sqrt(weights[pivotRow]));
637#endif
638  // could re-initialize here (could be expensive)
639  norm /= model_->alpha() * model_->alpha();
640  assert(model_->alpha());
641  assert(norm);
642  double multiplier = 2.0 / model_->alpha();
643#ifdef PRINT_WEIGHTS_FREQUENCY
644  xweights1++;
645  if ((xweights1 + xweights2 + xweights3) % 100 == 0) {
646    printf("ZZ iteration %d sum weights %d - %d %d %d\n",
647      model_->numberIterations(), xweights1 + xweights2 + xweights3, xweights1, xweights2, xweights3);
648    assert(abs(model_->numberIterations() - (xweights1 + xweights2 + xweights3)) <= 1);
649  }
650#endif
651  // look at updated column
652  double *COIN_RESTRICT work = updateColumn.denseVector();
653  numberNonZero = updateColumn.getNumElements();
654  which = updateColumn.getIndices();
655
656  // pivot element
657  alpha = work[pivotRow];
658  for (i = 0; i < numberNonZero; i++) {
659    int iRow = which[i];
660    double theta = work[iRow];
661    double devex = weights[iRow];
662    //work3[iRow] = devex; // save old
663    //which3[nSave++] = iRow;
664    double value = work2[iRow];
665    devex += theta * (theta * norm + value * multiplier);
666    if (devex < DEVEX_TRY_NORM)
667      devex = DEVEX_TRY_NORM;
668    weights[iRow] = devex;
669  }
670  if (norm < DEVEX_TRY_NORM)
671    norm = DEVEX_TRY_NORM;
672  // Try this to make less likely will happen again and stop cycling
673  //norm *= 1.02;
674  weights[pivotRow] = norm;
675  input.clear();
676#ifdef CLP_DEBUG
677  input.checkClear();
678#endif
679#ifdef CLP_DEBUG
680  input.checkClear();
681#endif
682  return alpha;
683}
684double
685AbcDualRowSteepest::updateWeights1(CoinIndexedVector &input, CoinIndexedVector &updateColumn)
686{
687  if (mode_ == 2) {
688    abort();
689    // no update
690    // Do FT update
691    model_->factorization()->updateColumnFT(updateColumn);
692    // pivot element
693    double alpha = 0.0;
694    // look at updated column
695    double *work = updateColumn.denseVector();
696    int pivotRow = model_->lastPivotRow();
697    assert(pivotRow == model_->pivotRow());
698
699    assert(!updateColumn.packedMode());
700    alpha = work[pivotRow];
701    return alpha;
702  }
703#if CLP_DEBUG > 2
704  // Very expensive debug
705  {
706    double *weights = weights_->denseVector();
707    int numberRows = model_->numberRows();
708    const int *pivotVariable = model_->pivotVariable();
709    CoinIndexedVector temp;
710    temp.reserve(numberRows);
711    double *array = temp.denseVector();
712    int *which = temp.getIndices();
713    for (int i = 0; i < numberRows; i++) {
714      double value = 0.0;
715      array[i] = 1.0;
716      which[0] = i;
717      temp.setNumElements(1);
718      model_->factorization()->updateColumnTranspose(temp);
719      int number = temp.getNumElements();
720      for (int j = 0; j < number; j++) {
721        int iRow = which[j];
722        value += array[iRow] * array[iRow];
723        array[iRow] = 0.0;
724      }
725      temp.setNumElements(0);
726      double w = CoinMax(weights[i], value) * .1;
727      if (fabs(weights[i] - value) > w) {
728        printf("XX row %d (variable %d) old %g, true %g\n", i, pivotVariable[i],
729          weights[i], value);
730        weights[i] = value; // to reduce printout
731      }
732    }
733  }
734#endif
735  norm_ = 0.0;
736  double *COIN_RESTRICT work2 = input.denseVector();
737  int numberNonZero = input.getNumElements();
738  int *COIN_RESTRICT which = input.getIndices();
739
740  //compute norm
741  for (int i = 0; i < numberNonZero; i++) {
742    int iRow = which[i];
743    double value = work2[iRow];
744    norm_ += value * value;
745  }
746  // Do FT update
747  if (true) {
748    model_->factorization()->updateTwoColumnsFT(updateColumn, input);
749#if CLP_DEBUG > 3
750    printf("REGION after %d els\n", input.getNumElements());
751    input.print();
752#endif
753  } else {
754    // Leave as old way
755    model_->factorization()->updateColumnFT(updateColumn);
756    model_->factorization()->updateWeights(input);
757  }
758  // pivot element
759  assert(model_->pivotRow() == model_->lastPivotRow());
760  return updateColumn.denseVector()[model_->lastPivotRow()];
761}
762void AbcDualRowSteepest::updateWeightsOnly(CoinIndexedVector &input)
763{
764  if (mode_ == 2)
765    return;
766  norm_ = 0.0;
767  double *COIN_RESTRICT work2 = input.denseVector();
768  int numberNonZero = input.getNumElements();
769  int *COIN_RESTRICT which = input.getIndices();
770
771  //compute norm
772  for (int i = 0; i < numberNonZero; i++) {
773    int iRow = which[i];
774    double value = work2[iRow];
775    norm_ += value * value;
776  }
777  model_->factorization()->updateWeights(input);
778}
779// Actually updates weights
780void AbcDualRowSteepest::updateWeights2(CoinIndexedVector &input, CoinIndexedVector &updateColumn)
781{
782  if (mode_ == 2)
783    return;
784  double *COIN_RESTRICT weights = weights_->denseVector();
785  const double *COIN_RESTRICT work2 = input.denseVector();
786  int pivotRow = model_->lastPivotRow();
787  assert(model_->pivotRow() == model_->lastPivotRow());
788#ifdef CLP_DEBUG
789  if (model_->logLevel() > 4 && fabs(norm_ - weights[pivotRow]) > 1.0e-2 * (1.0 + norm_))
790    printf("on row %d, true weight %g, old %g\n",
791      pivotRow, sqrt(norm_), sqrt(weights[pivotRow]));
792#endif
793  // pivot element
794  double alpha = model_->alpha();
795  // could re-initialize here (could be expensive)
796  norm_ /= alpha * alpha;
797  assert(alpha);
798  assert(norm_);
799  double multiplier = 2.0 / alpha;
800  // look at updated column
801  const double *COIN_RESTRICT work = updateColumn.denseVector();
802  int numberNonZero = updateColumn.getNumElements();
803  const int *which = updateColumn.getIndices();
804  double multiplier2 = 1.0;
805#ifdef PRINT_WEIGHTS_FREQUENCY
806  xweights2++;
807  if ((xweights1 + xweights2 + xweights3) % 100 == 0) {
808    printf("ZZ iteration %d sum weights %d - %d %d %d\n",
809      model_->numberIterations(), xweights1 + xweights2 + xweights3, xweights1, xweights2, xweights3);
810    assert(abs(model_->numberIterations() - (xweights1 + xweights2 + xweights3)) <= 1);
811  }
812#endif
813  if (model_->directionOut() < 0) {
814#ifdef ABC_DEBUG
815    printf("Minus out in %d %d, out %d %d\n", model_->sequenceIn(), model_->directionIn(),
816      model_->sequenceOut(), model_->directionOut());
817#endif
818    multiplier2 = -1.0;
819  }
820  for (int i = 0; i < numberNonZero; i++) {
821    int iRow = which[i];
822    double theta = multiplier2 * work[iRow];
823    double devex = weights[iRow];
824    double value = work2[iRow];
825    devex += theta * (theta * norm_ + value * multiplier);
826    if (devex < DEVEX_TRY_NORM)
827      devex = DEVEX_TRY_NORM;
828    weights[iRow] = devex;
829  }
830  if (norm_ < DEVEX_TRY_NORM)
831    norm_ = DEVEX_TRY_NORM;
832  // Try this to make less likely will happen again and stop cycling
833  //norm_ *= 1.02;
834  weights[pivotRow] = norm_;
835  input.clear();
836#ifdef CLP_DEBUG
837  input.checkClear();
838#endif
839}
840
841/* Updates primal solution (and maybe list of candidates)
842   Uses input vector which it deletes
843   Computes change in objective function
844*/
845void AbcDualRowSteepest::updatePrimalSolution(
846  CoinIndexedVector &primalUpdate,
847  double primalRatio)
848{
849  double *COIN_RESTRICT work = primalUpdate.denseVector();
850  int number = primalUpdate.getNumElements();
851  int *COIN_RESTRICT which = primalUpdate.getIndices();
852  double tolerance = model_->currentPrimalTolerance();
853  double *COIN_RESTRICT infeas = infeasible_->denseVector();
854  double *COIN_RESTRICT solutionBasic = model_->solutionBasic();
855  const double *COIN_RESTRICT lowerBasic = model_->lowerBasic();
856  const double *COIN_RESTRICT upperBasic = model_->upperBasic();
857  assert(!primalUpdate.packedMode());
858  for (int i = 0; i < number; i++) {
859    int iRow = which[i];
860    double value = solutionBasic[iRow];
861    double change = primalRatio * work[iRow];
862    // no need to keep for steepest edge
863    work[iRow] = 0.0;
864    value -= change;
865    double lower = lowerBasic[iRow];
866    double upper = upperBasic[iRow];
867    solutionBasic[iRow] = value;
868    if (value < lower - tolerance) {
869      value -= lower;
870      value *= value;
871#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
872      if (lower == upper)
873        value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
874#endif
875      // store square in list
876      if (infeas[iRow])
877        infeas[iRow] = value; // already there
878      else
879        infeasible_->quickAdd(iRow, value);
880    } else if (value > upper + tolerance) {
881      value -= upper;
882      value *= value;
883#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
884      if (lower == upper)
885        value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
886#endif
887      // store square in list
888      if (infeas[iRow])
889        infeas[iRow] = value; // already there
890      else
891        infeasible_->quickAdd(iRow, value);
892    } else {
893      // feasible - was it infeasible - if so set tiny
894      if (infeas[iRow])
895        infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
896    }
897  }
898  // Do pivot row
899  int iRow = model_->lastPivotRow();
900  assert(model_->pivotRow() == model_->lastPivotRow());
901  // feasible - was it infeasible - if so set tiny
902  //assert (infeas[iRow]);
903  if (infeas[iRow])
904    infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
905  primalUpdate.setNumElements(0);
906}
907#if ABC_PARALLEL == 2
908static void update(int first, int last,
909  const int *COIN_RESTRICT which, double *COIN_RESTRICT work,
910  const double *COIN_RESTRICT work2, double *COIN_RESTRICT weights,
911  const double *COIN_RESTRICT lowerBasic, double *COIN_RESTRICT solutionBasic,
912  const double *COIN_RESTRICT upperBasic,
913  double multiplier, double multiplier2,
914  double norm, double theta, double tolerance)
915{
916  if (last - first > 256) {
917    int mid = (last + first) >> 1;
918    cilk_spawn update(first, mid, which, work, work2, weights, lowerBasic, solutionBasic,
919      upperBasic, multiplier, multiplier2, norm, theta, tolerance);
920    update(mid, last, which, work, work2, weights, lowerBasic, solutionBasic,
921      upperBasic, multiplier, multiplier2, norm, theta, tolerance);
922    cilk_sync;
923  } else {
924    for (int i = first; i < last; i++) {
925      int iRow = which[i];
926      double updateValue = work[iRow];
927      double thetaDevex = multiplier2 * updateValue;
928      double devex = weights[iRow];
929      double valueDevex = work2[iRow];
930      devex += thetaDevex * (thetaDevex * norm + valueDevex * multiplier);
931      if (devex < DEVEX_TRY_NORM)
932        devex = DEVEX_TRY_NORM;
933      weights[iRow] = devex;
934      double value = solutionBasic[iRow];
935      double change = theta * updateValue;
936      value -= change;
937      double lower = lowerBasic[iRow];
938      double upper = upperBasic[iRow];
939      solutionBasic[iRow] = value;
940      if (value < lower - tolerance) {
941        value -= lower;
942        value *= value;
943#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
944        if (lower == upper)
945          value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
946#endif
947      } else if (value > upper + tolerance) {
948        value -= upper;
949        value *= value;
950#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
951        if (lower == upper)
952          value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
953#endif
954      } else {
955        // feasible
956        value = 0.0;
957      }
958      // store square
959      work[iRow] = value;
960    }
961  }
962}
963#endif
964void AbcDualRowSteepest::updatePrimalSolutionAndWeights(CoinIndexedVector &weightsVector,
965  CoinIndexedVector &primalUpdate,
966  double theta)
967{
968  //printf("updatePrimal iteration %d - weightsVector has %d, primalUpdate has %d\n",
969  //     model_->numberIterations(),weightsVector.getNumElements(),primalUpdate.getNumElements());
970  double *COIN_RESTRICT weights = weights_->denseVector();
971  double *COIN_RESTRICT work2 = weightsVector.denseVector();
972  int pivotRow = model_->lastPivotRow();
973  assert(model_->pivotRow() == model_->lastPivotRow());
974  double *COIN_RESTRICT work = primalUpdate.denseVector();
975  int numberNonZero = primalUpdate.getNumElements();
976  int *COIN_RESTRICT which = primalUpdate.getIndices();
977  double tolerance = model_->currentPrimalTolerance();
978  double *COIN_RESTRICT infeas = infeasible_->denseVector();
979  double *COIN_RESTRICT solutionBasic = model_->solutionBasic();
980  const double *COIN_RESTRICT lowerBasic = model_->lowerBasic();
981  const double *COIN_RESTRICT upperBasic = model_->upperBasic();
982  assert(!primalUpdate.packedMode());
983#ifdef CLP_DEBUG
984  if (model_->logLevel() > 4 && fabs(norm_ - weights[pivotRow]) > 1.0e-2 * (1.0 + norm_))
985    printf("on row %d, true weight %g, old %g\n",
986      pivotRow, sqrt(norm_), sqrt(weights[pivotRow]));
987#endif
988  // pivot element
989  double alpha = model_->alpha();
990  // could re-initialize here (could be expensive)
991  norm_ /= alpha * alpha;
992  assert(alpha);
993  assert(norm_);
994  double multiplier = 2.0 / alpha;
995#ifdef PRINT_WEIGHTS_FREQUENCY
996  // only weights3 is used? - slightly faster to just update weights if going to invert
997  xweights3++;
998  if ((xweights1 + xweights2 + xweights3) % 1000 == 0) {
999    printf("ZZ iteration %d sum weights %d - %d %d %d\n",
1000      model_->numberIterations(), xweights1 + xweights2 + xweights3, xweights1, xweights2, xweights3);
1001    assert(abs(model_->numberIterations() - (xweights1 + xweights2 + xweights3)) <= 1);
1002  }
1003#endif
1004  // look at updated column
1005  double multiplier2 = 1.0;
1006  if (model_->directionOut() < 0) {
1007#ifdef ABC_DEBUG
1008    printf("Minus out in %d %d, out %d %d\n", model_->sequenceIn(), model_->directionIn(),
1009      model_->sequenceOut(), model_->directionOut());
1010#endif
1011    multiplier2 = -1.0;
1012  }
1013#if ABC_PARALLEL == 2
1014  update(0, numberNonZero, which, work, work2, weights,
1015    lowerBasic, solutionBasic, upperBasic,
1016    multiplier, multiplier2, norm_, theta, tolerance);
1017  for (int i = 0; i < numberNonZero; i++) {
1018    int iRow = which[i];
1019    double infeasibility = work[iRow];
1020    work[iRow] = 0.0;
1021    if (infeasibility) {
1022      if (infeas[iRow])
1023        infeas[iRow] = infeasibility; // already there
1024      else
1025        infeasible_->quickAdd(iRow, infeasibility);
1026    } else {
1027      // feasible - was it infeasible - if so set tiny
1028      if (infeas[iRow])
1029        infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
1030    }
1031  }
1032#else
1033  for (int i = 0; i < numberNonZero; i++) {
1034    int iRow = which[i];
1035    double updateValue = work[iRow];
1036    work[iRow] = 0.0;
1037    double thetaDevex = multiplier2 * updateValue;
1038    double devex = weights[iRow];
1039    double valueDevex = work2[iRow];
1040    devex += thetaDevex * (thetaDevex * norm_ + valueDevex * multiplier);
1041    if (devex < DEVEX_TRY_NORM)
1042      devex = DEVEX_TRY_NORM;
1043    weights[iRow] = devex;
1044    double value = solutionBasic[iRow];
1045    double change = theta * updateValue;
1046    value -= change;
1047    double lower = lowerBasic[iRow];
1048    double upper = upperBasic[iRow];
1049    solutionBasic[iRow] = value;
1050    if (value < lower - tolerance) {
1051      value -= lower;
1052      value *= value;
1053#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
1054      if (lower == upper)
1055        value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
1056#endif
1057      // store square in list
1058      if (infeas[iRow])
1059        infeas[iRow] = value; // already there
1060      else
1061        infeasible_->quickAdd(iRow, value);
1062    } else if (value > upper + tolerance) {
1063      value -= upper;
1064      value *= value;
1065#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
1066      if (lower == upper)
1067        value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
1068#endif
1069      // store square in list
1070      if (infeas[iRow])
1071        infeas[iRow] = value; // already there
1072      else
1073        infeasible_->quickAdd(iRow, value);
1074    } else {
1075      // feasible - was it infeasible - if so set tiny
1076      if (infeas[iRow])
1077        infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
1078    }
1079  }
1080#endif
1081  weightsVector.clear();
1082  primalUpdate.setNumElements(0);
1083#ifdef CLP_DEBUG
1084  weightsVector.checkClear();
1085#endif
1086  // Do pivot row
1087  if (norm_ < DEVEX_TRY_NORM)
1088    norm_ = DEVEX_TRY_NORM;
1089  // Try this to make less likely will happen again and stop cycling
1090  //norm_ *= 1.02;
1091  weights[pivotRow] = norm_;
1092  // feasible - was it infeasible - if so set tiny
1093  //assert (infeas[iRow]);
1094  if (infeas[pivotRow])
1095    infeas[pivotRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
1096  // see 1/2 above primalUpdate.setNumElements(0);
1097}
1098extern void parallelDual5(AbcSimplexFactorization *factorization,
1099  CoinIndexedVector **whichVector,
1100  int numberCpu,
1101  int whichCpu,
1102  double *weights);
1103/* Saves any weights round factorization as pivot rows may change
1104   1) before factorization
1105   2) after factorization
1106   3) just redo infeasibilities
1107   4) restore weights
1108*/
1109void AbcDualRowSteepest::saveWeights(AbcSimplex *model, int mode)
1110{
1111  model_ = model;
1112  int numberRows = model_->numberRows();
1113  // faster to save indices in array numberTotal in length
1114  const int *pivotVariable = model_->pivotVariable();
1115  double *weights = weights_ ? weights_->denseVector() : NULL;
1116  // see if we need to switch
1117  //if (mode_==2&&mode==2) {
1118  // switch
1119  //mode=5;
1120  //}
1121  if (mode == 1) {
1122    if (weights_) {
1123      // Check if size has changed
1124      if (infeasible_->capacity() == numberRows) {
1125        // change from row numbers to sequence numbers
1126        CoinAbcMemcpy(weights_->getIndices(), pivotVariable, numberRows);
1127        state_ = 1;
1128      } else {
1129        // size has changed - clear everything
1130        delete weights_;
1131        weights_ = NULL;
1132        delete infeasible_;
1133        infeasible_ = NULL;
1134        delete savedWeights_;
1135        savedWeights_ = NULL;
1136        state_ = -1;
1137      }
1138    }
1139  } else if (mode != 3) {
1140    // restore
1141    if (!weights_ || state_ == -1 || mode == 5) {
1142      // see if crash basis???
1143      int numberBasicSlacks = 0;
1144      bool initializeNorms = false;
1145      for (int i = 0; i < numberRows; i++) {
1146        if (model_->getInternalStatus(i) == AbcSimplex::basic)
1147          numberBasicSlacks++;
1148      }
1149      if (numberBasicSlacks < numberRows && !model_->numberIterations()) {
1150        char line[100];
1151        if (numberBasicSlacks > numberRows - (numberRows >> 1) && mode_ == 3 && (model_->moreSpecialOptions() & 256) == 0) {
1152          sprintf(line, "Steep initial basis has %d structurals out of %d - initializing norms\n", numberRows - numberBasicSlacks, numberRows);
1153          initializeNorms = true;
1154        } else {
1155          sprintf(line, "Steep initial basis has %d structurals out of %d - too many\n", numberRows - numberBasicSlacks, numberRows);
1156        }
1157        model_->messageHandler()->message(CLP_GENERAL2, *model_->messagesPointer())
1158          << line << CoinMessageEol;
1159      }
1160      // initialize weights
1161      delete weights_;
1162      delete savedWeights_;
1163      delete infeasible_;
1164      weights_ = new CoinIndexedVector();
1165      weights_->reserve(numberRows);
1166      savedWeights_ = new CoinIndexedVector();
1167      // allow for dense multiple of 8!
1168      savedWeights_->reserve(numberRows + ((initializeNorms) ? 8 : 0));
1169      infeasible_ = new CoinIndexedVector();
1170      infeasible_->reserve(numberRows);
1171      weights = weights_->denseVector();
1172      if (mode == 5 || (mode_ != 1 && !initializeNorms)) {
1173        // initialize to 1.0 (can we do better?)
1174        for (int i = 0; i < numberRows; i++) {
1175          weights[i] = 1.0;
1176        }
1177      } else {
1178        double *COIN_RESTRICT array = savedWeights_->denseVector();
1179        int *COIN_RESTRICT which = savedWeights_->getIndices();
1180        // Later look at factorization to see which can be skipped
1181        // also use threads to parallelize
1182#if 0
1183        // use infeasible to collect row counts
1184        int * COIN_RESTRICT counts = infeasible_->getIndices();
1185        CoinAbcMemset0(counts,numberRows);
1186        // get matrix data pointers
1187        int maximumRows = model_->maximumAbcNumberRows();
1188        const CoinBigIndex * columnStart = model_->abcMatrix()->getPackedMatrix()->getVectorStarts()-maximumRows;
1189        //const int * columnLength = model_->abcMatrix()->getPackedMatrix()->getVectorLengths()-maximumRows;
1190        const int * row = model_->abcMatrix()->getPackedMatrix()->getIndices();
1191        const int * pivotVariable = model_->pivotVariable();
1192        for (int i = 0; i < numberRows; i++) {
1193          int iSequence = pivotVariable[i];
1194          if (iSequence>=numberRows) {
1195            for (CoinBigIndex j=columnStart[iSequence];j<columnStart[iSequence+1];j++) {
1196              int iRow=row[j];
1197              counts[iRow]++;
1198            }
1199          }
1200        }
1201        int nSlack=0;
1202#endif
1203        //#undef ABC_PARALLEL
1204        // for now just with cilk
1205#if ABC_PARALLEL == 2
1206        if (model_->parallelMode() == 0) {
1207#endif
1208          for (int i = 0; i < numberRows; i++) {
1209            double value = 0.0;
1210            array[i] = 1.0;
1211            which[0] = i;
1212            savedWeights_->setNumElements(1);
1213            model_->factorization()->updateColumnTranspose(*savedWeights_);
1214            int number = savedWeights_->getNumElements();
1215            for (int j = 0; j < number; j++) {
1216              int k = which[j];
1217              value += array[k] * array[k];
1218              array[k] = 0.0;
1219            }
1220#if 0
1221          if (!counts[i]) {
1222            assert (value==1.0);
1223            nSlack++;
1224          }
1225#endif
1226            savedWeights_->setNumElements(0);
1227            weights[i] = value;
1228          }
1229#if ABC_PARALLEL == 2
1230        } else {
1231          // parallel
1232          // get arrays
1233          int numberCpuMinusOne = CoinMin(model_->parallelMode(), 3);
1234          int which[3];
1235          CoinIndexedVector *whichVector[4];
1236          for (int i = 0; i < numberCpuMinusOne; i++) {
1237            which[i] = model_->getAvailableArray();
1238            whichVector[i] = model_->usefulArray(which[i]);
1239            assert(!whichVector[i]->packedMode());
1240          }
1241          whichVector[numberCpuMinusOne] = savedWeights_;
1242          parallelDual5(model_->factorization(), whichVector, numberCpuMinusOne + 1, numberCpuMinusOne, weights);
1243          for (int i = 0; i < numberCpuMinusOne; i++)
1244            model_->clearArrays(which[i]);
1245        }
1246#endif
1247#if 0
1248        printf("Steep %d can be skipped\n",nSlack);
1249#endif
1250      }
1251    } else {
1252      int useArray = model_->getAvailableArrayPublic();
1253      CoinIndexedVector *usefulArray = model_->usefulArray(useArray);
1254      int *COIN_RESTRICT which = usefulArray->getIndices();
1255      double *COIN_RESTRICT element = usefulArray->denseVector();
1256      if (mode != 4) {
1257        CoinAbcMemcpy(which, weights_->getIndices(), numberRows);
1258        CoinAbcScatterTo(weights_->denseVector(), element, which, numberRows);
1259      } else {
1260        // restore weights
1261        CoinAbcMemcpy(which, savedWeights_->getIndices(), numberRows);
1262        CoinAbcScatterTo(savedWeights_->denseVector(), element, which, numberRows);
1263      }
1264      for (int i = 0; i < numberRows; i++) {
1265        int iPivot = pivotVariable[i];
1266        double thisWeight = element[iPivot];
1267        if (thisWeight) {
1268          if (thisWeight < DEVEX_TRY_NORM)
1269            thisWeight = DEVEX_TRY_NORM; // may need to check more
1270          weights[i] = thisWeight;
1271        } else {
1272          // odd
1273          weights[i] = 1.0;
1274        }
1275      }
1276      if (0) {
1277        double w[1000];
1278        int s[1000];
1279        assert(numberRows <= 1000);
1280        memcpy(w, weights, numberRows * sizeof(double));
1281        memcpy(s, pivotVariable, numberRows * sizeof(int));
1282        CoinSort_2(s, s + numberRows, w);
1283        printf("Weights===========\n");
1284        for (int i = 0; i < numberRows; i++) {
1285          printf("%d seq %d weight %g\n", i, s[i], w[i]);
1286        }
1287        printf("End weights===========\n");
1288      }
1289      usefulArray->setNumElements(numberRows);
1290      usefulArray->clear();
1291      // won't be done in parallel
1292      model_->clearArraysPublic(useArray);
1293    }
1294    state_ = 0;
1295    CoinAbcMemcpy(savedWeights_->denseVector(), weights_->denseVector(), numberRows);
1296    CoinAbcMemcpy(savedWeights_->getIndices(), pivotVariable, numberRows);
1297  }
1298  if (mode >= 2) {
1299    infeasible_->clear();
1300    double tolerance = model_->currentPrimalTolerance();
1301    const double *COIN_RESTRICT solutionBasic = model_->solutionBasic();
1302    const double *COIN_RESTRICT lowerBasic = model_->lowerBasic();
1303    const double *COIN_RESTRICT upperBasic = model_->upperBasic();
1304    if ((model_->stateOfProblem() & VALUES_PASS2) == 0) {
1305      for (int iRow = 0; iRow < numberRows; iRow++) {
1306        double value = solutionBasic[iRow];
1307        double lower = lowerBasic[iRow];
1308        double upper = upperBasic[iRow];
1309        if (value < lower - tolerance) {
1310          value -= lower;
1311          value *= value;
1312#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
1313          if (lower == upper)
1314            value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
1315#endif
1316          // store square in list
1317          infeasible_->quickAdd(iRow, value);
1318        } else if (value > upper + tolerance) {
1319          value -= upper;
1320          value *= value;
1321#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
1322          if (lower == upper)
1323            value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
1324#endif
1325          // store square in list
1326          infeasible_->quickAdd(iRow, value);
1327        }
1328      }
1329    } else {
1330      // check values pass
1331      const double *fakeDjs = model_->fakeDjs();
1332      int numberFake = 0;
1333      for (int iRow = 0; iRow < numberRows; iRow++) {
1334        double value = solutionBasic[iRow];
1335        double lower = lowerBasic[iRow];
1336        double upper = upperBasic[iRow];
1337        if (value < lower - tolerance) {
1338          value -= lower;
1339        } else if (value > upper + tolerance) {
1340          value -= upper;
1341        } else {
1342          value = 0.0;
1343        }
1344        if (value) {
1345          value *= value;
1346#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
1347          if (lower == upper)
1348            value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
1349#endif
1350          // store square in list
1351          infeasible_->quickAdd(iRow, value);
1352          int iSequence = pivotVariable[iRow];
1353          double fakeDj = fakeDjs[iSequence];
1354          if (fabs(fakeDj) > 10.0 * tolerance)
1355            numberFake++;
1356        }
1357      }
1358      if (numberFake)
1359        printf("%d fake basic djs\n", numberFake);
1360      else
1361        model_->setStateOfProblem(model_->stateOfProblem() & ~VALUES_PASS2);
1362    }
1363  }
1364}
1365// Recompute infeasibilities
1366void AbcDualRowSteepest::recomputeInfeasibilities()
1367{
1368  int numberRows = model_->numberRows();
1369  infeasible_->clear();
1370  double tolerance = model_->currentPrimalTolerance();
1371  const double *COIN_RESTRICT solutionBasic = model_->solutionBasic();
1372  const double *COIN_RESTRICT lowerBasic = model_->lowerBasic();
1373  const double *COIN_RESTRICT upperBasic = model_->upperBasic();
1374  for (int iRow = 0; iRow < numberRows; iRow++) {
1375    double value = solutionBasic[iRow];
1376    double lower = lowerBasic[iRow];
1377    double upper = upperBasic[iRow];
1378    if (value < lower - tolerance) {
1379      value -= lower;
1380      value *= value;
1381#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
1382      if (lower == upper)
1383        value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
1384#endif
1385      // store square in list
1386      infeasible_->quickAdd(iRow, value);
1387    } else if (value > upper + tolerance) {
1388      value -= upper;
1389      value *= value;
1390#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
1391      if (lower == upper)
1392        value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
1393#endif
1394      // store square in list
1395      infeasible_->quickAdd(iRow, value);
1396    }
1397  }
1398}
1399//-------------------------------------------------------------------
1400// Clone
1401//-------------------------------------------------------------------
1402AbcDualRowPivot *AbcDualRowSteepest::clone(bool CopyData) const
1403{
1404  if (CopyData) {
1405    return new AbcDualRowSteepest(*this);
1406  } else {
1407    return new AbcDualRowSteepest();
1408  }
1409}
1410// Gets rid of all arrays
1411void AbcDualRowSteepest::clearArrays()
1412{
1413  if (persistence_ == normal) {
1414    delete weights_;
1415    weights_ = NULL;
1416    delete infeasible_;
1417    infeasible_ = NULL;
1418    delete savedWeights_;
1419    savedWeights_ = NULL;
1420  }
1421  state_ = -1;
1422}
1423// Returns true if would not find any row
1424bool AbcDualRowSteepest::looksOptimal() const
1425{
1426  int iRow;
1427  double tolerance = model_->currentPrimalTolerance();
1428  // we can't really trust infeasibilities if there is primal error
1429  // this coding has to mimic coding in checkPrimalSolution
1430  double error = CoinMin(1.0e-2, model_->largestPrimalError());
1431  // allow tolerance at least slightly bigger than standard
1432  tolerance = tolerance + error;
1433  // But cap
1434  tolerance = CoinMin(1000.0, tolerance);
1435  int numberRows = model_->numberRows();
1436  int numberInfeasible = 0;
1437  const double *COIN_RESTRICT solutionBasic = model_->solutionBasic();
1438  const double *COIN_RESTRICT lowerBasic = model_->lowerBasic();
1439  const double *COIN_RESTRICT upperBasic = model_->upperBasic();
1440  for (iRow = 0; iRow < numberRows; iRow++) {
1441    double value = solutionBasic[iRow];
1442    double lower = lowerBasic[iRow];
1443    double upper = upperBasic[iRow];
1444    if (value < lower - tolerance) {
1445      numberInfeasible++;
1446    } else if (value > upper + tolerance) {
1447      numberInfeasible++;
1448    }
1449  }
1450  return (numberInfeasible == 0);
1451}
1452
1453/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
1454*/
Note: See TracBrowser for help on using the repository browser.