source: trunk/ClpSolve.cpp @ 437

Last change on this file since 437 was 437, checked in by forrest, 15 years ago

improvements for dual

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