source: trunk/ClpDualRowSteepest.cpp @ 225

Last change on this file since 225 was 225, checked in by forrest, 16 years ago

This should break everything

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