source: branches/pre/include/ClpInterior.hpp @ 212

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

lots of stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 36.7 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4/*
5   Authors
6   
7   John Tomlin (with some help from John Forrest)
8
9 */
10#ifndef ClpInterior_H
11#define ClpInterior_H
12
13#include <iostream>
14#include <cfloat>
15#include "ClpModel.hpp"
16#include "ClpMatrixBase.hpp"
17#include "ClpSolve.hpp"
18class ClpDualRowPivot;
19class ClpPrimalColumnPivot;
20class ClpFactorization;
21class CoinIndexedVector;
22class ClpNonLinearCost;
23class ClpInteriorProgress;
24
25/** This solves LPs using the simplex method
26
27    It inherits from ClpModel and all its arrays are created at
28    algorithm time. Originally I tried to work with model arrays
29    but for simplicity of coding I changed to single arrays with
30    structural variables then row variables.  Some coding is still
31    based on old style and needs cleaning up.
32
33    For a description of algorithms:
34
35    for dual see ClpInteriorDual.hpp and at top of ClpInteriorDual.cpp
36    for primal see ClpInteriorPrimal.hpp and at top of ClpInteriorPrimal.cpp
37
38    There is an algorithm data member.  + for primal variations
39    and - for dual variations
40
41    This file also includes (at end) a very simple class ClpInteriorProgress
42    which is where anti-looping stuff should migrate to
43
44*/
45
46class ClpInterior : public ClpModel {
47   friend void ClpInteriorUnitTest(const std::string & mpsDir,
48                                  const std::string & netlibDir);
49
50public:
51  /** enums for status of various sorts.
52      First 4 match CoinWarmStartBasis,
53      isFixed means fixed at lower bound and out of basis
54  */
55  enum Status {
56    isFree = 0x00,
57    basic = 0x01,
58    atUpperBound = 0x02,
59    atLowerBound = 0x03,
60    superBasic = 0x04,
61    isFixed = 0x05
62  };
63
64  enum FakeBound {
65    noFake = 0x00,
66    bothFake = 0x01,
67    upperFake = 0x02,
68    lowerFake = 0x03
69  };
70
71  /**@name Constructors and destructor and copy */
72  //@{
73  /// Default constructor
74    ClpInterior (  );
75
76  /// Copy constructor.
77  ClpInterior(const ClpInterior &);
78  /// Copy constructor from model.
79  ClpInterior(const ClpModel &);
80  /** Subproblem constructor.  A subset of whole model is created from the
81      row and column lists given.  The new order is given by list order and
82      duplicates are allowed.  Name and integer information can be dropped
83  */
84  ClpInterior (const ClpModel * wholeModel,
85              int numberRows, const int * whichRows,
86              int numberColumns, const int * whichColumns,
87              bool dropNames=true, bool dropIntegers=true);
88  /// Assignment operator. This copies the data
89    ClpInterior & operator=(const ClpInterior & rhs);
90  /// Destructor
91   ~ClpInterior (  );
92  // Ones below are just ClpModel with some changes
93  /** Loads a problem (the constraints on the
94        rows are given by lower and upper bounds). If a pointer is 0 then the
95        following values are the default:
96        <ul>
97          <li> <code>colub</code>: all columns have upper bound infinity
98          <li> <code>collb</code>: all columns have lower bound 0
99          <li> <code>rowub</code>: all rows have upper bound infinity
100          <li> <code>rowlb</code>: all rows have lower bound -infinity
101          <li> <code>obj</code>: all variables have 0 objective coefficient
102        </ul>
103    */
104  void loadProblem (  const ClpMatrixBase& matrix,
105                     const double* collb, const double* colub,   
106                     const double* obj,
107                     const double* rowlb, const double* rowub,
108                      const double * rowObjective=NULL);
109  void loadProblem (  const CoinPackedMatrix& matrix,
110                     const double* collb, const double* colub,   
111                     const double* obj,
112                     const double* rowlb, const double* rowub,
113                      const double * rowObjective=NULL);
114
115  /** Just like the other loadProblem() method except that the matrix is
116        given in a standard column major ordered format (without gaps). */
117  void loadProblem (  const int numcols, const int numrows,
118                     const CoinBigIndex* start, const int* index,
119                     const double* value,
120                     const double* collb, const double* colub,   
121                     const double* obj,
122                      const double* rowlb, const double* rowub,
123                      const double * rowObjective=NULL);
124  /// This one is for after presolve to save memory
125  void loadProblem (  const int numcols, const int numrows,
126                     const CoinBigIndex* start, const int* index,
127                      const double* value,const int * length,
128                     const double* collb, const double* colub,   
129                     const double* obj,
130                      const double* rowlb, const double* rowub,
131                      const double * rowObjective=NULL);
132  /// Read an mps file from the given filename
133  int readMps(const char *filename,
134              bool keepNames=false,
135              bool ignoreErrors = false);
136  /** Borrow model.  This is so we dont have to copy large amounts
137      of data around.  It assumes a derived class wants to overwrite
138      an empty model with a real one - while it does an algorithm.
139      This is same as ClpModel one, but sets scaling on etc. */
140  void borrowModel(ClpModel & otherModel);
141  //@}
142
143  /**@name Functions most useful to user */
144  //@{
145  /** General solve algorithm which can do presolve.
146      See  ClpSolve.hpp for options
147   */
148  int initialSolve(ClpSolve & options);
149  /** Dual algorithm - see ClpInteriorDual.hpp for method */
150  int dual(int ifValuesPass=0);
151  /** Primal algorithm - see ClpInteriorPrimal.hpp for method */
152  int primal(int ifValuesPass=0);
153  /** Solves quadratic problem using SLP - may be used as crash
154      for other algorithms when number of iterations small.
155      Also exits if all problematical variables are changing
156      less than deltaTolerance
157  */
158  int quadraticSLP(int numberPasses,double deltaTolerance);
159  /// Solves quadratic using Dantzig's algorithm - primal
160  int quadraticPrimal(int phase=2);
161  /// Passes in factorization
162  void setFactorization( ClpFactorization & factorization);
163  /// Sets or unsets scaling, 0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
164  void scaling(int mode=1);
165  /// Gets scalingFlag
166  inline int scalingFlag() const {return scalingFlag_;};
167  /** Tightens primal bounds to make dual faster.  Unless
168      fixed, bounds are slightly looser than they could be.
169      This is to make dual go faster and is probably not needed
170      with a presolve.  Returns non-zero if problem infeasible.
171
172      Fudge for branch and bound - put bounds on columns of factor *
173      largest value (at continuous) - should improve stability
174      in branch and bound on infeasible branches (0.0 is off)
175  */
176  int tightenPrimalBounds(double factor=0.0);
177  /** Crash - at present just aimed at dual, returns
178      -2 if dual preferred and crash basis created
179      -1 if dual preferred and all slack basis preferred
180       0 if basis going in was not all slack
181       1 if primal preferred and all slack basis preferred
182       2 if primal preferred and crash basis created.
183       
184       if gap between bounds <="gap" variables can be flipped
185
186       If "pivot" is
187       0 No pivoting (so will just be choice of algorithm)
188       1 Simple pivoting e.g. gub
189       2 Mini iterations
190  */
191  int crash(double gap,int pivot);
192  /// Sets row pivot choice algorithm in dual
193  void setDualRowPivotAlgorithm(ClpDualRowPivot & choice);
194  /// Sets column pivot choice algorithm in primal
195  void setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice);
196  /** For strong branching.  On input lower and upper are new bounds
197      while on output they are change in objective function values
198      (>1.0e50 infeasible).
199      Return code is 0 if nothing interesting, -1 if infeasible both
200      ways and +1 if infeasible one way (check values to see which one(s))
201      Solutions are filled in as well - even down, odd up - also
202      status and number of iterations
203  */
204  int strongBranching(int numberVariables,const int * variables,
205                      double * newLower, double * newUpper,
206                      double ** outputSolution,
207                      int * outputStatus, int * outputIterations,
208                      bool stopOnFirstInfeasible=true,
209                      bool alwaysFinish=false);
210  //@}
211
212  /**@name Needed for functionality of OsiSimplexInterface */
213  //@{
214  /** Pivot in a variable and out a variable.  Returns 0 if okay,
215      1 if inaccuracy forced re-factorization, -1 if would be singular.
216      Also updates primal/dual infeasibilities.
217      Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
218  */
219  int pivot();
220
221  /** Pivot in a variable and choose an outgoing one.  Assumes primal
222      feasible - will not go through a bound.  Returns step length in theta
223      Returns ray in ray_ (or NULL if no pivot)
224      Return codes as before but -1 means no acceptable pivot
225  */
226  int primalPivotResult();
227 
228  /** Pivot out a variable and choose an incoing one.  Assumes dual
229      feasible - will not go through a reduced cost. 
230      Returns step length in theta
231      Returns ray in ray_ (or NULL if no pivot)
232      Return codes as before but -1 means no acceptable pivot
233  */
234  int dualPivotResult();
235
236  /** Common bits of coding for dual and primal.  Return s0 if okay,
237      1 if bad matrix, 2 if very bad factorization
238  */
239  int startup(int ifValuesPass);
240  void finish();
241 
242  /** Factorizes and returns true if optimal.  Used by user */
243  bool statusOfProblem();
244  //@}
245
246  /**@name most useful gets and sets */
247  //@{
248  /// If problem is primal feasible
249  inline bool primalFeasible() const
250         { return (numberPrimalInfeasibilities_==0);};
251  /// If problem is dual feasible
252  inline bool dualFeasible() const
253         { return (numberDualInfeasibilities_==0);};
254  /// factorization
255  inline ClpFactorization * factorization() const 
256          { return factorization_;};
257  /// Sparsity on or off
258  bool sparseFactorization() const;
259  void setSparseFactorization(bool value);
260  /// Factorization frequency
261  int factorizationFrequency() const;
262  void setFactorizationFrequency(int value);
263  /// Dual bound
264  inline double dualBound() const
265          { return dualBound_;};
266  void setDualBound(double value);
267  /// Infeasibility cost
268  inline double infeasibilityCost() const
269          { return infeasibilityCost_;};
270  void setInfeasibilityCost(double value);
271  /** Amount of print out:
272      0 - none
273      1 - just final
274      2 - just factorizations
275      3 - as 2 plus a bit more
276      4 - verbose
277      above that 8,16,32 etc just for selective debug
278  */
279  /** Perturbation:
280      -50 to +50 - perturb by this power of ten (-6 sounds good)
281      100 - auto perturb if takes too long (1.0e-6 largest nonzero)
282      101 - we are perturbed
283      102 - don't try perturbing again
284      default is 100
285  */
286  inline int perturbation() const
287    { return perturbation_;};
288  void setPerturbation(int value);
289  /// Current (or last) algorithm
290  inline int algorithm() const 
291  {return algorithm_; } ;
292  /// Set algorithm
293  inline void setAlgorithm(int value)
294  {algorithm_=value; } ;
295  /// Sum of dual infeasibilities
296  inline double sumDualInfeasibilities() const 
297          { return sumDualInfeasibilities_;} ;
298  /// Number of dual infeasibilities
299  inline int numberDualInfeasibilities() const 
300          { return numberDualInfeasibilities_;} ;
301  /// Sum of primal infeasibilities
302  inline double sumPrimalInfeasibilities() const 
303          { return sumPrimalInfeasibilities_;} ;
304  /// Number of primal infeasibilities
305  inline int numberPrimalInfeasibilities() const 
306          { return numberPrimalInfeasibilities_;} ;
307  /** Save model to file, returns 0 if success.  This is designed for
308      use outside algorithms so does not save iterating arrays etc.
309  It does not save any messaging information.
310  Does not save scaling values.
311  It does not know about all types of virtual functions.
312  */
313  int saveModel(const char * fileName);
314  /** Restore model from file, returns 0 if success,
315      deletes current model */
316  int restoreModel(const char * fileName);
317 
318  /** Just check solution (for external use) - sets sum of
319      infeasibilities etc */
320  void checkSolution();
321  /// Useful row length arrays (0,1,2,3,4,5)
322  inline CoinIndexedVector * rowArray(int index) const
323  { return rowArray_[index];};
324  /// Useful column length arrays (0,1,2,3,4,5)
325  inline CoinIndexedVector * columnArray(int index) const
326  { return columnArray_[index];};
327  //@}
328
329  /******************** End of most useful part **************/
330  /**@name Functions less likely to be useful to casual user */
331  //@{
332  /** Given an existing factorization computes and checks
333      primal and dual solutions.  Uses input arrays for variables at
334      bounds.  Returns feasibility states */
335  int getSolution (  const double * rowActivities,
336                     const double * columnActivities);
337  /** Given an existing factorization computes and checks
338      primal and dual solutions.  Uses current problem arrays for
339      bounds.  Returns feasibility states */
340  int getSolution ();
341  /** Constructs a non linear cost from list of non-linearities (columns only)
342      First lower of each column is taken as real lower
343      Last lower is taken as real upper and cost ignored
344
345      Returns nonzero if bad data e.g. lowers not monotonic
346  */
347  int createPiecewiseLinearCosts(const int * starts,
348                   const double * lower, const double * gradient);
349  /** Return model - updates any scalars */
350  void returnModel(ClpInterior & otherModel);
351  /** Factorizes using current basis. 
352      solveType - 1 iterating, 0 initial, -1 external
353      If 10 added then in primal values pass
354      Return codes are as from ClpFactorization unless initial factorization
355      when total number of singularities is returned
356  */
357    /// Save data
358    ClpDataSave saveData() ;
359    /// Restore data
360    void restoreData(ClpDataSave saved);
361  int internalFactorize(int solveType);
362  /// Clean up status
363  void cleanStatus();
364  /// Factorizes using current basis. For external use
365  int factorize();
366  /** Computes duals from scratch. If givenDjs then
367      allows for nonzero basic djs */
368  void computeDuals(double * givenDjs);
369  /// Computes primals from scratch
370  void computePrimals (  const double * rowActivities,
371                     const double * columnActivities);
372  /**
373     Unpacks one column of the matrix into indexed array
374     Uses sequenceIn_
375     Also applies scaling if needed
376  */
377  void unpack(CoinIndexedVector * rowArray) const ;
378  /**
379     Unpacks one column of the matrix into indexed array
380     Slack if sequence>= numberColumns
381     Also applies scaling if needed
382  */
383  void unpack(CoinIndexedVector * rowArray,int sequence) const;
384  /**
385     Unpacks one column of the matrix into indexed array
386     ** as packed vector
387     Uses sequenceIn_
388     Also applies scaling if needed
389  */
390  void unpackPacked(CoinIndexedVector * rowArray) ;
391  /**
392     Unpacks one column of the matrix into indexed array
393     ** as packed vector
394     Slack if sequence>= numberColumns
395     Also applies scaling if needed
396  */
397  void unpackPacked(CoinIndexedVector * rowArray,int sequence);
398 
399  /**
400      This does basis housekeeping and does values for in/out variables.
401      Can also decide to re-factorize
402  */
403  int housekeeping(double objectiveChange);
404  /** This sets largest infeasibility and most infeasible and sum
405      and number of infeasibilities (Primal) */
406  void checkPrimalSolution(const double * rowActivities=NULL,
407                           const double * columnActivies=NULL);
408  /** This sets largest infeasibility and most infeasible and sum
409      and number of infeasibilities (Dual) */
410  void checkDualSolution();
411  /** For advanced use.  When doing iterative solves things can get
412      nasty so on values pass if incoming solution has largest
413      infeasibility < incomingInfeasibility throw out variables
414      from basis until largest infeasibility < allowedInfeasibility
415      or incoming largest infeasibility.
416      If allowedInfeasibility>= incomingInfeasibility this is
417      always possible altough you may end up with an all slack basis.
418
419      Defaults are 1.0,10.0
420  */
421  void setValuesPassAction(float incomingInfeasibility,
422                           float allowedInfeasibility);
423  //@}
424  /**@name Matrix times vector methods
425     They can be faster if scalar is +- 1
426     These are covers so user need not worry about scaling
427     Also for simplex I am not using basic/non-basic split */
428  //@{
429    /** Return <code>y + A * x * scalar</code> in <code>y</code>.
430        @pre <code>x</code> must be of size <code>numColumns()</code>
431        @pre <code>y</code> must be of size <code>numRows()</code> */
432   void times(double scalar,
433                       const double * x, double * y) const;
434    /** Return <code>y + x * scalar * A</code> in <code>y</code>.
435        @pre <code>x</code> must be of size <code>numRows()</code>
436        @pre <code>y</code> must be of size <code>numColumns()</code> */
437    void transposeTimes(double scalar,
438                                const double * x, double * y) const ;
439  //@}
440
441  /**@name most useful gets and sets */
442  //@{
443  /// Worst column primal infeasibility
444  inline double columnPrimalInfeasibility() const 
445          { return columnPrimalInfeasibility_;} ;
446  /// Sequence of worst (-1 if feasible)
447  inline int columnPrimalSequence() const 
448          { return columnPrimalSequence_;} ;
449  /// Worst row primal infeasibility
450  inline double rowPrimalInfeasibility() const 
451          { return rowPrimalInfeasibility_;} ;
452  /// Sequence of worst (-1 if feasible)
453  inline int rowPrimalSequence() const 
454          { return rowPrimalSequence_;} ;
455  /** Worst column dual infeasibility (note - these may not be as meaningful
456      if the problem is primal infeasible */
457  inline double columnDualInfeasibility() const 
458          { return columnDualInfeasibility_;} ;
459  /// Sequence of worst (-1 if feasible)
460  inline int columnDualSequence() const 
461          { return columnDualSequence_;} ;
462  /// Worst row dual infeasibility
463  inline double rowDualInfeasibility() const 
464          { return rowDualInfeasibility_;} ;
465  /// Sequence of worst (-1 if feasible)
466  inline int rowDualSequence() const 
467          { return rowDualSequence_;} ;
468  /// Primal tolerance needed to make dual feasible (<largeTolerance)
469  inline double primalToleranceToGetOptimal() const 
470          { return primalToleranceToGetOptimal_;} ;
471  /// Remaining largest dual infeasibility
472  inline double remainingDualInfeasibility() const 
473          { return remainingDualInfeasibility_;} ;
474  /// Large bound value (for complementarity etc)
475  inline double largeValue() const 
476          { return largeValue_;} ;
477  void setLargeValue( double value) ;
478  /// Largest error on Ax-b
479  inline double largestPrimalError() const
480          { return largestPrimalError_;} ;
481  /// Largest error on basic duals
482  inline double largestDualError() const
483          { return largestDualError_;} ;
484  /// Largest difference between input primal solution and computed
485  inline double largestSolutionError() const
486          { return largestSolutionError_;} ;
487  /// Basic variables pivoting on which rows
488  inline int * pivotVariable() const
489          { return pivotVariable_;};
490  /// Scaling of objective 12345.0 (auto), 0.0 (off), other user
491  inline double objectiveScale() const 
492          { return objectiveScale_;} ;
493  inline void setObjectiveScale(double value)
494          { objectiveScale_ = value;} ;
495  /// Current dual tolerance
496  inline double currentDualTolerance() const 
497          { return dualTolerance_;} ;
498  inline void setCurrentDualTolerance(double value)
499          { dualTolerance_ = value;} ;
500  /// Current primal tolerance
501  inline double currentPrimalTolerance() const 
502          { return primalTolerance_;} ;
503  inline void setCurrentPrimalTolerance(double value)
504          { primalTolerance_ = value;} ;
505  /// How many iterative refinements to do
506  inline int numberRefinements() const 
507          { return numberRefinements_;} ;
508  void setNumberRefinements( int value) ;
509  /// Alpha (pivot element) for use by classes e.g. steepestedge
510  inline double alpha() const { return alpha_;};
511  /// Reduced cost of last incoming for use by classes e.g. steepestedge
512  inline double dualIn() const { return dualIn_;};
513  /// Pivot Row for use by classes e.g. steepestedge
514  inline int pivotRow() const{ return pivotRow_;};
515  /// value of incoming variable (in Dual)
516  double valueIncomingDual() const;
517  //@}
518
519  protected:
520  /**@name protected methods */
521  //@{
522  /** May change basis and then returns number changed.
523      Computation of solutions may be overriden by given pi and solution
524  */
525  int gutsOfSolution ( double * givenDuals,
526                       const double * givenPrimals,
527                       bool valuesPass=false);
528  /// Does most of deletion (0 = all, 1 = most, 2 most + factorization)
529  void gutsOfDelete(int type);
530  /// Does most of copying
531  void gutsOfCopy(const ClpInterior & rhs);
532  /** puts in format I like (rowLower,rowUpper) also see StandardMatrix
533      1 bit does rows, 2 bit does column bounds, 4 bit does objective(s).
534      8 bit does solution scaling in
535      16 bit does rowArray and columnArray indexed vectors
536      and makes row copy if wanted, also sets columnStart_ etc
537      Also creates scaling arrays if needed.  It does scaling if needed.
538      16 also moves solutions etc in to work arrays
539      On 16 returns false if problem "bad" i.e. matrix or bounds bad
540  */
541  bool createRim(int what,bool makeRowCopy=false);
542  /** releases above arrays and does solution scaling out.  May also
543      get rid of factorization data -
544      0 get rid of nothing, 1 get rid of arrays, 2 also factorization
545  */
546  void deleteRim(int getRidOfFactorizationData=2);
547  /// Sanity check on input rim data (after scaling) - returns true if okay
548  bool sanityCheck();
549  //@}
550  public:
551  /**@name public methods */
552  //@{
553  /** Return row or column sections - not as much needed as it
554      once was.  These just map into single arrays */
555  inline double * solutionRegion(int section)
556  { if (!section) return rowActivityWork_; else return columnActivityWork_;};
557  inline double * djRegion(int section)
558  { if (!section) return rowReducedCost_; else return reducedCostWork_;};
559  inline double * lowerRegion(int section)
560  { if (!section) return rowLowerWork_; else return columnLowerWork_;};
561  inline double * upperRegion(int section)
562  { if (!section) return rowUpperWork_; else return columnUpperWork_;};
563  inline double * costRegion(int section)
564  { if (!section) return rowObjectiveWork_; else return objectiveWork_;};
565  /// Return region as single array
566  inline double * solutionRegion()
567  { return solution_;};
568  inline double * djRegion()
569  { return dj_;};
570  inline double * lowerRegion()
571  { return lower_;};
572  inline double * upperRegion()
573  { return upper_;};
574  inline double * costRegion()
575  { return cost_;};
576  inline Status getStatus(int sequence) const
577  {return static_cast<Status> (status_[sequence]&7);};
578  inline void setStatus(int sequence, Status status)
579  {
580    unsigned char & st_byte = status_[sequence];
581    st_byte &= ~7;
582    st_byte |= status;
583  };
584  /** Normally the first factorization does sparse coding because
585      the factorization could be singular.  This allows initial dense
586      factorization when it is known to be safe
587  */
588  void setInitialDenseFactorization(bool onOff);
589  bool  initialDenseFactorization() const;
590  /** Return sequence In or Out */
591  inline int sequenceIn() const
592  {return sequenceIn_;};
593  inline int sequenceOut() const
594  {return sequenceOut_;};
595  /** Set sequenceIn or Out */
596  inline void  setSequenceIn(int sequence)
597  { sequenceIn_=sequence;};
598  inline void  setSequenceOut(int sequence)
599  { sequenceOut_=sequence;};
600  /** Return direction In or Out */
601  inline int directionIn() const
602  {return directionIn_;};
603  inline int directionOut() const
604  {return directionOut_;};
605  /** Set directionIn or Out */
606  inline void  setDirectionIn(int direction)
607  { directionIn_=direction;};
608  inline void  setDirectionOut(int direction)
609  { directionOut_=direction;};
610  /// Returns 1 if sequence indicates column
611  inline int isColumn(int sequence) const
612  { return sequence<numberColumns_ ? 1 : 0;};
613  /// Returns sequence number within section
614  inline int sequenceWithin(int sequence) const
615  { return sequence<numberColumns_ ? sequence : sequence-numberColumns_;};
616  /// Return row or column values
617  inline double solution(int sequence)
618  { return solution_[sequence];};
619  /// Return address of row or column values
620  inline double & solutionAddress(int sequence)
621  { return solution_[sequence];};
622  inline double reducedCost(int sequence)
623   { return dj_[sequence];};
624  inline double & reducedCostAddress(int sequence)
625   { return dj_[sequence];};
626  inline double lower(int sequence)
627  { return lower_[sequence];};
628  /// Return address of row or column lower bound
629  inline double & lowerAddress(int sequence)
630  { return lower_[sequence];};
631  inline double upper(int sequence)
632  { return upper_[sequence];};
633  /// Return address of row or column upper bound
634  inline double & upperAddress(int sequence)
635  { return upper_[sequence];};
636  inline double cost(int sequence)
637  { return cost_[sequence];};
638  /// Return address of row or column cost
639  inline double & costAddress(int sequence)
640  { return cost_[sequence];};
641  /// Return original lower bound
642  inline double originalLower(int iSequence) const
643  { if (iSequence<numberColumns_) return columnLower_[iSequence]; else
644    return rowLower_[iSequence-numberColumns_];};
645  /// Return original lower bound
646  inline double originalUpper(int iSequence) const
647  { if (iSequence<numberColumns_) return columnUpper_[iSequence]; else
648    return rowUpper_[iSequence-numberColumns_];};
649  /// Theta (pivot change)
650  inline double theta() const
651  { return theta_;};
652  /// Scaling
653  const double * rowScale() const {return rowScale_;};
654  const double * columnScale() const {return columnScale_;};
655  void setRowScale(double * scale) { rowScale_ = scale;};
656  void setColumnScale(double * scale) { columnScale_ = scale;};
657  /// Return pointer to details of costs
658  inline ClpNonLinearCost * nonLinearCost() const
659  { return nonLinearCost_;};
660  //@}
661  /**@name status methods */
662  //@{
663  inline void setFakeBound(int sequence, FakeBound fakeBound)
664  {
665    unsigned char & st_byte = status_[sequence];
666    st_byte &= ~24;
667    st_byte |= fakeBound<<3;
668  };
669  inline FakeBound getFakeBound(int sequence) const
670  {return static_cast<FakeBound> ((status_[sequence]>>3)&3);};
671  inline void setRowStatus(int sequence, Status status)
672  {
673    unsigned char & st_byte = status_[sequence+numberColumns_];
674    st_byte &= ~7;
675    st_byte |= status;
676  };
677  inline Status getRowStatus(int sequence) const
678  {return static_cast<Status> (status_[sequence+numberColumns_]&7);};
679  inline void setColumnStatus(int sequence, Status status)
680  {
681    unsigned char & st_byte = status_[sequence];
682    st_byte &= ~7;
683    st_byte |= status;
684  };
685  inline Status getColumnStatus(int sequence) const
686  {return static_cast<Status> (status_[sequence]&7);};
687  inline void setPivoted( int sequence)
688  { status_[sequence] |= 32;};
689  inline void clearPivoted( int sequence)
690  { status_[sequence] &= ~32; };
691  inline bool pivoted(int sequence) const
692  {return (((status_[sequence]>>5)&1)!=0);};
693  inline void setFlagged( int sequence)
694  {
695    status_[sequence] |= 64;
696  };
697  inline void clearFlagged( int sequence)
698  {
699    status_[sequence] &= ~64;
700  };
701  inline bool flagged(int sequence) const
702  {return (((status_[sequence]>>6)&1)!=0);};
703  /** Set up status array (can be used by OsiClp).
704      Also can be used to set up all slack basis */
705  void createStatus() ;
706  inline void allSlackBasis()
707  { createStatus();};
708   
709  /// So we know when to be cautious
710  inline int lastBadIteration() const
711  {return lastBadIteration_;};
712  /// Progress flag - at present 0 bit says artificials out
713  inline int progressFlag() const
714  {return progressFlag_;};
715  /// Force re-factorization early
716  inline void forceFactorization(int value)
717  { forceFactorization_ = value;};
718  /// Raw objective value (so always minimize in primal)
719  inline double rawObjectiveValue() const
720  { return objectiveValue_;};
721  //@}
722
723////////////////// data //////////////////
724protected:
725
726  /**@name data.  Many arrays have a row part and a column part.
727   There is a single array with both - columns then rows and
728   then normally two arrays pointing to rows and columns.  The
729   single array is the owner of memory
730  */
731  //@{
732  /// Worst column primal infeasibility
733  double columnPrimalInfeasibility_;
734  /// Worst row primal infeasibility
735  double rowPrimalInfeasibility_;
736  /// Sequence of worst (-1 if feasible)
737  int columnPrimalSequence_;
738  /// Sequence of worst (-1 if feasible)
739  int rowPrimalSequence_;
740  /// Worst column dual infeasibility
741  double columnDualInfeasibility_;
742  /// Worst row dual infeasibility
743  double rowDualInfeasibility_;
744  /// Sequence of worst (-1 if feasible)
745  int columnDualSequence_;
746  /// Sequence of worst (-1 if feasible)
747  int rowDualSequence_;
748  /// Primal tolerance needed to make dual feasible (<largeTolerance)
749  double primalToleranceToGetOptimal_;
750  /// Remaining largest dual infeasibility
751  double remainingDualInfeasibility_;
752  /// Large bound value (for complementarity etc)
753  double largeValue_;
754  /// Scaling of objective
755  double objectiveScale_;
756  /// Largest error on Ax-b
757  double largestPrimalError_;
758  /// Largest error on basic duals
759  double largestDualError_;
760  /// Largest difference between input primal solution and computed
761  double largestSolutionError_;
762  /// Dual bound
763  double dualBound_;
764  /// Alpha (pivot element)
765  double alpha_;
766  /// Theta (pivot change)
767  double theta_;
768  /// Lower Bound on In variable
769  double lowerIn_;
770  /// Value of In variable
771  double valueIn_;
772  /// Upper Bound on In variable
773  double upperIn_;
774  /// Reduced cost of In variable
775  double dualIn_;
776  /// Lower Bound on Out variable
777  double lowerOut_;
778  /// Value of Out variable
779  double valueOut_;
780  /// Upper Bound on Out variable
781  double upperOut_;
782  /// Infeasibility (dual) or ? (primal) of Out variable
783  double dualOut_;
784  /// Current dual tolerance for algorithm
785  double dualTolerance_;
786  /// Current primal tolerance for algorithm
787  double primalTolerance_;
788  /// Sum of dual infeasibilities
789  double sumDualInfeasibilities_;
790  /// Sum of primal infeasibilities
791  double sumPrimalInfeasibilities_;
792  /// Weight assigned to being infeasible in primal
793  double infeasibilityCost_;
794  /// Sum of Dual infeasibilities using tolerance based on error in duals
795  double sumOfRelaxedDualInfeasibilities_;
796  /// Sum of Primal infeasibilities using tolerance based on error in primals
797  double sumOfRelaxedPrimalInfeasibilities_;
798  /// Working copy of lower bounds (Owner of arrays below)
799  double * lower_;
800  /// Row lower bounds - working copy
801  double * rowLowerWork_;
802  /// Column lower bounds - working copy
803  double * columnLowerWork_;
804  /// Working copy of upper bounds (Owner of arrays below)
805  double * upper_;
806  /// Row upper bounds - working copy
807  double * rowUpperWork_;
808  /// Column upper bounds - working copy
809  double * columnUpperWork_;
810  /// Working copy of objective (Owner of arrays below)
811  double * cost_;
812  /// Row objective - working copy
813  double * rowObjectiveWork_;
814  /// Column objective - working copy
815  double * objectiveWork_;
816  /// Useful row length arrays
817  CoinIndexedVector * rowArray_[6];
818  /// Useful column length arrays
819  CoinIndexedVector * columnArray_[6];
820  /// Sequence of In variable
821  int sequenceIn_;
822  /// Direction of In, 1 going up, -1 going down, 0 not a clude
823  int directionIn_;
824  /// Sequence of Out variable
825  int sequenceOut_;
826  /// Direction of Out, 1 to upper bound, -1 to lower bound, 0 - superbasic
827  int directionOut_;
828  /// Pivot Row
829  int pivotRow_;
830  /// Last good iteration (immediately after a re-factorization)
831  int lastGoodIteration_;
832  /// Working copy of reduced costs (Owner of arrays below)
833  double * dj_;
834  /// Reduced costs of slacks not same as duals (or - duals)
835  double * rowReducedCost_;
836  /// Possible scaled reduced costs
837  double * reducedCostWork_;
838  /// Working copy of primal solution (Owner of arrays below)
839  double * solution_;
840  /// Row activities - working copy
841  double * rowActivityWork_;
842  /// Column activities - working copy
843  double * columnActivityWork_;
844  /// Number of dual infeasibilities
845  int numberDualInfeasibilities_;
846  /// Number of dual infeasibilities (without free)
847  int numberDualInfeasibilitiesWithoutFree_;
848  /// Number of primal infeasibilities
849  int numberPrimalInfeasibilities_;
850  /// How many iterative refinements to do
851  int numberRefinements_;
852  /// dual row pivot choice
853  ClpDualRowPivot * dualRowPivot_;
854  /// primal column pivot choice
855  ClpPrimalColumnPivot * primalColumnPivot_;
856  /// Basic variables pivoting on which rows
857  int * pivotVariable_;
858  /// factorization
859  ClpFactorization * factorization_;
860  /// Row scale factors for matrix
861  // ****** get working simply then make coding more efficient
862  // on full matrix operations
863  double * rowScale_;
864  /// Saved version of solution
865  double * savedSolution_;
866  /// Column scale factors
867  double * columnScale_;
868  /// Scale flag, 0 none, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic
869  int scalingFlag_;
870  /// Number of times code has tentatively thought optimal
871  int numberTimesOptimal_;
872  /// If change has been made (first attempt at stopping looping)
873  int changeMade_;
874  /// Algorithm >0 == Primal, <0 == Dual
875  int algorithm_;
876  /** Now for some reliability aids
877      This forces re-factorization early */
878  int forceFactorization_;
879  /** Perturbation:
880      -50 to +50 - perturb by this power of ten (-6 sounds good)
881      100 - auto perturb if takes too long (1.0e-6 largest nonzero)
882      101 - we are perturbed
883      102 - don't try perturbing again
884      default is 100
885  */
886  int perturbation_;
887  /// Saved status regions
888  unsigned char * saveStatus_;
889  /** Very wasteful way of dealing with infeasibilities in primal.
890      However it will allow non-linearities and use of dual
891      analysis.  If it doesn't work it can easily be replaced.
892  */
893  ClpNonLinearCost * nonLinearCost_;
894  /** For advanced options
895      1 - Don't keep changing infeasibility weight
896      2 - Keep nonLinearCost round solves
897      4 - Force outgoing variables to exact bound (primal)
898      8 - Safe to use dense initial factorization
899  */
900  unsigned int specialOptions_;
901  /// So we know when to be cautious
902  int lastBadIteration_;
903  /// Can be used for count of fake bounds (dual) or fake costs (primal)
904  int numberFake_;
905  /// Progress flag - at present 0 bit says artificials out, 1 free in
906  int progressFlag_;
907  /// First free/super-basic variable (-1 if none)
908  int firstFree_;
909  /** For advanced use.  When doing iterative solves things can get
910      nasty so on values pass if incoming solution has largest
911      infeasibility < incomingInfeasibility throw out variables
912      from basis until largest infeasibility < allowedInfeasibility.
913      if allowedInfeasibility>= incomingInfeasibility this is
914      always possible altough you may end up with an all slack basis.
915
916      Defaults are 1.0,10.0
917  */
918  float incomingInfeasibility_;
919  float allowedInfeasibility_;
920  /// For dealing with all issues of cycling etc
921  ClpInteriorProgress * progress_;
922  //@}
923};
924//#############################################################################
925/** A function that tests the methods in the ClpInterior class. The
926    only reason for it not to be a member method is that this way it doesn't
927    have to be compiled into the library. And that's a gain, because the
928    library should be compiled with optimization on, but this method should be
929    compiled with debugging.
930
931    It also does some testing of ClpFactorization class
932 */
933void
934ClpInteriorUnitTest(const std::string & mpsDir,
935                   const std::string & netlibDir);
936
937
938/// For saving extra information to see if looping.
939class ClpInteriorProgress {
940
941public:
942
943
944  /**@name Constructors and destructor and copy */
945  //@{
946  /// Default constructor
947    ClpInteriorProgress (  );
948
949  /// Constructor from model
950    ClpInteriorProgress ( ClpInterior * model );
951
952  /// Copy constructor.
953  ClpInteriorProgress(const ClpInteriorProgress &);
954
955  /// Assignment operator. This copies the data
956    ClpInteriorProgress & operator=(const ClpInteriorProgress & rhs);
957  /// Destructor
958   ~ClpInteriorProgress (  );
959  //@}
960
961  /**@name Check progress */
962  //@{
963  /** Returns -1 if okay, -n+1 (n number of times bad) if bad but action taken,
964      >=0 if give up and use as problem status
965  */
966    int looping (  );
967  /// Start check at beginning of whileIterating
968  void startCheck();
969  /// Returns cycle length in whileIterating
970  int cycle(int in, int out,int wayIn,int wayOut); 
971
972  /// Returns previous objective (if -1) - current if (0)
973  double lastObjective(int back=1) const;
974  /// Modify objective e.g. if dual infeasible in dual
975  void modifyObjective(double value);
976  /// Returns previous iteration number (if -1) - current if (0)
977  int lastIterationNumber(int back=1) const;
978
979  //@}
980  /**@name Data  */
981#define CLP_PROGRESS 5
982  //@{
983  /// Objective values
984  double objective_[CLP_PROGRESS];
985  /// Sum of infeasibilities for algorithm
986  double infeasibility_[CLP_PROGRESS];
987  /// Pointer back to model so we can get information
988  ClpInterior * model_;
989  /// Number of infeasibilities
990  int numberInfeasibilities_[CLP_PROGRESS];
991  /// Iteration number at which occurred
992  int iterationNumber_[CLP_PROGRESS];
993  /// Number of times checked (so won't stop too early)
994  int numberTimes_;
995  /// Number of times it looked like loop
996  int numberBadTimes_;
997#define CLP_CYCLE 12
998  /// For cycle checking
999  int in_[CLP_CYCLE];
1000  int out_[CLP_CYCLE];
1001  char way_[CLP_CYCLE];
1002  //@}
1003};
1004#endif
Note: See TracBrowser for help on using the repository browser.