source: trunk/Clp/src/ClpDualRowSteepest.cpp @ 1304

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

changes for cbc event handler and multiple factorizations

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