source: trunk/ClpDualRowSteepest.cpp @ 437

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

to CoinMax? and CoinMin?

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