source: trunk/ClpSolve.cpp @ 461

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

Trying to make faster

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