source: stable/1.6/Clp/src/ClpDualRowSteepest.cpp @ 1609

Last change on this file since 1609 was 1034, checked in by forrest, 13 years ago

moving branches/devel to trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 28.2 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-2,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  if (!updatedColumn->packedMode()) {
367    // I think this means empty
368#ifdef COIN_DEVELOP
369    printf("updatedColumn not packed mode ClpDualRowSteepest::updateWeights\n");
370#endif
371    return 0.0;
372  }
373  double alpha=0.0;
374  if (!model_->factorization()->networkBasis()) {
375    // clear other region
376    alternateWeights_->clear();
377    double norm = 0.0;
378    int i;
379    double * work = input->denseVector();
380    int numberNonZero = input->getNumElements();
381    int * which = input->getIndices();
382    double * work2 = spare->denseVector();
383    int * which2 = spare->getIndices();
384    // ftran
385    //permute and move indices into index array
386    //also compute norm
387    //int *regionIndex = alternateWeights_->getIndices (  );
388    const int *permute = model_->factorization()->permute();
389    //double * region = alternateWeights_->denseVector();
390    for ( i = 0; i < numberNonZero; i ++ ) {
391      int iRow = which[i];
392      double value = work[i];
393      norm += value*value;
394      iRow = permute[iRow];
395      work2[iRow] = value;
396      which2[i] = iRow;
397    }
398    spare->setNumElements ( numberNonZero );
399    // Only one array active as already permuted
400    model_->factorization()->updateColumn(spare,spare,true);
401    // permute back
402    numberNonZero = spare->getNumElements();
403    // alternateWeights_ should still be empty
404    int pivotRow = model_->pivotRow();
405#ifdef CLP_DEBUG
406    if ( model_->logLevel (  ) >4  && 
407         fabs(norm-weights_[pivotRow])>1.0e-3*(1.0+norm)) 
408      printf("on row %d, true weight %g, old %g\n",
409             pivotRow,sqrt(norm),sqrt(weights_[pivotRow]));
410#endif
411    // could re-initialize here (could be expensive)
412    norm /= model_->alpha() * model_->alpha();
413   
414    assert(norm);
415    // pivot element
416    alpha=0.0;
417    double multiplier = 2.0 / model_->alpha();
418    // look at updated column
419    work = updatedColumn->denseVector();
420    numberNonZero = updatedColumn->getNumElements();
421    which = updatedColumn->getIndices();
422   
423    int nSave=0;
424    double * work3 = alternateWeights_->denseVector();
425    int * which3 = alternateWeights_->getIndices();
426    const int * pivotColumn = model_->factorization()->pivotColumn();
427    for (i =0; i < numberNonZero; i++) {
428      int iRow = which[i];
429      double theta = work[i];
430      if (iRow==pivotRow)
431        alpha = theta;
432      double devex = weights_[iRow];
433      work3[nSave]=devex; // save old
434      which3[nSave++]=iRow;
435      // transform to match spare
436      int jRow = pivotColumn[iRow];
437      double value = work2[jRow];
438      devex +=  theta * (theta*norm+value * multiplier);
439      if (devex < DEVEX_TRY_NORM) 
440        devex = DEVEX_TRY_NORM;
441      weights_[iRow]=devex;
442    }
443    alternateWeights_->setPackedMode(true);
444    alternateWeights_->setNumElements(nSave);
445    if (norm < DEVEX_TRY_NORM) 
446      norm = DEVEX_TRY_NORM;
447    // Try this to make less likely will happen again and stop cycling
448    //norm *= 1.02;
449    weights_[pivotRow] = norm;
450    spare->clear();
451#ifdef CLP_DEBUG
452    spare->checkClear();
453#endif
454  } else {
455    // clear other region
456    alternateWeights_->clear();
457    double norm = 0.0;
458    int i;
459    double * work = input->denseVector();
460    int number = input->getNumElements();
461    int * which = input->getIndices();
462    double * work2 = spare->denseVector();
463    int * which2 = spare->getIndices();
464    for (i=0;i<number;i++) {
465      int iRow = which[i];
466      double value = work[i];
467      norm += value*value;
468      work2[iRow]=value;
469      which2[i]=iRow;
470    }
471    spare->setNumElements(number);
472    // ftran
473    model_->factorization()->updateColumn(alternateWeights_,spare);
474    // alternateWeights_ should still be empty
475    int pivotRow = model_->pivotRow();
476#ifdef CLP_DEBUG
477    if ( model_->logLevel (  ) >4  && 
478         fabs(norm-weights_[pivotRow])>1.0e-3*(1.0+norm)) 
479      printf("on row %d, true weight %g, old %g\n",
480             pivotRow,sqrt(norm),sqrt(weights_[pivotRow]));
481#endif
482    // could re-initialize here (could be expensive)
483    norm /= model_->alpha() * model_->alpha();
484   
485    assert(norm);
486    //if (norm < DEVEX_TRY_NORM)
487    //norm = DEVEX_TRY_NORM;
488    // pivot element
489    alpha=0.0;
490    double multiplier = 2.0 / model_->alpha();
491    // look at updated column
492    work = updatedColumn->denseVector();
493    number = updatedColumn->getNumElements();
494    which = updatedColumn->getIndices();
495   
496    int nSave=0;
497    double * work3 = alternateWeights_->denseVector();
498    int * which3 = alternateWeights_->getIndices();
499    for (i =0; i < number; i++) {
500      int iRow = which[i];
501      double theta = work[i];
502      if (iRow==pivotRow)
503        alpha = theta;
504      double devex = weights_[iRow];
505      work3[nSave]=devex; // save old
506      which3[nSave++]=iRow;
507      double value = work2[iRow];
508      devex +=  theta * (theta*norm+value * multiplier);
509      if (devex < DEVEX_TRY_NORM) 
510        devex = DEVEX_TRY_NORM;
511      weights_[iRow]=devex;
512    }
513    if (!alpha) {
514      // error - but carry on
515      alpha =1.0e-50;
516    }
517    alternateWeights_->setPackedMode(true);
518    alternateWeights_->setNumElements(nSave);
519    if (norm < DEVEX_TRY_NORM) 
520      norm = DEVEX_TRY_NORM;
521    weights_[pivotRow] = norm;
522    spare->clear();
523  }
524#ifdef CLP_DEBUG
525  spare->checkClear();
526#endif
527  return alpha;
528}
529 
530/* Updates primal solution (and maybe list of candidates)
531   Uses input vector which it deletes
532   Computes change in objective function
533*/
534void 
535ClpDualRowSteepest::updatePrimalSolution(
536                                        CoinIndexedVector * primalUpdate,
537                                        double primalRatio,
538                                        double & objectiveChange)
539{
540  double * work = primalUpdate->denseVector();
541  int number = primalUpdate->getNumElements();
542  int * which = primalUpdate->getIndices();
543  int i;
544  double changeObj=0.0;
545  double tolerance=model_->currentPrimalTolerance();
546  const int * pivotVariable = model_->pivotVariable();
547  double * infeas = infeasible_->denseVector();
548  int pivotRow = model_->pivotRow();
549  double * solution = model_->solutionRegion();
550#ifdef COLUMN_BIAS
551  int numberColumns = model_->numberColumns();
552#endif
553  if (primalUpdate->packedMode()) {
554    for (i=0;i<number;i++) {
555      int iRow=which[i];
556      int iPivot=pivotVariable[iRow];
557      double value = solution[iPivot];
558      double cost = model_->cost(iPivot);
559      double change = primalRatio*work[i];
560      work[i]=0.0;
561      value -= change;
562      changeObj -= change*cost;
563      solution[iPivot] = value;
564      double lower = model_->lower(iPivot);
565      double upper = model_->upper(iPivot);
566      // But if pivot row then use value of incoming
567      // Although it is safer to recompute before next selection
568      // in case something odd happens
569      if (iRow==pivotRow) {
570        iPivot = model_->sequenceIn();
571        lower = model_->lower(iPivot);
572        upper = model_->upper(iPivot);
573        value = model_->valueIncomingDual();
574      }
575      if (value<lower-tolerance) {
576        value -= lower;
577        value *= value;
578#ifdef COLUMN_BIAS
579        if (iPivot<numberColumns)
580          value *= COLUMN_BIAS; // bias towards columns
581#endif
582#ifdef FIXED_BIAS
583        if (lower==upper)
584          value *= FIXED_BIAS; // bias towards taking out fixed variables
585#endif
586        // store square in list
587        if (infeas[iRow])
588          infeas[iRow] = value; // already there
589        else
590          infeasible_->quickAdd(iRow,value);
591      } else if (value>upper+tolerance) {
592        value -= upper;
593        value *= value;
594#ifdef COLUMN_BIAS
595        if (iPivot<numberColumns)
596          value *= COLUMN_BIAS; // bias towards columns
597#endif
598#ifdef FIXED_BIAS
599        if (lower==upper)
600          value *= FIXED_BIAS; // bias towards taking out fixed variables
601#endif
602        // store square in list
603        if (infeas[iRow])
604          infeas[iRow] = value; // already there
605        else
606          infeasible_->quickAdd(iRow,value);
607      } else {
608        // feasible - was it infeasible - if so set tiny
609        if (infeas[iRow])
610          infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
611      }
612    }
613  } else {
614    for (i=0;i<number;i++) {
615      int iRow=which[i];
616      int iPivot=pivotVariable[iRow];
617      double value = solution[iPivot];
618      double cost = model_->cost(iPivot);
619      double change = primalRatio*work[iRow];
620      value -= change;
621      changeObj -= change*cost;
622      solution[iPivot] = value;
623      double lower = model_->lower(iPivot);
624      double upper = model_->upper(iPivot);
625      // But if pivot row then use value of incoming
626      // Although it is safer to recompute before next selection
627      // in case something odd happens
628      if (iRow==pivotRow) {
629        iPivot = model_->sequenceIn();
630        lower = model_->lower(iPivot);
631        upper = model_->upper(iPivot);
632        value = model_->valueIncomingDual();
633      }
634      if (value<lower-tolerance) {
635        value -= lower;
636        value *= value;
637#ifdef COLUMN_BIAS
638        if (iPivot<numberColumns)
639          value *= COLUMN_BIAS; // bias towards columns
640#endif
641#ifdef FIXED_BIAS
642        if (lower==upper)
643          value *= FIXED_BIAS; // bias towards taking out fixed variables
644#endif
645        // store square in list
646        if (infeas[iRow])
647          infeas[iRow] = value; // already there
648        else
649          infeasible_->quickAdd(iRow,value);
650      } else if (value>upper+tolerance) {
651        value -= upper;
652        value *= value;
653#ifdef COLUMN_BIAS
654        if (iPivot<numberColumns)
655          value *= COLUMN_BIAS; // bias towards columns
656#endif
657#ifdef FIXED_BIAS
658        if (lower==upper)
659          value *= FIXED_BIAS; // bias towards taking out fixed variables
660#endif
661        // store square in list
662        if (infeas[iRow])
663          infeas[iRow] = value; // already there
664        else
665          infeasible_->quickAdd(iRow,value);
666      } else {
667        // feasible - was it infeasible - if so set tiny
668        if (infeas[iRow])
669          infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
670      }
671      work[iRow]=0.0;
672    }
673  }
674  primalUpdate->setNumElements(0);
675  objectiveChange += changeObj;
676}
677/* Saves any weights round factorization as pivot rows may change
678   1) before factorization
679   2) after factorization
680   3) just redo infeasibilities
681   4) restore weights
682*/
683void 
684ClpDualRowSteepest::saveWeights(ClpSimplex * model,int mode)
685{
686  // alternateWeights_ is defined as indexed but is treated oddly
687  model_ = model;
688  int numberRows = model_->numberRows();
689  int numberColumns = model_->numberColumns();
690  const int * pivotVariable = model_->pivotVariable();
691  int i;
692  if (mode==1) {
693    if(weights_) {
694      // Check if size has changed
695      if (infeasible_->capacity()==numberRows) {
696        alternateWeights_->clear();
697        // change from row numbers to sequence numbers
698        int * which = alternateWeights_->getIndices();
699        for (i=0;i<numberRows;i++) {
700          int iPivot=pivotVariable[i];
701          which[i]=iPivot;
702        }
703        state_=1;
704      } else {
705        // size has changed - clear everything
706        delete [] weights_;
707        weights_=NULL;
708        delete [] dubiousWeights_;
709        dubiousWeights_=NULL;
710        delete infeasible_;
711        infeasible_=NULL;
712        delete alternateWeights_;
713        alternateWeights_=NULL;
714        delete savedWeights_;
715        savedWeights_=NULL;
716        state_=-1;
717      }
718    }
719  } else if (mode==2||mode==4||mode==5) {
720    // restore
721    if (!weights_||state_==-1||mode==5) {
722      // initialize weights
723      delete [] weights_;
724      delete alternateWeights_;
725      weights_ = new double[numberRows];
726      alternateWeights_ = new CoinIndexedVector();
727      // enough space so can use it for factorization
728      alternateWeights_->reserve(numberRows+
729                                 model_->factorization()->maximumPivots());
730      if (mode_!=1||mode==5) {
731        // initialize to 1.0 (can we do better?)
732        for (i=0;i<numberRows;i++) {
733          weights_[i]=1.0;
734        }
735      } else {
736        CoinIndexedVector * temp = new CoinIndexedVector();
737        temp->reserve(numberRows+
738                      model_->factorization()->maximumPivots());
739        double * array = alternateWeights_->denseVector();
740        int * which = alternateWeights_->getIndices();
741        for (i=0;i<numberRows;i++) {
742          double value=0.0;
743          array[0] = 1.0;
744          which[0] = i;
745          alternateWeights_->setNumElements(1);
746          alternateWeights_->setPackedMode(true);
747          model_->factorization()->updateColumnTranspose(temp,
748                                                         alternateWeights_);
749          int number = alternateWeights_->getNumElements();
750          int j;
751          for (j=0;j<number;j++) {
752            value+=array[j]*array[j];
753            array[j]=0.0;
754          }
755          alternateWeights_->setNumElements(0);
756          weights_[i] = value;
757        }
758        delete temp;
759      }
760      // create saved weights (not really indexedvector)
761      savedWeights_ = new CoinIndexedVector();
762      savedWeights_->reserve(numberRows);
763     
764      double * array = savedWeights_->denseVector();
765      int * which = savedWeights_->getIndices();
766      for (i=0;i<numberRows;i++) {
767        array[i]=weights_[i];
768        which[i]=pivotVariable[i];
769      }
770    } else {
771      int * which = alternateWeights_->getIndices();
772      CoinIndexedVector * rowArray3 = model_->rowArray(3);
773      rowArray3->clear();
774      int * back = rowArray3->getIndices();
775      // In case something went wrong
776      for (i=0;i<numberRows+numberColumns;i++)
777        back[i]=-1;
778      if (mode!=4) {
779        // save
780        memcpy(savedWeights_->getIndices(),which,
781               numberRows*sizeof(int));
782        memcpy(savedWeights_->denseVector(),weights_,
783               numberRows*sizeof(double));
784      } else {
785        // restore
786        //memcpy(which,savedWeights_->getIndices(),
787        //     numberRows*sizeof(int));
788        //memcpy(weights_,savedWeights_->denseVector(),
789        //     numberRows*sizeof(double));
790        which = savedWeights_->getIndices();
791      }
792      // restore (a bit slow - but only every re-factorization)
793      double * array = savedWeights_->denseVector();
794      for (i=0;i<numberRows;i++) {
795        int iSeq=which[i];
796        back[iSeq]=i;
797      }
798      for (i=0;i<numberRows;i++) {
799        int iPivot=pivotVariable[i];
800        iPivot = back[iPivot];
801        if (iPivot>=0) {
802          weights_[i]=array[iPivot];
803          if (weights_[i]<DEVEX_TRY_NORM)
804            weights_[i] = DEVEX_TRY_NORM; // may need to check more
805        } else {
806          // odd
807          weights_[i]=1.0;
808        }
809      }
810    }
811    state_=0;
812    // set up infeasibilities
813    if (!infeasible_) {
814      infeasible_ = new CoinIndexedVector();
815      infeasible_->reserve(numberRows);
816    }
817  }
818  if (mode>=2) {
819    // Get dubious weights
820    //if (!dubiousWeights_)
821    //dubiousWeights_=new int[numberRows];
822    //model_->factorization()->getWeights(dubiousWeights_);
823    infeasible_->clear();
824    int iRow;
825    const int * pivotVariable = model_->pivotVariable();
826    double tolerance=model_->currentPrimalTolerance();
827    for (iRow=0;iRow<numberRows;iRow++) {
828      int iPivot=pivotVariable[iRow];
829      double value = model_->solution(iPivot);
830      double lower = model_->lower(iPivot);
831      double upper = model_->upper(iPivot);
832      if (value<lower-tolerance) {
833        value -= lower;
834        value *= value;
835#ifdef COLUMN_BIAS
836        if (iPivot<numberColumns)
837          value *= COLUMN_BIAS; // bias towards columns
838#endif
839#ifdef FIXED_BIAS
840        if (lower==upper)
841          value *= FIXED_BIAS; // bias towards taking out fixed variables
842#endif
843        // store square in list
844        infeasible_->quickAdd(iRow,value);
845      } else if (value>upper+tolerance) {
846        value -= upper;
847        value *= value;
848#ifdef COLUMN_BIAS
849        if (iPivot<numberColumns)
850          value *= COLUMN_BIAS; // bias towards columns
851#endif
852#ifdef FIXED_BIAS
853        if (lower==upper)
854          value *= FIXED_BIAS; // bias towards taking out fixed variables
855#endif
856        // store square in list
857        infeasible_->quickAdd(iRow,value);
858      }
859    }
860  }
861}
862// Gets rid of last update
863void 
864ClpDualRowSteepest::unrollWeights()
865{
866  double * saved = alternateWeights_->denseVector();
867  int number = alternateWeights_->getNumElements();
868  int * which = alternateWeights_->getIndices();
869  int i;
870  if (alternateWeights_->packedMode()) {
871    for (i=0;i<number;i++) {
872      int iRow = which[i];
873      weights_[iRow]=saved[i];
874      saved[i]=0.0;
875    }
876  } else {
877    for (i=0;i<number;i++) {
878      int iRow = which[i];
879      weights_[iRow]=saved[iRow];
880      saved[iRow]=0.0;
881    }
882  }
883  alternateWeights_->setNumElements(0);
884}
885//-------------------------------------------------------------------
886// Clone
887//-------------------------------------------------------------------
888ClpDualRowPivot * ClpDualRowSteepest::clone(bool CopyData) const
889{
890  if (CopyData) {
891    return new ClpDualRowSteepest(*this);
892  } else {
893    return new ClpDualRowSteepest();
894  }
895}
896// Gets rid of all arrays
897void 
898ClpDualRowSteepest::clearArrays()
899{
900  if (persistence_==normal) {
901    delete [] weights_;
902    weights_=NULL;
903    delete [] dubiousWeights_;
904    dubiousWeights_=NULL;
905    delete infeasible_;
906    infeasible_ = NULL;
907    delete alternateWeights_;
908    alternateWeights_ = NULL;
909    delete savedWeights_;
910    savedWeights_ = NULL;
911  }
912  state_ =-1;
913}
914// Returns true if would not find any row
915bool 
916ClpDualRowSteepest::looksOptimal() const
917{
918  int iRow;
919  const int * pivotVariable = model_->pivotVariable();
920  double tolerance=model_->currentPrimalTolerance();
921  // we can't really trust infeasibilities if there is primal error
922  // this coding has to mimic coding in checkPrimalSolution
923  double error = CoinMin(1.0e-2,model_->largestPrimalError());
924  // allow tolerance at least slightly bigger than standard
925  tolerance = tolerance +  error;
926  // But cap
927  tolerance = CoinMin(1000.0,tolerance);
928  int numberRows = model_->numberRows();
929  int numberInfeasible=0;
930  for (iRow=0;iRow<numberRows;iRow++) {
931    int iPivot=pivotVariable[iRow];
932    double value = model_->solution(iPivot);
933    double lower = model_->lower(iPivot);
934    double upper = model_->upper(iPivot);
935    if (value<lower-tolerance) {
936      numberInfeasible++;
937    } else if (value>upper+tolerance) {
938      numberInfeasible++;
939    }
940  }
941  return (numberInfeasible==0);
942}
943// Called when maximum pivots changes
944void 
945ClpDualRowSteepest::maximumPivotsChanged()
946{
947  if (alternateWeights_&&
948      alternateWeights_->capacity()!=model_->numberRows()+
949      model_->factorization()->maximumPivots()) {
950    delete alternateWeights_;
951    alternateWeights_ = new CoinIndexedVector();
952    // enough space so can use it for factorization
953    alternateWeights_->reserve(model_->numberRows()+
954                                   model_->factorization()->maximumPivots());
955  }
956}
957
Note: See TracBrowser for help on using the repository browser.