source: branches/pre/ClpDualRowSteepest.cpp @ 180

Last change on this file since 180 was 180, checked in by forrest, 18 years ago

new stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.5 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//#############################################################################
13// Constructors / Destructor / Assignment
14//#############################################################################
15
16//-------------------------------------------------------------------
17// Default Constructor
18//-------------------------------------------------------------------
19ClpDualRowSteepest::ClpDualRowSteepest (int mode) 
20  : ClpDualRowPivot(),
21    state_(-1),
22    mode_(mode),
23    weights_(NULL),
24    infeasible_(NULL),
25    alternateWeights_(NULL),
26    savedWeights_(NULL)
27{
28  type_=2+64*mode;
29}
30
31//-------------------------------------------------------------------
32// Copy constructor
33//-------------------------------------------------------------------
34ClpDualRowSteepest::ClpDualRowSteepest (const ClpDualRowSteepest & rhs) 
35: ClpDualRowPivot(rhs)
36{ 
37  state_=rhs.state_;
38  mode_ = rhs.mode_;
39  model_ = rhs.model_;
40  if (rhs.infeasible_) {
41    infeasible_= new CoinIndexedVector(rhs.infeasible_);
42  } else {
43    infeasible_=NULL;
44  }
45  if (rhs.weights_) {
46    assert(model_);
47    int number = model_->numberRows();
48    weights_= new double[number];
49    ClpDisjointCopyN(rhs.weights_,number,weights_);
50  } else {
51    weights_=NULL;
52  }
53  if (rhs.alternateWeights_) {
54    alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
55  } else {
56    alternateWeights_=NULL;
57  }
58  if (rhs.savedWeights_) {
59    savedWeights_= new CoinIndexedVector(rhs.savedWeights_);
60  } else {
61    savedWeights_=NULL;
62  }
63}
64
65//-------------------------------------------------------------------
66// Destructor
67//-------------------------------------------------------------------
68ClpDualRowSteepest::~ClpDualRowSteepest ()
69{
70  delete [] weights_;
71  delete infeasible_;
72  delete alternateWeights_;
73  delete savedWeights_;
74}
75
76//----------------------------------------------------------------
77// Assignment operator
78//-------------------------------------------------------------------
79ClpDualRowSteepest &
80ClpDualRowSteepest::operator=(const ClpDualRowSteepest& rhs)
81{
82  if (this != &rhs) {
83    ClpDualRowPivot::operator=(rhs);
84    state_=rhs.state_;
85    mode_ = rhs.mode_;
86    model_ = rhs.model_;
87    delete [] weights_;
88    delete infeasible_;
89    delete alternateWeights_;
90    delete savedWeights_;
91    if (rhs.infeasible_!=NULL) {
92      infeasible_= new CoinIndexedVector(rhs.infeasible_);
93    } else {
94      infeasible_=NULL;
95    }
96    if (rhs.weights_!=NULL) {
97      assert(model_);
98      int number = model_->numberRows();
99      weights_= new double[number];
100      ClpDisjointCopyN(rhs.weights_,number,weights_);
101    } else {
102      weights_=NULL;
103    }
104    if (rhs.alternateWeights_!=NULL) {
105      alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
106    } else {
107      alternateWeights_=NULL;
108    }
109    if (rhs.savedWeights_!=NULL) {
110      savedWeights_= new CoinIndexedVector(rhs.savedWeights_);
111    } else {
112      savedWeights_=NULL;
113    }
114  }
115  return *this;
116}
117
118// Returns pivot row, -1 if none
119int 
120ClpDualRowSteepest::pivotRow()
121{
122  assert(model_);
123  int i,iRow;
124  double * infeas = infeasible_->denseVector();
125  double largest=1.0e-50;
126  int * index = infeasible_->getIndices();
127  int number = infeasible_->getNumElements();
128  const int * pivotVariable =model_->pivotVariable();
129  int chosenRow=-1;
130  int lastPivotRow = model_->pivotRow();
131  double tolerance=model_->currentPrimalTolerance();
132  // we can't really trust infeasibilities if there is primal error
133  // this coding has to mimic coding in checkPrimalSolution
134  double error = min(1.0e-3,model_->largestPrimalError());
135  // allow tolerance at least slightly bigger than standard
136  tolerance = tolerance +  error;
137  // But cap
138  tolerance = min(1000.0,tolerance);
139  tolerance *= tolerance; // as we are using squares
140  double * solution = model_->solutionRegion();
141  double * lower = model_->lowerRegion();
142  double * upper = model_->upperRegion();
143  // do last pivot row one here
144  //#define COLUMN_BIAS 4.0
145  //#define FIXED_BIAS 10.0
146  if (lastPivotRow>=0) {
147#ifdef COLUMN_BIAS
148    int numberColumns = model_->numberColumns();
149#endif
150    int iPivot=pivotVariable[lastPivotRow];
151    double value = solution[iPivot];
152    double lower = model_->lower(iPivot);
153    double upper = model_->upper(iPivot);
154    if (value>upper+tolerance) {
155      value -= upper;
156      value *= value;
157#ifdef COLUMN_BIAS
158      if (iPivot<numberColumns)
159        value *= COLUMN_BIAS; // bias towards columns
160k
161#endif
162      // store square in list
163      if (infeas[lastPivotRow])
164        infeas[lastPivotRow] = value; // already there
165      else
166        infeasible_->quickAdd(lastPivotRow,value);
167    } else if (value<lower-tolerance) {
168      value -= lower;
169      value *= value;
170#ifdef COLUMN_BIAS
171      if (iPivot<numberColumns)
172        value *= COLUMN_BIAS; // bias towards columns
173#endif
174      // store square in list
175      if (infeas[lastPivotRow])
176        infeas[lastPivotRow] = value; // already there
177      else
178        infeasible_->add(lastPivotRow,value);
179    } else {
180      // feasible - was it infeasible - if so set tiny
181      if (infeas[lastPivotRow])
182        infeas[lastPivotRow] = 1.0e-100;
183    }
184    number = infeasible_->getNumElements();
185  }
186  if(model_->numberIterations()<model_->lastBadIteration()+200) {
187    // we can't really trust infeasibilities if there is dual error
188    double checkTolerance = 1.0e-8;
189    if (!model_->factorization()->pivots())
190      checkTolerance = 1.0e-6;
191    if (model_->largestPrimalError()>checkTolerance)
192      tolerance *= model_->largestPrimalError()/checkTolerance;
193  }
194  int numberWanted;
195  if (mode_<2)
196    numberWanted = number+1;
197  else
198    numberWanted = max(2000,number/8);
199  int iPass;
200  // Setup two passes
201  int start[4];
202  start[1]=number;
203  start[2]=0;
204  double dstart = ((double) number) * CoinDrand48();
205  start[0]=(int) dstart;
206  start[3]=start[0];
207  for (iPass=0;iPass<2;iPass++) {
208    int end = start[2*iPass+1];
209    for (i=start[2*iPass];i<end;i++) {
210      iRow = index[i];
211      double value = infeas[iRow];
212      if (value>tolerance) {
213        if (value>largest*weights_[iRow]) {
214          // make last pivot row last resort choice
215          if (iRow==lastPivotRow) {
216            if (value*1.0e-10<largest*weights_[iRow]) 
217              continue;
218            else 
219              value *= 1.0e-10;
220          }
221          int iSequence = pivotVariable[iRow];
222          if (!model_->flagged(iSequence)) {
223            //#define CLP_DEBUG 1
224#ifdef CLP_DEBUG
225            double value2=0.0;
226            if (solution[iSequence]>upper[iSequence]+tolerance) 
227              value2=solution[iSequence]-upper[iSequence];
228            else if (solution[iSequence]<lower[iSequence]-tolerance) 
229              value2=solution[iSequence]-lower[iSequence];
230            assert(fabs(value2*value2-infeas[iRow])<1.0e-8*min(value2*value2,infeas[iRow]));
231#endif
232            if (solution[iSequence]>upper[iSequence]+tolerance||
233                solution[iSequence]<lower[iSequence]-tolerance) {
234              chosenRow=iRow;
235              largest=value/weights_[iRow];
236            }
237          }
238        }
239        numberWanted--;
240        if (!numberWanted)
241          break;
242      }
243    }
244    if (!numberWanted)
245      break;
246  }
247  return chosenRow;
248}
249#define TRY_NORM 1.0e-4
250// Updates weights and returns pivot alpha
251double
252ClpDualRowSteepest::updateWeights(CoinIndexedVector * input,
253                                  CoinIndexedVector * spare,
254                                  CoinIndexedVector * updatedColumn)
255{
256#if CLP_DEBUG>2
257  // Very expensive debug
258  {
259    int numberRows = model_->numberRows();
260    CoinIndexedVector * temp = new CoinIndexedVector();
261    temp->reserve(numberRows+
262                  model_->factorization()->maximumPivots());
263    double * array = alternateWeights_->denseVector();
264    int * which = alternateWeights_->getIndices();
265    for (i=0;i<numberRows;i++) {
266      double value=0.0;
267      array[i] = 1.0;
268      which[0] = i;
269      alternateWeights_->setNumElements(1);
270      model_->factorization()->updateColumnTranspose(temp,
271                                                     alternateWeights_);
272      int number = alternateWeights_->getNumElements();
273      int j;
274      for (j=0;j<number;j++) {
275        int iRow=which[j];
276        value+=array[iRow]*array[iRow];
277        array[iRow]=0.0;
278      }
279      alternateWeights_->setNumElements(0);
280      if (fabs(weights_[i]-value)>1.0e-4)
281        printf("%d old %g, true %g\n",i,weights_[i],value);
282    }
283    delete temp;
284  }
285#endif
286  assert (input->packedMode());
287  assert (updatedColumn->packedMode());
288  double alpha=0.0;
289  if (!model_->factorization()->networkBasis()) {
290    // clear other region
291    alternateWeights_->clear();
292    double norm = 0.0;
293    int i;
294    double * work = input->denseVector();
295    int numberNonZero = input->getNumElements();
296    int * which = input->getIndices();
297    double * work2 = spare->denseVector();
298    int * which2 = spare->getIndices();
299    // ftran
300    //permute and move indices into index array
301    //also compute norm
302    //int *regionIndex = alternateWeights_->getIndices (  );
303    const int *permute = model_->factorization()->permute();
304    //double * region = alternateWeights_->denseVector();
305    for ( i = 0; i < numberNonZero; i ++ ) {
306      int iRow = which[i];
307      double value = work[i];
308      norm += value*value;
309      iRow = permute[iRow];
310      work2[iRow] = value;
311      which2[i] = iRow;
312    }
313    spare->setNumElements ( numberNonZero );
314    model_->factorization()->updateColumn(spare,spare,true);
315    // permute back
316    numberNonZero = spare->getNumElements();
317    // alternateWeights_ should still be empty
318    int pivotRow = model_->pivotRow();
319    // could re-initialize here (could be expensive)
320    norm /= model_->alpha() * model_->alpha();
321   
322    assert(norm);
323    // pivot element
324    alpha=0.0;
325    double multiplier = 2.0 / model_->alpha();
326    // look at updated column
327    work = updatedColumn->denseVector();
328    numberNonZero = updatedColumn->getNumElements();
329    which = updatedColumn->getIndices();
330   
331    int nSave=0;
332    double * work3 = alternateWeights_->denseVector();
333    int * which3 = alternateWeights_->getIndices();
334    const int * pivotColumn = model_->factorization()->pivotColumn();
335    for (i =0; i < numberNonZero; i++) {
336      int iRow = which[i];
337      double theta = work[i];
338      if (iRow==pivotRow)
339        alpha = theta;
340      double devex = weights_[iRow];
341      work3[nSave]=devex; // save old
342      which3[nSave++]=iRow;
343      // transform to match spare
344      int jRow = pivotColumn[iRow];
345      double value = work2[jRow];
346      devex +=  theta * (theta*norm+value * multiplier);
347      if (devex < TRY_NORM) 
348        devex = TRY_NORM;
349      weights_[iRow]=devex;
350    }
351    assert (alpha);
352    alternateWeights_->setPackedMode(true);
353    alternateWeights_->setNumElements(nSave);
354    if (norm < TRY_NORM) 
355      norm = TRY_NORM;
356    weights_[pivotRow] = norm;
357    spare->clear();
358#ifdef CLP_DEBUG
359    spare->checkClear();
360#endif
361  } else {
362    // clear other region
363    alternateWeights_->clear();
364    double norm = 0.0;
365    int i;
366    double * work = input->denseVector();
367    int number = input->getNumElements();
368    int * which = input->getIndices();
369    double * work2 = spare->denseVector();
370    int * which2 = spare->getIndices();
371    for (i=0;i<number;i++) {
372      int iRow = which[i];
373      double value = work[i];
374      norm += value*value;
375      work2[iRow]=value;
376      which2[i]=iRow;
377    }
378    spare->setNumElements(number);
379    // ftran
380    model_->factorization()->updateColumn(alternateWeights_,spare);
381    // alternateWeights_ should still be empty
382    int pivotRow = model_->pivotRow();
383#ifdef CLP_DEBUG
384    if ( model_->logLevel (  ) >4  && 
385         fabs(norm-weights_[pivotRow])>1.0e-3*(1.0+norm)) 
386      printf("on row %d, true weight %g, old %g\n",
387             pivotRow,sqrt(norm),sqrt(weights_[pivotRow]));
388#endif
389    // could re-initialize here (could be expensive)
390    norm /= model_->alpha() * model_->alpha();
391   
392    assert(norm);
393    //if (norm < TRY_NORM)
394    //norm = TRY_NORM;
395    // pivot element
396    alpha=0.0;
397    double multiplier = 2.0 / model_->alpha();
398    // look at updated column
399    work = updatedColumn->denseVector();
400    number = updatedColumn->getNumElements();
401    which = updatedColumn->getIndices();
402   
403    int nSave=0;
404    double * work3 = alternateWeights_->denseVector();
405    int * which3 = alternateWeights_->getIndices();
406    for (i =0; i < number; i++) {
407      int iRow = which[i];
408      double theta = work[i];
409      if (iRow==pivotRow)
410        alpha = theta;
411      double devex = weights_[iRow];
412      work3[nSave]=devex; // save old
413      which3[nSave++]=iRow;
414      double value = work2[iRow];
415      devex +=  theta * (theta*norm+value * multiplier);
416      if (devex < TRY_NORM) 
417        devex = TRY_NORM;
418      weights_[iRow]=devex;
419    }
420    assert (alpha);
421    alternateWeights_->setPackedMode(true);
422    alternateWeights_->setNumElements(nSave);
423    if (norm < TRY_NORM) 
424      norm = TRY_NORM;
425    weights_[pivotRow] = norm;
426    spare->clear();
427  }
428#ifdef CLP_DEBUG
429  spare->checkClear();
430#endif
431  return alpha;
432}
433 
434/* Updates primal solution (and maybe list of candidates)
435   Uses input vector which it deletes
436   Computes change in objective function
437*/
438void 
439ClpDualRowSteepest::updatePrimalSolution(
440                                        CoinIndexedVector * primalUpdate,
441                                        double primalRatio,
442                                        double & objectiveChange)
443{
444  double * work = primalUpdate->denseVector();
445  int number = primalUpdate->getNumElements();
446  int * which = primalUpdate->getIndices();
447  int i;
448  double changeObj=0.0;
449  double tolerance=model_->currentPrimalTolerance();
450  const int * pivotVariable = model_->pivotVariable();
451  double * infeas = infeasible_->denseVector();
452  int pivotRow = model_->pivotRow();
453  double * solution = model_->solutionRegion();
454#ifdef COLUMN_BIAS
455  int numberColumns = model_->numberColumns();
456#endif
457  if (primalUpdate->packedMode()) {
458    for (i=0;i<number;i++) {
459      int iRow=which[i];
460      int iPivot=pivotVariable[iRow];
461      double value = solution[iPivot];
462      double cost = model_->cost(iPivot);
463      double change = primalRatio*work[i];
464      work[i]=0.0;
465      value -= change;
466      changeObj -= change*cost;
467      solution[iPivot] = value;
468      double lower = model_->lower(iPivot);
469      double upper = model_->upper(iPivot);
470      // But if pivot row then use value of incoming
471      // Although it is safer to recompute before next selection
472      // in case something odd happens
473      if (iRow==pivotRow) {
474        iPivot = model_->sequenceIn();
475        lower = model_->lower(iPivot);
476        upper = model_->upper(iPivot);
477        value = model_->valueIncomingDual();
478      }
479      if (value<lower-tolerance) {
480        value -= lower;
481        value *= value;
482#ifdef COLUMN_BIAS
483        if (iPivot<numberColumns)
484          value *= COLUMN_BIAS; // bias towards columns
485#endif
486#ifdef FIXED_BIAS
487        if (lower==upper)
488          value *= FIXED_BIAS; // bias towards taking out fixed variables
489#endif
490        // store square in list
491        if (infeas[iRow])
492          infeas[iRow] = value; // already there
493        else
494          infeasible_->quickAdd(iRow,value);
495      } else if (value>upper+tolerance) {
496        value -= upper;
497        value *= value;
498#ifdef COLUMN_BIAS
499        if (iPivot<numberColumns)
500          value *= COLUMN_BIAS; // bias towards columns
501#endif
502#ifdef FIXED_BIAS
503        if (lower==upper)
504          value *= FIXED_BIAS; // bias towards taking out fixed variables
505#endif
506        // store square in list
507        if (infeas[iRow])
508          infeas[iRow] = value; // already there
509        else
510          infeasible_->quickAdd(iRow,value);
511      } else {
512        // feasible - was it infeasible - if so set tiny
513        if (infeas[iRow])
514          infeas[iRow] = 1.0e-100;
515      }
516    }
517  } else {
518    for (i=0;i<number;i++) {
519      int iRow=which[i];
520      int iPivot=pivotVariable[iRow];
521      double value = solution[iPivot];
522      double cost = model_->cost(iPivot);
523      double change = primalRatio*work[iRow];
524      value -= change;
525      changeObj -= change*cost;
526      solution[iPivot] = value;
527      double lower = model_->lower(iPivot);
528      double upper = model_->upper(iPivot);
529      // But if pivot row then use value of incoming
530      // Although it is safer to recompute before next selection
531      // in case something odd happens
532      if (iRow==pivotRow) {
533        iPivot = model_->sequenceIn();
534        lower = model_->lower(iPivot);
535        upper = model_->upper(iPivot);
536        value = model_->valueIncomingDual();
537      }
538      if (value<lower-tolerance) {
539        value -= lower;
540        value *= value;
541#ifdef COLUMN_BIAS
542        if (iPivot<numberColumns)
543          value *= COLUMN_BIAS; // bias towards columns
544#endif
545#ifdef FIXED_BIAS
546        if (lower==upper)
547          value *= FIXED_BIAS; // bias towards taking out fixed variables
548#endif
549        // store square in list
550        if (infeas[iRow])
551          infeas[iRow] = value; // already there
552        else
553          infeasible_->quickAdd(iRow,value);
554      } else if (value>upper+tolerance) {
555        value -= upper;
556        value *= value;
557#ifdef COLUMN_BIAS
558        if (iPivot<numberColumns)
559          value *= COLUMN_BIAS; // bias towards columns
560#endif
561#ifdef FIXED_BIAS
562        if (lower==upper)
563          value *= FIXED_BIAS; // bias towards taking out fixed variables
564#endif
565        // store square in list
566        if (infeas[iRow])
567          infeas[iRow] = value; // already there
568        else
569          infeasible_->quickAdd(iRow,value);
570      } else {
571        // feasible - was it infeasible - if so set tiny
572        if (infeas[iRow])
573          infeas[iRow] = 1.0e-100;
574      }
575      work[iRow]=0.0;
576    }
577  }
578  primalUpdate->setNumElements(0);
579  objectiveChange += changeObj;
580}
581/* Saves any weights round factorization as pivot rows may change
582   1) before factorization
583   2) after factorization
584   3) just redo infeasibilities
585   4) restore weights
586*/
587void 
588ClpDualRowSteepest::saveWeights(ClpSimplex * model,int mode)
589{
590  // alternateWeights_ is defined as indexed but is treated oddly
591  model_ = model;
592  int numberRows = model_->numberRows();
593  int numberColumns = model_->numberColumns();
594  const int * pivotVariable = model_->pivotVariable();
595  int i;
596  if (mode==1&&weights_) {
597    alternateWeights_->clear();
598    // change from row numbers to sequence numbers
599    int * which = alternateWeights_->getIndices();
600    for (i=0;i<numberRows;i++) {
601      int iPivot=pivotVariable[i];
602      which[i]=iPivot;
603    }
604    state_=1;
605  } else if (mode==2||mode==4||mode==5) {
606    // restore
607    if (!weights_||state_==-1||mode==5) {
608      // initialize weights
609      delete [] weights_;
610      delete alternateWeights_;
611      weights_ = new double[numberRows];
612      alternateWeights_ = new CoinIndexedVector();
613      // enough space so can use it for factorization
614      alternateWeights_->reserve(numberRows+
615                                 model_->factorization()->maximumPivots());
616      if (!mode_||mode_==2||mode==5) {
617        // initialize to 1.0 (can we do better?)
618        for (i=0;i<numberRows;i++) {
619          weights_[i]=1.0;
620        }
621      } else {
622        CoinIndexedVector * temp = new CoinIndexedVector();
623        temp->reserve(numberRows+
624                      model_->factorization()->maximumPivots());
625        double * array = alternateWeights_->denseVector();
626        int * which = alternateWeights_->getIndices();
627        for (i=0;i<numberRows;i++) {
628          double value=0.0;
629          array[0] = 1.0;
630          which[0] = i;
631          alternateWeights_->setNumElements(1);
632          alternateWeights_->setPackedMode(true);
633          model_->factorization()->updateColumnTranspose(temp,
634                                                         alternateWeights_);
635          int number = alternateWeights_->getNumElements();
636          int j;
637          for (j=0;j<number;j++) {
638            value+=array[j]*array[j];
639            array[j]=0.0;
640          }
641          alternateWeights_->setNumElements(0);
642          weights_[i] = value;
643        }
644        delete temp;
645      }
646      // create saved weights (not really indexedvector)
647      savedWeights_ = new CoinIndexedVector();
648      savedWeights_->reserve(numberRows);
649     
650      double * array = savedWeights_->denseVector();
651      int * which = savedWeights_->getIndices();
652      for (i=0;i<numberRows;i++) {
653        array[i]=weights_[i];
654        which[i]=pivotVariable[i];
655      }
656    } else {
657      int * which = alternateWeights_->getIndices();
658      if (mode!=4) {
659        // save
660        memcpy(savedWeights_->getIndices(),which,
661               numberRows*sizeof(int));
662        memcpy(savedWeights_->denseVector(),weights_,
663               numberRows*sizeof(double));
664      } else {
665        // restore
666        memcpy(which,savedWeights_->getIndices(),
667               numberRows*sizeof(int));
668        memcpy(weights_,savedWeights_->denseVector(),
669               numberRows*sizeof(double));
670      }
671      // restore (a bit slow - but only every re-factorization)
672      double * array = new double[numberRows+numberColumns];
673      for (i=0;i<numberRows;i++) {
674        int iSeq=which[i];
675        array[iSeq]=weights_[i];
676      }
677      for (i=0;i<numberRows;i++) {
678        int iPivot=pivotVariable[i];
679        weights_[i]=array[iPivot];
680        if (weights_[i]<TRY_NORM)
681          weights_[i] = TRY_NORM; // may need to check more
682      }
683      delete [] array;
684    }
685    state_=0;
686    // set up infeasibilities
687    if (!infeasible_) {
688      infeasible_ = new CoinIndexedVector();
689      infeasible_->reserve(numberRows);
690    }
691  }
692  if (mode>=2) {
693    infeasible_->clear();
694    int iRow;
695    const int * pivotVariable = model_->pivotVariable();
696    double tolerance=model_->currentPrimalTolerance();
697    for (iRow=0;iRow<numberRows;iRow++) {
698      int iPivot=pivotVariable[iRow];
699      double value = model_->solution(iPivot);
700      double lower = model_->lower(iPivot);
701      double upper = model_->upper(iPivot);
702      if (value<lower-tolerance) {
703        value -= lower;
704        value *= value;
705#ifdef COLUMN_BIAS
706        if (iPivot<numberColumns)
707          value *= COLUMN_BIAS; // bias towards columns
708#endif
709#ifdef FIXED_BIAS
710        if (lower==upper)
711          value *= FIXED_BIAS; // bias towards taking out fixed variables
712#endif
713        // store square in list
714        infeasible_->quickAdd(iRow,value);
715      } else if (value>upper+tolerance) {
716        value -= upper;
717        value *= value;
718#ifdef COLUMN_BIAS
719        if (iPivot<numberColumns)
720          value *= COLUMN_BIAS; // bias towards columns
721#endif
722#ifdef FIXED_BIAS
723        if (lower==upper)
724          value *= FIXED_BIAS; // bias towards taking out fixed variables
725#endif
726        // store square in list
727        infeasible_->quickAdd(iRow,value);
728      }
729    }
730  }
731}
732// Gets rid of last update
733void 
734ClpDualRowSteepest::unrollWeights()
735{
736  double * saved = alternateWeights_->denseVector();
737  int number = alternateWeights_->getNumElements();
738  int * which = alternateWeights_->getIndices();
739  int i;
740  if (alternateWeights_->packedMode()) {
741    for (i=0;i<number;i++) {
742      int iRow = which[i];
743      weights_[iRow]=saved[i];
744      saved[i]=0.0;
745    }
746  } else {
747    for (i=0;i<number;i++) {
748      int iRow = which[i];
749      weights_[iRow]=saved[iRow];
750      saved[iRow]=0.0;
751    }
752  }
753  alternateWeights_->setNumElements(0);
754}
755//-------------------------------------------------------------------
756// Clone
757//-------------------------------------------------------------------
758ClpDualRowPivot * ClpDualRowSteepest::clone(bool CopyData) const
759{
760  if (CopyData) {
761    return new ClpDualRowSteepest(*this);
762  } else {
763    return new ClpDualRowSteepest();
764  }
765}
766// Gets rid of all arrays
767void 
768ClpDualRowSteepest::clearArrays()
769{
770  delete [] weights_;
771  weights_=NULL;
772  delete infeasible_;
773  infeasible_ = NULL;
774  delete alternateWeights_;
775  alternateWeights_ = NULL;
776  delete savedWeights_;
777  savedWeights_ = NULL;
778  state_ =-1;
779}
780// Returns true if would not find any row
781bool 
782ClpDualRowSteepest::looksOptimal() const
783{
784  int iRow;
785  const int * pivotVariable = model_->pivotVariable();
786  double tolerance=model_->currentPrimalTolerance();
787  // we can't really trust infeasibilities if there is primal error
788  // this coding has to mimic coding in checkPrimalSolution
789  double error = min(1.0e-3,model_->largestPrimalError());
790  // allow tolerance at least slightly bigger than standard
791  tolerance = tolerance +  error;
792  // But cap
793  tolerance = min(1000.0,tolerance);
794  int numberRows = model_->numberRows();
795  int numberInfeasible=0;
796  for (iRow=0;iRow<numberRows;iRow++) {
797    int iPivot=pivotVariable[iRow];
798    double value = model_->solution(iPivot);
799    double lower = model_->lower(iPivot);
800    double upper = model_->upper(iPivot);
801    if (value<lower-tolerance) {
802      numberInfeasible++;
803    } else if (value>upper+tolerance) {
804      numberInfeasible++;
805    }
806  }
807  return (numberInfeasible==0);
808}
809
Note: See TracBrowser for help on using the repository browser.