source: branches/pre/ClpSolve.cpp @ 222

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

Start of mini sprint

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 25.1 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4// This file has higher level solve functions
5
6
7#include "CoinPragma.hpp"
8
9#include <math.h>
10
11#include "CoinHelperFunctions.hpp"
12#include "CoinSort.hpp"
13#include "ClpFactorization.hpp"
14#include "ClpSimplex.hpp"
15#include "ClpSolve.hpp"
16#include "ClpPackedMatrix.hpp"
17#include "ClpPlusMinusOneMatrix.hpp"
18#include "ClpNetworkMatrix.hpp"
19#include "ClpMessage.hpp"
20#include "CoinTime.hpp"
21
22#include "ClpPresolve.hpp"
23#ifdef CLP_IDIOT
24#include "Idiot.hpp"
25#endif
26//#############################################################################
27// Allow for interrupts
28// But is this threadsafe ? (so switched off by option
29#include <signal.h>
30static ClpSimplex * currentModel = NULL;
31static void signal_handler(int whichSignal)
32{
33  if (currentModel!=NULL) 
34    currentModel->setMaximumIterations(0); // stop at next iterations
35  return;
36}
37
38/** General solve algorithm which can do presolve
39    special options (bits)
40    1 - do not perturb
41    2 - do not scale
42    4 - use crash (default allslack in dual, idiot in primal)
43    8 - all slack basis in primal
44    16 - switch off interrupt handling
45    32 - do not try and make plus minus one matrix
46 */
47int 
48ClpSimplex::initialSolve(ClpSolve & options)
49{
50  ClpSolve::SolveType method=options.getSolveType();
51  ClpSolve::PresolveType presolve = options.getPresolveType();
52  int saveMaxIterations = maximumIterations();
53  int finalStatus=-1;
54  int numberIterations=0;
55  double time1 = CoinCpuTime();
56  double timeX = time1;
57  double time2;
58  ClpMatrixBase * saveMatrix=NULL;
59  ClpSimplex * model2 = this;
60  bool interrupt = (options.getSpecialOption(2)==0);
61  sighandler_t saveSignal=SIG_DFL;
62  if (interrupt) {
63    currentModel = model2;
64    // register signal handler
65    saveSignal = signal(SIGINT,signal_handler);
66  }
67  ClpPresolve pinfo;
68  double timePresolve=0.0;
69  double timeIdiot=0.0;
70  double timeCore=0.0;
71  double timeSimplex=0.0;
72  int savePerturbation=perturbation_;
73  int saveScaling = scalingFlag_;
74  if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
75    // network - switch off stuff
76    presolve = ClpSolve::presolveOff;
77  }
78
79  if (presolve!=ClpSolve::presolveOff) {
80    int numberPasses=5;
81    if (presolve==ClpSolve::presolveNumber) {
82      numberPasses=options.getPresolvePasses();
83      presolve = ClpSolve::presolveOn;
84    }
85    model2 = pinfo.presolvedModel(*this,1.0e-8,
86                                  false,numberPasses,true);
87    time2 = CoinCpuTime();
88    timePresolve = time2-timeX;
89    handler_->message(CLP_INTERVAL_TIMING,messages_)
90      <<"Presolve"<<timePresolve<<time2-time1
91      <<CoinMessageEol;
92    timeX=time2;
93    if (model2) {
94      if (method==ClpSolve::useDual) {
95        int numberInfeasibilities = model2->tightenPrimalBounds();
96        if (numberInfeasibilities) {
97          handler_->message(CLP_INFEASIBLE,messages_)
98            <<CoinMessageEol;
99          model2 = this;
100          presolve=ClpSolve::presolveOff;
101        }
102      }
103    } else {
104      handler_->message(CLP_INFEASIBLE,messages_)
105        <<CoinMessageEol;
106      model2 = this;
107      presolve=ClpSolve::presolveOff;
108    }
109    // We may be better off using original (but if dual leave because of bounds)
110    if (presolve!=ClpSolve::presolveOff&&
111        numberRows_<1.01*model2->numberRows_&&numberColumns_<1.01*model2->numberColumns_
112        &&method!=ClpSolve::useDual) {
113      delete model2;
114      model2 = this;
115      presolve=ClpSolve::presolveOff;
116    }
117  }
118  if (interrupt)
119    currentModel = model2;
120  // See if worth trying +- one matrix
121  bool plusMinus=false;
122  int numberElements=model2->getNumElements();
123  // For below >0 overrides
124  // 0 means no, -1 means maybe
125  int doIdiot=0;
126  int doCrash=0;
127  int doSprint=0;
128  switch (options.getSpecialOption(1)) {
129  case 0:
130    doIdiot=-1;
131    doCrash=-1;
132    doSprint=-1;
133    break;
134  case 1:
135    doIdiot=0;
136    doCrash=1;
137    doSprint=0;
138    break;
139  case 2:
140    doIdiot=1;
141    doCrash=0;
142    doSprint=0;
143    break;
144  case 3:
145    doIdiot=0;
146    doCrash=0;
147    doSprint=1;
148    break;
149  case 4:
150    doIdiot=0;
151    doCrash=0;
152    doSprint=0;
153    break;
154  case 5:
155    doIdiot=0;
156    doCrash=-1;
157    doSprint=-1;
158    break;
159  case 6:
160    doIdiot=-1;
161    doCrash=-1;
162    doSprint=0;
163    break;
164  case 7:
165    doIdiot=-1;
166    doCrash=0;
167    doSprint=-1;
168    break;
169  case 8:
170    doIdiot=-1;
171    doCrash=0;
172    doSprint=0;
173    break;
174  case 9:
175    doIdiot=0;
176    doCrash=0;
177    doSprint=-1;
178    break;
179  default:
180    abort();
181  }
182  if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
183    // network - switch off stuff
184    doIdiot=0;
185    doSprint=0;
186  }
187  int numberColumns = model2->numberColumns();
188  int numberRows = model2->numberRows();
189  // If not all slack basis - switch off all
190  int number=0;
191  int iRow;
192  for (iRow=0;iRow<numberRows;iRow++)
193    if (model2->getRowStatus(iRow)==basic)
194      number++;
195  if (number<numberRows) {
196    doIdiot=0;
197    doCrash=0;
198    doSprint=0;
199  }
200  if (options.getSpecialOption(3)==0) {
201    if(numberElements>100000)
202      plusMinus=true;
203    if(numberElements>10000&&((doIdiot||doSprint)&&method==ClpSolve::usePrimal)) 
204      plusMinus=true;
205  }
206  if (plusMinus) {
207    saveMatrix = model2->clpMatrix();
208    ClpPackedMatrix* clpMatrix =
209      dynamic_cast< ClpPackedMatrix*>(saveMatrix);
210    if (clpMatrix) {
211      ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
212      if (newMatrix->getIndices()) {
213        model2->replaceMatrix(newMatrix);
214      } else {
215        handler_->message(CLP_MATRIX_CHANGE,messages_)
216          <<"+- 1"
217          <<CoinMessageEol;
218        saveMatrix=NULL;
219        plusMinus=false;
220        delete newMatrix;
221      }
222    } else {
223      saveMatrix=NULL;
224      plusMinus=false;
225    }
226  }
227  if (model2->factorizationFrequency()==200) {
228    // User did not touch preset
229    model2->setFactorizationFrequency(100+model2->numberRows()/200);
230  }
231  if (method==ClpSolve::usePrimalorSprint) {
232    if (doSprint<0) { 
233      if(numberRows*10>numberColumns||numberColumns<6000
234         ||(numberRows*20>numberColumns&&!plusMinus))
235        method=ClpSolve::usePrimal; // switch off sprint
236    } else if (doSprint==0) {
237      method=ClpSolve::usePrimal; // switch off sprint
238    }
239  }
240  if (method==ClpSolve::useDual) {
241    if (options.getSpecialOption(0)!=0)
242      model2->crash(1000,1);
243    model2->dual();
244    time2 = CoinCpuTime();
245    timeCore = time2-timeX;
246    handler_->message(CLP_INTERVAL_TIMING,messages_)
247      <<"Dual"<<timeCore<<time2-time1
248      <<CoinMessageEol;
249    timeX=time2;
250  } else if (method==ClpSolve::usePrimal) {
251#ifdef CLP_IDIOT
252    if (doIdiot) {
253      int nPasses=0;
254      Idiot info(*model2);
255      // Get average number of elements per column
256      double ratio  = ((double) numberElements/(double) numberColumns);
257      // look at rhs
258      int iRow;
259      double largest=0.0;
260      double smallest = 1.0e30;
261      double largestGap=0.0;
262      int numberNotE=0;
263      for (iRow=0;iRow<numberRows;iRow++) {
264        double value1 = model2->rowLower_[iRow];
265        if (value1&&value1>-1.0e31) {
266          largest = max(largest,fabs(value1));
267          smallest=min(smallest,fabs(value1));
268        }
269        double value2 = model2->rowUpper_[iRow];
270        if (value2&&value2<1.0e31) {
271          largest = max(largest,fabs(value2));
272          smallest=min(smallest,fabs(value2));
273        }
274        if (value2>value1) {
275          numberNotE++;
276          if (value2>1.0e31||value1<-1.0e31)
277            largestGap = COIN_DBL_MAX;
278          else
279            largestGap = value2-value1;
280        }
281      }
282      if (numberRows>200&&numberColumns>5000&&numberColumns>2*numberRows) {
283        if (plusMinus) {
284          if (largest/smallest>2.0) {
285            nPasses = 10+numberColumns/100000;
286            nPasses = min(nPasses,50);
287            nPasses = max(nPasses,15);
288            if (numberRows>25000&&nPasses>5) {
289              // Might as well go for it
290              nPasses = max(nPasses,71);
291            } else if (numberElements<3*numberColumns) {
292              nPasses=min(nPasses,10); // probably not worh it
293            }
294            if (doIdiot<0)
295              info.setLightweight(1); // say lightweight idiot
296          } else if (largest/smallest>1.01||numberElements<=3*numberColumns) {
297            nPasses = 10+numberColumns/1000;
298            nPasses = min(nPasses,100);
299            nPasses = max(nPasses,30);
300            if (numberRows>25000) {
301              // Might as well go for it
302              nPasses = max(nPasses,71);
303            }
304            if (!largestGap)
305              nPasses *= 2;
306          } else {
307            nPasses = 10+numberColumns/1000;
308            nPasses = min(nPasses,200);
309            nPasses = max(nPasses,100);
310            info.setStartingWeight(1.0e-1);
311            info.setReduceIterations(6);
312            if (!largestGap)
313              nPasses *= 2;
314            //info.setFeasibilityTolerance(1.0e-7);
315          }
316          // If few passes - don't bother
317          if (nPasses<=5)
318            nPasses=0;
319        } else {
320          if (doIdiot<0)
321            info.setLightweight(1); // say lightweight idiot
322          if (largest/smallest>1.01||numberNotE) {
323            if (numberRows>25000||numberColumns>5*numberRows) {
324              nPasses = 50;
325            } else if (numberColumns>4*numberRows) {
326              nPasses = 20;
327            } else {
328              nPasses=5;
329            }
330          } else {
331            if (numberRows>25000||numberColumns>5*numberRows) {
332              nPasses = 50;
333              info.setLightweight(0); // say not lightweight idiot
334            } else if (numberColumns>4*numberRows) {
335              nPasses = 20;
336            } else {
337              nPasses=15;
338            }
339          }
340          if (numberElements<3*numberColumns) { 
341            nPasses=(int) ((2.0*(double) nPasses)/ratio); // probably not worh it
342          } else {
343            nPasses = max(nPasses,5);
344            nPasses = (int) (((double) nPasses)*5.0/ratio); // reduce if lots of elements per column
345          }
346          if (numberRows>25000&&nPasses>5) {
347            // Might as well go for it
348            nPasses = max(nPasses,71);
349          } else if (plusMinus) {
350            nPasses *= 2;
351            nPasses=min(nPasses,71);
352          }
353          if (nPasses<=5)
354            nPasses=0;
355          //info.setStartingWeight(1.0e-1);
356        }
357      }
358      if (doIdiot>0) {
359        // pick up number passes
360        nPasses=options.getExtraInfo(1);
361        if (nPasses>70) {
362          info.setStartingWeight(1.0e3);
363          info.setReduceIterations(6);
364        }
365      }
366      if (nPasses) {
367        doCrash=0;
368        info.crash(nPasses,model2->messageHandler(),model2->messagesPointer());
369        time2 = CoinCpuTime();
370        timeIdiot = time2-timeX;
371        handler_->message(CLP_INTERVAL_TIMING,messages_)
372          <<"Idiot Crash"<<timeIdiot<<time2-time1
373          <<CoinMessageEol;
374        timeX=time2;
375      }
376    }
377#endif
378    // ?
379    if (doCrash)
380      model2->crash(1000,1);
381    model2->primal(2);
382    time2 = CoinCpuTime();
383    timeCore = time2-timeX;
384    timeSimplex = timeCore;
385    handler_->message(CLP_INTERVAL_TIMING,messages_)
386      <<"Primal"<<timeCore<<time2-time1
387      <<CoinMessageEol;
388    timeX=time2;
389  } else if (method==ClpSolve::usePrimalorSprint) {
390    // Sprint
391    /*
392      This driver implements what I called Sprint when I introduced the idea
393      many years ago.  Cplex calls it "sifting" which I think is just as silly.
394      When I thought of this trivial idea
395      it reminded me of an LP code of the 60's called sprint which after
396      every factorization took a subset of the matrix into memory (all
397      64K words!) and then iterated very fast on that subset.  On the
398      problems of those days it did not work very well, but it worked very
399      well on aircrew scheduling problems where there were very large numbers
400      of columns all with the same flavor.
401    */
402   
403    /* The idea works best if you can get feasible easily.  To make it
404       more general we can add in costed slacks */
405   
406    int originalNumberColumns = model2->numberColumns();
407    int numberRows = model2->numberRows();
408   
409    // We will need arrays to choose variables.  These are too big but ..
410    double * weight = new double [numberRows+originalNumberColumns];
411    int * sort = new int [numberRows+originalNumberColumns];
412    int numberSort=0;
413    // We are going to add slacks to get feasible.
414    // initial list will just be artificials
415    // first we will set all variables as close to zero as possible
416    int iColumn;
417    const double * columnLower = model2->columnLower();
418    const double * columnUpper = model2->columnUpper();
419    double * columnSolution = model2->primalColumnSolution();
420   
421    for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
422      double value =0.0;
423      if (columnLower[iColumn]>0.0)
424        value = columnLower[iColumn];
425      else if (columnUpper[iColumn]<0.0)
426        value = columnUpper[iColumn];
427      columnSolution[iColumn]=value;
428    }
429    // now see what that does to row solution
430    double * rowSolution = model2->primalRowSolution();
431    memset (rowSolution,0,numberRows*sizeof(double));
432    model2->times(1.0,columnSolution,rowSolution);
433   
434    int * addStarts = new int [numberRows+1];
435    int * addRow = new int[numberRows];
436    double * addElement = new double[numberRows];
437    const double * lower = model2->rowLower();
438    const double * upper = model2->rowUpper();
439    addStarts[0]=0;
440    int numberArtificials=0;
441    double * addCost = new double [numberRows];
442    const double penalty=1.0e8;
443    int iRow;
444    for (iRow=0;iRow<numberRows;iRow++) {
445      if (lower[iRow]>rowSolution[iRow]) {
446        addRow[numberArtificials]=iRow;
447        addElement[numberArtificials]=1.0;
448        addCost[numberArtificials]=penalty;
449        numberArtificials++;
450        addStarts[numberArtificials]=numberArtificials;
451      } else if (upper[iRow]<rowSolution[iRow]) {
452        addRow[numberArtificials]=iRow;
453        addElement[numberArtificials]=-1.0;
454        addCost[numberArtificials]=penalty;
455        numberArtificials++;
456        addStarts[numberArtificials]=numberArtificials;
457      }
458    }
459    model2->addColumns(numberArtificials,NULL,NULL,addCost,
460                       addStarts,addRow,addElement);
461    delete [] addStarts;
462    delete [] addRow;
463    delete [] addElement;
464    delete [] addCost;
465    // look at rhs to see if to perturb
466    double largest=0.0;
467    double smallest = 1.0e30;
468    for (iRow=0;iRow<numberRows;iRow++) {
469      double value;
470      value = fabs(model2->rowLower_[iRow]);
471      if (value&&value<1.0e30) {
472        largest = max(largest,value);
473        smallest=min(smallest,value);
474      }
475      value = fabs(model2->rowUpper_[iRow]);
476      if (value&&value<1.0e30) {
477        largest = max(largest,value);
478        smallest=min(smallest,value);
479      }
480    }
481    double * saveLower = NULL;
482    double * saveUpper = NULL;
483    if (largest<2.01*smallest) {
484      // perturb - so switch off standard
485      model2->setPerturbation(100);
486      saveLower = new double[numberRows];
487      memcpy(saveLower,model2->rowLower_,numberRows*sizeof(double));
488      saveUpper = new double[numberRows];
489      memcpy(saveUpper,model2->rowUpper_,numberRows*sizeof(double));
490      double * lower = model2->rowLower();
491      double * upper = model2->rowUpper();
492      for (iRow=0;iRow<numberRows;iRow++) {
493        double lowerValue=lower[iRow], upperValue=upper[iRow];
494        double value = CoinDrand48();
495        if (upperValue>lowerValue+primalTolerance_) {
496          if (lowerValue>-1.0e20&&lowerValue)
497            lowerValue -= value * 1.0e-4*fabs(lowerValue); 
498          if (upperValue<1.0e20&&upperValue)
499            upperValue += value * 1.0e-4*fabs(upperValue); 
500        } else if (upperValue>0.0) {
501          upperValue -= value * 1.0e-4*fabs(lowerValue); 
502          lowerValue -= value * 1.0e-4*fabs(lowerValue); 
503        } else if (upperValue<0.0) {
504          upperValue += value * 1.0e-4*fabs(lowerValue); 
505          lowerValue += value * 1.0e-4*fabs(lowerValue); 
506        } else {
507        }
508        lower[iRow]=lowerValue;
509        upper[iRow]=upperValue;
510      }
511    }
512    int i;
513    // Set up initial list
514    if (numberArtificials) {
515      numberSort=numberArtificials;
516      for (i=0;i<numberSort;i++)
517        sort[i] = i+originalNumberColumns;
518    } else {
519      numberSort = min(numberRows_,numberColumns_);
520      for (i=0;i<numberSort;i++)
521        sort[i] = i;
522    }
523   
524    // redo as will have changed
525    columnLower = model2->columnLower();
526    columnUpper = model2->columnUpper();
527    int numberColumns = model2->numberColumns();
528    double * fullSolution = model2->primalColumnSolution();
529   
530    // Just do this number of passes
531    int maxPass=100;
532    if (doSprint>0)
533      maxPass=options.getExtraInfo(1);
534    int iPass;
535    double lastObjective=1.0e31;
536    // It will be safe to allow dense
537    model2->setInitialDenseFactorization(true);
538   
539    // Just take this number of columns in small problem
540    int smallNumberColumns = min(3*numberRows,numberColumns);
541    smallNumberColumns = max(smallNumberColumns,3000);
542    //int smallNumberColumns = min(12*numberRows/10,numberColumns);
543    //smallNumberColumns = max(smallNumberColumns,3000);
544    //smallNumberColumns = max(smallNumberColumns,numberRows+1000);
545    // We will be using all rows
546    int * whichRows = new int [numberRows];
547    for (int iRow=0;iRow<numberRows;iRow++)
548      whichRows[iRow]=iRow;
549    double originalOffset;
550    model2->getDblParam(ClpObjOffset,originalOffset);
551    int totalIterations=0;
552    for (iPass=0;iPass<maxPass;iPass++) {
553      //printf("Bug until submodel new version\n");
554      //CoinSort_2(sort,sort+numberSort,weight);
555      // Create small problem
556      ClpSimplex small(model2,numberRows,whichRows,numberSort,sort);
557      small.setPerturbation(model2->perturbation());
558      // now see what variables left out do to row solution
559      double * rowSolution = model2->primalRowSolution();
560      double * sumFixed = new double[numberRows];
561      memset (sumFixed,0,numberRows*sizeof(double));
562      int iRow,iColumn;
563      // zero out ones in small problem
564      for (iColumn=0;iColumn<numberSort;iColumn++) {
565        int kColumn = sort[iColumn];
566        fullSolution[kColumn]=0.0;
567      }
568      // Get objective offset
569      double offset=0.0;
570      const double * objective = model2->objective();
571      for (iColumn=0;iColumn<numberColumns;iColumn++) 
572        offset += fullSolution[iColumn]*objective[iColumn];
573      small.setDblParam(ClpObjOffset,originalOffset-offset);
574      model2->times(1.0,fullSolution,sumFixed);
575     
576      double * lower = small.rowLower();
577      double * upper = small.rowUpper();
578      for (iRow=0;iRow<numberRows;iRow++) {
579        if (lower[iRow]>-1.0e50) 
580          lower[iRow] -= sumFixed[iRow];
581        if (upper[iRow]<1.0e50)
582          upper[iRow] -= sumFixed[iRow];
583        rowSolution[iRow] -= sumFixed[iRow];
584      }
585      delete [] sumFixed;
586      // Solve
587      if (interrupt)
588        currentModel = &small;
589      small.primal();
590      totalIterations += small.numberIterations();
591      // move solution back
592      const double * solution = small.primalColumnSolution();
593      for (iColumn=0;iColumn<numberSort;iColumn++) {
594        int kColumn = sort[iColumn];
595        model2->setColumnStatus(kColumn,small.getColumnStatus(iColumn));
596        fullSolution[kColumn]=solution[iColumn];
597      }
598      for (iRow=0;iRow<numberRows;iRow++) 
599        model2->setRowStatus(iRow,small.getRowStatus(iRow));
600      memcpy(model2->primalRowSolution(),small.primalRowSolution(),
601             numberRows*sizeof(double));
602      // get reduced cost for large problem
603      memcpy(weight,model2->objective(),numberColumns*sizeof(double));
604      model2->transposeTimes(-1.0,small.dualRowSolution(),weight);
605      int numberNegative=0;
606      double sumNegative = 0.0;
607      // now massage weight so all basic in plus good djs
608      for (iColumn=0;iColumn<numberColumns;iColumn++) {
609        double dj = weight[iColumn]*optimizationDirection_;
610        double value = fullSolution[iColumn];
611        if (model2->getColumnStatus(iColumn)==ClpSimplex::basic) 
612          dj = -1.0e50;
613        else if (dj<0.0&&value<columnUpper[iColumn])
614          dj = dj;
615        else if (dj>0.0&&value>columnLower[iColumn])
616          dj = -dj;
617        else if (columnUpper[iColumn]>columnLower[iColumn])
618          dj = fabs(dj);
619        else
620          dj = 1.0e50;
621        weight[iColumn] = dj;
622        if (dj<-dualTolerance_&&dj>-1.0e50) {
623          numberNegative++;
624          sumNegative -= dj;
625        }
626        sort[iColumn] = iColumn;
627      }
628      handler_->message(CLP_SPRINT,messages_)
629        <<iPass+1<<small.numberIterations()<<small.objectiveValue()<<sumNegative
630        <<numberNegative
631        <<CoinMessageEol;
632      if ((small.objectiveValue()*optimizationDirection_>lastObjective-1.0e-7&&iPass>5)||
633          !small.numberIterations()||
634          iPass==maxPass-1||small.status()==3) {
635       
636        break; // finished
637      } else {
638        lastObjective = small.objectiveValue()*optimizationDirection_;
639        // sort
640        CoinSort_2(weight,weight+numberColumns,sort);
641        numberSort = smallNumberColumns;
642      }
643    }
644    if (interrupt) 
645      currentModel = model2;
646    for (i=0;i<numberArtificials;i++)
647      sort[i] = i + originalNumberColumns;
648    model2->deleteColumns(numberArtificials,sort);
649    delete [] weight;
650    delete [] sort;
651    delete [] whichRows;
652    if (saveLower) {
653      // unperturb and clean
654      for (iRow=0;iRow<numberRows;iRow++) {
655        double diffLower = saveLower[iRow]-model2->rowLower_[iRow];
656        double diffUpper = saveUpper[iRow]-model2->rowUpper_[iRow];
657        model2->rowLower_[iRow]=saveLower[iRow];
658        model2->rowUpper_[iRow]=saveUpper[iRow];
659        if (diffLower) 
660          assert (!diffUpper||fabs(diffLower-diffUpper)<1.0e-5);
661        else
662          diffLower = diffUpper;
663        model2->rowActivity_[iRow] += diffLower;
664      }
665      delete [] saveLower;
666      delete [] saveUpper;
667    }
668    model2->primal(0);
669    model2->setPerturbation(savePerturbation);
670    time2 = CoinCpuTime();
671    timeCore = time2-timeX;
672    timeSimplex = timeCore;
673    handler_->message(CLP_INTERVAL_TIMING,messages_)
674      <<"Sprint"<<timeCore<<time2-time1
675      <<CoinMessageEol;
676    timeX=time2;
677    model2->setNumberIterations(model2->numberIterations()+totalIterations);
678  } else {
679    assert (method!=ClpSolve::automatic); // later
680    time2=0.0;
681  }
682  if (saveMatrix) {
683    // delete and replace
684    delete model2->clpMatrix();
685    model2->replaceMatrix(saveMatrix);
686  }
687  numberIterations = model2->numberIterations();
688  finalStatus=model2->status();
689  if (presolve==ClpSolve::presolveOn) {
690    int saveLevel = logLevel();
691    setLogLevel(1);
692    pinfo.postsolve(true);
693    time2 = CoinCpuTime();
694    timePresolve += time2-timeX;
695    handler_->message(CLP_INTERVAL_TIMING,messages_)
696      <<"Postsolve"<<time2-timeX<<time2-time1
697      <<CoinMessageEol;
698    timeX=time2;
699    delete model2;
700    if (interrupt)
701      currentModel = this;
702    checkSolution();
703    setLogLevel(saveLevel);
704    if (finalStatus!=3&&(finalStatus||status()==-1)) {
705      int savePerturbation = perturbation();
706      setPerturbation(100);
707      primal(1);
708      setPerturbation(savePerturbation);
709      numberIterations += this->numberIterations();
710      finalStatus=status();
711      time2 = CoinCpuTime();
712      timeSimplex += time2-timeX;
713      handler_->message(CLP_INTERVAL_TIMING,messages_)
714      <<"Cleanup"<<time2-timeX<<time2-time1
715      <<CoinMessageEol;
716      timeX=time2;
717    }
718  }
719  setMaximumIterations(saveMaxIterations);
720  std::string statusMessage[]={"Unknown","Optimal","PrimalInfeasible","DualInfeasible","Stopped"};
721  assert (finalStatus>=-1&&finalStatus<=3);
722  handler_->message(CLP_TIMING,messages_)
723    <<statusMessage[finalStatus+1]<<objectiveValue()<<numberIterations<<time2-time1;
724  handler_->printing(presolve==ClpSolve::presolveOn)
725    <<timePresolve;
726  handler_->printing(timeIdiot)
727    <<timeIdiot;
728  handler_->message()<<CoinMessageEol;
729  if (interrupt) 
730    signal(SIGINT,saveSignal);
731  perturbation_=savePerturbation;
732  scalingFlag_=saveScaling;
733  return finalStatus;
734}
735// General solve
736int 
737ClpSimplex::initialSolve()
738{
739  // Default so use dual
740  ClpSolve options;
741  return initialSolve(options);
742}
743// General dual solve
744int 
745ClpSimplex::initialDualSolve()
746{
747  ClpSolve options;
748  // Use dual
749  options.setSolveType(ClpSolve::useDual);
750  return initialSolve(options);
751}
752// General dual solve
753int 
754ClpSimplex::initialPrimalSolve()
755{
756  ClpSolve options;
757  // Use primal
758  options.setSolveType(ClpSolve::usePrimal);
759  return initialSolve(options);
760}
761// Default constructor
762ClpSolve::ClpSolve (  )
763{
764  method_ = useDual;
765  presolveType_=presolveOn;
766  numberPasses_=5;
767  for (int i=0;i<4;i++)
768    options_[i]=0;
769  for (int i=0;i<4;i++)
770    extraInfo_[i]=-1;
771}
772
773// Copy constructor.
774ClpSolve::ClpSolve(const ClpSolve & rhs)
775{
776  method_ = rhs.method_;
777  presolveType_=rhs.presolveType_;
778  numberPasses_=rhs.numberPasses_;
779  for (int i=0;i<4;i++)
780    options_[i]=rhs.options_[i];
781  for (int i=0;i<4;i++)
782    extraInfo_[i]=rhs.extraInfo_[i];
783}
784// Assignment operator. This copies the data
785ClpSolve & 
786ClpSolve::operator=(const ClpSolve & rhs)
787{
788  if (this != &rhs) {
789    method_ = rhs.method_;
790    presolveType_=rhs.presolveType_;
791    numberPasses_=rhs.numberPasses_;
792    for (int i=0;i<4;i++)
793      options_[i]=rhs.options_[i];
794    for (int i=0;i<4;i++)
795      extraInfo_[i]=rhs.extraInfo_[i];
796  }
797  return *this;
798
799}
800// Destructor
801ClpSolve::~ClpSolve (  )
802{
803}
804/*   which translation is:
805     which:
806     0 - startup in Dual  (nothing if basis exists).:
807             0 - no basis, 1 crash
808     1 - startup in Primal (nothing if basis exists):
809        0 - use initiative
810        1 - use crash
811        2 - use idiot and look at further info
812        3 - use sprint and look at further info
813        4 - use all slack
814        5 - use initiative but no idiot
815        6 - use initiative but no sprint
816        7 - use initiative but no crash
817        8 - do allslack or idiot
818        9 - do allslack or sprint
819     2 - interrupt handling - 0 yes, 1 no (for threadsafe)
820     3 - whether to make +- 1matrix - 0 yes, 1 no
821*/
822void 
823ClpSolve::setSpecialOption(int which,int value,int extraInfo)
824{
825  options_[which]=value;
826  extraInfo_[which]=extraInfo;
827}
828int 
829ClpSolve::getSpecialOption(int which) const
830{
831  return options_[which];
832}
833
834// Solve types
835void 
836ClpSolve::setSolveType(SolveType method, int extraInfo)
837{
838  method_=method;
839}
840
841ClpSolve::SolveType
842ClpSolve::getSolveType()
843{
844  return method_;
845}
846
847// Presolve types
848void 
849ClpSolve::setPresolveType(PresolveType amount, int extraInfo)
850{
851  presolveType_=amount;
852  numberPasses_=extraInfo;
853}
854ClpSolve::PresolveType
855ClpSolve::getPresolveType()
856{
857  return presolveType_;
858}
859// Extra info for idiot (or sprint)
860int 
861ClpSolve::getExtraInfo(int which) const
862{
863  return extraInfo_[which];
864}
865int 
866ClpSolve::getPresolvePasses() const
867{
868  return numberPasses_;
869}
Note: See TracBrowser for help on using the repository browser.