source: trunk/ClpSolve.cpp @ 340

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

static_cast

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 37.7 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 "ClpHelperFunctions.hpp"
13#include "CoinSort.hpp"
14#include "ClpFactorization.hpp"
15#include "ClpSimplex.hpp"
16#include "ClpInterior.hpp"
17#ifdef REAL_BARRIER
18#include "ClpCholeskyWssmp.hpp"
19//#include "ClpCholeskyWssmpKKT.hpp"
20//#include "ClpCholeskyTaucs.hpp"
21#endif
22#include "ClpSolve.hpp"
23#include "ClpPackedMatrix.hpp"
24#include "ClpPlusMinusOneMatrix.hpp"
25#include "ClpNetworkMatrix.hpp"
26#include "ClpMessage.hpp"
27#include "CoinTime.hpp"
28
29#include "ClpPresolve.hpp"
30#include "Idiot.hpp"
31
32//#############################################################################
33// Allow for interrupts
34// But is this threadsafe ? (so switched off by option)
35
36#include "CoinSignal.hpp"
37static ClpSimplex * currentModel = NULL;
38static ClpInterior * currentModel2 = NULL;
39
40extern "C" {
41   static void signal_handler(int whichSignal)
42   {
43      if (currentModel!=NULL) 
44         currentModel->setMaximumIterations(0); // stop at next iterations
45      if (currentModel2!=NULL) 
46         currentModel2->setMaximumBarrierIterations(0); // stop at next iterations
47      return;
48   }
49}
50
51/** General solve algorithm which can do presolve
52    special options (bits)
53    1 - do not perturb
54    2 - do not scale
55    4 - use crash (default allslack in dual, idiot in primal)
56    8 - all slack basis in primal
57    16 - switch off interrupt handling
58    32 - do not try and make plus minus one matrix
59 */
60int 
61ClpSimplex::initialSolve(ClpSolve & options)
62{
63  ClpSolve::SolveType method=options.getSolveType();
64  ClpSolve::PresolveType presolve = options.getPresolveType();
65  int saveMaxIterations = maximumIterations();
66  int finalStatus=-1;
67  int numberIterations=0;
68  double time1 = CoinCpuTime();
69  double timeX = time1;
70  double time2;
71  ClpMatrixBase * saveMatrix=NULL;
72  ClpSimplex * model2 = this;
73  bool interrupt = (options.getSpecialOption(2)==0);
74  CoinSighandler_t saveSignal=SIG_DFL;
75  if (interrupt) {
76    currentModel = model2;
77    // register signal handler
78    saveSignal = signal(SIGINT,signal_handler);
79  }
80  ClpPresolve pinfo;
81  double timePresolve=0.0;
82  double timeIdiot=0.0;
83  double timeCore=0.0;
84  int savePerturbation=perturbation_;
85  int saveScaling = scalingFlag_;
86#ifndef NO_RTTI
87  if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
88    // network - switch off stuff
89    presolve = ClpSolve::presolveOff;
90  }
91#else
92  if (matrix_->type()==11) {
93    // network - switch off stuff
94    presolve = ClpSolve::presolveOff;
95  }
96#endif
97  // For below >0 overrides
98  // 0 means no, -1 means maybe
99  int doIdiot=0;
100  int doCrash=0;
101  int doSprint=0;
102  if (method!=ClpSolve::useDual) {
103    switch (options.getSpecialOption(1)) {
104    case 0:
105      doIdiot=-1;
106      doCrash=-1;
107      doSprint=-1;
108      break;
109    case 1:
110      doIdiot=0;
111      doCrash=1;
112      doSprint=0;
113      break;
114    case 2:
115      doIdiot=1;
116      doCrash=0;
117      doSprint=0;
118      break;
119    case 3:
120      doIdiot=0;
121      doCrash=0;
122      doSprint=1;
123      break;
124    case 4:
125      doIdiot=0;
126      doCrash=0;
127      doSprint=0;
128      break;
129    case 5:
130      doIdiot=0;
131      doCrash=-1;
132      doSprint=-1;
133      break;
134    case 6:
135      doIdiot=-1;
136      doCrash=-1;
137      doSprint=0;
138      break;
139    case 7:
140      doIdiot=-1;
141      doCrash=0;
142      doSprint=-1;
143      break;
144    case 8:
145      doIdiot=-1;
146      doCrash=0;
147      doSprint=0;
148      break;
149    case 9:
150      doIdiot=0;
151      doCrash=0;
152      doSprint=-1;
153      break;
154    default:
155      abort();
156    }
157  } else {
158    // Dual
159    switch (options.getSpecialOption(0)) {
160    case 0:
161      doIdiot=0;
162      doCrash=0;
163      doSprint=0;
164      break;
165    case 1:
166      doIdiot=0;
167      doCrash=1;
168      doSprint=0;
169      break;
170    case 2:
171      doIdiot=-1;
172      if (options.getExtraInfo(0)>0)
173        doIdiot = options.getExtraInfo(0);
174      doCrash=0;
175      doSprint=0;
176      break;
177    default:
178      abort();
179    }
180  }
181  // Just do this number of passes in Sprint
182  int maxSprintPass=100;
183  // PreSolve to file - not fully tested
184  bool presolveToFile=false;
185
186  if (presolve!=ClpSolve::presolveOff) {
187    int numberPasses=5;
188    if (presolve==ClpSolve::presolveNumber) {
189      numberPasses=options.getPresolvePasses();
190      presolve = ClpSolve::presolveOn;
191    }
192    if (presolveToFile) {
193      printf("***** temp test\n");
194      pinfo.presolvedModelToFile(*this,"ss.sav",1.0e-8,
195                           false,numberPasses);
196      model2=this;
197    } else {
198      model2 = pinfo.presolvedModel(*this,1.0e-8,
199                                    false,numberPasses,true);
200    }
201    time2 = CoinCpuTime();
202    timePresolve = time2-timeX;
203    handler_->message(CLP_INTERVAL_TIMING,messages_)
204      <<"Presolve"<<timePresolve<<time2-time1
205      <<CoinMessageEol;
206    timeX=time2;
207    if (model2) {
208    } else {
209      handler_->message(CLP_INFEASIBLE,messages_)
210        <<CoinMessageEol;
211      model2 = this;
212      presolve=ClpSolve::presolveOff;
213    }
214    // We may be better off using original (but if dual leave because of bounds)
215    if (presolve!=ClpSolve::presolveOff&&
216        numberRows_<1.01*model2->numberRows_&&numberColumns_<1.01*model2->numberColumns_
217        &&method!=ClpSolve::useDual&&model2!=this) {
218      delete model2;
219      model2 = this;
220      presolve=ClpSolve::presolveOff;
221    }
222  }
223  if (interrupt)
224    currentModel = model2;
225  // See if worth trying +- one matrix
226  bool plusMinus=false;
227  int numberElements=model2->getNumElements();
228#ifndef NO_RTTI
229  if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
230    // network - switch off stuff
231    doIdiot=0;
232    doSprint=0;
233  }
234#else
235  if (matrix_->type()==11) {
236    // network - switch off stuff
237    doIdiot=0;
238    doSprint=0;
239  }
240#endif
241  int numberColumns = model2->numberColumns();
242  int numberRows = model2->numberRows();
243  // If not all slack basis - switch off all
244  int number=0;
245  int iRow;
246  for (iRow=0;iRow<numberRows;iRow++)
247    if (model2->getRowStatus(iRow)==basic)
248      number++;
249  if (number<numberRows) {
250    doIdiot=0;
251    doCrash=0;
252    doSprint=0;
253  }
254  if (options.getSpecialOption(3)==0) {
255    if(numberElements>100000)
256      plusMinus=true;
257    if(numberElements>10000&&(doIdiot||doSprint)) 
258      plusMinus=true;
259  }
260  if (plusMinus) {
261    saveMatrix = model2->clpMatrix();
262#ifndef NO_RTTI
263    ClpPackedMatrix* clpMatrix =
264      dynamic_cast< ClpPackedMatrix*>(saveMatrix);
265#else
266    ClpPackedMatrix* clpMatrix = NULL;
267    if (saveMatrix->type()==1)
268      clpMatrix =
269        static_cast< ClpPackedMatrix*>(saveMatrix);
270#endif
271    if (clpMatrix) {
272      ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
273      if (newMatrix->getIndices()) {
274        model2->replaceMatrix(newMatrix);
275      } else {
276        handler_->message(CLP_MATRIX_CHANGE,messages_)
277          <<"+- 1"
278          <<CoinMessageEol;
279        saveMatrix=NULL;
280        plusMinus=false;
281        delete newMatrix;
282      }
283    } else {
284      saveMatrix=NULL;
285      plusMinus=false;
286    }
287  }
288  if (model2->factorizationFrequency()==200) {
289    // User did not touch preset
290    model2->setFactorizationFrequency(min(2000,100+model2->numberRows()/200));
291  }
292  if (method==ClpSolve::usePrimalorSprint) {
293    if (doSprint<0) { 
294      if (numberElements<500000) {
295        // Small problem
296        if(numberRows*10>numberColumns||numberColumns<6000
297           ||(numberRows*20>numberColumns&&!plusMinus))
298          method=ClpSolve::usePrimal; // switch off sprint
299      } else {
300        // larger problem
301        if(numberRows*8>numberColumns)
302          method=ClpSolve::usePrimal; // switch off sprint
303        // but make lightweight
304        if(numberRows*10>numberColumns||numberColumns<6000
305           ||(numberRows*20>numberColumns&&!plusMinus))
306          maxSprintPass=5;
307      }
308    } else if (doSprint==0) {
309      method=ClpSolve::usePrimal; // switch off sprint
310    }
311  }
312  if (method==ClpSolve::useDual) {
313    // pick up number passes
314    int nPasses=0;
315    int numberNotE=0;
316    if ((doIdiot<0&&plusMinus)||doIdiot>0) {
317      // See if candidate for idiot
318      nPasses=0;
319      Idiot info(*model2);
320      // Get average number of elements per column
321      double ratio  = ((double) numberElements/(double) numberColumns);
322      // look at rhs
323      int iRow;
324      double largest=0.0;
325      double smallest = 1.0e30;
326      double largestGap=0.0;
327      for (iRow=0;iRow<numberRows;iRow++) {
328        double value1 = model2->rowLower_[iRow];
329        if (value1&&value1>-1.0e31) {
330          largest = max(largest,fabs(value1));
331          smallest=min(smallest,fabs(value1));
332        }
333        double value2 = model2->rowUpper_[iRow];
334        if (value2&&value2<1.0e31) {
335          largest = max(largest,fabs(value2));
336          smallest=min(smallest,fabs(value2));
337        }
338        if (value2>value1) {
339          numberNotE++;
340          if (value2>1.0e31||value1<-1.0e31)
341            largestGap = COIN_DBL_MAX;
342          else
343            largestGap = value2-value1;
344        }
345      }
346      if (doIdiot<0) {
347        if (numberRows>200&&numberColumns>5000&&ratio>=3.0&&
348            largest/smallest<1.1&&!numberNotE) {
349          nPasses = 71;
350        }
351      } 
352      if (doIdiot>0) {
353        nPasses=max(nPasses,doIdiot);
354        if (nPasses>70) 
355          info.setStartingWeight(1.0e3);
356      }
357      if (nPasses) {
358        info.setReduceIterations(5);
359        doCrash=0;
360        info.crash(nPasses,model2->messageHandler(),model2->messagesPointer());
361        time2 = CoinCpuTime();
362        timeIdiot = time2-timeX;
363        handler_->message(CLP_INTERVAL_TIMING,messages_)
364          <<"Idiot Crash"<<timeIdiot<<time2-time1
365          <<CoinMessageEol;
366        timeX=time2;
367      }
368    }
369    if (presolve==ClpSolve::presolveOn) {
370      int numberInfeasibilities = model2->tightenPrimalBounds();
371      if (numberInfeasibilities) {
372        handler_->message(CLP_INFEASIBLE,messages_)
373          <<CoinMessageEol;
374        model2 = this;
375        presolve=ClpSolve::presolveOff;
376      }
377    }
378    if (options.getSpecialOption(0)==1)
379      model2->crash(1000,1);
380    if (!nPasses) {
381      model2->dual(0);
382    } else if (!numberNotE&&0) {
383      // E so we can do in another way
384      double * pi = model2->dualRowSolution();
385      int i;
386      int numberColumns = model2->numberColumns();
387      int numberRows = model2->numberRows();
388      double * saveObj = new double[numberColumns];
389      memcpy(saveObj,model2->objective(),numberColumns*sizeof(double));
390      memcpy(model2->dualColumnSolution(),model2->objective(),
391             numberColumns*sizeof(double));
392      model2->clpMatrix()->transposeTimes(-1.0,pi,model2->dualColumnSolution());
393      memcpy(model2->objective(),model2->dualColumnSolution(),
394             numberColumns*sizeof(double));
395      const double * rowsol = model2->primalRowSolution();
396      double offset=0.0;
397      for (i=0;i<numberRows;i++) {
398        offset += pi[i]*rowsol[i];
399      }
400      double value2;
401      model2->getDblParam(ClpObjOffset,value2);
402      printf("Offset %g %g\n",offset,value2);
403      model2->setRowObjective(pi);
404      // zero out pi
405      memset(pi,0,numberRows*sizeof(double));
406      // Could put some in basis - only partially tested
407      model2->allSlackBasis(); 
408      model2->factorization()->maximumPivots(200);
409      //model2->setLogLevel(63);
410      // solve
411      model2->dual(1);
412      memcpy(model2->objective(),saveObj,numberColumns*sizeof(double));
413      // zero out pi
414      memset(pi,0,numberRows*sizeof(double));
415      model2->setRowObjective(pi);
416      delete [] saveObj;
417      model2->primal();
418    } else {
419      // solve
420      model2->dual(1);
421    }
422    time2 = CoinCpuTime();
423    timeCore = time2-timeX;
424    handler_->message(CLP_INTERVAL_TIMING,messages_)
425      <<"Dual"<<timeCore<<time2-time1
426      <<CoinMessageEol;
427    timeX=time2;
428  } else if (method==ClpSolve::usePrimal) {
429    if (doIdiot) {
430      int nPasses=0;
431      Idiot info(*model2);
432      // Get average number of elements per column
433      double ratio  = ((double) numberElements/(double) numberColumns);
434      // look at rhs
435      int iRow;
436      double largest=0.0;
437      double smallest = 1.0e30;
438      double largestGap=0.0;
439      int numberNotE=0;
440      for (iRow=0;iRow<numberRows;iRow++) {
441        double value1 = model2->rowLower_[iRow];
442        if (value1&&value1>-1.0e31) {
443          largest = max(largest,fabs(value1));
444          smallest=min(smallest,fabs(value1));
445        }
446        double value2 = model2->rowUpper_[iRow];
447        if (value2&&value2<1.0e31) {
448          largest = max(largest,fabs(value2));
449          smallest=min(smallest,fabs(value2));
450        }
451        if (value2>value1) {
452          numberNotE++;
453          if (value2>1.0e31||value1<-1.0e31)
454            largestGap = COIN_DBL_MAX;
455          else
456            largestGap = value2-value1;
457        }
458      }
459      if (numberRows>200&&numberColumns>5000&&numberColumns>2*numberRows) {
460        if (plusMinus) {
461          if (largest/smallest>2.0) {
462            nPasses = 10+numberColumns/100000;
463            nPasses = min(nPasses,50);
464            nPasses = max(nPasses,15);
465            if (numberRows>25000&&nPasses>5) {
466              // Might as well go for it
467              nPasses = max(nPasses,71);
468            } else if (numberElements<3*numberColumns) {
469              nPasses=min(nPasses,10); // probably not worh it
470            }
471            if (doIdiot<0)
472              info.setLightweight(1); // say lightweight idiot
473          } else if (largest/smallest>1.01||numberElements<=3*numberColumns) {
474            nPasses = 10+numberColumns/1000;
475            nPasses = min(nPasses,100);
476            nPasses = max(nPasses,30);
477            if (numberRows>25000) {
478              // Might as well go for it
479              nPasses = max(nPasses,71);
480            }
481            if (!largestGap)
482              nPasses *= 2;
483          } else {
484            nPasses = 10+numberColumns/1000;
485            nPasses = min(nPasses,200);
486            nPasses = max(nPasses,100);
487            info.setStartingWeight(1.0e-1);
488            info.setReduceIterations(6);
489            if (!largestGap)
490              nPasses *= 2;
491            //info.setFeasibilityTolerance(1.0e-7);
492          }
493          // If few passes - don't bother
494          if (nPasses<=5)
495            nPasses=0;
496        } else {
497          if (doIdiot<0)
498            info.setLightweight(1); // say lightweight idiot
499          if (largest/smallest>1.01||numberNotE) {
500            if (numberRows>25000||numberColumns>5*numberRows) {
501              nPasses = 50;
502            } else if (numberColumns>4*numberRows) {
503              nPasses = 20;
504            } else {
505              nPasses=5;
506            }
507          } else {
508            if (numberRows>25000||numberColumns>5*numberRows) {
509              nPasses = 50;
510              info.setLightweight(0); // say not lightweight idiot
511            } else if (numberColumns>4*numberRows) {
512              nPasses = 20;
513            } else {
514              nPasses=15;
515            }
516          }
517          if (numberElements<3*numberColumns) { 
518            nPasses=(int) ((2.0*(double) nPasses)/ratio); // probably not worh it
519          } else {
520            nPasses = max(nPasses,5);
521            nPasses = (int) (((double) nPasses)*5.0/ratio); // reduce if lots of elements per column
522          }
523          if (numberRows>25000&&nPasses>5) {
524            // Might as well go for it
525            nPasses = max(nPasses,71);
526          } else if (plusMinus) {
527            nPasses *= 2;
528            nPasses=min(nPasses,71);
529          }
530          if (nPasses<=5)
531            nPasses=0;
532          //info.setStartingWeight(1.0e-1);
533        }
534      }
535      if (doIdiot>0) {
536        // pick up number passes
537        nPasses=options.getExtraInfo(1);
538        if (nPasses>70) {
539          info.setStartingWeight(1.0e3);
540          info.setReduceIterations(6);
541          // also be more lenient on infeasibilities
542          info.setDropEnoughFeasibility(0.5*info.getDropEnoughFeasibility());
543          info.setDropEnoughWeighted(-2.0);
544        } else if (nPasses>=50) {
545          info.setStartingWeight(1.0e3);
546          //info.setReduceIterations(6);
547        } 
548        // For experimenting
549        if (nPasses<70&&(nPasses%10)>0&&(nPasses%10)<4) {
550          info.setStartingWeight(1.0e3);
551          info.setLightweight(nPasses%10); // special testing
552          //info.setReduceIterations(6);
553        }
554      }
555      if (nPasses) {
556        doCrash=0;
557#if 0
558        double * solution = model2->primalColumnSolution();
559        int iColumn;
560        double * saveLower = new double[numberColumns];
561        memcpy(saveLower,model2->columnLower(),numberColumns*sizeof(double));
562        double * saveUpper = new double[numberColumns];
563        memcpy(saveUpper,model2->columnUpper(),numberColumns*sizeof(double));
564        printf("doing tighten before idiot\n");
565        model2->tightenPrimalBounds();
566        // Move solution
567        double * columnLower = model2->columnLower();
568        double * columnUpper = model2->columnUpper();
569        for (iColumn=0;iColumn<numberColumns;iColumn++) {
570          if (columnLower[iColumn]>0.0)
571            solution[iColumn]=columnLower[iColumn];
572          else if (columnUpper[iColumn]<0.0)
573            solution[iColumn]=columnUpper[iColumn];
574          else
575            solution[iColumn]=0.0;
576        }
577        memcpy(columnLower,saveLower,numberColumns*sizeof(double));
578        memcpy(columnUpper,saveUpper,numberColumns*sizeof(double));
579        delete [] saveLower;
580        delete [] saveUpper;
581#else
582        // Allow for crossover
583        info.setStrategy(512|info.getStrategy());
584        // Allow for scaling
585        info.setStrategy(32|info.getStrategy());
586        info.crash(nPasses,model2->messageHandler(),model2->messagesPointer());
587#endif
588        time2 = CoinCpuTime();
589        timeIdiot = time2-timeX;
590        handler_->message(CLP_INTERVAL_TIMING,messages_)
591          <<"Idiot Crash"<<timeIdiot<<time2-time1
592          <<CoinMessageEol;
593        timeX=time2;
594      }
595    }
596    // ?
597    if (doCrash)
598      model2->crash(1000,1);
599    model2->primal(1);
600    time2 = CoinCpuTime();
601    timeCore = time2-timeX;
602    handler_->message(CLP_INTERVAL_TIMING,messages_)
603      <<"Primal"<<timeCore<<time2-time1
604      <<CoinMessageEol;
605    timeX=time2;
606  } else if (method==ClpSolve::usePrimalorSprint) {
607    // Sprint
608    /*
609      This driver implements what I called Sprint when I introduced the idea
610      many years ago.  Cplex calls it "sifting" which I think is just as silly.
611      When I thought of this trivial idea
612      it reminded me of an LP code of the 60's called sprint which after
613      every factorization took a subset of the matrix into memory (all
614      64K words!) and then iterated very fast on that subset.  On the
615      problems of those days it did not work very well, but it worked very
616      well on aircrew scheduling problems where there were very large numbers
617      of columns all with the same flavor.
618    */
619   
620    /* The idea works best if you can get feasible easily.  To make it
621       more general we can add in costed slacks */
622   
623    int originalNumberColumns = model2->numberColumns();
624    int numberRows = model2->numberRows();
625   
626    // We will need arrays to choose variables.  These are too big but ..
627    double * weight = new double [numberRows+originalNumberColumns];
628    int * sort = new int [numberRows+originalNumberColumns];
629    int numberSort=0;
630    // We are going to add slacks to get feasible.
631    // initial list will just be artificials
632    // first we will set all variables as close to zero as possible
633    int iColumn;
634    const double * columnLower = model2->columnLower();
635    const double * columnUpper = model2->columnUpper();
636    double * columnSolution = model2->primalColumnSolution();
637   
638    for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
639      double value =0.0;
640      if (columnLower[iColumn]>0.0)
641        value = columnLower[iColumn];
642      else if (columnUpper[iColumn]<0.0)
643        value = columnUpper[iColumn];
644      columnSolution[iColumn]=value;
645    }
646    // now see what that does to row solution
647    double * rowSolution = model2->primalRowSolution();
648    memset (rowSolution,0,numberRows*sizeof(double));
649    model2->times(1.0,columnSolution,rowSolution);
650   
651    int * addStarts = new int [numberRows+1];
652    int * addRow = new int[numberRows];
653    double * addElement = new double[numberRows];
654    const double * lower = model2->rowLower();
655    const double * upper = model2->rowUpper();
656    addStarts[0]=0;
657    int numberArtificials=0;
658    double * addCost = new double [numberRows];
659    const double penalty=1.0e8;
660    int iRow;
661    for (iRow=0;iRow<numberRows;iRow++) {
662      if (lower[iRow]>rowSolution[iRow]) {
663        addRow[numberArtificials]=iRow;
664        addElement[numberArtificials]=1.0;
665        addCost[numberArtificials]=penalty;
666        numberArtificials++;
667        addStarts[numberArtificials]=numberArtificials;
668      } else if (upper[iRow]<rowSolution[iRow]) {
669        addRow[numberArtificials]=iRow;
670        addElement[numberArtificials]=-1.0;
671        addCost[numberArtificials]=penalty;
672        numberArtificials++;
673        addStarts[numberArtificials]=numberArtificials;
674      }
675    }
676    model2->addColumns(numberArtificials,NULL,NULL,addCost,
677                       addStarts,addRow,addElement);
678    delete [] addStarts;
679    delete [] addRow;
680    delete [] addElement;
681    delete [] addCost;
682    // look at rhs to see if to perturb
683    double largest=0.0;
684    double smallest = 1.0e30;
685    for (iRow=0;iRow<numberRows;iRow++) {
686      double value;
687      value = fabs(model2->rowLower_[iRow]);
688      if (value&&value<1.0e30) {
689        largest = max(largest,value);
690        smallest=min(smallest,value);
691      }
692      value = fabs(model2->rowUpper_[iRow]);
693      if (value&&value<1.0e30) {
694        largest = max(largest,value);
695        smallest=min(smallest,value);
696      }
697    }
698    double * saveLower = NULL;
699    double * saveUpper = NULL;
700    if (largest<2.01*smallest) {
701      // perturb - so switch off standard
702      model2->setPerturbation(100);
703      saveLower = new double[numberRows];
704      memcpy(saveLower,model2->rowLower_,numberRows*sizeof(double));
705      saveUpper = new double[numberRows];
706      memcpy(saveUpper,model2->rowUpper_,numberRows*sizeof(double));
707      double * lower = model2->rowLower();
708      double * upper = model2->rowUpper();
709      for (iRow=0;iRow<numberRows;iRow++) {
710        double lowerValue=lower[iRow], upperValue=upper[iRow];
711        double value = CoinDrand48();
712        if (upperValue>lowerValue+primalTolerance_) {
713          if (lowerValue>-1.0e20&&lowerValue)
714            lowerValue -= value * 1.0e-4*fabs(lowerValue); 
715          if (upperValue<1.0e20&&upperValue)
716            upperValue += value * 1.0e-4*fabs(upperValue); 
717        } else if (upperValue>0.0) {
718          upperValue -= value * 1.0e-4*fabs(lowerValue); 
719          lowerValue -= value * 1.0e-4*fabs(lowerValue); 
720        } else if (upperValue<0.0) {
721          upperValue += value * 1.0e-4*fabs(lowerValue); 
722          lowerValue += value * 1.0e-4*fabs(lowerValue); 
723        } else {
724        }
725        lower[iRow]=lowerValue;
726        upper[iRow]=upperValue;
727      }
728    }
729    int i;
730    // Set up initial list
731    if (numberArtificials) {
732      numberSort=numberArtificials;
733      for (i=0;i<numberSort;i++)
734        sort[i] = i+originalNumberColumns;
735    } else {
736      numberSort = min(numberRows_,numberColumns_);
737      for (i=0;i<numberSort;i++)
738        sort[i] = i;
739    }
740   
741    // redo as will have changed
742    columnLower = model2->columnLower();
743    columnUpper = model2->columnUpper();
744    int numberColumns = model2->numberColumns();
745    double * fullSolution = model2->primalColumnSolution();
746   
747    // Just do this number of passes in Sprint
748    if (doSprint>0)
749      maxSprintPass=options.getExtraInfo(1);
750    int iPass;
751    double lastObjective=1.0e31;
752    // It will be safe to allow dense
753    model2->setInitialDenseFactorization(true);
754   
755    // Just take this number of columns in small problem
756    int smallNumberColumns = min(3*numberRows,numberColumns);
757    smallNumberColumns = max(smallNumberColumns,3000);
758    //int smallNumberColumns = min(12*numberRows/10,numberColumns);
759    //smallNumberColumns = max(smallNumberColumns,3000);
760    //smallNumberColumns = max(smallNumberColumns,numberRows+1000);
761    // We will be using all rows
762    int * whichRows = new int [numberRows];
763    for (iRow=0;iRow<numberRows;iRow++)
764      whichRows[iRow]=iRow;
765    double originalOffset;
766    model2->getDblParam(ClpObjOffset,originalOffset);
767    int totalIterations=0;
768    for (iPass=0;iPass<maxSprintPass;iPass++) {
769      //printf("Bug until submodel new version\n");
770      //CoinSort_2(sort,sort+numberSort,weight);
771      // Create small problem
772      ClpSimplex small(model2,numberRows,whichRows,numberSort,sort);
773      small.setPerturbation(model2->perturbation());
774      // now see what variables left out do to row solution
775      double * rowSolution = model2->primalRowSolution();
776      double * sumFixed = new double[numberRows];
777      memset (sumFixed,0,numberRows*sizeof(double));
778      int iRow,iColumn;
779      // zero out ones in small problem
780      for (iColumn=0;iColumn<numberSort;iColumn++) {
781        int kColumn = sort[iColumn];
782        fullSolution[kColumn]=0.0;
783      }
784      // Get objective offset
785      double offset=0.0;
786      const double * objective = model2->objective();
787      for (iColumn=0;iColumn<numberColumns;iColumn++) 
788        offset += fullSolution[iColumn]*objective[iColumn];
789      small.setDblParam(ClpObjOffset,originalOffset-offset);
790      model2->times(1.0,fullSolution,sumFixed);
791     
792      double * lower = small.rowLower();
793      double * upper = small.rowUpper();
794      for (iRow=0;iRow<numberRows;iRow++) {
795        if (lower[iRow]>-1.0e50) 
796          lower[iRow] -= sumFixed[iRow];
797        if (upper[iRow]<1.0e50)
798          upper[iRow] -= sumFixed[iRow];
799        rowSolution[iRow] -= sumFixed[iRow];
800      }
801      delete [] sumFixed;
802      // Solve
803      if (interrupt)
804        currentModel = &small;
805      small.primal();
806      totalIterations += small.numberIterations();
807      // move solution back
808      const double * solution = small.primalColumnSolution();
809      for (iColumn=0;iColumn<numberSort;iColumn++) {
810        int kColumn = sort[iColumn];
811        model2->setColumnStatus(kColumn,small.getColumnStatus(iColumn));
812        fullSolution[kColumn]=solution[iColumn];
813      }
814      for (iRow=0;iRow<numberRows;iRow++) 
815        model2->setRowStatus(iRow,small.getRowStatus(iRow));
816      memcpy(model2->primalRowSolution(),small.primalRowSolution(),
817             numberRows*sizeof(double));
818      // get reduced cost for large problem
819      memcpy(weight,model2->objective(),numberColumns*sizeof(double));
820      model2->transposeTimes(-1.0,small.dualRowSolution(),weight);
821      int numberNegative=0;
822      double sumNegative = 0.0;
823      // now massage weight so all basic in plus good djs
824      for (iColumn=0;iColumn<numberColumns;iColumn++) {
825        double dj = weight[iColumn]*optimizationDirection_;
826        double value = fullSolution[iColumn];
827        if (model2->getColumnStatus(iColumn)==ClpSimplex::basic) 
828          dj = -1.0e50;
829        else if (dj<0.0&&value<columnUpper[iColumn])
830          dj = dj;
831        else if (dj>0.0&&value>columnLower[iColumn])
832          dj = -dj;
833        else if (columnUpper[iColumn]>columnLower[iColumn])
834          dj = fabs(dj);
835        else
836          dj = 1.0e50;
837        weight[iColumn] = dj;
838        if (dj<-dualTolerance_&&dj>-1.0e50) {
839          numberNegative++;
840          sumNegative -= dj;
841        }
842        sort[iColumn] = iColumn;
843      }
844      handler_->message(CLP_SPRINT,messages_)
845        <<iPass+1<<small.numberIterations()<<small.objectiveValue()<<sumNegative
846        <<numberNegative
847        <<CoinMessageEol;
848      if ((small.objectiveValue()*optimizationDirection_>lastObjective-1.0e-7&&iPass>5)||
849          !small.numberIterations()||
850          iPass==maxSprintPass-1||small.status()==3) {
851       
852        break; // finished
853      } else {
854        lastObjective = small.objectiveValue()*optimizationDirection_;
855        // sort
856        CoinSort_2(weight,weight+numberColumns,sort);
857        numberSort = smallNumberColumns;
858      }
859    }
860    if (interrupt) 
861      currentModel = model2;
862    for (i=0;i<numberArtificials;i++)
863      sort[i] = i + originalNumberColumns;
864    model2->deleteColumns(numberArtificials,sort);
865    delete [] weight;
866    delete [] sort;
867    delete [] whichRows;
868    if (saveLower) {
869      // unperturb and clean
870      for (iRow=0;iRow<numberRows;iRow++) {
871        double diffLower = saveLower[iRow]-model2->rowLower_[iRow];
872        double diffUpper = saveUpper[iRow]-model2->rowUpper_[iRow];
873        model2->rowLower_[iRow]=saveLower[iRow];
874        model2->rowUpper_[iRow]=saveUpper[iRow];
875        if (diffLower) 
876          assert (!diffUpper||fabs(diffLower-diffUpper)<1.0e-5);
877        else
878          diffLower = diffUpper;
879        model2->rowActivity_[iRow] += diffLower;
880      }
881      delete [] saveLower;
882      delete [] saveUpper;
883    }
884    model2->primal(0);
885    model2->setPerturbation(savePerturbation);
886    time2 = CoinCpuTime();
887    timeCore = time2-timeX;
888    handler_->message(CLP_INTERVAL_TIMING,messages_)
889      <<"Sprint"<<timeCore<<time2-time1
890      <<CoinMessageEol;
891    timeX=time2;
892    model2->setNumberIterations(model2->numberIterations()+totalIterations);
893  } else if (method==ClpSolve::useBarrier) {
894    printf("***** experimental pretty crude barrier\n");
895    ClpInterior barrier(*model2);
896    if (interrupt) 
897      currentModel2 = &barrier;
898#ifdef REAL_BARRIER
899    // uncomment this if you have Anshul Gupta's wsmp package
900    //ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(max(100,model2->numberRows()/10));
901    ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
902    barrier.setCholesky(cholesky);
903    //ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(max(100,model2->numberRows()/10));
904    //barrier.setCholesky(cholesky);
905    // uncomment this if you have Sivan Toledo's Taucs package
906    //ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
907    //barrier.setCholesky(cholesky);
908#endif
909    //#define SAVEIT 2
910#ifndef SAVEIT
911    barrier.primalDual();
912#elif SAVEIT==1
913    barrier.primalDual();
914#else
915    model2->restoreModel("xx.save");
916    // move solutions
917    CoinMemcpyN(model2->primalRowSolution(),
918                barrier.numberRows(),barrier.primalRowSolution());
919    CoinMemcpyN(model2->dualRowSolution(),
920                barrier.numberRows(),barrier.dualRowSolution());
921    CoinMemcpyN(model2->primalColumnSolution(),
922                barrier.numberColumns(),barrier.primalColumnSolution());
923    CoinMemcpyN(model2->dualColumnSolution(),
924                barrier.numberColumns(),barrier.dualColumnSolution());
925#endif
926    time2 = CoinCpuTime();
927    timeCore = time2-timeX;
928    handler_->message(CLP_INTERVAL_TIMING,messages_)
929      <<"Barrier"<<timeCore<<time2-time1
930      <<CoinMessageEol;
931    timeX=time2;
932    if (barrier.maximumBarrierIterations()&&barrier.status()!=4) {
933      printf("***** crossover - needs more thought on difficult models\n");
934      // move solutions
935      CoinMemcpyN(barrier.primalRowSolution(),
936                  model2->numberRows(),model2->primalRowSolution());
937      CoinMemcpyN(barrier.dualRowSolution(),
938                  model2->numberRows(),model2->dualRowSolution());
939      CoinMemcpyN(barrier.primalColumnSolution(),
940                  model2->numberColumns(),model2->primalColumnSolution());
941      CoinMemcpyN(barrier.dualColumnSolution(),
942                  model2->numberColumns(),model2->dualColumnSolution());
943#if SAVEIT==1
944      model2->ClpSimplex::saveModel("xx.save");
945#endif
946      // make sure no status left
947      model2->createStatus();
948      // solve
949      model2->setPerturbation(100);
950#if 1
951      // throw some into basis
952      {
953        int numberColumns = model2->numberColumns();
954        int numberRows = model2->numberRows();
955        double * dsort = new double[numberColumns];
956        int * sort = new int[numberColumns];
957        int n=0;
958        const double * columnLower = model2->columnLower();
959        const double * columnUpper = model2->columnUpper();
960        const double * primalSolution = model2->primalColumnSolution();
961        double tolerance = 10.0*primalTolerance_;
962        int i;
963        for ( i=0;i<numberRows;i++) 
964          model2->setRowStatus(i,superBasic);
965        for ( i=0;i<numberColumns;i++) {
966          double distance = min(columnUpper[i]-primalSolution[i],
967                                primalSolution[i]-columnLower[i]);
968          if (distance>tolerance) {
969            dsort[n]=-distance;
970            sort[n++]=i;
971            model2->setStatus(i,superBasic);
972          } else if (distance>primalTolerance_) {
973            model2->setStatus(i,superBasic);
974          } else if (primalSolution[i]<=columnLower[i]+primalTolerance_) {
975            model2->setStatus(i,atLowerBound);
976          } else {
977            model2->setStatus(i,atUpperBound);
978          }
979        }
980        CoinSort_2(dsort,dsort+n,sort);
981        n = min(numberRows,n);
982        for ( i=0;i<n;i++) {
983          int iColumn = sort[i];
984          model2->setStatus(iColumn,basic);
985        }
986        delete [] sort;
987        delete [] dsort;
988      }
989      // just primal values pass
990      model2->primal(2);
991      // move solutions
992      CoinMemcpyN(model2->primalRowSolution(),
993                  model2->numberRows(),barrier.primalRowSolution());
994      CoinMemcpyN(barrier.dualRowSolution(),
995                  model2->numberRows(),model2->dualRowSolution());
996      CoinMemcpyN(model2->primalColumnSolution(),
997                  model2->numberColumns(),barrier.primalColumnSolution());
998      CoinMemcpyN(barrier.dualColumnSolution(),
999                  model2->numberColumns(),model2->dualColumnSolution());
1000      //model2->primal(1);
1001      // clean up reduced costs and flag variables
1002      {
1003        int numberColumns = model2->numberColumns();
1004        double * dj = model2->dualColumnSolution();
1005        double * cost = model2->objective();
1006        double * saveCost = new double[numberColumns];
1007        memcpy(saveCost,cost,numberColumns*sizeof(double));
1008        double * saveLower = new double[numberColumns];
1009        double * lower = model2->columnLower();
1010        memcpy(saveLower,lower,numberColumns*sizeof(double));
1011        double * saveUpper = new double[numberColumns];
1012        double * upper = model2->columnUpper();
1013        memcpy(saveUpper,upper,numberColumns*sizeof(double));
1014        int i;
1015        double tolerance = 10.0*dualTolerance_;
1016        for ( i=0;i<numberColumns;i++) {
1017          if (model2->getStatus(i)==basic) {
1018            dj[i]=0.0;
1019          } else if (model2->getStatus(i)==atLowerBound) {
1020            if (optimizationDirection_*dj[i]<tolerance) {
1021              if (optimizationDirection_*dj[i]<0.0) {
1022                //if (dj[i]<-1.0e-3)
1023                //printf("bad dj at lb %d %g\n",i,dj[i]);
1024                cost[i] -= dj[i];
1025                dj[i]=0.0;
1026              }
1027            } else {
1028              upper[i]=lower[i];
1029            }
1030          } else if (model2->getStatus(i)==atUpperBound) {
1031            if (optimizationDirection_*dj[i]>tolerance) {
1032              if (optimizationDirection_*dj[i]>0.0) {
1033                //if (dj[i]>1.0e-3)
1034                //printf("bad dj at ub %d %g\n",i,dj[i]);
1035                cost[i] -= dj[i];
1036                dj[i]=0.0;
1037              }
1038            } else {
1039              lower[i]=upper[i];
1040            }
1041          }
1042        }
1043        // just dual values pass
1044        //model2->setLogLevel(63);
1045        //model2->setFactorizationFrequency(1);
1046        model2->dual(2);
1047        memcpy(cost,saveCost,numberColumns*sizeof(double));
1048        delete [] saveCost;
1049        memcpy(lower,saveLower,numberColumns*sizeof(double));
1050        delete [] saveLower;
1051        memcpy(upper,saveUpper,numberColumns*sizeof(double));
1052        delete [] saveUpper;
1053      }
1054      // and finish
1055      // move solutions
1056      CoinMemcpyN(barrier.primalRowSolution(),
1057                  model2->numberRows(),model2->primalRowSolution());
1058      CoinMemcpyN(barrier.primalColumnSolution(),
1059                  model2->numberColumns(),model2->primalColumnSolution());
1060      model2->primal(1);
1061#else
1062      // just primal
1063      model2->primal(1);
1064#endif
1065    } else if (barrier.status()==4) {
1066      // memory problems
1067      model2->setPerturbation(savePerturbation);
1068      model2->createStatus();
1069      model2->dual();
1070    }
1071    model2->setPerturbation(savePerturbation);
1072    time2 = CoinCpuTime();
1073    timeCore = time2-timeX;
1074    handler_->message(CLP_INTERVAL_TIMING,messages_)
1075      <<"Crossover"<<timeCore<<time2-time1
1076      <<CoinMessageEol;
1077    timeX=time2;
1078  } else {
1079    assert (method!=ClpSolve::automatic); // later
1080    time2=0.0;
1081  }
1082  if (saveMatrix) {
1083    // delete and replace
1084    delete model2->clpMatrix();
1085    model2->replaceMatrix(saveMatrix);
1086  }
1087  numberIterations = model2->numberIterations();
1088  finalStatus=model2->status();
1089  if (presolve==ClpSolve::presolveOn) {
1090    int saveLevel = logLevel();
1091    setLogLevel(1);
1092    pinfo.postsolve(true);
1093    time2 = CoinCpuTime();
1094    timePresolve += time2-timeX;
1095    handler_->message(CLP_INTERVAL_TIMING,messages_)
1096      <<"Postsolve"<<time2-timeX<<time2-time1
1097      <<CoinMessageEol;
1098    timeX=time2;
1099    if (!presolveToFile)
1100      delete model2;
1101    if (interrupt)
1102      currentModel = this;
1103    checkSolution();
1104    setLogLevel(saveLevel);
1105    if (finalStatus!=3&&(finalStatus||status()==-1)) {
1106      int savePerturbation = perturbation();
1107      setPerturbation(100);
1108      primal(1);
1109      setPerturbation(savePerturbation);
1110      numberIterations += this->numberIterations();
1111      finalStatus=status();
1112      time2 = CoinCpuTime();
1113      handler_->message(CLP_INTERVAL_TIMING,messages_)
1114      <<"Cleanup"<<time2-timeX<<time2-time1
1115      <<CoinMessageEol;
1116      timeX=time2;
1117    }
1118  }
1119  setMaximumIterations(saveMaxIterations);
1120  std::string statusMessage[]={"Unknown","Optimal","PrimalInfeasible","DualInfeasible","Stopped",
1121                               "Errors"};
1122  assert (finalStatus>=-1&&finalStatus<=4);
1123  handler_->message(CLP_TIMING,messages_)
1124    <<statusMessage[finalStatus+1]<<objectiveValue()<<numberIterations<<time2-time1;
1125  handler_->printing(presolve==ClpSolve::presolveOn)
1126    <<timePresolve;
1127  handler_->printing(timeIdiot)
1128    <<timeIdiot;
1129  handler_->message()<<CoinMessageEol;
1130  if (interrupt) 
1131    signal(SIGINT,saveSignal);
1132  perturbation_=savePerturbation;
1133  scalingFlag_=saveScaling;
1134  return finalStatus;
1135}
1136// General solve
1137int 
1138ClpSimplex::initialSolve()
1139{
1140  // Default so use dual
1141  ClpSolve options;
1142  return initialSolve(options);
1143}
1144// General dual solve
1145int 
1146ClpSimplex::initialDualSolve()
1147{
1148  ClpSolve options;
1149  // Use dual
1150  options.setSolveType(ClpSolve::useDual);
1151  return initialSolve(options);
1152}
1153// General dual solve
1154int 
1155ClpSimplex::initialPrimalSolve()
1156{
1157  ClpSolve options;
1158  // Use primal
1159  options.setSolveType(ClpSolve::usePrimal);
1160  return initialSolve(options);
1161}
1162// Default constructor
1163ClpSolve::ClpSolve (  )
1164{
1165  method_ = useDual;
1166  presolveType_=presolveOn;
1167  numberPasses_=5;
1168  int i;
1169  for (i=0;i<4;i++)
1170    options_[i]=0;
1171  for (i=0;i<4;i++)
1172    extraInfo_[i]=-1;
1173}
1174
1175// Copy constructor.
1176ClpSolve::ClpSolve(const ClpSolve & rhs)
1177{
1178  method_ = rhs.method_;
1179  presolveType_=rhs.presolveType_;
1180  numberPasses_=rhs.numberPasses_;
1181  int i;
1182  for ( i=0;i<4;i++)
1183    options_[i]=rhs.options_[i];
1184  for ( i=0;i<4;i++)
1185    extraInfo_[i]=rhs.extraInfo_[i];
1186}
1187// Assignment operator. This copies the data
1188ClpSolve & 
1189ClpSolve::operator=(const ClpSolve & rhs)
1190{
1191  if (this != &rhs) {
1192    method_ = rhs.method_;
1193    presolveType_=rhs.presolveType_;
1194    numberPasses_=rhs.numberPasses_;
1195    int i;
1196    for (i=0;i<4;i++)
1197      options_[i]=rhs.options_[i];
1198    for (i=0;i<4;i++)
1199      extraInfo_[i]=rhs.extraInfo_[i];
1200  }
1201  return *this;
1202
1203}
1204// Destructor
1205ClpSolve::~ClpSolve (  )
1206{
1207}
1208/*   which translation is:
1209     which:
1210     0 - startup in Dual  (nothing if basis exists).:
1211             0 - no basis, 1 crash
1212     1 - startup in Primal (nothing if basis exists):
1213        0 - use initiative
1214        1 - use crash
1215        2 - use idiot and look at further info
1216        3 - use sprint and look at further info
1217        4 - use all slack
1218        5 - use initiative but no idiot
1219        6 - use initiative but no sprint
1220        7 - use initiative but no crash
1221        8 - do allslack or idiot
1222        9 - do allslack or sprint
1223     2 - interrupt handling - 0 yes, 1 no (for threadsafe)
1224     3 - whether to make +- 1matrix - 0 yes, 1 no
1225*/
1226void 
1227ClpSolve::setSpecialOption(int which,int value,int extraInfo)
1228{
1229  options_[which]=value;
1230  extraInfo_[which]=extraInfo;
1231}
1232int 
1233ClpSolve::getSpecialOption(int which) const
1234{
1235  return options_[which];
1236}
1237
1238// Solve types
1239void 
1240ClpSolve::setSolveType(SolveType method, int extraInfo)
1241{
1242  method_=method;
1243}
1244
1245ClpSolve::SolveType
1246ClpSolve::getSolveType()
1247{
1248  return method_;
1249}
1250
1251// Presolve types
1252void 
1253ClpSolve::setPresolveType(PresolveType amount, int extraInfo)
1254{
1255  presolveType_=amount;
1256  numberPasses_=extraInfo;
1257}
1258ClpSolve::PresolveType
1259ClpSolve::getPresolveType()
1260{
1261  return presolveType_;
1262}
1263// Extra info for idiot (or sprint)
1264int 
1265ClpSolve::getExtraInfo(int which) const
1266{
1267  return extraInfo_[which];
1268}
1269int 
1270ClpSolve::getPresolvePasses() const
1271{
1272  return numberPasses_;
1273}
Note: See TracBrowser for help on using the repository browser.