source: branches/pre/ClpSolve.cpp @ 221

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

Still trying to go faster

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.4 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    // We will be using all rows
543    int * whichRows = new int [numberRows];
544    for (int iRow=0;iRow<numberRows;iRow++)
545      whichRows[iRow]=iRow;
546    double originalOffset;
547    model2->getDblParam(ClpObjOffset,originalOffset);
548    int totalIterations=0;
549    for (iPass=0;iPass<maxPass;iPass++) {
550      //printf("Bug until submodel new version\n");
551      //CoinSort_2(sort,sort+numberSort,weight);
552      // Create small problem
553      ClpSimplex small(model2,numberRows,whichRows,numberSort,sort);
554      small.setPerturbation(model2->perturbation());
555      // now see what variables left out do to row solution
556      double * rowSolution = model2->primalRowSolution();
557      double * sumFixed = new double[numberRows];
558      memset (sumFixed,0,numberRows*sizeof(double));
559      int iRow,iColumn;
560      // zero out ones in small problem
561      for (iColumn=0;iColumn<numberSort;iColumn++) {
562        int kColumn = sort[iColumn];
563        fullSolution[kColumn]=0.0;
564      }
565      // Get objective offset
566      double offset=0.0;
567      const double * objective = model2->objective();
568      for (iColumn=0;iColumn<numberColumns;iColumn++) 
569        offset += fullSolution[iColumn]*objective[iColumn];
570      small.setDblParam(ClpObjOffset,originalOffset-offset);
571      model2->times(1.0,fullSolution,sumFixed);
572     
573      double * lower = small.rowLower();
574      double * upper = small.rowUpper();
575      for (iRow=0;iRow<numberRows;iRow++) {
576        if (lower[iRow]>-1.0e50) 
577          lower[iRow] -= sumFixed[iRow];
578        if (upper[iRow]<1.0e50)
579          upper[iRow] -= sumFixed[iRow];
580        rowSolution[iRow] -= sumFixed[iRow];
581      }
582      delete [] sumFixed;
583      // Solve
584      if (interrupt)
585        currentModel = &small;
586      small.primal();
587      totalIterations += small.numberIterations();
588      // move solution back
589      const double * solution = small.primalColumnSolution();
590      for (iColumn=0;iColumn<numberSort;iColumn++) {
591        int kColumn = sort[iColumn];
592        model2->setColumnStatus(kColumn,small.getColumnStatus(iColumn));
593        fullSolution[kColumn]=solution[iColumn];
594      }
595      for (iRow=0;iRow<numberRows;iRow++) 
596        model2->setRowStatus(iRow,small.getRowStatus(iRow));
597      memcpy(model2->primalRowSolution(),small.primalRowSolution(),
598             numberRows*sizeof(double));
599      // get reduced cost for large problem
600      memcpy(weight,model2->objective(),numberColumns*sizeof(double));
601      model2->transposeTimes(-1.0,small.dualRowSolution(),weight);
602      int numberNegative=0;
603      double sumNegative = 0.0;
604      // now massage weight so all basic in plus good djs
605      for (iColumn=0;iColumn<numberColumns;iColumn++) {
606        double dj = weight[iColumn]*optimizationDirection_;
607        double value = fullSolution[iColumn];
608        if (model2->getColumnStatus(iColumn)==ClpSimplex::basic) 
609          dj = -1.0e50;
610        else if (dj<0.0&&value<columnUpper[iColumn])
611          dj = dj;
612        else if (dj>0.0&&value>columnLower[iColumn])
613          dj = -dj;
614        else if (columnUpper[iColumn]>columnLower[iColumn])
615          dj = fabs(dj);
616        else
617          dj = 1.0e50;
618        weight[iColumn] = dj;
619        if (dj<-dualTolerance_&&dj>-1.0e50) {
620          numberNegative++;
621          sumNegative -= dj;
622        }
623        sort[iColumn] = iColumn;
624      }
625      handler_->message(CLP_SPRINT,messages_)
626        <<iPass+1<<small.numberIterations()<<small.objectiveValue()<<sumNegative
627        <<numberNegative
628        <<CoinMessageEol;
629      if ((small.objectiveValue()*optimizationDirection_>lastObjective-1.0e-7&&iPass>5)||
630          !small.numberIterations()||
631          iPass==maxPass-1||small.status()==3) {
632       
633        break; // finished
634      } else {
635        lastObjective = small.objectiveValue()*optimizationDirection_;
636        // sort
637        CoinSort_2(weight,weight+numberColumns,sort);
638        numberSort = smallNumberColumns;
639      }
640    }
641    if (interrupt) 
642      currentModel = model2;
643    for (i=0;i<numberArtificials;i++)
644      sort[i] = i + originalNumberColumns;
645    model2->deleteColumns(numberArtificials,sort);
646    delete [] weight;
647    delete [] sort;
648    delete [] whichRows;
649    if (saveLower) {
650      // unperturb and clean
651      for (iRow=0;iRow<numberRows;iRow++) {
652        double diffLower = saveLower[iRow]-model2->rowLower_[iRow];
653        double diffUpper = saveUpper[iRow]-model2->rowUpper_[iRow];
654        model2->rowLower_[iRow]=saveLower[iRow];
655        model2->rowUpper_[iRow]=saveUpper[iRow];
656        if (diffLower) 
657          assert (!diffUpper||fabs(diffLower-diffUpper)<1.0e-5);
658        else
659          diffLower = diffUpper;
660        model2->rowActivity_[iRow] += diffLower;
661      }
662      delete [] saveLower;
663      delete [] saveUpper;
664    }
665    model2->primal(0);
666    model2->setPerturbation(savePerturbation);
667    time2 = CoinCpuTime();
668    timeCore = time2-timeX;
669    timeSimplex = timeCore;
670    handler_->message(CLP_INTERVAL_TIMING,messages_)
671      <<"Sprint"<<timeCore<<time2-time1
672      <<CoinMessageEol;
673    timeX=time2;
674    model2->setNumberIterations(model2->numberIterations()+totalIterations);
675  } else {
676    assert (method!=ClpSolve::automatic); // later
677    time2=0.0;
678  }
679  if (saveMatrix) {
680    // delete and replace
681    delete model2->clpMatrix();
682    model2->replaceMatrix(saveMatrix);
683  }
684  numberIterations = model2->numberIterations();
685  finalStatus=model2->status();
686  if (presolve==ClpSolve::presolveOn) {
687    int saveLevel = logLevel();
688    setLogLevel(1);
689    pinfo.postsolve(true);
690    time2 = CoinCpuTime();
691    timePresolve += time2-timeX;
692    handler_->message(CLP_INTERVAL_TIMING,messages_)
693      <<"Postsolve"<<time2-timeX<<time2-time1
694      <<CoinMessageEol;
695    timeX=time2;
696    delete model2;
697    if (interrupt)
698      currentModel = this;
699    checkSolution();
700    setLogLevel(saveLevel);
701    if (finalStatus!=3&&(finalStatus||status()==-1)) {
702      int savePerturbation = perturbation();
703      setPerturbation(100);
704      primal(1);
705      setPerturbation(savePerturbation);
706      numberIterations += this->numberIterations();
707      finalStatus=status();
708      time2 = CoinCpuTime();
709      timeSimplex += time2-timeX;
710      handler_->message(CLP_INTERVAL_TIMING,messages_)
711      <<"Cleanup"<<time2-timeX<<time2-time1
712      <<CoinMessageEol;
713      timeX=time2;
714    }
715  }
716  setMaximumIterations(saveMaxIterations);
717  std::string statusMessage[]={"Unknown","Optimal","PrimalInfeasible","DualInfeasible","Stopped"};
718  assert (finalStatus>=-1&&finalStatus<=3);
719  handler_->message(CLP_TIMING,messages_)
720    <<statusMessage[finalStatus+1]<<objectiveValue()<<numberIterations<<time2-time1;
721  handler_->printing(presolve==ClpSolve::presolveOn)
722    <<timePresolve;
723  handler_->printing(timeIdiot)
724    <<timeIdiot;
725  handler_->message()<<CoinMessageEol;
726  if (interrupt) 
727    signal(SIGINT,saveSignal);
728  perturbation_=savePerturbation;
729  scalingFlag_=saveScaling;
730  return finalStatus;
731}
732// Default constructor
733ClpSolve::ClpSolve (  )
734{
735  method_ = useDual;
736  presolveType_=presolveOn;
737  numberPasses_=5;
738  for (int i=0;i<4;i++)
739    options_[i]=0;
740  for (int i=0;i<4;i++)
741    extraInfo_[i]=-1;
742}
743
744// Copy constructor.
745ClpSolve::ClpSolve(const ClpSolve & rhs)
746{
747  method_ = rhs.method_;
748  presolveType_=rhs.presolveType_;
749  numberPasses_=rhs.numberPasses_;
750  for (int i=0;i<4;i++)
751    options_[i]=rhs.options_[i];
752  for (int i=0;i<4;i++)
753    extraInfo_[i]=rhs.extraInfo_[i];
754}
755// Assignment operator. This copies the data
756ClpSolve & 
757ClpSolve::operator=(const ClpSolve & rhs)
758{
759  if (this != &rhs) {
760    method_ = rhs.method_;
761    presolveType_=rhs.presolveType_;
762    numberPasses_=rhs.numberPasses_;
763    for (int i=0;i<4;i++)
764      options_[i]=rhs.options_[i];
765    for (int i=0;i<4;i++)
766      extraInfo_[i]=rhs.extraInfo_[i];
767  }
768  return *this;
769
770}
771// Destructor
772ClpSolve::~ClpSolve (  )
773{
774}
775/*   which translation is:
776     which:
777     0 - startup in Dual  (nothing if basis exists).:
778             0 - no basis, 1 crash
779     1 - startup in Primal (nothing if basis exists):
780        0 - use initiative
781        1 - use crash
782        2 - use idiot and look at further info
783        3 - use sprint and look at further info
784        4 - use all slack
785        5 - use initiative but no idiot
786        6 - use initiative but no sprint
787        7 - use initiative but no crash
788        8 - do allslack or idiot
789        9 - do allslack or sprint
790     2 - interrupt handling - 0 yes, 1 no (for threadsafe)
791     3 - whether to make +- 1matrix - 0 yes, 1 no
792*/
793void 
794ClpSolve::setSpecialOption(int which,int value,int extraInfo)
795{
796  options_[which]=value;
797  extraInfo_[which]=extraInfo;
798}
799int 
800ClpSolve::getSpecialOption(int which) const
801{
802  return options_[which];
803}
804
805// Solve types
806void 
807ClpSolve::setSolveType(SolveType method, int extraInfo)
808{
809  method_=method;
810}
811
812ClpSolve::SolveType
813ClpSolve::getSolveType()
814{
815  return method_;
816}
817
818// Presolve types
819void 
820ClpSolve::setPresolveType(PresolveType amount, int extraInfo)
821{
822  presolveType_=amount;
823  numberPasses_=extraInfo;
824}
825ClpSolve::PresolveType
826ClpSolve::getPresolveType()
827{
828  return presolveType_;
829}
830// Extra info for idiot (or sprint)
831int 
832ClpSolve::getExtraInfo(int which) const
833{
834  return extraInfo_[which];
835}
836int 
837ClpSolve::getPresolvePasses() const
838{
839  return numberPasses_;
840}
Note: See TracBrowser for help on using the repository browser.