source: trunk/ClpSolve.cpp @ 225

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

This should break everything

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