source: tags/move-to-subversion/ClpDualRowSteepest.cpp @ 753

Last change on this file since 753 was 691, checked in by forrest, 15 years ago

mainly for Francois

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