source: stable/1.8/Clp/src/ClpDualRowSteepest.cpp @ 1260

Last change on this file since 1260 was 1260, checked in by forrest, 12 years ago

fix tolerance bug?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 31.8 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5#include "ClpSimplex.hpp"
6#include "ClpDualRowSteepest.hpp"
7#include "CoinIndexedVector.hpp"
8#include "ClpFactorization.hpp"
9#include "CoinHelperFunctions.hpp"
10#include <cstdio>
11//#############################################################################
12// Constructors / Destructor / Assignment
13//#############################################################################
14//#define CLP_DEBUG 4
15//-------------------------------------------------------------------
16// Default Constructor
17//-------------------------------------------------------------------
18ClpDualRowSteepest::ClpDualRowSteepest (int mode) 
19  : ClpDualRowPivot(),
20    state_(-1),
21    mode_(mode),
22    persistence_(normal),
23    weights_(NULL),
24    infeasible_(NULL),
25    alternateWeights_(NULL),
26    savedWeights_(NULL),
27    dubiousWeights_(NULL)
28{
29  type_=2+64*mode;
30}
31
32//-------------------------------------------------------------------
33// Copy constructor
34//-------------------------------------------------------------------
35ClpDualRowSteepest::ClpDualRowSteepest (const ClpDualRowSteepest & rhs) 
36: ClpDualRowPivot(rhs)
37{ 
38  state_=rhs.state_;
39  mode_ = rhs.mode_;
40  persistence_ = rhs.persistence_;
41  model_ = rhs.model_;
42  if ((model_&&model_->whatsChanged()&1)!=0) {
43    int number = model_->numberRows();
44    if (rhs.savedWeights_) 
45      number = CoinMin(number,rhs.savedWeights_->capacity());
46    if (rhs.infeasible_) {
47      infeasible_= new CoinIndexedVector(rhs.infeasible_);
48    } else {
49      infeasible_=NULL;
50    }
51    if (rhs.weights_) {
52      weights_= new double[number];
53      ClpDisjointCopyN(rhs.weights_,number,weights_);
54    } else {
55      weights_=NULL;
56    }
57    if (rhs.alternateWeights_) {
58      alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
59    } else {
60      alternateWeights_=NULL;
61    }
62    if (rhs.savedWeights_) {
63      savedWeights_= new CoinIndexedVector(rhs.savedWeights_);
64    } else {
65      savedWeights_=NULL;
66    }
67    if (rhs.dubiousWeights_) {
68      assert(model_);
69      int number = model_->numberRows();
70      dubiousWeights_= new int[number];
71      ClpDisjointCopyN(rhs.dubiousWeights_,number,dubiousWeights_);
72    } else {
73      dubiousWeights_=NULL;
74    }
75  } else {
76    infeasible_=NULL;
77    weights_=NULL;
78    alternateWeights_=NULL;
79    savedWeights_=NULL;
80    dubiousWeights_=NULL;
81  }
82}
83
84//-------------------------------------------------------------------
85// Destructor
86//-------------------------------------------------------------------
87ClpDualRowSteepest::~ClpDualRowSteepest ()
88{
89  delete [] weights_;
90  delete [] dubiousWeights_;
91  delete infeasible_;
92  delete alternateWeights_;
93  delete savedWeights_;
94}
95
96//----------------------------------------------------------------
97// Assignment operator
98//-------------------------------------------------------------------
99ClpDualRowSteepest &
100ClpDualRowSteepest::operator=(const ClpDualRowSteepest& rhs)
101{
102  if (this != &rhs) {
103    ClpDualRowPivot::operator=(rhs);
104    state_=rhs.state_;
105    mode_ = rhs.mode_;
106    persistence_ = rhs.persistence_;
107    model_ = rhs.model_;
108    delete [] weights_;
109    delete [] dubiousWeights_;
110    delete infeasible_;
111    delete alternateWeights_;
112    delete savedWeights_;
113    assert(model_);
114    int number = model_->numberRows();
115    if (rhs.savedWeights_) 
116      number = CoinMin(number,rhs.savedWeights_->capacity());
117    if (rhs.infeasible_!=NULL) {
118      infeasible_= new CoinIndexedVector(rhs.infeasible_);
119    } else {
120      infeasible_=NULL;
121    }
122    if (rhs.weights_!=NULL) {
123      weights_= new double[number];
124      ClpDisjointCopyN(rhs.weights_,number,weights_);
125    } else {
126      weights_=NULL;
127    }
128    if (rhs.alternateWeights_!=NULL) {
129      alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
130    } else {
131      alternateWeights_=NULL;
132    }
133    if (rhs.savedWeights_!=NULL) {
134      savedWeights_= new CoinIndexedVector(rhs.savedWeights_);
135    } else {
136      savedWeights_=NULL;
137    }
138    if (rhs.dubiousWeights_) {
139      assert(model_);
140      int number = model_->numberRows();
141      dubiousWeights_= new int[number];
142      ClpDisjointCopyN(rhs.dubiousWeights_,number,dubiousWeights_);
143    } else {
144      dubiousWeights_=NULL;
145    }
146  }
147  return *this;
148}
149// Fill most values
150void 
151ClpDualRowSteepest::fill(const ClpDualRowSteepest& rhs)
152{
153  state_=rhs.state_;
154  mode_ = rhs.mode_;
155  persistence_ = rhs.persistence_;
156  assert (model_->numberRows()==rhs.model_->numberRows());
157  model_ = rhs.model_;
158  assert(model_);
159  int number = model_->numberRows();
160  if (rhs.savedWeights_) 
161    number = CoinMin(number,rhs.savedWeights_->capacity());
162  if (rhs.infeasible_!=NULL) {
163    if (!infeasible_)
164      infeasible_= new CoinIndexedVector(rhs.infeasible_);
165    else
166      *infeasible_ = *rhs.infeasible_;
167  } else {
168    delete infeasible_;
169    infeasible_=NULL;
170  }
171  if (rhs.weights_!=NULL) {
172    if (!weights_)
173      weights_= new double[number];
174    ClpDisjointCopyN(rhs.weights_,number,weights_);
175  } else {
176    delete [] weights_;
177    weights_=NULL;
178  }
179  if (rhs.alternateWeights_!=NULL) {
180    if (!alternateWeights_)
181      alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
182    else
183      *alternateWeights_ = *rhs.alternateWeights_;
184  } else {
185    delete alternateWeights_;
186    alternateWeights_=NULL;
187  }
188  if (rhs.savedWeights_!=NULL) {
189    if (!savedWeights_)
190      savedWeights_= new CoinIndexedVector(rhs.savedWeights_);
191    else
192      *savedWeights_ = *rhs.savedWeights_;
193  } else {
194    delete savedWeights_;
195    savedWeights_=NULL;
196  }
197  if (rhs.dubiousWeights_) {
198    assert(model_);
199    int number = model_->numberRows();
200    if (!dubiousWeights_)
201      dubiousWeights_= new int[number];
202    ClpDisjointCopyN(rhs.dubiousWeights_,number,dubiousWeights_);
203  } else {
204    delete [] dubiousWeights_;
205    dubiousWeights_=NULL;
206  }
207}
208// Returns pivot row, -1 if none
209int 
210ClpDualRowSteepest::pivotRow()
211{
212  assert(model_);
213  int i,iRow;
214  double * infeas = infeasible_->denseVector();
215  double largest=0.0;
216  int * index = infeasible_->getIndices();
217  int number = infeasible_->getNumElements();
218  const int * pivotVariable =model_->pivotVariable();
219  int chosenRow=-1;
220  int lastPivotRow = model_->pivotRow();
221  assert (lastPivotRow<model_->numberRows());
222  double tolerance=model_->currentPrimalTolerance();
223  // we can't really trust infeasibilities if there is primal error
224  // this coding has to mimic coding in checkPrimalSolution
225  double error = CoinMin(1.0e-2,model_->largestPrimalError());
226  // allow tolerance at least slightly bigger than standard
227  tolerance = tolerance +  error;
228  // But cap
229  tolerance = CoinMin(1000.0,tolerance);
230  tolerance *= tolerance; // as we are using squares
231  bool toleranceChanged = false;
232  double * solution = model_->solutionRegion();
233  double * lower = model_->lowerRegion();
234  double * upper = model_->upperRegion();
235  // do last pivot row one here
236  //#define COLUMN_BIAS 4.0
237  //#define FIXED_BIAS 10.0
238  if (lastPivotRow>=0&&lastPivotRow<model_->numberRows()) {
239#if defined (__MINGW32__) || defined(__CYGWIN32__)
240    if (model_->numberIterations()<0)
241      printf("aab_p it %d\n",model_->numberIterations());
242#endif
243#if defined (__MINGW32__) || defined(__CYGWIN32__)
244    if (model_->numberIterations()<0)
245      printf("aab_p it %d\n",model_->numberIterations());
246#endif
247#ifdef COLUMN_BIAS
248    int numberColumns = model_->numberColumns();
249#endif
250    int iPivot=pivotVariable[lastPivotRow];
251    double value = solution[iPivot];
252    double lower = model_->lower(iPivot);
253    double upper = model_->upper(iPivot);
254    if (value>upper+tolerance) {
255      value -= upper;
256      value *= value;
257#ifdef COLUMN_BIAS
258      if (iPivot<numberColumns)
259        value *= COLUMN_BIAS; // bias towards columns
260k
261#endif
262      // store square in list
263      if (infeas[lastPivotRow])
264        infeas[lastPivotRow] = value; // already there
265      else
266        infeasible_->quickAdd(lastPivotRow,value);
267    } else if (value<lower-tolerance) {
268      value -= lower;
269      value *= value;
270#ifdef COLUMN_BIAS
271      if (iPivot<numberColumns)
272        value *= COLUMN_BIAS; // bias towards columns
273#endif
274      // store square in list
275      if (infeas[lastPivotRow])
276        infeas[lastPivotRow] = value; // already there
277      else
278        infeasible_->add(lastPivotRow,value);
279    } else {
280      // feasible - was it infeasible - if so set tiny
281      if (infeas[lastPivotRow])
282        infeas[lastPivotRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
283    }
284    number = infeasible_->getNumElements();
285  }
286#if defined(__MINGW32__) || defined(__CYGWIN32__)
287  if (model_->numberIterations()<0)
288    printf("aac_p it %d\n",model_->numberIterations());
289#endif
290  if(model_->numberIterations()<model_->lastBadIteration()+200) {
291    // we can't really trust infeasibilities if there is dual error
292    if (model_->largestDualError()>model_->largestPrimalError()) {
293      tolerance *= CoinMin(model_->largestDualError()/model_->largestPrimalError(),1000.0);
294      toleranceChanged=true;
295    }
296  }
297  int numberWanted;
298  if (mode_<2 ) {
299    numberWanted = number+1;
300  } else if (mode_==2) {
301    numberWanted = CoinMax(2000,number/8);
302  } else {
303    int numberElements = model_->factorization()->numberElements();
304    double ratio = (double) numberElements/(double) model_->numberRows();
305    numberWanted = CoinMax(2000,number/8);
306    if (ratio<1.0) {
307      numberWanted = CoinMax(2000,number/20);
308    } else if (ratio>10.0) {
309      ratio = number * (ratio/80.0);
310      if (ratio>number)
311        numberWanted=number+1;
312      else
313        numberWanted = CoinMax(2000,(int) ratio);
314    }
315  }
316  if (model_->largestPrimalError()>1.0e-3)
317    numberWanted = number+1; // be safe
318  int iPass;
319  // Setup two passes
320  int start[4];
321  start[1]=number;
322  start[2]=0;
323  double dstart = ((double) number) * model_->randomNumberGenerator()->randomDouble();
324  start[0]=(int) dstart;
325  start[3]=start[0];
326  //double largestWeight=0.0;
327  //double smallestWeight=1.0e100;
328  for (iPass=0;iPass<2;iPass++) {
329    int end = start[2*iPass+1];
330    for (i=start[2*iPass];i<end;i++) {
331      iRow = index[i];
332      double value = infeas[iRow];
333      if (value>tolerance) {
334        //#define OUT_EQ
335#ifdef OUT_EQ
336        {
337          int iSequence = pivotVariable[iRow];
338          if (upper[iSequence]==lower[iSequence])
339            value *= 2.0;
340        }
341#endif
342        double weight = CoinMin(weights_[iRow],1.0e50);
343        //largestWeight = CoinMax(largestWeight,weight);
344        //smallestWeight = CoinMin(smallestWeight,weight);
345        //double dubious = dubiousWeights_[iRow];
346        //weight *= dubious;
347        //if (value>2.0*largest*weight||(value>0.5*largest*weight&&value*largestWeight>dubious*largest*weight)) {
348        if (value>largest*weight) {
349          // make last pivot row last resort choice
350          if (iRow==lastPivotRow) {
351            if (value*1.0e-10<largest*weight) 
352              continue;
353            else 
354              value *= 1.0e-10;
355          }
356          int iSequence = pivotVariable[iRow];
357          if (!model_->flagged(iSequence)) {
358            //#define CLP_DEBUG 3
359#ifdef CLP_DEBUG
360            double value2=0.0;
361            if (solution[iSequence]>upper[iSequence]+tolerance) 
362              value2=solution[iSequence]-upper[iSequence];
363            else if (solution[iSequence]<lower[iSequence]-tolerance) 
364              value2=solution[iSequence]-lower[iSequence];
365            assert(fabs(value2*value2-infeas[iRow])<1.0e-8*CoinMin(value2*value2,infeas[iRow]));
366#endif
367            if (solution[iSequence]>upper[iSequence]+tolerance||
368                solution[iSequence]<lower[iSequence]-tolerance) {
369              chosenRow=iRow;
370              largest=value/weight;
371              //largestWeight = dubious;
372            }
373          } else {
374            // just to make sure we don't exit before got something
375            numberWanted++;
376          }
377        }
378        numberWanted--;
379        if (!numberWanted)
380          break;
381      }
382    }
383    if (!numberWanted)
384      break;
385  }
386  //printf("smallest %g largest %g\n",smallestWeight,largestWeight);
387  if (chosenRow<0 && toleranceChanged) {
388    // won't line up with checkPrimalSolution - do again
389    double saveError = model_->largestDualError();
390    model_->setLargestDualError(0.0);
391    // can't loop
392    chosenRow=pivotRow();
393    model_->setLargestDualError(saveError);
394  }
395  return chosenRow;
396}
397#if 0
398static double ft_count=0.0;
399static double up_count=0.0;
400static double ft_count_in=0.0;
401static double up_count_in=0.0;
402static int xx_count=0;
403#endif
404/* Updates weights and returns pivot alpha.
405   Also does FT update */
406double
407ClpDualRowSteepest::updateWeights(CoinIndexedVector * input,
408                                  CoinIndexedVector * spare,
409                                  CoinIndexedVector * spare2,
410                                  CoinIndexedVector * updatedColumn)
411{
412  //#define CLP_DEBUG 3
413#if CLP_DEBUG>2
414  // Very expensive debug
415  {
416    int numberRows = model_->numberRows();
417    CoinIndexedVector * temp = new CoinIndexedVector();
418    temp->reserve(numberRows+
419                  model_->factorization()->maximumPivots());
420    double * array = alternateWeights_->denseVector();
421    int * which = alternateWeights_->getIndices();
422    alternateWeights_->clear();
423    int i;
424    for (i=0;i<numberRows;i++) {
425      double value=0.0;
426      array[i] = 1.0;
427      which[0] = i;
428      alternateWeights_->setNumElements(1);
429      model_->factorization()->updateColumnTranspose(temp,
430                                                     alternateWeights_);
431      int number = alternateWeights_->getNumElements();
432      int j;
433      for (j=0;j<number;j++) {
434        int iRow=which[j];
435        value+=array[iRow]*array[iRow];
436        array[iRow]=0.0;
437      }
438      alternateWeights_->setNumElements(0);
439      double w = CoinMax(weights_[i],value)*.1;
440      if (fabs(weights_[i]-value)>w) {
441        printf("%d old %g, true %g\n",i,weights_[i],value);
442        weights_[i]=value; // to reduce printout
443      }
444      //else
445      //printf("%d matches %g\n",i,value);
446    }
447    delete temp;
448  }
449#endif
450  assert (input->packedMode());
451  if (!updatedColumn->packedMode()) {
452    // I think this means empty
453#ifdef COIN_DEVELOP
454    printf("updatedColumn not packed mode ClpDualRowSteepest::updateWeights\n");
455#endif
456    return 0.0;
457  }
458  double alpha=0.0;
459  if (!model_->factorization()->networkBasis()) {
460    // clear other region
461    alternateWeights_->clear();
462    double norm = 0.0;
463    int i;
464    double * work = input->denseVector();
465    int numberNonZero = input->getNumElements();
466    int * which = input->getIndices();
467    double * work2 = spare->denseVector();
468    int * which2 = spare->getIndices();
469    // ftran
470    //permute and move indices into index array
471    //also compute norm
472    //int *regionIndex = alternateWeights_->getIndices (  );
473    const int *permute = model_->factorization()->permute();
474    //double * region = alternateWeights_->denseVector();
475    if (permute) {
476      for ( i = 0; i < numberNonZero; i ++ ) {
477        int iRow = which[i];
478        double value = work[i];
479        norm += value*value;
480        iRow = permute[iRow];
481        work2[iRow] = value;
482        which2[i] = iRow;
483      }
484    } else {
485      for ( i = 0; i < numberNonZero; i ++ ) {
486        int iRow = which[i];
487        double value = work[i];
488        norm += value*value;
489        //iRow = permute[iRow];
490        work2[iRow] = value;
491        which2[i] = iRow;
492      }
493    }
494    spare->setNumElements ( numberNonZero );
495    // Do FT update
496#if 0
497    ft_count_in+= updatedColumn->getNumElements();
498    up_count_in+= spare->getNumElements();
499#endif
500    if (permute) {
501#if CLP_DEBUG>2
502      printf("REGION before %d els\n",spare->getNumElements());
503      spare->print();
504#endif
505      model_->factorization()->updateTwoColumnsFT(spare2,updatedColumn,
506                                                  spare,true);
507#if CLP_DEBUG>2
508      printf("REGION after %d els\n",spare->getNumElements());
509      spare->print();
510#endif
511   } else {
512      // Leave as old way
513      model_->factorization()->updateColumnFT(spare2,updatedColumn);
514      model_->factorization()->updateColumn(spare2,spare,false);
515    }
516#undef CLP_DEBUG
517#if 0
518    ft_count+= updatedColumn->getNumElements();
519    up_count+= spare->getNumElements();
520    xx_count++;
521    if ((xx_count%1000)==0)
522      printf("zz %d ft %g up %g (in %g %g)\n",xx_count,ft_count,up_count,
523             ft_count_in,up_count_in);
524#endif
525    numberNonZero = spare->getNumElements();
526    // alternateWeights_ should still be empty
527    int pivotRow = model_->pivotRow();
528#ifdef CLP_DEBUG
529    if ( model_->logLevel (  ) >4  && 
530         fabs(norm-weights_[pivotRow])>1.0e-3*(1.0+norm)) 
531      printf("on row %d, true weight %g, old %g\n",
532             pivotRow,sqrt(norm),sqrt(weights_[pivotRow]));
533#endif
534    // could re-initialize here (could be expensive)
535    norm /= model_->alpha() * model_->alpha();
536   
537    assert(norm);
538    // pivot element
539    alpha=0.0;
540    double multiplier = 2.0 / model_->alpha();
541    // look at updated column
542    work = updatedColumn->denseVector();
543    numberNonZero = updatedColumn->getNumElements();
544    which = updatedColumn->getIndices();
545   
546    int nSave=0;
547    double * work3 = alternateWeights_->denseVector();
548    int * which3 = alternateWeights_->getIndices();
549    const int * pivotColumn = model_->factorization()->pivotColumn();
550    for (i =0; i < numberNonZero; i++) {
551      int iRow = which[i];
552      double theta = work[i];
553      if (iRow==pivotRow)
554        alpha = theta;
555      double devex = weights_[iRow];
556      work3[nSave]=devex; // save old
557      which3[nSave++]=iRow;
558      // transform to match spare
559      int jRow = permute ? pivotColumn[iRow] : iRow;
560      double value = work2[jRow];
561      devex +=  theta * (theta*norm+value * multiplier);
562      if (devex < DEVEX_TRY_NORM) 
563        devex = DEVEX_TRY_NORM;
564      weights_[iRow]=devex;
565    }
566    alternateWeights_->setPackedMode(true);
567    alternateWeights_->setNumElements(nSave);
568    if (norm < DEVEX_TRY_NORM) 
569      norm = DEVEX_TRY_NORM;
570    // Try this to make less likely will happen again and stop cycling
571    //norm *= 1.02;
572    weights_[pivotRow] = norm;
573    spare->clear();
574#ifdef CLP_DEBUG
575    spare->checkClear();
576#endif
577  } else {
578    // Do FT update
579    model_->factorization()->updateColumnFT(spare,updatedColumn);
580    // clear other region
581    alternateWeights_->clear();
582    double norm = 0.0;
583    int i;
584    double * work = input->denseVector();
585    int number = input->getNumElements();
586    int * which = input->getIndices();
587    double * work2 = spare->denseVector();
588    int * which2 = spare->getIndices();
589    for (i=0;i<number;i++) {
590      int iRow = which[i];
591      double value = work[i];
592      norm += value*value;
593      work2[iRow]=value;
594      which2[i]=iRow;
595    }
596    spare->setNumElements(number);
597    // ftran
598#ifndef NDEBUG
599    alternateWeights_->checkClear();
600#endif
601    model_->factorization()->updateColumn(alternateWeights_,spare);
602    // alternateWeights_ should still be empty
603#ifndef NDEBUG
604    alternateWeights_->checkClear();
605#endif
606    int pivotRow = model_->pivotRow();
607#ifdef CLP_DEBUG
608    if ( model_->logLevel (  ) >4  && 
609         fabs(norm-weights_[pivotRow])>1.0e-3*(1.0+norm)) 
610      printf("on row %d, true weight %g, old %g\n",
611             pivotRow,sqrt(norm),sqrt(weights_[pivotRow]));
612#endif
613    // could re-initialize here (could be expensive)
614    norm /= model_->alpha() * model_->alpha();
615   
616    assert(norm);
617    //if (norm < DEVEX_TRY_NORM)
618    //norm = DEVEX_TRY_NORM;
619    // pivot element
620    alpha=0.0;
621    double multiplier = 2.0 / model_->alpha();
622    // look at updated column
623    work = updatedColumn->denseVector();
624    number = updatedColumn->getNumElements();
625    which = updatedColumn->getIndices();
626   
627    int nSave=0;
628    double * work3 = alternateWeights_->denseVector();
629    int * which3 = alternateWeights_->getIndices();
630    for (i =0; i < number; i++) {
631      int iRow = which[i];
632      double theta = work[i];
633      if (iRow==pivotRow)
634        alpha = theta;
635      double devex = weights_[iRow];
636      work3[nSave]=devex; // save old
637      which3[nSave++]=iRow;
638      double value = work2[iRow];
639      devex +=  theta * (theta*norm+value * multiplier);
640      if (devex < DEVEX_TRY_NORM) 
641        devex = DEVEX_TRY_NORM;
642      weights_[iRow]=devex;
643    }
644    if (!alpha) {
645      // error - but carry on
646      alpha =1.0e-50;
647    }
648    alternateWeights_->setPackedMode(true);
649    alternateWeights_->setNumElements(nSave);
650    if (norm < DEVEX_TRY_NORM) 
651      norm = DEVEX_TRY_NORM;
652    weights_[pivotRow] = norm;
653    spare->clear();
654  }
655#ifdef CLP_DEBUG
656  spare->checkClear();
657#endif
658  return alpha;
659}
660 
661/* Updates primal solution (and maybe list of candidates)
662   Uses input vector which it deletes
663   Computes change in objective function
664*/
665void 
666ClpDualRowSteepest::updatePrimalSolution(
667                                        CoinIndexedVector * primalUpdate,
668                                        double primalRatio,
669                                        double & objectiveChange)
670{
671  double * work = primalUpdate->denseVector();
672  int number = primalUpdate->getNumElements();
673  int * which = primalUpdate->getIndices();
674  int i;
675  double changeObj=0.0;
676  double tolerance=model_->currentPrimalTolerance();
677  const int * pivotVariable = model_->pivotVariable();
678  double * infeas = infeasible_->denseVector();
679  int pivotRow = model_->pivotRow();
680  double * solution = model_->solutionRegion();
681#ifdef COLUMN_BIAS
682  int numberColumns = model_->numberColumns();
683#endif
684  if (primalUpdate->packedMode()) {
685    for (i=0;i<number;i++) {
686      int iRow=which[i];
687      int iPivot=pivotVariable[iRow];
688      double value = solution[iPivot];
689      double cost = model_->cost(iPivot);
690      double change = primalRatio*work[i];
691      work[i]=0.0;
692      value -= change;
693      changeObj -= change*cost;
694      solution[iPivot] = value;
695      double lower = model_->lower(iPivot);
696      double upper = model_->upper(iPivot);
697      // But if pivot row then use value of incoming
698      // Although it is safer to recompute before next selection
699      // in case something odd happens
700      if (iRow==pivotRow) {
701        iPivot = model_->sequenceIn();
702        lower = model_->lower(iPivot);
703        upper = model_->upper(iPivot);
704        value = model_->valueIncomingDual();
705      }
706      if (value<lower-tolerance) {
707        value -= lower;
708        value *= value;
709#ifdef COLUMN_BIAS
710        if (iPivot<numberColumns)
711          value *= COLUMN_BIAS; // bias towards columns
712#endif
713#ifdef FIXED_BIAS
714        if (lower==upper)
715          value *= FIXED_BIAS; // bias towards taking out fixed variables
716#endif
717        // store square in list
718        if (infeas[iRow])
719          infeas[iRow] = value; // already there
720        else
721          infeasible_->quickAdd(iRow,value);
722      } else if (value>upper+tolerance) {
723        value -= upper;
724        value *= value;
725#ifdef COLUMN_BIAS
726        if (iPivot<numberColumns)
727          value *= COLUMN_BIAS; // bias towards columns
728#endif
729#ifdef FIXED_BIAS
730        if (lower==upper)
731          value *= FIXED_BIAS; // bias towards taking out fixed variables
732#endif
733        // store square in list
734        if (infeas[iRow])
735          infeas[iRow] = value; // already there
736        else
737          infeasible_->quickAdd(iRow,value);
738      } else {
739        // feasible - was it infeasible - if so set tiny
740        if (infeas[iRow])
741          infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
742      }
743    }
744  } else {
745    for (i=0;i<number;i++) {
746      int iRow=which[i];
747      int iPivot=pivotVariable[iRow];
748      double value = solution[iPivot];
749      double cost = model_->cost(iPivot);
750      double change = primalRatio*work[iRow];
751      value -= change;
752      changeObj -= change*cost;
753      solution[iPivot] = value;
754      double lower = model_->lower(iPivot);
755      double upper = model_->upper(iPivot);
756      // But if pivot row then use value of incoming
757      // Although it is safer to recompute before next selection
758      // in case something odd happens
759      if (iRow==pivotRow) {
760        iPivot = model_->sequenceIn();
761        lower = model_->lower(iPivot);
762        upper = model_->upper(iPivot);
763        value = model_->valueIncomingDual();
764      }
765      if (value<lower-tolerance) {
766        value -= lower;
767        value *= value;
768#ifdef COLUMN_BIAS
769        if (iPivot<numberColumns)
770          value *= COLUMN_BIAS; // bias towards columns
771#endif
772#ifdef FIXED_BIAS
773        if (lower==upper)
774          value *= FIXED_BIAS; // bias towards taking out fixed variables
775#endif
776        // store square in list
777        if (infeas[iRow])
778          infeas[iRow] = value; // already there
779        else
780          infeasible_->quickAdd(iRow,value);
781      } else if (value>upper+tolerance) {
782        value -= upper;
783        value *= value;
784#ifdef COLUMN_BIAS
785        if (iPivot<numberColumns)
786          value *= COLUMN_BIAS; // bias towards columns
787#endif
788#ifdef FIXED_BIAS
789        if (lower==upper)
790          value *= FIXED_BIAS; // bias towards taking out fixed variables
791#endif
792        // store square in list
793        if (infeas[iRow])
794          infeas[iRow] = value; // already there
795        else
796          infeasible_->quickAdd(iRow,value);
797      } else {
798        // feasible - was it infeasible - if so set tiny
799        if (infeas[iRow])
800          infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
801      }
802      work[iRow]=0.0;
803    }
804  }
805  primalUpdate->setNumElements(0);
806  objectiveChange += changeObj;
807}
808/* Saves any weights round factorization as pivot rows may change
809   1) before factorization
810   2) after factorization
811   3) just redo infeasibilities
812   4) restore weights
813*/
814void 
815ClpDualRowSteepest::saveWeights(ClpSimplex * model,int mode)
816{
817  // alternateWeights_ is defined as indexed but is treated oddly
818  model_ = model;
819  int numberRows = model_->numberRows();
820  int numberColumns = model_->numberColumns();
821  const int * pivotVariable = model_->pivotVariable();
822  int i;
823  if (mode==1) {
824    if(weights_) {
825      // Check if size has changed
826      if (infeasible_->capacity()==numberRows) {
827        alternateWeights_->clear();
828        // change from row numbers to sequence numbers
829        int * which = alternateWeights_->getIndices();
830        for (i=0;i<numberRows;i++) {
831          int iPivot=pivotVariable[i];
832          which[i]=iPivot;
833        }
834        state_=1;
835      } else {
836        // size has changed - clear everything
837        delete [] weights_;
838        weights_=NULL;
839        delete [] dubiousWeights_;
840        dubiousWeights_=NULL;
841        delete infeasible_;
842        infeasible_=NULL;
843        delete alternateWeights_;
844        alternateWeights_=NULL;
845        delete savedWeights_;
846        savedWeights_=NULL;
847        state_=-1;
848      }
849    }
850  } else if (mode==2||mode==4||mode==5) {
851    // restore
852    if (!weights_||state_==-1||mode==5) {
853      // initialize weights
854      delete [] weights_;
855      delete alternateWeights_;
856      weights_ = new double[numberRows];
857      alternateWeights_ = new CoinIndexedVector();
858      // enough space so can use it for factorization
859      alternateWeights_->reserve(numberRows+
860                                 model_->factorization()->maximumPivots());
861      if (mode_!=1||mode==5) {
862        // initialize to 1.0 (can we do better?)
863        for (i=0;i<numberRows;i++) {
864          weights_[i]=1.0;
865        }
866      } else {
867        CoinIndexedVector * temp = new CoinIndexedVector();
868        temp->reserve(numberRows+
869                      model_->factorization()->maximumPivots());
870        double * array = alternateWeights_->denseVector();
871        int * which = alternateWeights_->getIndices();
872        for (i=0;i<numberRows;i++) {
873          double value=0.0;
874          array[0] = 1.0;
875          which[0] = i;
876          alternateWeights_->setNumElements(1);
877          alternateWeights_->setPackedMode(true);
878          model_->factorization()->updateColumnTranspose(temp,
879                                                         alternateWeights_);
880          int number = alternateWeights_->getNumElements();
881          int j;
882          for (j=0;j<number;j++) {
883            value+=array[j]*array[j];
884            array[j]=0.0;
885          }
886          alternateWeights_->setNumElements(0);
887          weights_[i] = value;
888        }
889        delete temp;
890      }
891      // create saved weights (not really indexedvector)
892      savedWeights_ = new CoinIndexedVector();
893      savedWeights_->reserve(numberRows);
894     
895      double * array = savedWeights_->denseVector();
896      int * which = savedWeights_->getIndices();
897      for (i=0;i<numberRows;i++) {
898        array[i]=weights_[i];
899        which[i]=pivotVariable[i];
900      }
901    } else {
902      int * which = alternateWeights_->getIndices();
903      CoinIndexedVector * rowArray3 = model_->rowArray(3);
904      rowArray3->clear();
905      int * back = rowArray3->getIndices();
906      // In case something went wrong
907      for (i=0;i<numberRows+numberColumns;i++)
908        back[i]=-1;
909      if (mode!=4) {
910        // save
911 CoinMemcpyN(which,     numberRows,savedWeights_->getIndices());
912 CoinMemcpyN(weights_,  numberRows,savedWeights_->denseVector());
913      } else {
914        // restore
915        //memcpy(which,savedWeights_->getIndices(),
916        //     numberRows*sizeof(int));
917        //memcpy(weights_,savedWeights_->denseVector(),
918        //     numberRows*sizeof(double));
919        which = savedWeights_->getIndices();
920      }
921      // restore (a bit slow - but only every re-factorization)
922      double * array = savedWeights_->denseVector();
923      for (i=0;i<numberRows;i++) {
924        int iSeq=which[i];
925        back[iSeq]=i;
926      }
927      for (i=0;i<numberRows;i++) {
928        int iPivot=pivotVariable[i];
929        iPivot = back[iPivot];
930        if (iPivot>=0) {
931          weights_[i]=array[iPivot];
932          if (weights_[i]<DEVEX_TRY_NORM)
933            weights_[i] = DEVEX_TRY_NORM; // may need to check more
934        } else {
935          // odd
936          weights_[i]=1.0;
937        }
938      }
939    }
940    state_=0;
941    // set up infeasibilities
942    if (!infeasible_) {
943      infeasible_ = new CoinIndexedVector();
944      infeasible_->reserve(numberRows);
945    }
946  }
947  if (mode>=2) {
948    // Get dubious weights
949    //if (!dubiousWeights_)
950    //dubiousWeights_=new int[numberRows];
951    //model_->factorization()->getWeights(dubiousWeights_);
952    infeasible_->clear();
953    int iRow;
954    const int * pivotVariable = model_->pivotVariable();
955    double tolerance=model_->currentPrimalTolerance();
956    for (iRow=0;iRow<numberRows;iRow++) {
957      int iPivot=pivotVariable[iRow];
958      double value = model_->solution(iPivot);
959      double lower = model_->lower(iPivot);
960      double upper = model_->upper(iPivot);
961      if (value<lower-tolerance) {
962        value -= lower;
963        value *= value;
964#ifdef COLUMN_BIAS
965        if (iPivot<numberColumns)
966          value *= COLUMN_BIAS; // bias towards columns
967#endif
968#ifdef FIXED_BIAS
969        if (lower==upper)
970          value *= FIXED_BIAS; // bias towards taking out fixed variables
971#endif
972        // store square in list
973        infeasible_->quickAdd(iRow,value);
974      } else if (value>upper+tolerance) {
975        value -= upper;
976        value *= value;
977#ifdef COLUMN_BIAS
978        if (iPivot<numberColumns)
979          value *= COLUMN_BIAS; // bias towards columns
980#endif
981#ifdef FIXED_BIAS
982        if (lower==upper)
983          value *= FIXED_BIAS; // bias towards taking out fixed variables
984#endif
985        // store square in list
986        infeasible_->quickAdd(iRow,value);
987      }
988    }
989  }
990}
991// Gets rid of last update
992void 
993ClpDualRowSteepest::unrollWeights()
994{
995  double * saved = alternateWeights_->denseVector();
996  int number = alternateWeights_->getNumElements();
997  int * which = alternateWeights_->getIndices();
998  int i;
999  if (alternateWeights_->packedMode()) {
1000    for (i=0;i<number;i++) {
1001      int iRow = which[i];
1002      weights_[iRow]=saved[i];
1003      saved[i]=0.0;
1004    }
1005  } else {
1006    for (i=0;i<number;i++) {
1007      int iRow = which[i];
1008      weights_[iRow]=saved[iRow];
1009      saved[iRow]=0.0;
1010    }
1011  }
1012  alternateWeights_->setNumElements(0);
1013}
1014//-------------------------------------------------------------------
1015// Clone
1016//-------------------------------------------------------------------
1017ClpDualRowPivot * ClpDualRowSteepest::clone(bool CopyData) const
1018{
1019  if (CopyData) {
1020    return new ClpDualRowSteepest(*this);
1021  } else {
1022    return new ClpDualRowSteepest();
1023  }
1024}
1025// Gets rid of all arrays
1026void 
1027ClpDualRowSteepest::clearArrays()
1028{
1029  if (persistence_==normal) {
1030    delete [] weights_;
1031    weights_=NULL;
1032    delete [] dubiousWeights_;
1033    dubiousWeights_=NULL;
1034    delete infeasible_;
1035    infeasible_ = NULL;
1036    delete alternateWeights_;
1037    alternateWeights_ = NULL;
1038    delete savedWeights_;
1039    savedWeights_ = NULL;
1040  }
1041  state_ =-1;
1042}
1043// Returns true if would not find any row
1044bool 
1045ClpDualRowSteepest::looksOptimal() const
1046{
1047  int iRow;
1048  const int * pivotVariable = model_->pivotVariable();
1049  double tolerance=model_->currentPrimalTolerance();
1050  // we can't really trust infeasibilities if there is primal error
1051  // this coding has to mimic coding in checkPrimalSolution
1052  double error = CoinMin(1.0e-2,model_->largestPrimalError());
1053  // allow tolerance at least slightly bigger than standard
1054  tolerance = tolerance +  error;
1055  // But cap
1056  tolerance = CoinMin(1000.0,tolerance);
1057  int numberRows = model_->numberRows();
1058  int numberInfeasible=0;
1059  for (iRow=0;iRow<numberRows;iRow++) {
1060    int iPivot=pivotVariable[iRow];
1061    double value = model_->solution(iPivot);
1062    double lower = model_->lower(iPivot);
1063    double upper = model_->upper(iPivot);
1064    if (value<lower-tolerance) {
1065      numberInfeasible++;
1066    } else if (value>upper+tolerance) {
1067      numberInfeasible++;
1068    }
1069  }
1070  return (numberInfeasible==0);
1071}
1072// Called when maximum pivots changes
1073void 
1074ClpDualRowSteepest::maximumPivotsChanged()
1075{
1076  if (alternateWeights_&&
1077      alternateWeights_->capacity()!=model_->numberRows()+
1078      model_->factorization()->maximumPivots()) {
1079    delete alternateWeights_;
1080    alternateWeights_ = new CoinIndexedVector();
1081    // enough space so can use it for factorization
1082    alternateWeights_->reserve(model_->numberRows()+
1083                                   model_->factorization()->maximumPivots());
1084  }
1085}
1086
Note: See TracBrowser for help on using the repository browser.