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

Last change on this file since 1879 was 1878, checked in by forrest, 7 years ago

minor changes to implement Aboca

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