source: trunk/ClpSolve.cpp @ 460

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

a bit faster plus compiler bug in ClpsimplexPrimal?

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