source: trunk/ClpSolve.cpp @ 369

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

improving interior point code

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 37.8 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    //printf("********** Stopping as this is debug run\n");
913    //exit(99);
914#elif SAVEIT==1
915    barrier.primalDual();
916#else
917    model2->restoreModel("xx.save");
918    // move solutions
919    CoinMemcpyN(model2->primalRowSolution(),
920                barrier.numberRows(),barrier.primalRowSolution());
921    CoinMemcpyN(model2->dualRowSolution(),
922                barrier.numberRows(),barrier.dualRowSolution());
923    CoinMemcpyN(model2->primalColumnSolution(),
924                barrier.numberColumns(),barrier.primalColumnSolution());
925    CoinMemcpyN(model2->dualColumnSolution(),
926                barrier.numberColumns(),barrier.dualColumnSolution());
927#endif
928    time2 = CoinCpuTime();
929    timeCore = time2-timeX;
930    handler_->message(CLP_INTERVAL_TIMING,messages_)
931      <<"Barrier"<<timeCore<<time2-time1
932      <<CoinMessageEol;
933    timeX=time2;
934    if (barrier.maximumBarrierIterations()&&barrier.status()<4) {
935      printf("***** crossover - needs more thought on difficult models\n");
936      // move solutions
937      CoinMemcpyN(barrier.primalRowSolution(),
938                  model2->numberRows(),model2->primalRowSolution());
939      CoinMemcpyN(barrier.dualRowSolution(),
940                  model2->numberRows(),model2->dualRowSolution());
941      CoinMemcpyN(barrier.primalColumnSolution(),
942                  model2->numberColumns(),model2->primalColumnSolution());
943      CoinMemcpyN(barrier.dualColumnSolution(),
944                  model2->numberColumns(),model2->dualColumnSolution());
945#if SAVEIT==1
946      model2->ClpSimplex::saveModel("xx.save");
947#endif
948      // make sure no status left
949      model2->createStatus();
950      // solve
951      model2->setPerturbation(100);
952#if 1
953      // throw some into basis
954      {
955        int numberColumns = model2->numberColumns();
956        int numberRows = model2->numberRows();
957        double * dsort = new double[numberColumns];
958        int * sort = new int[numberColumns];
959        int n=0;
960        const double * columnLower = model2->columnLower();
961        const double * columnUpper = model2->columnUpper();
962        const double * primalSolution = model2->primalColumnSolution();
963        double tolerance = 10.0*primalTolerance_;
964        int i;
965        for ( i=0;i<numberRows;i++) 
966          model2->setRowStatus(i,superBasic);
967        for ( i=0;i<numberColumns;i++) {
968          double distance = min(columnUpper[i]-primalSolution[i],
969                                primalSolution[i]-columnLower[i]);
970          if (distance>tolerance) {
971            dsort[n]=-distance;
972            sort[n++]=i;
973            model2->setStatus(i,superBasic);
974          } else if (distance>primalTolerance_) {
975            model2->setStatus(i,superBasic);
976          } else if (primalSolution[i]<=columnLower[i]+primalTolerance_) {
977            model2->setStatus(i,atLowerBound);
978          } else {
979            model2->setStatus(i,atUpperBound);
980          }
981        }
982        CoinSort_2(dsort,dsort+n,sort);
983        n = min(numberRows,n);
984        for ( i=0;i<n;i++) {
985          int iColumn = sort[i];
986          model2->setStatus(iColumn,basic);
987        }
988        delete [] sort;
989        delete [] dsort;
990      }
991      // just primal values pass
992      model2->primal(2);
993      // move solutions
994      CoinMemcpyN(model2->primalRowSolution(),
995                  model2->numberRows(),barrier.primalRowSolution());
996      CoinMemcpyN(barrier.dualRowSolution(),
997                  model2->numberRows(),model2->dualRowSolution());
998      CoinMemcpyN(model2->primalColumnSolution(),
999                  model2->numberColumns(),barrier.primalColumnSolution());
1000      CoinMemcpyN(barrier.dualColumnSolution(),
1001                  model2->numberColumns(),model2->dualColumnSolution());
1002      //model2->primal(1);
1003      // clean up reduced costs and flag variables
1004      {
1005        int numberColumns = model2->numberColumns();
1006        double * dj = model2->dualColumnSolution();
1007        double * cost = model2->objective();
1008        double * saveCost = new double[numberColumns];
1009        memcpy(saveCost,cost,numberColumns*sizeof(double));
1010        double * saveLower = new double[numberColumns];
1011        double * lower = model2->columnLower();
1012        memcpy(saveLower,lower,numberColumns*sizeof(double));
1013        double * saveUpper = new double[numberColumns];
1014        double * upper = model2->columnUpper();
1015        memcpy(saveUpper,upper,numberColumns*sizeof(double));
1016        int i;
1017        double tolerance = 10.0*dualTolerance_;
1018        for ( i=0;i<numberColumns;i++) {
1019          if (model2->getStatus(i)==basic) {
1020            dj[i]=0.0;
1021          } else if (model2->getStatus(i)==atLowerBound) {
1022            if (optimizationDirection_*dj[i]<tolerance) {
1023              if (optimizationDirection_*dj[i]<0.0) {
1024                //if (dj[i]<-1.0e-3)
1025                //printf("bad dj at lb %d %g\n",i,dj[i]);
1026                cost[i] -= dj[i];
1027                dj[i]=0.0;
1028              }
1029            } else {
1030              upper[i]=lower[i];
1031            }
1032          } else if (model2->getStatus(i)==atUpperBound) {
1033            if (optimizationDirection_*dj[i]>tolerance) {
1034              if (optimizationDirection_*dj[i]>0.0) {
1035                //if (dj[i]>1.0e-3)
1036                //printf("bad dj at ub %d %g\n",i,dj[i]);
1037                cost[i] -= dj[i];
1038                dj[i]=0.0;
1039              }
1040            } else {
1041              lower[i]=upper[i];
1042            }
1043          }
1044        }
1045        // just dual values pass
1046        //model2->setLogLevel(63);
1047        //model2->setFactorizationFrequency(1);
1048        model2->dual(2);
1049        memcpy(cost,saveCost,numberColumns*sizeof(double));
1050        delete [] saveCost;
1051        memcpy(lower,saveLower,numberColumns*sizeof(double));
1052        delete [] saveLower;
1053        memcpy(upper,saveUpper,numberColumns*sizeof(double));
1054        delete [] saveUpper;
1055      }
1056      // and finish
1057      // move solutions
1058      CoinMemcpyN(barrier.primalRowSolution(),
1059                  model2->numberRows(),model2->primalRowSolution());
1060      CoinMemcpyN(barrier.primalColumnSolution(),
1061                  model2->numberColumns(),model2->primalColumnSolution());
1062      model2->primal(1);
1063#else
1064      // just primal
1065      model2->primal(1);
1066#endif
1067    } else if (barrier.status()==4) {
1068      // memory problems
1069      model2->setPerturbation(savePerturbation);
1070      model2->createStatus();
1071      model2->dual();
1072    }
1073    model2->setPerturbation(savePerturbation);
1074    time2 = CoinCpuTime();
1075    timeCore = time2-timeX;
1076    handler_->message(CLP_INTERVAL_TIMING,messages_)
1077      <<"Crossover"<<timeCore<<time2-time1
1078      <<CoinMessageEol;
1079    timeX=time2;
1080  } else {
1081    assert (method!=ClpSolve::automatic); // later
1082    time2=0.0;
1083  }
1084  if (saveMatrix) {
1085    // delete and replace
1086    delete model2->clpMatrix();
1087    model2->replaceMatrix(saveMatrix);
1088  }
1089  numberIterations = model2->numberIterations();
1090  finalStatus=model2->status();
1091  if (presolve==ClpSolve::presolveOn) {
1092    int saveLevel = logLevel();
1093    setLogLevel(1);
1094    pinfo.postsolve(true);
1095    time2 = CoinCpuTime();
1096    timePresolve += time2-timeX;
1097    handler_->message(CLP_INTERVAL_TIMING,messages_)
1098      <<"Postsolve"<<time2-timeX<<time2-time1
1099      <<CoinMessageEol;
1100    timeX=time2;
1101    if (!presolveToFile)
1102      delete model2;
1103    if (interrupt)
1104      currentModel = this;
1105    checkSolution();
1106    setLogLevel(saveLevel);
1107    if (finalStatus!=3&&(finalStatus||status()==-1)) {
1108      int savePerturbation = perturbation();
1109      setPerturbation(100);
1110      primal(1);
1111      setPerturbation(savePerturbation);
1112      numberIterations += this->numberIterations();
1113      finalStatus=status();
1114      time2 = CoinCpuTime();
1115      handler_->message(CLP_INTERVAL_TIMING,messages_)
1116      <<"Cleanup"<<time2-timeX<<time2-time1
1117      <<CoinMessageEol;
1118      timeX=time2;
1119    }
1120  }
1121  setMaximumIterations(saveMaxIterations);
1122  std::string statusMessage[]={"Unknown","Optimal","PrimalInfeasible","DualInfeasible","Stopped",
1123                               "Errors"};
1124  assert (finalStatus>=-1&&finalStatus<=4);
1125  handler_->message(CLP_TIMING,messages_)
1126    <<statusMessage[finalStatus+1]<<objectiveValue()<<numberIterations<<time2-time1;
1127  handler_->printing(presolve==ClpSolve::presolveOn)
1128    <<timePresolve;
1129  handler_->printing(timeIdiot)
1130    <<timeIdiot;
1131  handler_->message()<<CoinMessageEol;
1132  if (interrupt) 
1133    signal(SIGINT,saveSignal);
1134  perturbation_=savePerturbation;
1135  scalingFlag_=saveScaling;
1136  return finalStatus;
1137}
1138// General solve
1139int 
1140ClpSimplex::initialSolve()
1141{
1142  // Default so use dual
1143  ClpSolve options;
1144  return initialSolve(options);
1145}
1146// General dual solve
1147int 
1148ClpSimplex::initialDualSolve()
1149{
1150  ClpSolve options;
1151  // Use dual
1152  options.setSolveType(ClpSolve::useDual);
1153  return initialSolve(options);
1154}
1155// General dual solve
1156int 
1157ClpSimplex::initialPrimalSolve()
1158{
1159  ClpSolve options;
1160  // Use primal
1161  options.setSolveType(ClpSolve::usePrimal);
1162  return initialSolve(options);
1163}
1164// Default constructor
1165ClpSolve::ClpSolve (  )
1166{
1167  method_ = useDual;
1168  presolveType_=presolveOn;
1169  numberPasses_=5;
1170  int i;
1171  for (i=0;i<4;i++)
1172    options_[i]=0;
1173  for (i=0;i<4;i++)
1174    extraInfo_[i]=-1;
1175}
1176
1177// Copy constructor.
1178ClpSolve::ClpSolve(const ClpSolve & rhs)
1179{
1180  method_ = rhs.method_;
1181  presolveType_=rhs.presolveType_;
1182  numberPasses_=rhs.numberPasses_;
1183  int i;
1184  for ( i=0;i<4;i++)
1185    options_[i]=rhs.options_[i];
1186  for ( i=0;i<4;i++)
1187    extraInfo_[i]=rhs.extraInfo_[i];
1188}
1189// Assignment operator. This copies the data
1190ClpSolve & 
1191ClpSolve::operator=(const ClpSolve & rhs)
1192{
1193  if (this != &rhs) {
1194    method_ = rhs.method_;
1195    presolveType_=rhs.presolveType_;
1196    numberPasses_=rhs.numberPasses_;
1197    int i;
1198    for (i=0;i<4;i++)
1199      options_[i]=rhs.options_[i];
1200    for (i=0;i<4;i++)
1201      extraInfo_[i]=rhs.extraInfo_[i];
1202  }
1203  return *this;
1204
1205}
1206// Destructor
1207ClpSolve::~ClpSolve (  )
1208{
1209}
1210/*   which translation is:
1211     which:
1212     0 - startup in Dual  (nothing if basis exists).:
1213             0 - no basis, 1 crash
1214     1 - startup in Primal (nothing if basis exists):
1215        0 - use initiative
1216        1 - use crash
1217        2 - use idiot and look at further info
1218        3 - use sprint and look at further info
1219        4 - use all slack
1220        5 - use initiative but no idiot
1221        6 - use initiative but no sprint
1222        7 - use initiative but no crash
1223        8 - do allslack or idiot
1224        9 - do allslack or sprint
1225     2 - interrupt handling - 0 yes, 1 no (for threadsafe)
1226     3 - whether to make +- 1matrix - 0 yes, 1 no
1227*/
1228void 
1229ClpSolve::setSpecialOption(int which,int value,int extraInfo)
1230{
1231  options_[which]=value;
1232  extraInfo_[which]=extraInfo;
1233}
1234int 
1235ClpSolve::getSpecialOption(int which) const
1236{
1237  return options_[which];
1238}
1239
1240// Solve types
1241void 
1242ClpSolve::setSolveType(SolveType method, int extraInfo)
1243{
1244  method_=method;
1245}
1246
1247ClpSolve::SolveType
1248ClpSolve::getSolveType()
1249{
1250  return method_;
1251}
1252
1253// Presolve types
1254void 
1255ClpSolve::setPresolveType(PresolveType amount, int extraInfo)
1256{
1257  presolveType_=amount;
1258  numberPasses_=extraInfo;
1259}
1260ClpSolve::PresolveType
1261ClpSolve::getPresolveType()
1262{
1263  return presolveType_;
1264}
1265// Extra info for idiot (or sprint)
1266int 
1267ClpSolve::getExtraInfo(int which) const
1268{
1269  return extraInfo_[which];
1270}
1271int 
1272ClpSolve::getPresolvePasses() const
1273{
1274  return numberPasses_;
1275}
Note: See TracBrowser for help on using the repository browser.