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

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