Changeset 214 for branches


Ignore:
Timestamp:
Oct 2, 2003 2:55:27 PM (16 years ago)
Author:
forrest
Message:

More pdco

Location:
branches/pre
Files:
1 added
6 edited

Legend:

Unmodified
Added
Removed
  • branches/pre/ClpInterior.cpp

    r213 r214  
    324324  return true;
    325325}
     326/* Loads a problem (the constraints on the
     327   rows are given by lower and upper bounds). If a pointer is 0 then the
     328   following values are the default:
     329   <ul>
     330   <li> <code>colub</code>: all columns have upper bound infinity
     331   <li> <code>collb</code>: all columns have lower bound 0
     332   <li> <code>rowub</code>: all rows have upper bound infinity
     333   <li> <code>rowlb</code>: all rows have lower bound -infinity
     334   <li> <code>obj</code>: all variables have 0 objective coefficient
     335   </ul>
     336*/
     337void
     338ClpInterior::loadProblem (  const ClpMatrixBase& matrix,
     339                    const double* collb, const double* colub,   
     340                    const double* obj,
     341                    const double* rowlb, const double* rowub,
     342                    const double * rowObjective)
     343{
     344  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
     345                        rowObjective);
     346}
     347void
     348ClpInterior::loadProblem (  const CoinPackedMatrix& matrix,
     349                    const double* collb, const double* colub,   
     350                    const double* obj,
     351                    const double* rowlb, const double* rowub,
     352                    const double * rowObjective)
     353{
     354  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
     355                        rowObjective);
     356}
     357
     358/* Just like the other loadProblem() method except that the matrix is
     359   given in a standard column major ordered format (without gaps). */
     360void
     361ClpInterior::loadProblem (  const int numcols, const int numrows,
     362                    const CoinBigIndex* start, const int* index,
     363                    const double* value,
     364                    const double* collb, const double* colub,   
     365                    const double* obj,
     366                    const double* rowlb, const double* rowub,
     367                    const double * rowObjective)
     368{
     369  ClpModel::loadProblem(numcols, numrows, start, index, value,
     370                          collb, colub, obj, rowlb, rowub,
     371                          rowObjective);
     372}
     373void
     374ClpInterior::loadProblem (  const int numcols, const int numrows,
     375                           const CoinBigIndex* start, const int* index,
     376                           const double* value,const int * length,
     377                           const double* collb, const double* colub,   
     378                           const double* obj,
     379                           const double* rowlb, const double* rowub,
     380                           const double * rowObjective)
     381{
     382  ClpModel::loadProblem(numcols, numrows, start, index, value, length,
     383                          collb, colub, obj, rowlb, rowub,
     384                          rowObjective);
     385}
     386// Read an mps file from the given filename
     387int
     388ClpInterior::readMps(const char *filename,
     389            bool keepNames,
     390            bool ignoreErrors)
     391{
     392  int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
     393  return status;
     394}
     395#include "ClpPdco.hpp"
     396/* Pdco algorithm - see ClpPdco.hpp for method */
     397int
     398ClpInterior::pdco()
     399{
     400  return ((ClpPdco *) this)->pdco();
     401}
     402// ** Temporary version
     403int 
     404ClpInterior::pdco( Lsqr *lsqr, Options &options, Info &info, Outfo &outfo)
     405{
     406  return ((ClpPdco *) this)->pdco(lsqr,options,info,outfo);
     407}
  • branches/pre/ClpPdco.cpp

    r212 r214  
    22// Corporation and others.  All Rights Reserved.
    33
    4 
    5 /* Notes on implementation of dual simplex algorithm.
    6 
    7    When dual feasible:
    8 
    9    If primal feasible, we are optimal.  Otherwise choose an infeasible
    10    basic variable to leave basis (normally going to nearest bound) (B).  We
    11    now need to find an incoming variable which will leave problem
    12    dual feasible so we get the row of the tableau corresponding to
    13    the basic variable (with the correct sign depending if basic variable
    14    above or below feasibility region - as that affects whether reduced
    15    cost on outgoing variable has to be positive or negative).
    16 
    17    We now perform a ratio test to determine which incoming variable will
    18    preserve dual feasibility (C).  If no variable found then problem
    19    is infeasible (in primal sense).  If there is a variable, we then
    20    perform pivot and repeat.  Trivial?
    21 
    22    -------------------------------------------
    23 
    24    A) How do we get dual feasible?  If all variables have bounds then
    25    it is trivial to get feasible by putting non-basic variables to
    26    correct bounds.  OSL did not have a phase 1/phase 2 approach but
    27    instead effectively put fake bounds on variables and this is the
    28    approach here, although I had hoped to make it cleaner.
    29 
    30    If there is a weight of X on getting dual feasible:
    31      Non-basic variables with negative reduced costs are put to
    32      lesser of their upper bound and their lower bound + X.
    33      Similarly, mutatis mutandis, for positive reduced costs.
    34 
    35    Free variables should normally be in basis, otherwise I have
    36    coding which may be able to come out (and may not be correct).
    37 
    38    In OSL, this weight was changed heuristically, here at present
    39    it is only increased if problem looks finished.  If problem is
    40    feasible I check for unboundedness.  If not unbounded we
    41    could play with going into primal.  As long as weights increase
    42    any algorithm would be finite.
     4/* Pdco algorithm
    435   
    44    B) Which outgoing variable to choose is a virtual base class.
    45    For difficult problems steepest edge is preferred while for
    46    very easy (large) problems we will need partial scan.
    47 
    48    C) Sounds easy, but this is hardest part of algorithm.
    49       1) Instead of stopping at first choice, we may be able
    50       to flip that variable to other bound and if objective
    51       still improving choose again.  These mini iterations can
    52       increase speed by orders of magnitude but we may need to
    53       go to more of a bucket choice of variable rather than looking
    54       at them one by one (for speed).
    55       2) Accuracy.  Reduced costs may be of wrong sign but less than
    56       tolerance.  Pivoting on these makes objective go backwards.
    57       OSL modified cost so a zero move was made, Gill et al
    58       (in primal analogue) modified so a strictly positive move was
    59       made.  It is not quite as neat in dual but that is what we
    60       try and do.  The two problems are that re-factorizations can
    61       change reduced costs above and below tolerances and that when
    62       finished we need to reset costs and try again.
    63       3) Degeneracy.  Gill et al helps but may not be enough.  We
    64       may need more.  Also it can improve speed a lot if we perturb
    65       the costs significantly. 
    66 
    67   References:
    68      Forrest and Goldfarb, Steepest-edge simplex algorithms for
    69        linear programming - Mathematical Programming 1992
    70      Forrest and Tomlin, Implementing the simplex method for
    71        the Optimization Subroutine Library - IBM Systems Journal 1992
    72      Gill, Murray, Saunders, Wright A Practical Anti-Cycling
    73        Procedure for Linear and Nonlinear Programming SOL report 1988
    74 
    75 
    76   TODO:
    77  
    78   a) Better recovery procedures.  At present I never check on forward
    79      progress.  There is checkpoint/restart with reducing
    80      re-factorization frequency, but this is only on singular
    81      factorizations.
    82   b) Fast methods for large easy problems (and also the option for
    83      the code to automatically choose which method).
    84   c) We need to be able to stop in various ways for OSI - this
    85      is fairly easy.
    86 
    87  */
     6Method
     7
     8
     9*/
     10
    8811
    8912
     
    9417#include "CoinHelperFunctions.hpp"
    9518#include "ClpPdco.hpp"
    96 #include "ClpFactorization.hpp"
     19#include "CoinDenseVector.hpp"
     20#include "ClpHelperFunctions.hpp"
    9721#include "CoinPackedMatrix.hpp"
    98 #include "CoinIndexedVector.hpp"
    99 #include "CoinWarmStartBasis.hpp"
    100 #include "ClpDualRowDantzig.hpp"
    101 #include "ClpPlusMinusOneMatrix.hpp"
    10222#include "ClpMessage.hpp"
    10323#include <cfloat>
     
    10626#include <stdio.h>
    10727#include <iostream>
    108 //#define CLP_DEBUG 1
    109 #define CHECK_DJ
    110 // dual
    111 int ClpPdco::dual (int ifValuesPass )
     28
     29int
     30ClpPdco::pdco()
    11231{
    113 
    114   /* *** Method
    115      This is a vanilla version of dual simplex.
    116 
    117      It tries to be a single phase approach with a weight of 1.0 being
    118      given to getting optimal and a weight of dualBound_ being
    119      given to getting dual feasible.  In this version I have used the
    120      idea that this weight can be thought of as a fake bound.  If the
    121      distance between the lower and upper bounds on a variable is less
    122      than the feasibility weight then we are always better off flipping
    123      to other bound to make dual feasible.  If the distance is greater
    124      then we make up a fake bound dualBound_ away from one bound.
    125      If we end up optimal or primal infeasible, we check to see if
    126      bounds okay.  If so we have finished, if not we increase dualBound_
    127      and continue (after checking if unbounded). I am undecided about
    128      free variables - there is coding but I am not sure about it.  At
    129      present I put them in basis anyway.
    130 
    131      The code is designed to take advantage of sparsity so arrays are
    132      seldom zeroed out from scratch or gone over in their entirety.
    133      The only exception is a full scan to find outgoing variable.  This
    134      will be changed to keep an updated list of infeasibilities (or squares
    135      if steepest edge).  Also on easy problems we don't need full scan - just
    136      pick first reasonable.
    137 
    138      One problem is how to tackle degeneracy and accuracy.  At present
    139      I am using the modification of costs which I put in OSL and which was
    140      extended by Gill et al.  I am still not sure of the exact details.
    141 
    142      The flow of dual is three while loops as follows:
    143 
    144      while (not finished) {
    145 
    146        while (not clean solution) {
    147 
    148           Factorize and/or clean up solution by flipping variables so
    149           dual feasible.  If looks finished check fake dual bounds.
    150           Repeat until status is iterating (-1) or finished (0,1,2)
    151 
    152        }
    153 
    154        while (status==-1) {
    155 
    156          Iterate until no pivot in or out or time to re-factorize.
    157 
    158          Flow is:
    159 
    160          choose pivot row (outgoing variable).  if none then
    161          we are primal feasible so looks as if done but we need to
    162          break and check bounds etc.
    163 
    164          Get pivot row in tableau
    165 
    166          Choose incoming column.  If we don't find one then we look
    167          primal infeasible so break and check bounds etc.  (Also the
    168          pivot tolerance is larger after any iterations so that may be
    169          reason)
    170 
    171          If we do find incoming column, we may have to adjust costs to
    172          keep going forwards (anti-degeneracy).  Check pivot will be stable
    173          and if unstable throw away iteration (we will need to implement
    174          flagging of basic variables sometime) and break to re-factorize.
    175          If minor error re-factorize after iteration.
    176 
    177          Update everything (this may involve flipping variables to stay
    178          dual feasible.
    179 
    180        }
    181 
    182      }
    183 
    184      At present we never check we are going forwards.  I overdid that in
    185      OSL so will try and make a last resort.
    186 
    187      Needs partial scan pivot out option.
    188      Needs dantzig, uninitialized and full steepest edge options (can still
    189      use partial scan)
    190 
    191      May need other anti-degeneracy measures, especially if we try and use
    192      loose tolerances as a way to solve in fewer iterations.
    193 
    194      I like idea of dynamic scaling.  This gives opportunity to decouple
    195      different implications of scaling for accuracy, iteration count and
    196      feasibility tolerance.
    197 
     32  printf("Dummy pdco solve\n");
     33  return 0;
     34}
     35// ** Temporary version
     36int 
     37ClpPdco::pdco( Lsqr *lsqr, Options &options, Info &info, Outfo &outfo)
     38{
     39//    D1, D2 are positive-definite diagonal matrices defined from d1, d2.
     40//           In particular, d2 indicates the accuracy required for
     41//           satisfying each row of Ax = b.
     42//
     43// D1 and D2 (via d1 and d2) provide primal and dual regularization
     44// respectively.  They ensure that the primal and dual solutions
     45// (x,r) and (y,z) are unique and bounded.
     46//
     47// A scalar d1 is equivalent to d1 = ones(n,1), D1 = diag(d1).
     48// A scalar d2 is equivalent to d2 = ones(m,1), D2 = diag(d2).
     49// Typically, d1 = d2 = 1e-4.
     50// These values perturb phi(x) only slightly  (by about 1e-8) and request
     51// that A*x = b be satisfied quite accurately (to about 1e-8).
     52// Set d1 = 1e-4, d2 = 1 for least-squares problems with bound constraints.
     53// The problem is then
     54//
     55//    minimize    phi(x) + 1/2 norm(d1*x)^2 + 1/2 norm(A*x - b)^2
     56//    subject to  bl <= x <= bu.
     57//
     58// More generally, d1 and d2 may be n and m vectors containing any positive
     59// values (preferably not too small, and typically no larger than 1).
     60// Bigger elements of d1 and d2 improve the stability of the solver.
     61//
     62// At an optimal solution, if x(j) is on its lower or upper bound,
     63// the corresponding z(j) is positive or negative respectively.
     64// If x(j) is between its bounds, z(j) = 0.
     65// If bl(j) = bu(j), x(j) is fixed at that value and z(j) may have
     66// either sign.
     67//
     68// Also, r and y satisfy r = D2 y, so that Ax + D2^2 y = b.
     69// Thus if d2(i) = 1e-4, the i-th row of Ax = b will be satisfied to
     70// approximately 1e-8.  This determines how large d2(i) can safely be.
     71//
     72//
     73// EXTERNAL FUNCTIONS:
     74// options         = pdcoSet;                  provided with pdco.m
     75// [obj,grad,hess] = pdObj( x );               provided by user
     76//               y = pdMat( name,mode,m,n,x ); provided by user if pdMat
     77//                                             is a string, not a matrix
     78//
     79// INPUT ARGUMENTS:
     80// pdObj      is a string containing the name of a function pdObj.m
     81//            or a function_handle for such a function
     82//            such that [obj,grad,hess] = pdObj(x) defines
     83//            obj  = phi(x)              : a scalar,
     84//            grad = gradient of phi(x)  : an n-vector,
     85//            hess = diag(Hessian of phi): an n-vector.
     86//         Examples:
     87//            If phi(x) is the linear function c"x, pdObj should return
     88//               [obj,grad,hess] = [c"*x, c, zeros(n,1)].
     89//            If phi(x) is the entropy function E(x) = sum x(j) log x(j),
     90//               [obj,grad,hess] = [E(x), log(x)+1, 1./x].
     91// pdMat      may be an ifexplicit m x n matrix A (preferably sparse!),
     92//            or a string containing the name of a function pdMat.m
     93//            or a function_handle for such a function
     94//            such that y = pdMat( name,mode,m,n,x )
     95//            returns   y = A*x (mode=1)  or  y = A"*x (mode=2).
     96//            The input parameter "name" will be the string pdMat.
     97// b          is an m-vector.
     98// bl         is an n-vector of lower bounds.  Non-existent bounds
     99//            may be represented by bl(j) = -Inf or bl(j) <= -1e+20.
     100// bu         is an n-vector of upper bounds.  Non-existent bounds
     101//            may be represented by bu(j) =  Inf or bu(j) >=  1e+20.
     102// d1, d2     may be positive scalars or positive vectors (see above).
     103// options    is a structure that may be set and altered by pdcoSet
     104//            (type help pdcoSet).
     105// x0, y0, z0 provide an initial solution.
     106// xsize, zsize are estimates of the biggest x and z at the solution.
     107//            They are used to scale (x,y,z).  Good estimates
     108//            should improve the performance of the barrier method.
     109//
     110//
     111// OUTPUT ARGUMENTS:
     112// x          is the primal solution.
     113// y          is the dual solution associated with Ax + D2 r = b.
     114// z          is the dual solution associated with bl <= x <= bu.
     115// inform = 0 if a solution is found;
     116//        = 1 if too many iterations were required;
     117//        = 2 if the linesearch failed too often.
     118// PDitns     is the number of Primal-Dual Barrier iterations required.
     119// CGitns     is the number of Conjugate-Gradient  iterations required
     120//            if an iterative solver is used (LSQR).
     121// time       is the cpu time used.
     122//----------------------------------------------------------------------
     123
     124// PRIVATE FUNCTIONS:
     125//    pdxxxbounds
     126//    pdxxxdistrib
     127//    pdxxxlsqr
     128//    pdxxxlsqrmat
     129//    pdxxxmat
     130//    pdxxxmerit
     131//    pdxxxresid1
     132//    pdxxxresid2
     133//    pdxxxstep
     134//
     135// GLOBAL VARIABLES:
     136//    global pdDDD1 pdDDD2 pdDDD3
     137//
     138//
     139// NOTES:
     140// The matrix A should be reasonably well scaled: norm(A,inf) =~ 1.
     141// The vector b and objective phi(x) may be of any size, but ensure that
     142// xsize and zsize are reasonably close to norm(x,inf) and norm(z,inf)
     143// at the solution.
     144//
     145// The files defining pdObj  and pdMat
     146// must not be called Fname.m or Aname.m!!
     147//
     148//
     149// AUTHOR:
     150//    Michael Saunders, Systems Optimization Laboratory (SOL),
     151//    Stanford University, Stanford, California, USA.
     152//    saunders@stanford.edu
     153//
     154// CONTRIBUTORS:
     155//    Byunggyoo Kim, SOL, Stanford University.
     156//    bgkim@stanford.edu
     157//
     158// DEVELOPMENT:
     159// 20 Jun 1997: Original version of pdsco.m derived from pdlp0.m.
     160// 29 Sep 2002: Original version of pdco.m  derived from pdsco.m.
     161//              Introduced D1, D2 in place of gamma*I, delta*I
     162//              and allowed for general bounds bl <= x <= bu.
     163// 06 Oct 2002: Allowed for fixed variabes: bl(j) = bu(j) for any j.
     164// 15 Oct 2002: Eliminated some work vectors (since m, n might be LARGE).
     165//              Modularized residuals, linesearch
     166// 16 Oct 2002: pdxxx..., pdDDD... names rationalized.
     167//              pdAAA eliminated (global copy of A).
     168//              Aname is now used directly as an ifexplicit A or a function.
     169//              NOTE: If Aname is a function, it now has an extra parameter.
     170// 23 Oct 2002: Fname and Aname can now be function handles.
     171// 01 Nov 2002: Bug fixed in feval in pdxxxmat.
     172//-----------------------------------------------------------------------
     173
     174//  global pdDDD1 pdDDD2 pdDDD3
     175  real inf = 1.0e30;
     176  real eps = 1.0e-15;
     177  real atolold, r3ratio, Pinf, Dinf, Cinf, Cinf0;
     178 
     179  printf("\n   --------------------------------------------------------");
     180  printf("\n   pdco.m                            Version of 01 Nov 2002");
     181  printf("\n   Primal-dual barrier method to minimize a convex function");
     182  printf("\n   subject to linear constraints Ax + r = b,  bl <= x <= bu");
     183  printf("\n   --------------------------------------------------------\n");
     184
     185  int m = model->rowSize();
     186  int n = model->colSize();
     187  bool ifexplicit = true;
     188
     189  CoinDenseVector b(m, model->rhs); 
     190  CoinDenseVector x(n, model->x);
     191  CoinDenseVector y(m, model->y);
     192  CoinDenseVector z(n, model->dj);
     193  //delete old arrays
     194  delete model->rhs;
     195  delete model->x;
     196  delete model->y;
     197  delete model->dj;
     198
     199  real normb  = b.infNorm();
     200  real normx0 = x.infNorm();
     201  real normy0 = y.infNorm();
     202  real normz0 = z.infNorm();
     203
     204  printf("\nmax |b | = %8g     max |x0| = %8g", normb , normx0);
     205  printf(                "      xsize   = %8g", model->xsize);
     206  printf("\nmax |y0| = %8g     max |z0| = %8g", normy0, normz0);
     207  printf(                "      zsize   = %8g", model->zsize);
     208
     209  //---------------------------------------------------------------------
     210  // Initialize.
     211  //---------------------------------------------------------------------
     212  //true   = 1;
     213  //false  = 0;
     214  //zn     = zeros(n,1);
     215  int nb     = n + m;
     216  int nkkt   = nb;
     217  int CGitns = 0;
     218  int inform = 0;
     219  //--------------------------------------------------------------------- 
     220  //  Only allow scalar d1, d2 for now
     221  //---------------------------------------------------------------------
     222  /*
     223  if (model->d1->size()==1)
     224      d1->resize(n, d1->getElements()[0]);  // Allow scalar d1, d2
     225  if (model->d2->size()==1)
     226      d2->resize(m, d2->getElements()[0]);  // to mean dk * unit vector
     227 */
     228  real d1 = model->d1;
     229  real d2 = model->d2;
     230
     231  //---------------------------------------------------------------------
     232  // Grab input options.
     233  //---------------------------------------------------------------------
     234  int  maxitn    = options.MaxIter;
     235  real featol    = options.FeaTol;
     236  real opttol    = options.OptTol;
     237  real steptol   = options.StepTol;
     238  int  stepSame  = 1;  /* options.StepSame;   // 1 means stepx == stepz */
     239  real x0min     = options.x0min;
     240  real z0min     = options.z0min;
     241  real mu0       = options.mu0;
     242  int  LSproblem = options.LSproblem;  // See below
     243  int  LSmethod  = options.LSmethod;   // 1=Cholesky    2=QR    3=LSQR
     244  int  itnlim    = options.LSQRMaxIter * min(m,n);
     245  real atol1     = options.LSQRatol1;  // Initial  atol
     246  real atol2     = options.LSQRatol2;  // Smallest atol,unless atol1 is smaller
     247  real conlim    = options.LSQRconlim;
     248  int  wait      = options.wait;
     249
     250  // LSproblem:
     251  //  1 = dy          2 = dy shifted, DLS
     252  // 11 = s          12 =  s shifted, DLS    (dx = Ds)
     253  // 21 = dx
     254  // 31 = 3x3 system, symmetrized by Z^{1/2}
     255  // 32 = 2x2 system, symmetrized by X^{1/2}
     256
     257  //---------------------------------------------------------------------
     258  // Set other parameters.
     259  //---------------------------------------------------------------------
     260  int  kminor    = 0;      // 1 stops after each iteration
     261  real eta       = 1e-4;   // Linesearch tolerance for "sufficient descent"
     262  real maxf      = 10;     // Linesearch backtrack limit (function evaluations)
     263  real maxfail   = 1;      // Linesearch failure limit (consecutive iterations)
     264  real bigcenter = 1e+3;   // mu is reduced if center < bigcenter.
     265
     266  // Parameters for LSQR.
     267  real atolmin   = eps;    // Smallest atol if linesearch back-tracks
     268  real btol      = 0;      // Should be small (zero is ok)
     269  real show      = false;  // Controls lsqr iteration log
     270  /*
     271  real gamma     = d1->infNorm();
     272  real delta     = d2->infNorm();
    198273  */
    199   algorithm_ = -1;
    200 
    201   // save data
    202   ClpDataSave data = saveData();
    203 
    204   // Save so can see if doing after primal
    205   int initialStatus=problemStatus_;
    206 
    207   // If values pass then save given duals round check solution
    208   double * saveDuals = NULL;
    209   if (ifValuesPass) {
    210     saveDuals = new double [numberRows_+numberColumns_];
    211     memcpy(saveDuals,dual_,numberRows_*sizeof(double));
     274  real gamma = d1;
     275  real delta = d2;
     276
     277  printf("\n\nx0min    = %8g     featol   = %8.1e", x0min, featol);
     278  printf(                  "      d1max   = %8.1e", gamma);
     279  printf(  "\nz0min    = %8g     opttol   = %8.1e", z0min, opttol);
     280  printf(                  "      d2max   = %8.1e", delta);
     281  printf(  "\nmu0      = %8.1e     steptol  = %8g", mu0  , steptol);
     282  printf(                  "     bigcenter= %8g"  , bigcenter);
     283
     284  printf("\n\nLSQR:");
     285  printf("\natol1    = %8.1e     atol2    = %8.1e", atol1 , atol2 );
     286  printf(                  "      btol    = %8.1e", btol );
     287  printf("\nconlim   = %8.1e     itnlim   = %8d"  , conlim, itnlim);
     288  printf(                  "      show    = %8g"  , show );
     289
     290// LSmethod  = 3;  ////// Hardwire LSQR
     291// LSproblem = 1;  ////// and LS problem defining "dy".
     292/*
     293  if wait
     294    printf("\n\nReview parameters... then type "return"\n")
     295    keyboard
     296  end
     297*/
     298  if (eta < 0)
     299    printf("\n\nLinesearch disabled by eta < 0");
     300
     301  //---------------------------------------------------------------------
     302  // All parameters have now been set.
     303  //---------------------------------------------------------------------
     304  //real time    = cputime();
     305  double time = now();
     306  bool useChol = (LSmethod == 1);
     307  bool useQR   = (LSmethod == 2);
     308  bool direct  = (LSmethod <= 2 && ifexplicit);
     309  char solver[6];
     310  strcpy(solver,"  LSQR");
     311
     312
     313  //---------------------------------------------------------------------
     314  // Categorize bounds and allow for fixed variables by modifying b.
     315  //---------------------------------------------------------------------
     316
     317  int nlow, nupp, nfix;
     318  int *bptrs[3] = {0};
     319  model->getBoundTypes(&nlow, &nupp, &nfix, bptrs );
     320  int *low = bptrs[0];
     321  int *upp = bptrs[1];
     322  int *fix = bptrs[2];
     323
     324  int nU = n;
     325  if (nupp ==0) nU = 1;   //Make dummy vectors if no Upper bounds
     326
     327  //---------------------------------------------------------------------
     328  //  Get pointers to local copy of model bounds
     329  //---------------------------------------------------------------------
     330
     331  CoinDenseVector bl(n, model->colLower);
     332  real *bl_elts = bl.getElements();
     333  CoinDenseVector bu(nU, model->colUpper);  // this is dummy if no UB
     334  real *bu_elts = bu.getElements();
     335  delete model->colLower;
     336  delete model->colUpper;
     337
     338  CoinDenseVector r1(m, 0.0);
     339  real *r1_elts = r1.getElements();
     340  CoinDenseVector x1(n, 0.0);
     341  real *x1_elts = x1.getElements();
     342
     343  if (nfix > 0){
     344    for (int k=0; k<nfix; k++)
     345      x1_elts[fix[k]] = bl[fix[k]];
     346    model->matVecMult(1, r1, x1);
     347    b = b - r1;
     348      // At some stage, might want to look at normfix = norm(r1,inf);
    212349  }
     350
     351  //---------------------------------------------------------------------
     352  // Scale the input data.
     353  // The scaled variables are
     354  //    xbar     = x/beta,
     355  //    ybar     = y/zeta,
     356  //    zbar     = z/zeta.
     357  // Define
     358  //    theta    = beta*zeta;
     359  // The scaled function is
     360  //    phibar   = ( 1   /theta) fbar(beta*xbar),
     361  //    gradient = (beta /theta) grad,
     362  //    Hessian  = (beta2/theta) hess.
     363  //---------------------------------------------------------------------
     364  real beta = model->xsize;   if (beta==0) beta = 1;   // beta scales b, x.
     365  real zeta = model->zsize;   if (zeta==0) zeta = 1;   // zeta scales y, z.
     366  real theta  = beta*zeta;                            // theta scales obj.
     367  // (theta could be anything, but theta = beta*zeta makes
     368  // scaled grad = grad/zeta = 1 approximately if zeta is chosen right.)
     369
     370  for (int k=0; k<nlow; k++)
     371    bl_elts[low[k]] = bl_elts[low[k]]/beta;
     372  for (int k=0; k<nupp; k++)
     373    bu_elts[upp[k]] = bu_elts[upp[k]]/beta;
     374  d1     = d1 * ( beta/sqrt(theta) );
     375  d2     = d2 * ( sqrt(theta)/beta );
     376
     377  real beta2  = beta*beta;
     378  b.scale( (1.0/beta) );
     379  y.scale( (1.0/zeta) );
     380  x.scale( (1.0/beta) );
     381  z.scale( (1.0/zeta) );
     382 
     383  //---------------------------------------------------------------------
     384  // Initialize vectors that are not fully used if bounds are missing.
     385  //---------------------------------------------------------------------
     386  CoinDenseVector rL(n, 0.0);
     387  CoinDenseVector cL(n, 0.0);
     388  CoinDenseVector z1(n, 0.0);
     389  CoinDenseVector dx1(n, 0.0);
     390  CoinDenseVector dz1(n, 0.0);
     391  CoinDenseVector r2(n, 0.0);
     392
     393  // Assign upper bd regions (dummy if no UBs)
     394
     395  CoinDenseVector rU(nU, 0.0);
     396  CoinDenseVector cU(nU, 0.0);
     397  CoinDenseVector x2(nU, 0.0);
     398  CoinDenseVector z2(nU, 0.0);
     399  CoinDenseVector dx2(nU, 0.0);   
     400  CoinDenseVector dz2(nU, 0.0);
     401
     402  //---------------------------------------------------------------------
     403  // Initialize x, y, z, objective, etc.
     404  //---------------------------------------------------------------------
     405  CoinDenseVector dx(n, 0.0);
     406  CoinDenseVector dy(m, 0.0);
     407  CoinDenseVector Pr(m);
     408  CoinDenseVector D(n);
     409  real *D_elts = D.getElements();
     410  CoinDenseVector w(n);
     411  real *w_elts = w.getElements();
     412  CoinDenseVector rhs(m+n);
     413   
     414
     415  //---------------------------------------------------------------------
     416  // Pull out the element array pointers for efficiency
     417  //---------------------------------------------------------------------
     418  real *x_elts  = x.getElements();
     419  real *x2_elts = x2.getElements();
     420  real *z_elts  = z.getElements();
     421  real *z1_elts = z1.getElements();
     422  real *z2_elts = z2.getElements();
     423
     424  for (int k=0; k<nlow; k++){
     425    x_elts[low[k]]  = max( x_elts[low[k]], bl[low[k]]);
     426    x1_elts[low[k]] = max( x_elts[low[k]] - bl[low[k]], x0min  );
     427    z1_elts[low[k]] = max( z_elts[low[k]], z0min  );
     428  }
     429  for (int k=0; k<nupp; k++){
     430    x_elts[upp[k]]  = min( x_elts[upp[k]], bu[upp[k]]);
     431    x2_elts[upp[k]] = max(bu[upp[k]] -  x_elts[upp[k]], x0min  );
     432    z2_elts[upp[k]] = max(-z_elts[upp[k]], z0min  );
     433  }
     434  //////////////////// Assume hessian is diagonal. //////////////////////
     435
     436//  [obj,grad,hess] = feval( Fname, (x*beta) );
     437  x.scale(beta);
     438  real obj = model->getObj(x);
     439  CoinDenseVector grad(n);
     440  model->getGrad(x, grad);
     441  CoinDenseVector H(n);
     442  model->getHessian(x ,H);
     443  x.scale((1.0/beta));
    213444 
    214   // sanity check
    215   // initialize - no values pass and algorithm_ is -1
    216   // put in standard form (and make row copy)
    217   // create modifiable copies of model rim and do optional scaling
    218   // If problem looks okay
    219   // Do initial factorization
    220   // If user asked for perturbation - do it
    221   if (!startup(0)) {
    222     // looks okay
    223    
    224     // If values pass then scale pi
    225     if (ifValuesPass) {
    226       if (problemStatus_&&perturbation_<100)
    227         perturb();
    228       int i;
    229       if (scalingFlag_>0) {
    230         for (i=0;i<numberRows_;i++) {
    231           dual_[i] = saveDuals[i]/rowScale_[i];
    232         }
    233       } else {
    234         memcpy(dual_,saveDuals,numberRows_*sizeof(double));
     445  real * g_elts=grad.getElements();
     446  real * H_elts=H.getElements();
     447
     448  obj /= theta;                       // Scaled obj.
     449  grad = grad*(beta/theta) + (d1*d1)*x; // grad includes x regularization.
     450  H  = H*(beta2/theta) + (d1*d1)      ;     // H    includes x regularization.
     451 
     452
     453  /*---------------------------------------------------------------------
     454  // Compute primal and dual residuals:
     455  // r1 =  b - Aprod(x) - d2*d2*y;
     456  // r2 =  grad - Atprod(y) + z2 - z1; 
     457  //  rL =  bl - x + x1;
     458  //  rU =  x + x2 - bu; */
     459  //---------------------------------------------------------------------
     460  //  [r1,r2,rL,rU,Pinf,Dinf] = ...
     461  //      pdxxxresid1( Aname,fix,low,upp, ...
     462  //                   b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 );
     463  pdxxxresid1( model, nlow, nupp, nfix, low, upp, fix,
     464               b,bl_elts,bu_elts,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2,
     465               r1, r2, &Pinf, &Dinf);
     466  //---------------------------------------------------------------------
     467  // Initialize mu and complementarity residuals:
     468  //    cL   = mu*e - X1*z1.
     469  //    cU   = mu*e - X2*z2.
     470  //
     471  // 25 Jan 2001: Now that b and obj are scaled (and hence x,y,z),
     472  //              we should be able to use mufirst = mu0 (absolute value).
     473  //              0.1 worked poorly on StarTest1 with x0min = z0min = 0.1.
     474  // 29 Jan 2001: We might as well use mu0 = x0min * z0min;
     475  //              so that most variables are centered after a warm start.
     476  // 29 Sep 2002: Use mufirst = mu0*(x0min * z0min),
     477  //              regarding mu0 as a scaling of the initial center.
     478  //---------------------------------------------------------------------
     479  //  real mufirst = mu0*(x0min * z0min);
     480  real mufirst = mu0;   // revert to absolute value
     481  real mulast  = 0.1 * opttol;
     482  mulast  = min( mulast, mufirst );
     483  real mu      = mufirst;
     484  real center,  fmerit;
     485  pdxxxresid2( mu, nlow, nupp, low,upp, cL, cU, x1, x2,
     486                        z1, z2, &center, &Cinf, &Cinf0 );
     487  fmerit = pdxxxmerit(nlow, nupp, low, upp, r1, r2, rL, rU, cL, cU );
     488
     489  // Initialize other things.
     490
     491  bool  precon   = true;       
     492  real PDitns    = 0;
     493  bool converged = false;
     494  real atol      = atol1;
     495  atol2     = max( atol2, atolmin );
     496  atolmin   = atol2;
     497  //  pdDDD2    = d2;    // Global vector for diagonal matrix D2
     498
     499  //  Iteration log.
     500
     501  real stepx   = 0;
     502  real stepz   = 0;
     503  int nf      = 0;
     504  int itncg   = 0;
     505  int nfail   = 0;
     506
     507  printf("\n\nItn   mu stepx stepz  Pinf  Dinf");
     508  printf("  Cinf   Objective    nf  center");
     509  if (direct) {
     510    printf("\n");
     511  }else{
     512    printf("  atol   solver   Inexact\n");
     513  }
     514
     515  real regx = (d1*x).twoNorm();
     516  real regy = (d2*y).twoNorm();
     517  //  regterm = twoNorm(d1.*x)^2  +  norm(d2.*y)^2;
     518  real regterm = regx*regx + regy*regy;
     519  real objreg  = obj  +  0.5*regterm;
     520  real objtrue = objreg * theta;
     521
     522  printf("\n%3g                 ", PDitns        );
     523  printf("%6.1f%6.1f" , log10(Pinf ), log10(Dinf));
     524  printf("%6.1f%15.7e", log10(Cinf0), objtrue    );
     525  printf("   %8.1f\n"   , center                   );
     526  /*
     527  if kminor
     528    printf("\n\nStart of first minor itn...\n");
     529    keyboard
     530  end
     531  */
     532  //---------------------------------------------------------------------
     533  // Main loop.
     534  //---------------------------------------------------------------------
     535
     536  //  while (converged) {
     537  while(PDitns < maxitn) {
     538    PDitns = PDitns + 1;
     539
     540    // 31 Jan 2001: Set atol according to progress, a la Inexact Newton.
     541    // 07 Feb 2001: 0.1 not small enough for Satellite problem.  Try 0.01.
     542    // 25 Apr 2001: 0.01 seems wasteful for Star problem.
     543    //              Now that starting conditions are better, go back to 0.1.
     544
     545    real r3norm = max(Pinf,   max(Dinf,  Cinf));
     546    atol   = min(atol,  r3norm*0.1);
     547    atol   = max(atol,  atolmin   );
     548    info.r3norm = r3norm;
     549
     550    //-------------------------------------------------------------------
     551    //  Define a damped Newton iteration for solving f = 0,
     552    //  keeping  x1, x2, z1, z2 > 0.  We eliminate dx1, dx2, dz1, dz2
     553    //  to obtain the system
     554    //
     555    //     [-H2  A"  ] [ dx ] = [ w ],   H2 = H + D1^2 + X1inv Z1 + X2inv Z2,
     556    //     [ A   D2^2] [ dy ] = [ r1]    w  = r2 - X1inv(cL + Z1 rL)
     557    //                                           + X2inv(cU + Z2 rU),
     558    //
     559    //  which is equivalent to the least-squares problem
     560    //
     561    //     min || [ D A"]dy  -  [  D w   ] ||,   D = H2^{-1/2}.         (*)
     562    //         || [  D2 ]       [D2inv r1] ||
     563    //-------------------------------------------------------------------
     564    for (int k=0; k<nlow; k++)
     565      H_elts[low[k]]  = H_elts[low[k]] + z1[low[k]]/x1[low[k]];
     566    for (int k=0; k<nupp; k++)
     567      H[upp[k]]  = H[upp[k]] + z2[upp[k]]/x2[upp[k]];
     568    w = r2;
     569    for (int k=0; k<nlow; k++)
     570      w[low[k]]  = w[low[k]] - (cL[low[k]] + z1[low[k]]*rL[low[k]])/x1[low[k]];
     571    for (int k=0; k<nupp; k++)
     572      w[upp[k]]  = w[upp[k]] + (cU[upp[k]] + z2[upp[k]]*rU[upp[k]])/x2[upp[k]];
     573
     574    if (LSproblem == 1){
     575      //-----------------------------------------------------------------
     576      //  Solve (*) for dy.
     577      //-----------------------------------------------------------------
     578      H      = 1.0/H;    // H is now Hinv (NOTE!)
     579      for (int k=0; k<nfix; k++) 
     580        H[fix[k]] = 0;
     581      for (int k=0; k<n; k++)
     582        D_elts[k]= sqrt(H_elts[k]);
     583
     584      lsqr->diag1 = D_elts;
     585      lsqr->diag2 = d2;
     586
     587      if (direct){
     588        // Omit direct option for now
     589      }else {// Iterative solve using LSQR.
     590        //rhs     = [ D.*w; r1./d2 ];
     591        for (int k=0; k<n; k++)
     592          rhs[k] = D_elts[k]*w_elts[k];
     593        for (int k=0; k<m; k++)
     594          rhs[n+k] = r1_elts[k]*(1.0/d2);
     595        real damp    = 0;
     596       
     597        if (precon){    // Construct diagonal preconditioner for LSQR
     598          model->matPrecon(d2, Pr, D);
     599        } 
     600            /*
     601        rw(7)        = precon;
     602        info.atolmin = atolmin;
     603        info.r3norm  = fmerit;  // Must be the 2-norm here.
     604       
     605        [ dy, istop, itncg, outfo ] = ...
     606            pdxxxlsqr( nb,m,"pdxxxlsqrmat",Aname,rw,rhs,damp, ...
     607                       atol,btol,conlim,itnlim,show,info );
     608       
     609
     610        lsqr->input->rhs_vec = &rhs;
     611        lsqr->input->sol_vec = &dy;
     612        lsqr->input->rel_mat_err = atol;       
     613        lsqr->do_lsqr(model);
     614        */
     615        //  New version of lsqr
     616
     617        int istop, itn;
     618        dy.clear();
     619        show = false;
     620        info.atolmin = atolmin;
     621        info.r3norm  = fmerit;  // Must be the 2-norm here.
     622
     623        lsqr->do_lsqr(model, rhs, damp, atol, btol, conlim, itnlim,
     624                      show, info, dy , &istop, &itncg, &outfo, precon, Pr);
     625        if (precon)
     626                  dy = dy*Pr;
     627
     628        if (!precon && itncg>999999)
     629          precon=true;
     630
     631        if (istop == 3  ||  istop == 7 )  // conlim or itnlim
     632          printf("\n    LSQR stopped early:  istop = //3d", istop);
     633       
     634
     635        atolold   = outfo.atolold;
     636        atol      = outfo.atolnew;
     637        r3ratio   = outfo.r3ratio;
     638      }// LSproblem 1
     639
     640      //      grad      = pdxxxmat( Aname,2,m,n,dy );   // grad = A"dy
     641      grad.clear();
     642      model->matVecMult(2, grad, dy);
     643      for (int k=0; k<nfix; k++)
     644        grad[fix[k]] = 0;                            // grad is a work vector
     645      dx = H * (grad - w);
     646   
     647    }else{
     648      perror( "This LSproblem not yet implemented\n" );
     649    }
     650    //-------------------------------------------------------------------
     651
     652    CGitns += itncg;
     653
     654    //-------------------------------------------------------------------
     655    // dx and dy are now known.  Get dx1, dx2, dz1, dz2.
     656    //-------------------------------------------------------------------
     657    for (int k=0; k<nlow; k++){   
     658      dx1[low[k]] = - rL[low[k]] + dx[low[k]];
     659      dz1[low[k]] =  (cL[low[k]] - z1[low[k]]*dx1[low[k]]) / x1[low[k]];
     660    }
     661    for (int k=0; k<nupp; k++){   
     662      dx2[upp[k]] = - rU[upp[k]] - dx[upp[k]];
     663      dz2[upp[k]] =  (cU[upp[k]] - z2[upp[k]]*dx2[upp[k]]) / x2[upp[k]];
     664    }
     665    //-------------------------------------------------------------------
     666    // Find the maximum step.
     667    //--------------------------------------------------------------------
     668    real stepx1 = pdxxxstep(nlow, low, x1, dx1 );
     669    real stepx2 = inf;
     670    if (nupp > 0)
     671      stepx2 = pdxxxstep(nupp, upp, x2, dx2 );
     672    real stepz1 = pdxxxstep( z1     , dz1      );
     673    real stepz2 = inf;
     674    if (nupp > 0)
     675      stepz2 = pdxxxstep( z2     , dz2      );
     676    real stepx  = min( stepx1, stepx2 );
     677    real stepz  = min( stepz1, stepz2 );
     678    stepx  = min( steptol*stepx, 1 );
     679    stepz  = min( steptol*stepz, 1 );
     680    if (stepSame){                   // For NLPs, force same step
     681      stepx = min( stepx, stepz );   // (true Newton method)
     682      stepz = stepx;
     683    }
     684   
     685    //-------------------------------------------------------------------
     686    // Backtracking linesearch.
     687    //-------------------------------------------------------------------
     688    bool fail     =  true;
     689    nf       =  0;
     690
     691    while (nf < maxf){
     692      nf      = nf + 1;
     693      x       = x        +  stepx * dx;
     694      y       = y        +  stepz * dy;
     695      for (int k=0; k<nlow; k++){   
     696        x1[low[k]] = x1[low[k]]  +  stepx * dx1[low[k]];
     697        z1[low[k]] = z1[low[k]]  +  stepz * dz1[low[k]];
    235698      }
    236       // now create my duals
    237       for (i=0;i<numberRows_;i++) {
    238         // slack
    239         double value = dual_[i];
    240         value += rowObjectiveWork_[i];
    241         saveDuals[i+numberColumns_]=value;
     699      for (int k=0; k<nupp; k++){   
     700        x2[upp[k]] = x2[upp[k]]  +  stepx * dx2[upp[k]];
     701        z2[upp[k]] = z2[upp[k]]  +  stepz * dz2[upp[k]];
    242702      }
    243       ClpDisjointCopyN(objectiveWork_,numberColumns_,saveDuals);
    244       transposeTimes(-1.0,dual_,saveDuals);
    245       memcpy(dj_,saveDuals,(numberColumns_+numberRows_)*sizeof(double));
    246       // set up possible ones
    247       for (i=0;i<numberRows_+numberColumns_;i++)
    248         clearPivoted(i);
    249       int iRow;
    250       for (iRow=0;iRow<numberRows_;iRow++) {
    251         int iPivot=pivotVariable_[iRow];
    252         if (fabs(saveDuals[iPivot])>dualTolerance_) {
    253           if (getStatus(iPivot)!=isFree)
    254             setPivoted(iPivot);
    255         }
     703      //      [obj,grad,hess] = feval( Fname, (x*beta) );
     704      x.scale(beta);
     705      obj = model->getObj(x);
     706      model->getGrad(x, grad);
     707      model->getHessian(x, H);
     708      x.scale((1.0/beta));
     709     
     710      obj        /= theta;
     711      grad       = grad*(beta /theta)  +  d1*d1*x;
     712      H          = H*(beta2/theta)  +  d1*d1;
     713     
     714      //      [r1,r2,rL,rU,Pinf,Dinf] = ...
     715      pdxxxresid1( model, nlow, nupp, nfix, low, upp, fix,
     716                   b,bl_elts,bu_elts,d1,d2,grad,rL,rU,x,x1,x2,
     717                   y,z1,z2, r1, r2, &Pinf, &Dinf );
     718      //real center, Cinf, Cinf0;
     719      //      [cL,cU,center,Cinf,Cinf0] = ...
     720                  pdxxxresid2( mu, nlow, nupp, low,upp,cL,cU,x1,x2,z1,z2,
     721                               &center, &Cinf, &Cinf0);
     722      real fmeritnew = pdxxxmerit(nlow, nupp, low,upp,r1,r2,rL,rU,cL,cU );
     723      real step      = min( stepx, stepz );
     724
     725      if (fmeritnew <= (1 - eta*step)*fmerit){
     726         fail = false;
     727         break;
    256728      }
    257     }
    258 
    259     double objectiveChange;
    260     numberFake_ =0; // Number of variables at fake bounds
    261     changeBounds(true,NULL,objectiveChange);
    262    
    263     int lastCleaned=0; // last time objective or bounds cleaned up
    264 
    265     if (!ifValuesPass) {
    266       // Check optimal
    267       if (!numberDualInfeasibilities_&&!numberPrimalInfeasibilities_)
    268         problemStatus_=0;
    269     }
    270     if (problemStatus_<0&&perturbation_<100) {
    271       perturb();
    272       // Can't get here if values pass
    273       gutsOfSolution(NULL,NULL);
    274       if (handler_->logLevel()>2) {
    275         handler_->message(CLP_SIMPLEX_STATUS,messages_)
    276           <<numberIterations_<<objectiveValue();
    277         handler_->printing(sumPrimalInfeasibilities_>0.0)
    278           <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
    279         handler_->printing(sumDualInfeasibilities_>0.0)
    280           <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
    281         handler_->printing(numberDualInfeasibilitiesWithoutFree_
    282                            <numberDualInfeasibilities_)
    283                              <<numberDualInfeasibilitiesWithoutFree_;
    284         handler_->message()<<CoinMessageEol;
     729
     730      // Merit function didn"t decrease.
     731      // Restore variables to previous values.
     732      // (This introduces a little error, but save lots of space.)
     733     
     734      x       = x        -  stepx * dx;
     735      y       = y        -  stepz * dy;
     736      for (int k=0; k<nlow; k++){   
     737        x1[low[k]] = x1[low[k]]  -  stepx * dx1[low[k]];
     738        z1[low[k]] = z1[low[k]]  -  stepz * dz1[low[k]];
    285739      }
    286     }
    287        
    288     // This says whether to restore things etc
    289     int factorType=0;
    290     /*
    291       Status of problem:
    292       0 - optimal
    293       1 - infeasible
    294       2 - unbounded
    295       -1 - iterating
    296       -2 - factorization wanted
    297       -3 - redo checking without factorization
    298       -4 - looks infeasible
    299     */
    300     while (problemStatus_<0) {
    301       int iRow, iColumn;
    302       // clear
    303       for (iRow=0;iRow<4;iRow++) {
    304         rowArray_[iRow]->clear();
    305       }   
    306      
    307       for (iColumn=0;iColumn<2;iColumn++) {
    308         columnArray_[iColumn]->clear();
    309       }   
    310      
    311       // give matrix (and model costs and bounds a chance to be
    312       // refreshed (normally null)
    313       matrix_->refresh(this);
    314       // If getting nowhere - why not give it a kick
    315       // does not seem to work too well - do some more work
    316       if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)
    317           &&initialStatus!=10) {
    318         perturb();
    319         // Can't get here if values pass
    320         gutsOfSolution(NULL,NULL);
    321         if (handler_->logLevel()>2) {
    322           handler_->message(CLP_SIMPLEX_STATUS,messages_)
    323             <<numberIterations_<<objectiveValue();
    324           handler_->printing(sumPrimalInfeasibilities_>0.0)
    325             <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
    326           handler_->printing(sumDualInfeasibilities_>0.0)
    327             <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
    328           handler_->printing(numberDualInfeasibilitiesWithoutFree_
    329                              <numberDualInfeasibilities_)
    330                                <<numberDualInfeasibilitiesWithoutFree_;
    331           handler_->message()<<CoinMessageEol;
    332         }
     740      for (int k=0; k<nupp; k++){   
     741        x2[upp[k]] = x2[upp[k]]  -  stepx * dx2[upp[k]];
     742        z2[upp[k]] = z2[upp[k]]  -  stepz * dz2[upp[k]];
    333743      }
    334       // may factorize, checks if problem finished
    335       statusOfProblemInDual(lastCleaned,factorType,progress_,saveDuals);
    336       // If values pass then do easy ones on first time
    337       if (ifValuesPass&&
    338           progress_->lastIterationNumber(0)<0) {
    339         doEasyOnesInValuesPass(saveDuals);
     744      // Back-track.
     745      // If it"s the first time,
     746      // make stepx and stepz the same.
     747
     748      if (nf == 1 && stepx != stepz){
     749         stepx = step;
     750      }else if (nf < maxf){
     751         stepx = stepx/2;
    340752      }
    341      
    342       // Say good factorization
    343       factorType=1;
    344       if (data.sparseThreshold_) {
    345         // use default at present
    346         factorization_->sparseThreshold(0);
    347         factorization_->goSparse();
    348       }
    349 
    350       // exit if victory declared
    351       if (problemStatus_>=0)
    352         break;
    353      
    354       // Do iterations
    355       whileIterating(saveDuals);
    356     }
     753      stepz = stepx;
     754    }
     755
     756    if (fail){
     757      printf("\n     Linesearch failed (nf too big)");
     758      nfail += 1;
     759    }else{
     760      nfail = 0;
     761    }
     762
     763    //-------------------------------------------------------------------
     764    // Set convergence measures.
     765    //--------------------------------------------------------------------
     766    regx = (d1*x).twoNorm();
     767    regy = (d2*y).twoNorm();
     768    regterm = regx*regx + regy*regy;
     769    objreg  = obj  +  0.5*regterm;
     770    objtrue = objreg * theta;
     771
     772    bool primalfeas    = Pinf  <=  featol;
     773    bool dualfeas      = Dinf  <=  featol;
     774    bool complementary = Cinf0 <=  opttol;
     775    bool enough        = PDitns>=       4;  // Prevent premature termination.
     776    bool converged     = primalfeas  &  dualfeas  &  complementary  &  enough;
     777
     778    //-------------------------------------------------------------------
     779    // Iteration log.
     780    //-------------------------------------------------------------------
     781    char str1[100],str2[100],str3[100],str4[100],str5[100];
     782    sprintf(str1, "\n%3g%5.1f" , PDitns      , log10(mu)   );
     783    sprintf(str2, "%6.3f%6.3f" , stepx       , stepz       );
     784    if (stepx < 0.001 || stepz < 0.001){
     785      sprintf(str2, "%6.1e%6.1e" , stepx       , stepz       );
     786    }
     787
     788    sprintf(str3, "%6.1f%6.1f" , log10(Pinf) , log10(Dinf));
     789    sprintf(str4, "%6.1f%15.7e", log10(Cinf0), objtrue     );
     790    sprintf(str5, "%3d%8.1f"   , nf          , center      );
     791    if (center > 99999){
     792      sprintf(str5, "%3d%8.1e"   , nf          , center      );
     793    }
     794    printf("%s%s%s%s%s", str1, str2, str3, str4, str5);
     795    if (direct){
     796      // relax
     797    }else{
     798      printf(" %5.1f%7d%7.3f", log10(atolold), itncg, r3ratio);
     799    }
     800
     801    //-------------------------------------------------------------------
     802    // Test for termination.
     803    //-------------------------------------------------------------------
     804    if (kminor){
     805      printf( "\nStart of next minor itn...\n");
     806      //      keyboard;
     807    }
     808
     809    if (converged){
     810      printf("\n   Converged");
     811      break;
     812    }else if (PDitns >= maxitn){
     813      printf("\n   Too many iterations");
     814      inform = 1;
     815      break;
     816    }else if (nfail  >= maxfail){
     817      printf("\n   Too many linesearch failures");
     818      inform = 2;
     819      break;
     820    }else{
     821
     822      // Reduce mu, and reset certain residuals.
     823
     824      real stepmu  = min( stepx , stepz   );
     825      stepmu  = min( stepmu, steptol );
     826      real muold   = mu;
     827      mu      = mu   -  stepmu * mu;
     828      if (center >= bigcenter)
     829        mu = muold; 
     830
     831      // mutrad = mu0*(sum(Xz)/n); // 24 May 1998: Traditional value, but
     832      // mu     = min(mu,mutrad ); // it seemed to decrease mu too much.
     833
     834      mu      = max(mu,mulast);  // 13 Jun 1998: No need for smaller mu.
     835      //      [cL,cU,center,Cinf,Cinf0] = ...
     836      pdxxxresid2( mu, nlow, nupp, low, upp, cL, cU, x1, x2, z1, z2,
     837                   &center, &Cinf, &Cinf0 );
     838      fmerit = pdxxxmerit( nlow, nupp, low,upp,r1,r2,rL,rU,cL,cU );
     839
     840      // Reduce atol for LSQR (and SYMMLQ).
     841      // NOW DONE AT TOP OF LOOP.
     842
     843      atolold = atol;
     844      // if atol > atol2
     845      //   atolfac = (mu/mufirst)^0.25;
     846      //   atol    = max( atol*atolfac, atol2 );
     847      // end
     848
     849      // atol = min( atol, mu );     // 22 Jan 2001: a la Inexact Newton.
     850      // atol = min( atol, 0.5*mu ); // 30 Jan 2001: A bit tighter
     851
     852      // If the linesearch took more than one function (nf > 1),
     853      // we assume the search direction needed more accuracy
     854      // (though this may be true only for LPs).
     855      // 12 Jun 1998: Ask for more accuracy if nf > 2.
     856      // 24 Nov 2000: Also if the steps are small.
     857      // 30 Jan 2001: Small steps might be ok with warm start.
     858      // 06 Feb 2001: Not necessarily.  Reinstated tests in next line.
     859
     860      if (nf > 2  ||  min( stepx, stepz ) <= 0.01)
     861        atol = atolold*0.1;
     862    }
     863  //---------------------------------------------------------------------
     864  // End of main loop.
     865  //---------------------------------------------------------------------
    357866  }
    358867
    359   assert (problemStatus_||!sumPrimalInfeasibilities_);
    360 
    361   // clean up
    362   finish();
    363   delete [] saveDuals;
    364 
    365   // Restore any saved stuff
    366   restoreData(data);
    367   return problemStatus_;
    368 }
    369 /* Reasons to come out:
    370    -1 iterations etc
    371    -2 inaccuracy
    372    -3 slight inaccuracy (and done iterations)
    373    +0 looks optimal (might be unbounded - but we will investigate)
    374    +1 looks infeasible
    375    +3 max iterations
    376  */
    377 int
    378 ClpPdco::whileIterating(double * & givenDuals)
    379 {
    380 #ifdef CLP_DEBUG
    381   int debugIteration=-1;
    382 #endif
    383   {
    384     int i;
    385     for (i=0;i<4;i++) {
    386       rowArray_[i]->clear();
    387     }   
    388     for (i=0;i<2;i++) {
    389       columnArray_[i]->clear();
    390     }   
    391   }     
    392   // if can't trust much and long way from optimal then relax
    393   if (largestPrimalError_>10.0)
    394     factorization_->relaxAccuracyCheck(min(1.0e2,largestPrimalError_/10.0));
    395   else
    396     factorization_->relaxAccuracyCheck(1.0);
    397   // status stays at -1 while iterating, >=0 finished, -2 to invert
    398   // status -3 to go to top without an invert
    399   int returnCode = -1;
    400 
    401   // compute average infeasibility for backward test
    402   double averagePrimalInfeasibility = sumPrimalInfeasibilities_/
    403     ((double ) (numberPrimalInfeasibilities_+1));
    404 
    405   // Get dubious weights
    406   factorization_->getWeights(rowArray_[0]->getIndices());
    407   CoinBigIndex * dubiousWeights = matrix_->dubiousWeights(this,rowArray_[0]->getIndices());
    408   // If values pass then get list of candidates
    409   int * candidateList = NULL;
    410   int numberCandidates = 0;
    411 #ifdef CLP_DEBUG
    412   bool wasInValuesPass= (givenDuals!=NULL);
    413 #endif
    414   int candidate=-1;
    415   if (givenDuals) {
    416     candidateList = new int[numberRows_];
    417     // move reduced costs across
    418     memcpy(dj_,givenDuals,(numberRows_+numberColumns_)*sizeof(double));
    419     int iRow;
    420     for (iRow=0;iRow<numberRows_;iRow++) {
    421       int iPivot=pivotVariable_[iRow];
    422       if (flagged(iPivot))
    423         continue;
    424       if (fabs(dj_[iPivot])>dualTolerance_) {
    425         if (pivoted(iPivot))
    426           candidateList[numberCandidates++]=iRow;
    427       } else {
    428         clearPivoted(iPivot);
    429       }
    430     }
    431     // and set first candidate
    432     if (!numberCandidates) {
    433       delete [] candidateList;
    434       delete [] givenDuals;
    435       givenDuals=NULL;
    436       candidateList=NULL;
    437       int iRow;
    438       for (iRow=0;iRow<numberRows_;iRow++) {
    439         int iPivot=pivotVariable_[iRow];
    440         clearPivoted(iPivot);
    441       }
    442     }
    443   }
    444   while (problemStatus_==-1) {
    445 #ifdef CLP_DEBUG
    446     if (givenDuals) {
    447       double value5=0.0;
    448       int i;
    449       for (i=0;i<numberRows_+numberColumns_;i++) {
    450         if (dj_[i]<-1.0e-6)
    451           if (upper_[i]<1.0e20)
    452             value5 += dj_[i]*upper_[i];
    453           else
    454             printf("bad dj %g on %d with large upper status %d\n",
    455                    dj_[i],i,status_[i]&7);
    456         else if (dj_[i] >1.0e-6)
    457           if (lower_[i]>-1.0e20)
    458             value5 += dj_[i]*lower_[i];
    459           else
    460             printf("bad dj %g on %d with large lower status %d\n",
    461                    dj_[i],i,status_[i]&7);
    462       }
    463       printf("Values objective Value %g\n",value5);
    464     }
    465     if ((handler_->logLevel()&32)&&wasInValuesPass) {
    466       double value5=0.0;
    467       int i;
    468       for (i=0;i<numberRows_+numberColumns_;i++) {
    469         if (dj_[i]<-1.0e-6)
    470           if (upper_[i]<1.0e20)
    471             value5 += dj_[i]*upper_[i];
    472         else if (dj_[i] >1.0e-6)
    473           if (lower_[i]>-1.0e20)
    474             value5 += dj_[i]*lower_[i];
    475       }
    476       printf("Values objective Value %g\n",value5);
    477       {
    478         int i;
    479         for (i=0;i<numberRows_+numberColumns_;i++) {
    480           int iSequence = i;
    481           double oldValue;
    482          
    483           switch(getStatus(iSequence)) {
    484            
    485           case basic:
    486           case ClpSimplex::isFixed:
    487             break;
    488           case isFree:
    489           case superBasic:
    490             abort();
    491             break;
    492           case atUpperBound:
    493             oldValue = dj_[iSequence];
    494             //assert (oldValue<=tolerance);
    495             assert (fabs(solution_[iSequence]-upper_[iSequence])<1.0e-7);
    496             break;
    497           case atLowerBound:
    498             oldValue = dj_[iSequence];
    499             //assert (oldValue>=-tolerance);
    500             assert (fabs(solution_[iSequence]-lower_[iSequence])<1.0e-7);
    501             break;
    502           }
    503         }
    504       }
    505     }
    506 #endif
    507 #ifdef CLP_DEBUG
    508     {
    509       int i;
    510       for (i=0;i<4;i++) {
    511         rowArray_[i]->checkClear();
    512       }   
    513       for (i=0;i<2;i++) {
    514         columnArray_[i]->checkClear();
    515       }   
    516     }     
    517 #endif
    518 #if CLP_DEBUG>2
    519     // very expensive
    520     if (numberIterations_>0&&numberIterations_<-801) {
    521       handler_->setLogLevel(63);
    522       double saveValue = objectiveValue_;
    523       double * saveRow1 = new double[numberRows_];
    524       double * saveRow2 = new double[numberRows_];
    525       memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
    526       memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
    527       double * saveColumn1 = new double[numberColumns_];
    528       double * saveColumn2 = new double[numberColumns_];
    529       memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
    530       memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
    531       gutsOfSolution(NULL,NULL);
    532       printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
    533              numberIterations_,
    534              saveValue,objectiveValue_,sumDualInfeasibilities_);
    535       if (saveValue>objectiveValue_+1.0e-2)
    536         printf("**bad**\n");
    537       memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
    538       memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
    539       memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
    540       memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
    541       delete [] saveRow1;
    542       delete [] saveRow2;
    543       delete [] saveColumn1;
    544       delete [] saveColumn2;
    545       objectiveValue_=saveValue;
    546     }
    547 #endif
    548 #if 0
    549     if (!factorization_->pivots()){
    550       int iPivot;
    551       double * array = rowArray_[3]->denseVector();
    552       int i;
    553       for (iPivot=0;iPivot<numberRows_;iPivot++) {
    554         int iSequence = pivotVariable_[iPivot];
    555         unpack(rowArray_[3],iSequence);
    556         factorization_->updateColumn(rowArray_[2],rowArray_[3]);
    557         assert (fabs(array[iPivot]-1.0)<1.0e-4);
    558         array[iPivot]=0.0;
    559         for (i=0;i<numberRows_;i++)
    560           assert (fabs(array[i])<1.0e-4);
    561         rowArray_[3]->clear();
    562       }
    563     }
    564 #endif
    565 #ifdef CLP_DEBUG
    566     {
    567       int iSequence, number=numberRows_+numberColumns_;
    568       for (iSequence=0;iSequence<number;iSequence++) {
    569         double lowerValue=lower_[iSequence];
    570         double upperValue=upper_[iSequence];
    571         double value=solution_[iSequence];
    572         if(getStatus(iSequence)!=basic) {
    573           assert(lowerValue>-1.0e20);
    574           assert(upperValue<1.0e20);
    575         }
    576         switch(getStatus(iSequence)) {
    577          
    578         case basic:
    579           break;
    580         case isFree:
    581         case superBasic:
    582           break;
    583         case atUpperBound:
    584           assert (fabs(value-upperValue)<=primalTolerance_) ;
    585           break;
    586         case atLowerBound:
    587         case ClpSimplex::isFixed:
    588           assert (fabs(value-lowerValue)<=primalTolerance_) ;
    589           break;
    590         }
    591       }
    592     }
    593     if(numberIterations_==debugIteration) {
    594       printf("dodgy iteration coming up\n");
    595     }
    596 #endif
    597     // choose row to go out
    598     // dualRow will go to virtual row pivot choice algorithm
    599     // make sure values pass off if it should be
    600     if (numberCandidates)
    601       candidate = candidateList[--numberCandidates];
    602     else
    603       candidate=-1;
    604     dualRow(candidate);
    605     if (pivotRow_>=0) {
    606       // we found a pivot row
    607       handler_->message(CLP_SIMPLEX_PIVOTROW,messages_)
    608         <<pivotRow_
    609         <<CoinMessageEol;
    610       // check accuracy of weights
    611       dualRowPivot_->checkAccuracy();
    612       // Get good size for pivot
    613       double acceptablePivot=1.0e-7;
    614       if (factorization_->pivots()>10)
    615         acceptablePivot=1.0e-5; // if we have iterated be more strict
    616       else if (factorization_->pivots()>5)
    617         acceptablePivot=1.0e-6; // if we have iterated be slightly more strict
    618       // get sign for finding row of tableau
    619       if (candidate<0) {
    620         // normal iteration
    621         // create as packed
    622         double direction=directionOut_;
    623         rowArray_[0]->createPacked(1,&pivotRow_,&direction);
    624         factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
    625         // put row of tableau in rowArray[0] and columnArray[0]
    626         matrix_->transposeTimes(this,-1.0,
    627                               rowArray_[0],rowArray_[3],columnArray_[0]);
    628         // do ratio test for normal iteration
    629         dualColumn(rowArray_[0],columnArray_[0],columnArray_[1],
    630                  rowArray_[3],acceptablePivot,dubiousWeights);
    631       } else {
    632         double direction=directionOut_;
    633         rowArray_[0]->createPacked(1,&pivotRow_,&direction);
    634         factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
    635         // put row of tableau in rowArray[0] and columnArray[0]
    636         matrix_->transposeTimes(this,-1.0,
    637                               rowArray_[0],rowArray_[3],columnArray_[0]);
    638         acceptablePivot *= 10.0;
    639         // do ratio test
    640         checkPossibleValuesMove(rowArray_[0],columnArray_[0],
    641                                             acceptablePivot,NULL);
    642 
    643         // recompute true dualOut_
    644         if (directionOut_<0) {
    645           dualOut_ = valueOut_ - upperOut_;
    646         } else {
    647           dualOut_ = lowerOut_ - valueOut_;
    648         }
    649         // check what happened if was values pass
    650         // may want to move part way i.e. movement
    651         bool normalIteration = (sequenceIn_!=sequenceOut_);
    652 
    653         clearPivoted(sequenceOut_);  // make sure won't be done again
    654         // see if end of values pass
    655         if (!numberCandidates) {
    656           int iRow;
    657           delete [] candidateList;
    658           delete [] givenDuals;
    659           candidate=-2; // -2 signals end
    660           givenDuals=NULL;
    661           handler_->message(CLP_END_VALUES_PASS,messages_)
    662             <<numberIterations_;
    663           candidateList=NULL;
    664           for (iRow=0;iRow<numberRows_;iRow++) {
    665             int iPivot=pivotVariable_[iRow];
    666             //assert (fabs(dj_[iPivot]),1.0e-5);
    667             clearPivoted(iPivot);
    668           }
    669         }
    670         if (!normalIteration) {
    671           //rowArray_[0]->cleanAndPackSafe(1.0e-60);
    672           //columnArray_[0]->cleanAndPackSafe(1.0e-60);
    673           updateDualsInValuesPass(rowArray_[0],columnArray_[0],theta_);
    674           if (candidate==-2)
    675             problemStatus_=-2;
    676           continue; // skip rest of iteration
    677         } else {
    678           // recompute dualOut_
    679           if (directionOut_<0) {
    680             dualOut_ = valueOut_ - upperOut_;
    681           } else {
    682             dualOut_ = lowerOut_ - valueOut_;
    683           }
    684         }
    685       }
    686       if (sequenceIn_>=0) {
    687         // normal iteration
    688         // update the incoming column
    689         double btranAlpha = -alpha_*directionOut_; // for check
    690         unpackPacked(rowArray_[1]);
    691         factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
    692         // and update dual weights (can do in parallel - with extra array)
    693         alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
    694                                               rowArray_[2],
    695                                               rowArray_[1]);
    696         // see if update stable
    697 #ifdef CLP_DEBUG
    698         if ((handler_->logLevel()&32))
    699           printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
    700 #endif
    701         double checkValue=1.0e-7;
    702         // if can't trust much and long way from optimal then relax
    703         if (largestPrimalError_>10.0)
    704           checkValue = min(1.0e-4,1.0e-8*largestPrimalError_);
    705         if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
    706             fabs(btranAlpha-alpha_)>checkValue*(1.0+fabs(alpha_))) {
    707           handler_->message(CLP_DUAL_CHECK,messages_)
    708             <<btranAlpha
    709             <<alpha_
    710             <<CoinMessageEol;
    711           if (factorization_->pivots()) {
    712             dualRowPivot_->unrollWeights();
    713             problemStatus_=-2; // factorize now
    714             rowArray_[0]->clear();
    715             rowArray_[1]->clear();
    716             columnArray_[0]->clear();
    717             returnCode=-2;
    718             break;
    719           } else {
    720             // take on more relaxed criterion
    721             if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
    722                 fabs(btranAlpha-alpha_)>1.0e-4*(1.0+fabs(alpha_))) {
    723               dualRowPivot_->unrollWeights();
    724               // need to reject something
    725               char x = isColumn(sequenceOut_) ? 'C' :'R';
    726               handler_->message(CLP_SIMPLEX_FLAG,messages_)
    727                 <<x<<sequenceWithin(sequenceOut_)
    728                 <<CoinMessageEol;
    729               setFlagged(sequenceOut_);
    730               lastBadIteration_ = numberIterations_; // say be more cautious
    731               rowArray_[0]->clear();
    732               rowArray_[1]->clear();
    733               columnArray_[0]->clear();
    734               continue;
    735             }
    736           }
    737         }
    738         // update duals BEFORE replaceColumn so can do updateColumn
    739         double objectiveChange=0.0;
    740         // do duals first as variables may flip bounds
    741         // rowArray_[0] and columnArray_[0] may have flips
    742         // so use rowArray_[3] for work array from here on
    743         int nswapped = 0;
    744         //rowArray_[0]->cleanAndPackSafe(1.0e-60);
    745         //columnArray_[0]->cleanAndPackSafe(1.0e-60);
    746         if (candidate==-1)
    747           nswapped = updateDualsInDual(rowArray_[0],columnArray_[0],
    748                                        rowArray_[2],theta_,
    749                                        objectiveChange,false);
    750         else
    751           updateDualsInValuesPass(rowArray_[0],columnArray_[0],theta_);
    752 
    753         // which will change basic solution
    754         if (nswapped) {
    755           factorization_->updateColumn(rowArray_[3],rowArray_[2]);
    756           dualRowPivot_->updatePrimalSolution(rowArray_[2],
    757                                               1.0,objectiveChange);
    758          
    759           // recompute dualOut_
    760           valueOut_ = solution_[sequenceOut_];
    761           if (directionOut_<0) {
    762             dualOut_ = valueOut_ - upperOut_;
    763           } else {
    764             dualOut_ = lowerOut_ - valueOut_;
    765           }
    766           if(dualOut_<-max(1.0e-12*averagePrimalInfeasibility,1.0e-8)
    767              &&factorization_->pivots()&&
    768              getStatus(sequenceIn_)!=isFree) {
    769             // going backwards - factorize
    770             dualRowPivot_->unrollWeights();
    771             problemStatus_=-2; // factorize now
    772             returnCode=-2;
    773             break;
    774           }
    775         }
    776         // amount primal will move
    777         double movement = -dualOut_*directionOut_/alpha_;
    778         // so objective should increase by fabs(dj)*movement
    779         // but we already have objective change - so check will be good
    780         if (objectiveChange+fabs(movement*dualIn_)<-1.0e-5) {
    781 #ifdef CLP_DEBUG
    782           if (handler_->logLevel()&32)
    783             printf("movement %g, swap change %g, rest %g  * %g\n",
    784                    objectiveChange+fabs(movement*dualIn_),
    785                    objectiveChange,movement,dualIn_);
    786 #endif
    787           if(factorization_->pivots()>5) {
    788             // going backwards - factorize
    789             dualRowPivot_->unrollWeights();
    790             problemStatus_=-2; // factorize now
    791             returnCode=-2;
    792             break;
    793           }
    794         }
    795         assert(fabs(dualOut_)<1.0e50);
    796         // if stable replace in basis
    797         int updateStatus = factorization_->replaceColumn(rowArray_[2],
    798                                                          pivotRow_,
    799                                                          alpha_);
    800         // if no pivots, bad update but reasonable alpha - take and invert
    801         if (updateStatus==2&&
    802                    !factorization_->pivots()&&fabs(alpha_)>1.0e-5)
    803           updateStatus=4;
    804         if (updateStatus==1||updateStatus==4) {
    805           // slight error
    806           if (factorization_->pivots()>5||updateStatus==4) {
    807             problemStatus_=-2; // factorize now
    808             returnCode=-3;
    809           }
    810         } else if (updateStatus==2) {
    811           // major error
    812           dualRowPivot_->unrollWeights();
    813           // later we may need to unwind more e.g. fake bounds
    814           if (factorization_->pivots()) {
    815             problemStatus_=-2; // factorize now
    816             returnCode=-2;
    817             break;
    818           } else {
    819             // need to reject something
    820             char x = isColumn(sequenceOut_) ? 'C' :'R';
    821             handler_->message(CLP_SIMPLEX_FLAG,messages_)
    822               <<x<<sequenceWithin(sequenceOut_)
    823               <<CoinMessageEol;
    824             setFlagged(sequenceOut_);
    825             lastBadIteration_ = numberIterations_; // say be more cautious
    826             rowArray_[0]->clear();
    827             rowArray_[1]->clear();
    828             columnArray_[0]->clear();
    829             // make sure dual feasible
    830             // look at all rows and columns
    831             double objectiveChange=0.0;
    832             updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
    833                               0.0,objectiveChange,true);
    834             continue;
    835           }
    836         } else if (updateStatus==3) {
    837           // out of memory
    838           // increase space if not many iterations
    839           if (factorization_->pivots()<
    840               0.5*factorization_->maximumPivots()&&
    841               factorization_->pivots()<200)
    842             factorization_->areaFactor(
    843                                        factorization_->areaFactor() * 1.1);
    844           problemStatus_=-2; // factorize now
    845         } else if (updateStatus==5) {
    846           problemStatus_=-2; // factorize now
    847         }
    848         // update primal solution
    849         if (theta_<0.0&&candidate==-1) {
    850 #ifdef CLP_DEBUG
    851           if (handler_->logLevel()&32)
    852             printf("negative theta %g\n",theta_);
    853 #endif
    854           theta_=0.0;
    855         }
    856         // do actual flips
    857         flipBounds(rowArray_[0],columnArray_[0],theta_);
    858         //rowArray_[1]->expand();
    859         dualRowPivot_->updatePrimalSolution(rowArray_[1],
    860                                             movement,
    861                                             objectiveChange);
    862 #ifdef CLP_DEBUG
    863         double oldobj=objectiveValue_;
    864 #endif
    865         // modify dualout
    866         dualOut_ /= alpha_;
    867         dualOut_ *= -directionOut_;
    868         //setStatus(sequenceIn_,basic);
    869         dj_[sequenceIn_]=0.0;
    870         double oldValue=valueIn_;
    871         if (directionIn_==-1) {
    872           // as if from upper bound
    873           valueIn_ = upperIn_+dualOut_;
    874         } else {
    875           // as if from lower bound
    876           valueIn_ = lowerIn_+dualOut_;
    877         }
    878         objectiveChange += cost_[sequenceIn_]*(valueIn_-oldValue);
    879         // outgoing
    880         // set dj to zero unless values pass
    881         if (directionOut_>0) {
    882           valueOut_ = lowerOut_;
    883           if (candidate==-1)
    884             dj_[sequenceOut_] = theta_;
    885         } else {
    886           valueOut_ = upperOut_;
    887           if (candidate==-1)
    888             dj_[sequenceOut_] = -theta_;
    889         }
    890         solution_[sequenceOut_]=valueOut_;
    891         int whatNext=housekeeping(objectiveChange);
    892         // and set bounds correctly
    893         originalBound(sequenceIn_);
    894         changeBound(sequenceOut_);
    895 #ifdef CLP_DEBUG
    896         if (objectiveValue_<oldobj-1.0e-5&&(handler_->logLevel()&16))
    897           printf("obj backwards %g %g\n",objectiveValue_,oldobj);
    898 #endif
    899         if (whatNext==1||candidate==-2) {
    900           problemStatus_ =-2; // refactorize
    901         } else if (whatNext==2) {
    902           // maximum iterations or equivalent
    903           problemStatus_= 3;
    904           returnCode=3;
    905           break;
    906         }
    907       } else {
    908         // no incoming column is valid
    909 #ifdef CLP_DEBUG
    910         if (handler_->logLevel()&32)
    911           printf("** no column pivot\n");
    912 #endif
    913         if (factorization_->pivots()<5) {
    914           problemStatus_=-4; //say looks infeasible
    915           // create ray anyway
    916           delete [] ray_;
    917           ray_ = new double [ numberRows_];
    918           rowArray_[0]->expand(); // in case packed
    919           ClpDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
    920         }
    921         rowArray_[0]->clear();
    922         columnArray_[0]->clear();
    923         returnCode=1;
     868
     869  for (int k=0; k<nfix; k++)   
     870    x[fix[k]] = bl[fix[k]];
     871  z      = z1;
     872  if (nupp > 0)
     873    z = z - z2;
     874  printf("\n\nmax |x| =%10.3f", x.infNorm() );
     875  printf("    max |y| =%10.3f", y.infNorm() );
     876  printf("    max |z| =%10.3f", z.infNorm() );
     877  printf("   scaled");
     878
     879  x.scale(beta);   y.scale(zeta);   z.scale(zeta);   // Unscale x, y, z.
     880
     881  printf(  "\nmax |x| =%10.3f", x.infNorm() );
     882  printf("    max |y| =%10.3f", y.infNorm() );
     883  printf("    max |z| =%10.3f", z.infNorm() );
     884  printf(" unscaled\n");
     885
     886  time   = now() - time;
     887  char str1[100],str2[100];
     888  sprintf(str1, "\nPDitns  =%10g", PDitns );
     889  sprintf(str2, "itns =%10g", CGitns );
     890  //  printf( [str1 " " solver str2] );
     891  printf("    time    =%10.1f\n", time);
     892  /*
     893  pdxxxdistrib( abs(x),abs(z) );   // Private function
     894 
     895  if (wait)
     896    keyboard;
     897  */
     898//-----------------------------------------------------------------------
     899// End function pdco.m
     900//-----------------------------------------------------------------------
     901/*  printf("Solution x values:\n\n");
     902  for (int k=0; k<n; k++)
     903    printf(" %d   %e\n", k, x[k]);
     904*/
     905// Print distribution
     906  float thresh[9]={ 0.00000001, 0.0000001, 0.000001,0.00001, 0.0001, 0.001, 0.01, 0.1, 1.00001};
     907  int counts[9]={0};
     908  for (int ij=0; ij<n; ij++){
     909    for (int j=0; j<9; j++){
     910      if(x[ij] < thresh[j]){
     911        counts[j] += 1;
    924912        break;
    925913      }
    926     } else {
    927       // no pivot row
    928 #ifdef CLP_DEBUG
    929       if (handler_->logLevel()&32)
    930         printf("** no row pivot\n");
    931 #endif
    932       if (!factorization_->pivots()) {
    933         // may have crept through - so may be optimal
    934         // check any flagged variables
    935         int iRow;
    936         for (iRow=0;iRow<numberRows_;iRow++) {
    937           int iPivot=pivotVariable_[iRow];
    938           if (flagged(iPivot))
    939             break;
    940         }
    941         if (numberFake_||numberDualInfeasibilities_) {
    942           // may be dual infeasible
    943           problemStatus_=-5;
    944         } else {
    945           if (iRow<numberRows_) {
    946 #ifdef CLP_DEBUG
    947             std::cerr<<"Flagged variables at end - infeasible?"<<std::endl;
    948             printf("Probably infeasible - pivot was %g\n",alpha_);
    949 #endif
    950             if (fabs(alpha_)<1.0e-4) {
    951               problemStatus_=1;
    952             } else {
    953 #ifdef CLP_DEBUG
    954               abort();
    955 #endif
    956             }
    957           } else {
    958             problemStatus_=0;
    959           }
    960         }
    961       } else {
    962         problemStatus_=-3;
    963       }
    964       returnCode=0;
    965       break;
    966914    }
    967915  }
    968   if (givenDuals) {
    969     memcpy(givenDuals,dj_,(numberRows_+numberColumns_)*sizeof(double));
    970     // get rid of any values pass array
    971     delete [] candidateList;
    972   }
    973   delete [] dubiousWeights;
    974   return returnCode;
     916  printf ("Distribution of Solution Values\n");
     917  for (int j=8; j>1; j--)
     918    printf(" %f  to  %f %d\n",thresh[j-1],thresh[j],counts[j]);
     919  printf("   Less than   %f %d\n",thresh[2],counts[0]);
     920
     921
     922return 0;
    975923}
    976 /* The duals are updated by the given arrays.
    977    Returns number of infeasibilities.
    978    rowArray and columnarray will have flipped
    979    The output vector has movement (row length array) */
    980 int
    981 ClpPdco::updateDualsInDual(CoinIndexedVector * rowArray,
    982                                   CoinIndexedVector * columnArray,
    983                                   CoinIndexedVector * outputArray,
    984                                   double theta,
    985                                   double & objectiveChange,
    986                                   bool fullRecompute)
     924// LSQR
     925void
     926ClpPdco::lsqr()
    987927{
    988  
    989   outputArray->clear();
    990  
    991  
    992   int numberInfeasibilities=0;
    993   int numberRowInfeasibilities=0;
    994  
    995   // use a tighter tolerance except for all being okay
    996   double tolerance = dualTolerance_;
    997  
    998   double changeObj=0.0;
    999 
    1000   // Coding is very similar but we can save a bit by splitting
    1001   // Do rows
    1002   if (!fullRecompute) {
    1003     int i;
    1004     double * reducedCost = djRegion(0);
    1005     const double * lower = lowerRegion(0);
    1006     const double * upper = upperRegion(0);
    1007     const double * cost = costRegion(0);
    1008     double * work;
    1009     int number;
    1010     int * which;
    1011     assert(rowArray->packedMode());
    1012     work = rowArray->denseVector();
    1013     number = rowArray->getNumElements();
    1014     which = rowArray->getIndices();
    1015     for (i=0;i<number;i++) {
    1016       int iSequence = which[i];
    1017       double alphaI = work[i];
    1018       work[i]=0.0;
    1019      
    1020       Status status = getStatus(iSequence+numberColumns_);
    1021       // more likely to be at upper bound ?
    1022       if (status==atUpperBound) {
    1023 
    1024         double value = reducedCost[iSequence]-theta*alphaI;
    1025         reducedCost[iSequence]=value;
    1026 
    1027         if (value>tolerance) {
    1028           // to lower bound (if swap)
    1029           which[numberInfeasibilities++]=iSequence;
    1030           double movement = lower[iSequence]-upper[iSequence];
    1031 #ifdef CLP_DEBUG
    1032           if ((handler_->logLevel()&32))
    1033             printf("%d %d, new dj %g, alpha %g, movement %g\n",
    1034                    0,iSequence,value,alphaI,movement);
    1035 #endif
    1036           changeObj += movement*cost[iSequence];
    1037           outputArray->quickAdd(iSequence,-movement);
    1038         }
    1039       } else if (status==atLowerBound) {
    1040 
    1041         double value = reducedCost[iSequence]-theta*alphaI;
    1042         reducedCost[iSequence]=value;
    1043 
    1044         if (value<-tolerance) {
    1045           // to upper bound
    1046           which[numberInfeasibilities++]=iSequence;
    1047           double movement = upper[iSequence] - lower[iSequence];
    1048 #ifdef CLP_DEBUG
    1049           if ((handler_->logLevel()&32))
    1050             printf("%d %d, new dj %g, alpha %g, movement %g\n",
    1051                    0,iSequence,value,alphaI,movement);
    1052 #endif
    1053           changeObj += movement*cost[iSequence];
    1054           outputArray->quickAdd(iSequence,-movement);
    1055         }
    1056       }
    1057     }
    1058   } else  {
    1059     double * solution = solutionRegion(0);
    1060     double * reducedCost = djRegion(0);
    1061     const double * lower = lowerRegion(0);
    1062     const double * upper = upperRegion(0);
    1063     const double * cost = costRegion(0);
    1064     int * which;
    1065     which = rowArray->getIndices();
    1066     int iSequence;
    1067     for (iSequence=0;iSequence<numberRows_;iSequence++) {
    1068       double value = reducedCost[iSequence];
    1069      
    1070       Status status = getStatus(iSequence+numberColumns_);
    1071       // more likely to be at upper bound ?
    1072       if (status==atUpperBound) {
    1073         double movement=0.0;
    1074 
    1075         if (value>tolerance) {
    1076           // to lower bound (if swap)
    1077           // put back alpha
    1078           which[numberInfeasibilities++]=iSequence;
    1079           movement = lower[iSequence]-upper[iSequence];
    1080           changeObj += movement*cost[iSequence];
    1081           outputArray->quickAdd(iSequence,-movement);
    1082         } else if (value>-tolerance) {
    1083           // at correct bound but may swap
    1084           FakeBound bound = getFakeBound(iSequence+numberColumns_);
    1085           if (bound==ClpPdco::upperFake) {
    1086             movement = lower[iSequence]-upper[iSequence];
    1087             setStatus(iSequence+numberColumns_,atLowerBound);
    1088             solution[iSequence] = lower[iSequence];
    1089             changeObj += movement*cost[iSequence];
    1090             numberFake_--;
    1091             setFakeBound(iSequence+numberColumns_,noFake);
    1092           }
    1093         }
    1094       } else if (status==atLowerBound) {
    1095         double movement=0.0;
    1096 
    1097         if (value<-tolerance) {
    1098           // to upper bound
    1099           // put back alpha
    1100           which[numberInfeasibilities++]=iSequence;
    1101           movement = upper[iSequence] - lower[iSequence];
    1102           changeObj += movement*cost[iSequence];
    1103           outputArray->quickAdd(iSequence,-movement);
    1104         } else if (value<tolerance) {
    1105           // at correct bound but may swap
    1106           FakeBound bound = getFakeBound(iSequence+numberColumns_);
    1107           if (bound==ClpPdco::lowerFake) {
    1108             movement = upper[iSequence]-lower[iSequence];
    1109             setStatus(iSequence+numberColumns_,atUpperBound);
    1110             solution[iSequence] = upper[iSequence];
    1111             changeObj += movement*cost[iSequence];
    1112             numberFake_--;
    1113             setFakeBound(iSequence+numberColumns_,noFake);
    1114           }
    1115         }
    1116       }
    1117     }
    1118   }
    1119 
    1120   // Do columns
    1121   if (!fullRecompute) {
    1122     int i;
    1123     double * reducedCost = djRegion(1);
    1124     const double * lower = lowerRegion(1);
    1125     const double * upper = upperRegion(1);
    1126     const double * cost = costRegion(1);
    1127     double * work;
    1128     int number;
    1129     int * which;
    1130     // set number of infeasibilities in row array
    1131     numberRowInfeasibilities=numberInfeasibilities;
    1132     rowArray->setNumElements(numberInfeasibilities);
    1133     numberInfeasibilities=0;
    1134     work = columnArray->denseVector();
    1135     number = columnArray->getNumElements();
    1136     which = columnArray->getIndices();
    1137    
    1138     for (i=0;i<number;i++) {
    1139       int iSequence = which[i];
    1140       double alphaI = work[i];
    1141       work[i]=0.0;
    1142      
    1143       Status status = getStatus(iSequence);
    1144       if (status==atLowerBound) {
    1145         double value = reducedCost[iSequence]-theta*alphaI;
    1146         reducedCost[iSequence]=value;
    1147         double movement=0.0;
    1148 
    1149         if (value<-tolerance) {
    1150           // to upper bound
    1151           which[numberInfeasibilities++]=iSequence;
    1152           movement = upper[iSequence] - lower[iSequence];
    1153 #ifdef CLP_DEBUG
    1154           if ((handler_->logLevel()&32))
    1155             printf("%d %d, new dj %g, alpha %g, movement %g\n",
    1156                    1,iSequence,value,alphaI,movement);
    1157 #endif
    1158           changeObj += movement*cost[iSequence];
    1159           matrix_->add(this,outputArray,iSequence,movement);
    1160         }
    1161       } else if (status==atUpperBound) {
    1162         double value = reducedCost[iSequence]-theta*alphaI;
    1163         reducedCost[iSequence]=value;
    1164         double movement=0.0;
    1165 
    1166         if (value>tolerance) {
    1167           // to lower bound (if swap)
    1168           which[numberInfeasibilities++]=iSequence;
    1169           movement = lower[iSequence]-upper[iSequence];
    1170 #ifdef CLP_DEBUG
    1171           if ((handler_->logLevel()&32))
    1172             printf("%d %d, new dj %g, alpha %g, movement %g\n",
    1173                    1,iSequence,value,alphaI,movement);
    1174 #endif
    1175           changeObj += movement*cost[iSequence];
    1176           matrix_->add(this,outputArray,iSequence,movement);
    1177         }
    1178       } else if (status==isFree) {
    1179         double value = reducedCost[iSequence]-theta*alphaI;
    1180         reducedCost[iSequence]=value;
    1181       }
    1182     }
    1183   } else {
    1184     double * solution = solutionRegion(1);
    1185     double * reducedCost = djRegion(1);
    1186     const double * lower = lowerRegion(1);
    1187     const double * upper = upperRegion(1);
    1188     const double * cost = costRegion(1);
    1189     int * which;
    1190     // set number of infeasibilities in row array
    1191     numberRowInfeasibilities=numberInfeasibilities;
    1192     rowArray->setNumElements(numberInfeasibilities);
    1193     numberInfeasibilities=0;
    1194     which = columnArray->getIndices();
    1195     int iSequence;
    1196     for (iSequence=0;iSequence<numberColumns_;iSequence++) {
    1197       double value = reducedCost[iSequence];
    1198      
    1199       Status status = getStatus(iSequence);
    1200       if (status==atLowerBound) {
    1201         double movement=0.0;
    1202 
    1203         if (value<-tolerance) {
    1204           // to upper bound
    1205           // put back alpha
    1206           which[numberInfeasibilities++]=iSequence;
    1207           movement = upper[iSequence] - lower[iSequence];
    1208           changeObj += movement*cost[iSequence];
    1209           matrix_->add(this,outputArray,iSequence,movement);
    1210         } else if (value<tolerance) {
    1211           // at correct bound but may swap
    1212           FakeBound bound = getFakeBound(iSequence);
    1213           if (bound==ClpPdco::lowerFake) {
    1214             movement = upper[iSequence]-lower[iSequence];
    1215             setStatus(iSequence,atUpperBound);
    1216             solution[iSequence] = upper[iSequence];
    1217             changeObj += movement*cost[iSequence];
    1218             numberFake_--;
    1219             setFakeBound(iSequence,noFake);
    1220           }
    1221         }
    1222       } else if (status==atUpperBound) {
    1223         double movement=0.0;
    1224 
    1225         if (value>tolerance) {
    1226           // to lower bound (if swap)
    1227           // put back alpha
    1228           which[numberInfeasibilities++]=iSequence;
    1229           movement = lower[iSequence]-upper[iSequence];
    1230           changeObj += movement*cost[iSequence];
    1231           matrix_->add(this,outputArray,iSequence,movement);
    1232         } else if (value>-tolerance) {
    1233           // at correct bound but may swap
    1234           FakeBound bound = getFakeBound(iSequence);
    1235           if (bound==ClpPdco::upperFake) {
    1236             movement = lower[iSequence]-upper[iSequence];
    1237             setStatus(iSequence,atLowerBound);
    1238             solution[iSequence] = lower[iSequence];
    1239             changeObj += movement*cost[iSequence];
    1240             numberFake_--;
    1241             setFakeBound(iSequence,noFake);
    1242           }
    1243         }
    1244       }
    1245     }
    1246   }
    1247 #ifdef CLP_DEBUG
    1248   if (fullRecompute&&numberFake_&&(handler_->logLevel()&16)!=0)
    1249     printf("%d fake after full update\n",numberFake_);
    1250 #endif
    1251   // set number of infeasibilities
    1252   columnArray->setNumElements(numberInfeasibilities);
    1253   numberInfeasibilities += numberRowInfeasibilities;
    1254   if (fullRecompute) {
    1255     // do actual flips
    1256     flipBounds(rowArray,columnArray,0.0);
    1257   }
    1258   objectiveChange += changeObj;
    1259   return numberInfeasibilities;
    1260928}
    1261 void
    1262 ClpPdco::updateDualsInValuesPass(CoinIndexedVector * rowArray,
    1263                                   CoinIndexedVector * columnArray,
    1264                                         double theta)
    1265 {
    1266  
    1267   // use a tighter tolerance except for all being okay
    1268   double tolerance = dualTolerance_;
    1269  
    1270   // Coding is very similar but we can save a bit by splitting
    1271   // Do rows
    1272   {
    1273     int i;
    1274     double * reducedCost = djRegion(0);
    1275     double * work;
    1276     int number;
    1277     int * which;
    1278     work = rowArray->denseVector();
    1279     number = rowArray->getNumElements();
    1280     which = rowArray->getIndices();
    1281     for (i=0;i<number;i++) {
    1282       int iSequence = which[i];
    1283       double alphaI = work[i];
    1284       double value = reducedCost[iSequence]-theta*alphaI;
    1285       work[i]=0.0;
    1286       reducedCost[iSequence]=value;
    1287      
    1288       Status status = getStatus(iSequence+numberColumns_);
    1289       // more likely to be at upper bound ?
    1290       if (status==atUpperBound) {
    1291 
    1292         if (value>tolerance)
    1293           reducedCost[iSequence]=0.0;
    1294       } else if (status==atLowerBound) {
    1295 
    1296         if (value<-tolerance) {
    1297           reducedCost[iSequence]=0.0;
    1298         }
    1299       }
    1300     }
    1301   }
    1302   rowArray->setNumElements(0);
    1303 
    1304   // Do columns
    1305   {
    1306     int i;
    1307     double * reducedCost = djRegion(1);
    1308     double * work;
    1309     int number;
    1310     int * which;
    1311     work = columnArray->denseVector();
    1312     number = columnArray->getNumElements();
    1313     which = columnArray->getIndices();
    1314    
    1315     for (i=0;i<number;i++) {
    1316       int iSequence = which[i];
    1317       double alphaI = work[i];
    1318       double value = reducedCost[iSequence]-theta*alphaI;
    1319       work[i]=0.0;
    1320       reducedCost[iSequence]=value;
    1321      
    1322       Status status = getStatus(iSequence);
    1323       if (status==atLowerBound) {
    1324         if (value<-tolerance)
    1325           reducedCost[iSequence]=0.0;
    1326       } else if (status==atUpperBound) {
    1327         if (value>tolerance)
    1328           reducedCost[iSequence]=0.0;
    1329       }
    1330     }
    1331   }
    1332   columnArray->setNumElements(0);
    1333 }
    1334 /*
    1335    Chooses dual pivot row
    1336    Would be faster with separate region to scan
    1337    and will have this (with square of infeasibility) when steepest
    1338    For easy problems we can just choose one of the first rows we look at
    1339 */
    1340 void
    1341 ClpPdco::dualRow(int alreadyChosen)
    1342 {
    1343   // get pivot row using whichever method it is
    1344   int chosenRow=-1;
    1345   if (alreadyChosen<0) {
    1346     // first see if any free variables and put them in basis
    1347     int nextFree = nextSuperBasic();
    1348     //nextFree=-1; //off
    1349     if (nextFree>=0) {
    1350       // unpack vector and find a good pivot
    1351       unpack(rowArray_[1],nextFree);
    1352       factorization_->updateColumn(rowArray_[2],rowArray_[1]);
    1353      
    1354       double * work=rowArray_[1]->denseVector();
    1355       int number=rowArray_[1]->getNumElements();
    1356       int * which=rowArray_[1]->getIndices();
    1357       double bestFeasibleAlpha=0.0;
    1358       int bestFeasibleRow=-1;
    1359       double bestInfeasibleAlpha=0.0;
    1360       int bestInfeasibleRow=-1;
    1361       int i;
    1362      
    1363       for (i=0;i<number;i++) {
    1364         int iRow = which[i];
    1365         double alpha = fabs(work[iRow]);
    1366         if (alpha>1.0e-3) {
    1367           int iSequence=pivotVariable_[iRow];
    1368           double value = solution_[iSequence];
    1369           double lower = lower_[iSequence];
    1370           double upper = upper_[iSequence];
    1371           double infeasibility=0.0;
    1372           if (value>upper)
    1373             infeasibility = value-upper;
    1374           else if (value<lower)
    1375             infeasibility = lower-value;
    1376           if (infeasibility*alpha>bestInfeasibleAlpha&&alpha>1.0e-1) {
    1377             if (!flagged(iSequence)) {
    1378               bestInfeasibleAlpha = infeasibility*alpha;
    1379               bestInfeasibleRow=iRow;
    1380             }
    1381           }
    1382           if (alpha>bestFeasibleAlpha&&(lower>-1.0e20||upper<1.0e20)) {
    1383             bestFeasibleAlpha = alpha;
    1384             bestFeasibleRow=iRow;
    1385           }
    1386         }
    1387       }
    1388       if (bestInfeasibleRow>=0)
    1389         chosenRow = bestInfeasibleRow;
    1390       else if (bestFeasibleAlpha>1.0e-2)
    1391         chosenRow = bestFeasibleRow;
    1392       if (chosenRow>=0)
    1393         pivotRow_=chosenRow;
    1394       rowArray_[1]->clear();
    1395     }
    1396   } else {
    1397     // in values pass
    1398     chosenRow=alreadyChosen;
    1399     pivotRow_=chosenRow;
    1400   }
    1401   if (chosenRow<0)
    1402     pivotRow_=dualRowPivot_->pivotRow();
    1403 
    1404   if (pivotRow_>=0) {
    1405     sequenceOut_ = pivotVariable_[pivotRow_];
    1406     valueOut_ = solution_[sequenceOut_];
    1407     lowerOut_ = lower_[sequenceOut_];
    1408     upperOut_ = upper_[sequenceOut_];
    1409     if (alreadyChosen<0) {
    1410       // if we have problems we could try other way and hope we get a
    1411       // zero pivot?
    1412       if (valueOut_>upperOut_) {
    1413         directionOut_ = -1;
    1414         dualOut_ = valueOut_ - upperOut_;
    1415       } else if (valueOut_<lowerOut_) {
    1416         directionOut_ = 1;
    1417         dualOut_ = lowerOut_ - valueOut_;
    1418       } else {
    1419         // odd (could be free) - it's feasible - go to nearest
    1420         if (valueOut_-lowerOut_<upperOut_-valueOut_) {
    1421           directionOut_ = 1;
    1422           dualOut_ = lowerOut_ - valueOut_;
    1423         } else {
    1424           directionOut_ = -1;
    1425           dualOut_ = valueOut_ - upperOut_;
    1426         }
    1427       }
    1428 #ifdef CLP_DEBUG
    1429       assert(dualOut_>=0.0);
    1430 #endif
    1431     } else {
    1432       // in values pass so just use sign of dj
    1433       // We don't want to go through any barriers so set dualOut low
    1434       // free variables will never be here
    1435       dualOut_ = 1.0e-6;
    1436       if (dj_[sequenceOut_]>0.0) {
    1437         // this will give a -1 in pivot row (as slacks are -1.0)
    1438         directionOut_ = 1;
    1439       } else {
    1440         directionOut_ = -1;
    1441       }
    1442     }
    1443   }
    1444   return ;
    1445 }
    1446 // Checks if any fake bounds active - if so returns number and modifies
    1447 // dualBound_ and everything.
    1448 // Free variables will be left as free
    1449 // Returns number of bounds changed if >=0
    1450 // Returns -1 if not initialize and no effect
    1451 // Fills in changeVector which can be used to see if unbounded
    1452 // and cost of change vector
    1453 int
    1454 ClpPdco::changeBounds(bool initialize,
    1455                                  CoinIndexedVector * outputArray,
    1456                                  double & changeCost)
    1457 {
    1458   numberFake_ = 0;
    1459   if (!initialize) {
    1460     int numberInfeasibilities;
    1461     double newBound;
    1462     newBound = 5.0*dualBound_;
    1463     numberInfeasibilities=0;
    1464     changeCost=0.0;
    1465     // put back original bounds and then check
    1466     createRim(3);
    1467     int iSequence;
    1468     // bounds will get bigger - just look at ones at bounds
    1469     for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
    1470       double lowerValue=lower_[iSequence];
    1471       double upperValue=upper_[iSequence];
    1472       double value=solution_[iSequence];
    1473       setFakeBound(iSequence,ClpPdco::noFake);
    1474       switch(getStatus(iSequence)) {
    1475 
    1476       case basic:
    1477       case ClpSimplex::isFixed:
    1478         break;
    1479       case isFree:
    1480       case superBasic:
    1481         break;
    1482       case atUpperBound:
    1483         if (fabs(value-upperValue)>primalTolerance_)
    1484           numberInfeasibilities++;
    1485         break;
    1486       case atLowerBound:
    1487         if (fabs(value-lowerValue)>primalTolerance_)
    1488           numberInfeasibilities++;
    1489         break;
    1490       }
    1491     }
    1492     // If dual infeasible then carry on
    1493     if (numberInfeasibilities) {
    1494       handler_->message(CLP_DUAL_CHECKB,messages_)
    1495         <<newBound
    1496         <<CoinMessageEol;
    1497       int iSequence;
    1498       for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
    1499         double lowerValue=lower_[iSequence];
    1500         double upperValue=upper_[iSequence];
    1501         double newLowerValue;
    1502         double newUpperValue;
    1503         Status status = getStatus(iSequence);
    1504         if (status==atUpperBound||
    1505             status==atLowerBound) {
    1506           double value = solution_[iSequence];
    1507           if (value-lowerValue<=upperValue-value) {
    1508             newLowerValue = max(lowerValue,value-0.666667*newBound);
    1509             newUpperValue = min(upperValue,newLowerValue+newBound);
    1510           } else {
    1511             newUpperValue = min(upperValue,value+0.666667*newBound);
    1512             newLowerValue = max(lowerValue,newUpperValue-newBound);
    1513           }
    1514           lower_[iSequence]=newLowerValue;
    1515           upper_[iSequence]=newUpperValue;
    1516           if (newLowerValue > lowerValue) {
    1517             if (newUpperValue < upperValue) {
    1518               setFakeBound(iSequence,ClpPdco::bothFake);
    1519               numberFake_++;
    1520             } else {
    1521               setFakeBound(iSequence,ClpPdco::lowerFake);
    1522               numberFake_++;
    1523             }
    1524           } else {
    1525             if (newUpperValue < upperValue) {
    1526               setFakeBound(iSequence,ClpPdco::upperFake);
    1527               numberFake_++;
    1528             }
    1529           }
    1530           if (status==atUpperBound)
    1531             solution_[iSequence] = newUpperValue;
    1532           else
    1533             solution_[iSequence] = newLowerValue;
    1534           double movement = solution_[iSequence] - value;
    1535           if (movement&&outputArray) {
    1536             if (iSequence>=numberColumns_) {
    1537               outputArray->quickAdd(iSequence,-movement);
    1538               changeCost += movement*cost_[iSequence];
    1539             } else {
    1540               matrix_->add(this,outputArray,iSequence,movement);
    1541               changeCost += movement*cost_[iSequence];
    1542             }
    1543           }
    1544         }
    1545       }
    1546       dualBound_ = newBound;
    1547     } else {
    1548       numberInfeasibilities=-1;
    1549     }
    1550     return numberInfeasibilities;
    1551   } else {
    1552     int iSequence;
    1553      
    1554     for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
    1555       Status status = getStatus(iSequence);
    1556       if (status==atUpperBound||
    1557           status==atLowerBound) {
    1558         double lowerValue=lower_[iSequence];
    1559         double upperValue=upper_[iSequence];
    1560         double value = solution_[iSequence];
    1561         if (lowerValue>-largeValue_||upperValue<largeValue_) {
    1562           if (lowerValue-value>-0.5*dualBound_||
    1563               upperValue-value<0.5*dualBound_) {
    1564             if (fabs(lowerValue-value)<=fabs(upperValue-value)) {
    1565               if (upperValue > lowerValue + dualBound_) {
    1566                 upper_[iSequence]=lowerValue+dualBound_;
    1567                 setFakeBound(iSequence,ClpPdco::upperFake);
    1568                 numberFake_++;
    1569               }
    1570             } else {
    1571               if (lowerValue < upperValue - dualBound_) {
    1572                 lower_[iSequence]=upperValue-dualBound_;
    1573                 setFakeBound(iSequence,ClpPdco::lowerFake);
    1574                 numberFake_++;
    1575               }
    1576             }
    1577           } else {
    1578             lower_[iSequence]=-0.5*dualBound_;
    1579             upper_[iSequence]= 0.5*dualBound_;
    1580             setFakeBound(iSequence,ClpPdco::bothFake);
    1581             numberFake_++;
    1582           }
    1583         } else {
    1584           // set non basic free variables to fake bounds
    1585           // I don't think we should ever get here
    1586           assert(!("should not be here"));
    1587           lower_[iSequence]=-0.5*dualBound_;
    1588           upper_[iSequence]= 0.5*dualBound_;
    1589           setFakeBound(iSequence,ClpPdco::bothFake);
    1590           numberFake_++;
    1591           setStatus(iSequence,atUpperBound);
    1592           solution_[iSequence]=0.5*dualBound_;
    1593         }
    1594       }
    1595     }
    1596 
    1597     return 1;
    1598   }
    1599 }
    1600 /*
    1601    Row array has row part of pivot row (as duals so sign may be switched)
    1602    Column array has column part.
    1603    This chooses pivot column.
    1604    Spare array will be needed when we start getting clever.
    1605    We will check for basic so spare array will never overflow.
    1606    If necessary will modify costs
    1607 */
    1608 void
    1609 ClpPdco::dualColumn(CoinIndexedVector * rowArray,
    1610                            CoinIndexedVector * columnArray,
    1611                            CoinIndexedVector * spareArray,
    1612                            CoinIndexedVector * spareArray2,
    1613                            double acceptablePivot,
    1614                            CoinBigIndex * dubiousWeights)
    1615 {
    1616   double * work;
    1617   int number;
    1618   int * which;
    1619   double * reducedCost;
    1620   int iSection;
    1621 
    1622   sequenceIn_=-1;
    1623   int numberPossiblySwapped=0;
    1624   int numberRemaining=0;
    1625  
    1626   double totalThru=0.0; // for when variables flip
    1627 
    1628   double bestEverPivot=acceptablePivot;
    1629   int lastSequence = -1;
    1630   double lastPivot=0.0;
    1631   double upperTheta;
    1632   double newTolerance = dualTolerance_;
    1633   //newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance];
    1634   // will we need to increase tolerance
    1635   bool thisIncrease=false;
    1636   // If we think we need to modify costs (not if something from broad sweep)
    1637   bool modifyCosts=false;
    1638   // Increase in objective due to swapping bounds (may be negative)
    1639   double increaseInObjective=0.0;
    1640 
    1641   // use spareArrays to put ones looked at in
    1642   // we are going to flip flop between
    1643   int iFlip = 0;
    1644   // Possible list of pivots
    1645   int interesting[2];
    1646   // where possible swapped ones are
    1647   int swapped[2];
    1648   // for zeroing out arrays after
    1649   int marker[2][2];
    1650   // pivot elements
    1651   double * array[2], * spare, * spare2;
    1652   // indices
    1653   int * indices[2], * index, * index2;
    1654   spareArray->clear();
    1655   spareArray2->clear();
    1656   array[0] = spareArray->denseVector();
    1657   indices[0] = spareArray->getIndices();
    1658   spare = array[0];
    1659   index = indices[0];
    1660   array[1] = spareArray2->denseVector();
    1661   indices[1] = spareArray2->getIndices();
    1662   int i;
    1663   const double * lower;
    1664   const double * upper;
    1665 
    1666   // initialize lists
    1667   for (i=0;i<2;i++) {
    1668     interesting[i]=0;
    1669     swapped[i]=numberColumns_;
    1670     marker[i][0]=0;
    1671     marker[i][1]=numberColumns_;
    1672   }
    1673 
    1674   /*
    1675     First we get a list of possible pivots.  We can also see if the
    1676     problem looks infeasible or whether we want to pivot in free variable.
    1677     This may make objective go backwards but can only happen a finite
    1678     number of times and I do want free variables basic.
    1679 
    1680     Then we flip back and forth.  At the start of each iteration
    1681     interesting[iFlip] should have possible candidates and swapped[iFlip]
    1682     will have pivots if we decide to take a previous pivot.
    1683     At end of each iteration interesting[1-iFlip] should have
    1684     candidates if we go through this theta and swapped[1-iFlip]
    1685     pivots if we don't go through.
    1686 
    1687     At first we increase theta and see what happens.  We start
    1688     theta at a reasonable guess.  If in right area then we do bit by bit.
    1689 
    1690    */
    1691 
    1692   // do first pass to get possibles
    1693   // We can also see if infeasible or pivoting on free
    1694   double tentativeTheta = 1.0e22;
    1695   upperTheta = 1.0e31;
    1696   double freePivot = acceptablePivot;
    1697 
    1698   for (iSection=0;iSection<2;iSection++) {
    1699 
    1700     int addSequence;
    1701 
    1702     if (!iSection) {
    1703       lower = rowLowerWork_;
    1704       upper = rowUpperWork_;
    1705       work = rowArray->denseVector();
    1706       number = rowArray->getNumElements();
    1707       which = rowArray->getIndices();
    1708       reducedCost = rowReducedCost_;
    1709       addSequence = numberColumns_;
    1710     } else {
    1711       lower = columnLowerWork_;
    1712       upper = columnUpperWork_;
    1713       work = columnArray->denseVector();
    1714       number = columnArray->getNumElements();
    1715       which = columnArray->getIndices();
    1716       reducedCost = reducedCostWork_;
    1717       addSequence = 0;
    1718     }
    1719    
    1720     for (i=0;i<number;i++) {
    1721       int iSequence = which[i];
    1722       double alpha;
    1723       double oldValue;
    1724       double value;
    1725       bool keep;
    1726 
    1727       switch(getStatus(iSequence+addSequence)) {
    1728          
    1729       case basic:
    1730       case ClpSimplex::isFixed:
    1731         break;
    1732       case isFree:
    1733       case superBasic:
    1734         alpha = work[i];
    1735         oldValue = reducedCost[iSequence];
    1736         if (oldValue>dualTolerance_) {
    1737           keep = true;
    1738         } else if (oldValue<-dualTolerance_) {
    1739           keep = true;
    1740         } else {
    1741           if (fabs(alpha)>10.0*acceptablePivot)
    1742             keep = true;
    1743           else
    1744             keep=false;
    1745         }
    1746         if (keep) {
    1747           // free - choose largest
    1748           if (fabs(alpha)>freePivot) {
    1749             freePivot=fabs(alpha);
    1750             sequenceIn_ = iSequence + addSequence;
    1751             theta_=oldValue/alpha;
    1752             alpha_=alpha;
    1753           }
    1754         }
    1755         break;
    1756       case atUpperBound:
    1757         alpha = work[i];
    1758         oldValue = reducedCost[iSequence];
    1759         value = oldValue-tentativeTheta*alpha;
    1760         assert (oldValue<=dualTolerance_*1.0001);
    1761         if (value>newTolerance) {
    1762           value = oldValue-upperTheta*alpha;
    1763           if (value>newTolerance && -alpha>=acceptablePivot)
    1764             upperTheta = (oldValue-newTolerance)/alpha;
    1765           // add to list
    1766           spare[numberRemaining]=alpha;
    1767           index[numberRemaining++]=iSequence+addSequence;
    1768         }
    1769         break;
    1770       case atLowerBound:
    1771         alpha = work[i];
    1772         oldValue = reducedCost[iSequence];
    1773         value = oldValue-tentativeTheta*alpha;
    1774         assert (oldValue>=-dualTolerance_*1.0001);
    1775         if (value<-newTolerance) {
    1776           value = oldValue-upperTheta*alpha;
    1777           if (value<-newTolerance && alpha>=acceptablePivot)
    1778             upperTheta = (oldValue+newTolerance)/alpha;
    1779           // add to list
    1780           spare[numberRemaining]=alpha;
    1781           index[numberRemaining++]=iSequence+addSequence;
    1782         }
    1783         break;
    1784       }
    1785     }
    1786   }
    1787   interesting[0]=numberRemaining;
    1788   marker[0][0] = numberRemaining;
    1789 
    1790   if (!numberRemaining&&sequenceIn_<0)
    1791     return ; // Looks infeasible
    1792 
    1793   if (sequenceIn_>=0) {
    1794     // free variable - always choose
    1795   } else {
    1796 
    1797     theta_=1.0e50;
    1798     // now flip flop between spare arrays until reasonable theta
    1799     tentativeTheta = max(10.0*upperTheta,1.0e-7);
    1800    
    1801     // loops increasing tentative theta until can't go through
    1802    
    1803     while (tentativeTheta < 1.0e22) {
    1804       double thruThis = 0.0;
    1805      
    1806       double bestPivot=acceptablePivot;
    1807       int bestSequence=-1;
    1808      
    1809       numberPossiblySwapped = numberColumns_;
    1810       numberRemaining = 0;
    1811 
    1812       upperTheta = 1.0e50;
    1813 
    1814       spare = array[iFlip];
    1815       index = indices[iFlip];
    1816       spare2 = array[1-iFlip];
    1817       index2 = indices[1-iFlip];
    1818 
    1819       // try 3 different ways
    1820       // 1 bias increase by ones with slightly wrong djs
    1821       // 2 bias by all
    1822       // 3 bias by all - tolerance (doesn't seem very good)
    1823 #define TRYBIAS 1
    1824 
    1825 
    1826       double increaseInThis=0.0; //objective increase in this loop
    1827      
    1828       for (i=0;i<interesting[iFlip];i++) {
    1829         int iSequence = index[i];
    1830         double alpha = spare[i];
    1831         double oldValue = dj_[iSequence];
    1832         double value = oldValue-tentativeTheta*alpha;
    1833 
    1834         if (alpha < 0.0) {
    1835           //at upper bound
    1836           if (value>newTolerance) {
    1837             double range = upper_[iSequence] - lower_[iSequence];
    1838             thruThis -= range*alpha;
    1839 #if TRYBIAS==1
    1840             if (oldValue>0.0)
    1841               increaseInThis -= oldValue*range;
    1842 #elif TRYBIAS==2
    1843             increaseInThis -= oldValue*range;
    1844 #else
    1845             increaseInThis -= (oldValue+dualTolerance_)*range;
    1846 #endif
    1847             // goes on swapped list (also means candidates if too many)
    1848             spare2[--numberPossiblySwapped]=alpha;
    1849             index2[numberPossiblySwapped]=iSequence;
    1850             if (fabs(alpha)>bestPivot) {
    1851               bestPivot=fabs(alpha);
    1852               bestSequence=numberPossiblySwapped;
    1853             }
    1854           } else {
    1855             value = oldValue-upperTheta*alpha;
    1856             if (value>newTolerance && -alpha>=acceptablePivot)
    1857               upperTheta = (oldValue-newTolerance)/alpha;
    1858             spare2[numberRemaining]=alpha;
    1859             index2[numberRemaining++]=iSequence;
    1860           }
    1861         } else {
    1862           // at lower bound
    1863           if (value<-newTolerance) {
    1864             double range = upper_[iSequence] - lower_[iSequence];
    1865             thruThis += range*alpha;
    1866             //?? is this correct - and should we look at good ones
    1867 #if TRYBIAS==1
    1868             if (oldValue<0.0)
    1869               increaseInThis += oldValue*range;
    1870 #elif TRYBIAS==2
    1871             increaseInThis += oldValue*range;
    1872 #else
    1873             increaseInThis += (oldValue-dualTolerance_)*range;
    1874 #endif
    1875             // goes on swapped list (also means candidates if too many)
    1876             spare2[--numberPossiblySwapped]=alpha;
    1877             index2[numberPossiblySwapped]=iSequence;
    1878             if (fabs(alpha)>bestPivot) {
    1879               bestPivot=fabs(alpha);
    1880               bestSequence=numberPossiblySwapped;
    1881             }
    1882           } else {
    1883             value = oldValue-upperTheta*alpha;
    1884             if (value<-newTolerance && alpha>=acceptablePivot)
    1885               upperTheta = (oldValue+newTolerance)/alpha;
    1886             spare2[numberRemaining]=alpha;
    1887             index2[numberRemaining++]=iSequence;
    1888           }
    1889         }
    1890       }
    1891       swapped[1-iFlip]=numberPossiblySwapped;
    1892       interesting[1-iFlip]=numberRemaining;
    1893       marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
    1894       marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
    1895      
    1896       if (totalThru+thruThis>=fabs(dualOut_)||
    1897           increaseInObjective+increaseInThis<0.0) {
    1898         // We should be pivoting in this batch
    1899         // so compress down to this lot
    1900         numberRemaining=0;
    1901         for (i=numberColumns_-1;i>=swapped[1-iFlip];i--) {
    1902           spare[numberRemaining]=spare2[i];
    1903           index[numberRemaining++]=index2[i];
    1904         }
    1905         interesting[iFlip]=numberRemaining;
    1906         int iTry;
    1907 #define MAXTRY 100
    1908         // first get ratio with tolerance
    1909         for (iTry=0;iTry<MAXTRY;iTry++) {
    1910          
    1911           upperTheta=1.0e50;
    1912           numberPossiblySwapped = numberColumns_;
    1913           numberRemaining = 0;
    1914 
    1915           increaseInThis=0.0; //objective increase in this loop
    1916 
    1917           thruThis=0.0;
    1918 
    1919           spare = array[iFlip];
    1920           index = indices[iFlip];
    1921           spare2 = array[1-iFlip];
    1922           index2 = indices[1-iFlip];
    1923      
    1924           for (i=0;i<interesting[iFlip];i++) {
    1925             int iSequence=index[i];
    1926             double alpha=spare[i];
    1927             double oldValue = dj_[iSequence];
    1928             double value = oldValue-upperTheta*alpha;
    1929            
    1930             if (alpha < 0.0) {
    1931               //at upper bound
    1932               if (value>newTolerance) {
    1933                 if (-alpha>=acceptablePivot) {
    1934                   upperTheta = (oldValue-newTolerance)/alpha;
    1935                   // recompute value and make sure works
    1936                   value = oldValue-upperTheta*alpha;
    1937                   if (value<0)
    1938                     upperTheta *= 1.0 +1.0e-11; // must be large
    1939                 }
    1940               }
    1941             } else {
    1942               // at lower bound
    1943               if (value<-newTolerance) {
    1944                 if (alpha>=acceptablePivot) {
    1945                   upperTheta = (oldValue+newTolerance)/alpha;
    1946                   // recompute value and make sure works
    1947                   value = oldValue-upperTheta*alpha;
    1948                   if (value>0)
    1949                     upperTheta *= 1.0 +1.0e-11; // must be large
    1950                 }
    1951               }
    1952             }
    1953           }
    1954           bestPivot=acceptablePivot;
    1955           sequenceIn_=-1;
    1956           double bestWeight=COIN_DBL_MAX;
    1957           double largestPivot=acceptablePivot;
    1958           // now choose largest and sum all ones which will go through
    1959           //printf("XX it %d number %d\n",numberIterations_,interesting[iFlip]);
    1960           for (i=0;i<interesting[iFlip];i++) {
    1961             int iSequence=index[i];
    1962             double alpha=spare[i];
    1963             double value = dj_[iSequence]-upperTheta*alpha;
    1964             double badDj=0.0;
    1965            
    1966             bool addToSwapped=false;
    1967            
    1968             if (alpha < 0.0) {
    1969               //at upper bound
    1970               if (value>=0.0) {
    1971                 addToSwapped=true;
    1972 #if TRYBIAS==1
    1973                 badDj = -max(dj_[iSequence],0.0);
    1974 #elif TRYBIAS==2
    1975                 badDj = -dj_[iSequence];
    1976 #else
    1977                 badDj = -dj_[iSequence]-dualTolerance_;
    1978 #endif
    1979               }
    1980             } else {
    1981               // at lower bound
    1982               if (value<=0.0) {
    1983                 addToSwapped=true;
    1984 #if TRYBIAS==1
    1985                 badDj = min(dj_[iSequence],0.0);
    1986 #elif TRYBIAS==2
    1987                 badDj = dj_[iSequence];
    1988 #else
    1989                 badDj = dj_[iSequence]-dualTolerance_;
    1990 #endif
    1991               }
    1992             }
    1993             if (!addToSwapped) {
    1994               // add to list of remaining
    1995               spare2[numberRemaining]=alpha;
    1996               index2[numberRemaining++]=iSequence;
    1997             } else {
    1998               // add to list of swapped
    1999               spare2[--numberPossiblySwapped]=alpha;
    2000               index2[numberPossiblySwapped]=iSequence;
    2001               // select if largest pivot
    2002               bool take=false;
    2003               double absAlpha = fabs(alpha);
    2004               // User could do anything to break ties here
    2005               double weight;
    2006               if (dubiousWeights)
    2007                 weight=dubiousWeights[iSequence];
    2008               else
    2009                 weight=1.0;
    2010               weight += CoinDrand48()*1.0e-2;
    2011               if (absAlpha>2.0*bestPivot) {
    2012                 take=true;
    2013               } else if (absAlpha>0.5*largestPivot) {
    2014                 // could multiply absAlpha and weight
    2015                 if (weight*bestPivot<bestWeight*absAlpha)
    2016                   take=true;
    2017               }
    2018               if (take) {
    2019                 sequenceIn_ = numberPossiblySwapped;
    2020                 bestPivot =  absAlpha;
    2021                 theta_ = dj_[iSequence]/alpha;
    2022                 largestPivot = max(largestPivot,bestPivot);
    2023                 bestWeight = weight;
    2024                 //printf(" taken seq %d alpha %g weight %d\n",
    2025                 //   iSequence,absAlpha,dubiousWeights[iSequence]);
    2026               } else {
    2027                 //printf(" not taken seq %d alpha %g weight %d\n",
    2028                 //   iSequence,absAlpha,dubiousWeights[iSequence]);
    2029               }
    2030               double range = upper[iSequence] - lower[iSequence];
    2031               thruThis += range*fabs(alpha);
    2032               increaseInThis += badDj*range;
    2033             }
    2034           }
    2035           swapped[1-iFlip]=numberPossiblySwapped;
    2036           interesting[1-iFlip]=numberRemaining;
    2037           marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
    2038           marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
    2039           // If we stop now this will be increase in objective (I think)
    2040           double increase = (fabs(dualOut_)-totalThru)*theta_;
    2041           increase += increaseInObjective;
    2042           if (theta_<0.0)
    2043             thruThis += fabs(dualOut_); // force using this one
    2044           if (increaseInObjective<0.0&&increase<0.0&&lastSequence>=0) {
    2045             // back
    2046             // We may need to be more careful - we could do by
    2047             // switch so we always do fine grained?
    2048             bestPivot=0.0;
    2049           } else {
    2050             // add in
    2051             totalThru += thruThis;
    2052             increaseInObjective += increaseInThis;
    2053           }
    2054           if (bestPivot<0.1*bestEverPivot&&
    2055               bestEverPivot>1.0e-6&&
    2056               (bestPivot<1.0e-3||totalThru*2.0>fabs(dualOut_))) {
    2057             // back to previous one
    2058             sequenceIn_=lastSequence;
    2059             // swap regions
    2060             iFlip = 1-iFlip;
    2061             break;
    2062           } else if (sequenceIn_==-1&&upperTheta>largeValue_) {
    2063             if (lastPivot>acceptablePivot) {
    2064               // back to previous one
    2065               sequenceIn_=lastSequence;
    2066               // swap regions
    2067               iFlip = 1-iFlip;
    2068             } else {
    2069               // can only get here if all pivots too small
    2070             }
    2071             break;
    2072           } else if (totalThru>=fabs(dualOut_)) {
    2073             modifyCosts=true; // fine grain - we can modify costs
    2074             break; // no point trying another loop
    2075           } else {
    2076             lastSequence=sequenceIn_;
    2077             if (bestPivot>bestEverPivot)
    2078               bestEverPivot=bestPivot;
    2079             iFlip = 1 -iFlip;
    2080             modifyCosts=true; // fine grain - we can modify costs
    2081           }
    2082         }
    2083         if (iTry==MAXTRY)
    2084           iFlip = 1-iFlip; // flip back
    2085         break;
    2086       } else {
    2087         // skip this lot
    2088         if (bestPivot>1.0e-3||bestPivot>bestEverPivot) {
    2089           bestEverPivot=bestPivot;
    2090           lastSequence=bestSequence;
    2091         } else {
    2092           // keep old swapped
    2093           memcpy(array[1-iFlip]+swapped[iFlip],
    2094                  array[iFlip]+swapped[iFlip],
    2095                  (numberColumns_-swapped[iFlip])*sizeof(double));
    2096           memcpy(indices[1-iFlip]+swapped[iFlip],
    2097                  indices[iFlip]+swapped[iFlip],
    2098                  (numberColumns_-swapped[iFlip])*sizeof(int));
    2099           marker[1-iFlip][1] = min(marker[1-iFlip][1],swapped[iFlip]);
    2100           swapped[1-iFlip]=swapped[iFlip];
    2101         }
    2102         increaseInObjective += increaseInThis;
    2103         iFlip = 1 - iFlip; // swap regions
    2104         tentativeTheta = 2.0*upperTheta;
    2105         totalThru += thruThis;
    2106       }
    2107     }
    2108    
    2109     // can get here without sequenceIn_ set but with lastSequence
    2110     if (sequenceIn_<0&&lastSequence>=0) {
    2111       // back to previous one
    2112       sequenceIn_=lastSequence;
    2113       // swap regions
    2114       iFlip = 1-iFlip;
    2115     }
    2116    
    2117 #define MINIMUMTHETA 1.0e-12
    2118     // Movement should be minimum for anti-degeneracy - unless
    2119     // fixed variable out
    2120     double minimumTheta;
    2121     if (upperOut_>lowerOut_)
    2122       minimumTheta=MINIMUMTHETA;
    2123     else
    2124       minimumTheta=0.0;
    2125     if (sequenceIn_>=0) {
    2126       // at this stage sequenceIn_ is just pointer into index array
    2127       // flip just so we can use iFlip
    2128       iFlip = 1 -iFlip;
    2129       spare = array[iFlip];
    2130       index = indices[iFlip];
    2131       double oldValue;
    2132       alpha_ = spare[sequenceIn_];
    2133       sequenceIn_ = indices[iFlip][sequenceIn_];
    2134       oldValue = dj_[sequenceIn_];
    2135       theta_ = oldValue/alpha_;
    2136       if (theta_<minimumTheta) {
    2137         // can't pivot to zero
    2138         if (oldValue-minimumTheta*alpha_>=-dualTolerance_) {
    2139           theta_=minimumTheta;
    2140         } else if (oldValue-minimumTheta*alpha_>=-newTolerance) {
    2141           theta_=minimumTheta;
    2142           thisIncrease=true;
    2143         } else {
    2144           theta_=(oldValue+newTolerance)/alpha_;
    2145           assert(theta_>=0.0);
    2146           thisIncrease=true;
    2147         }
    2148       }
    2149       // may need to adjust costs so all dual feasible AND pivoted is exactly 0
    2150       //int costOffset = numberRows_+numberColumns_;
    2151       if (modifyCosts) {
    2152         int i;
    2153         for (i=numberColumns_-1;i>=swapped[iFlip];i--) {
    2154           int iSequence=index[i];
    2155           double alpha=spare[i];
    2156           double value = dj_[iSequence]-theta_*alpha;
    2157            
    2158           // can't be free here
    2159          
    2160           if (alpha < 0.0) {
    2161             //at upper bound
    2162             if (value>dualTolerance_) {
    2163               thisIncrease=true;
    2164 #define MODIFYCOST 2
    2165 #if MODIFYCOST
    2166               // modify cost to hit new tolerance
    2167               double modification = alpha*theta_-dj_[iSequence]
    2168                 +newTolerance;
    2169               //modification = min(modification,dualTolerance_);
    2170               //assert (fabs(modification)<1.0e-7);
    2171               dj_[iSequence] += modification;
    2172               cost_[iSequence] +=  modification;
    2173               //cost_[iSequence+costOffset] += modification; // save change
    2174 #endif
    2175             }
    2176           } else {
    2177             // at lower bound
    2178             if (-value>dualTolerance_) {
    2179               thisIncrease=true;
    2180 #if MODIFYCOST
    2181               // modify cost to hit new tolerance
    2182               double modification = alpha*theta_-dj_[iSequence]
    2183                 -newTolerance;
    2184               //modification = max(modification,-dualTolerance_);
    2185               //assert (fabs(modification)<1.0e-7);
    2186               dj_[iSequence] += modification;
    2187               cost_[iSequence] +=  modification;
    2188               //cost_[iSequence+costOffset] += modification; // save change
    2189 #endif
    2190             }
    2191           }
    2192         }
    2193       }
    2194     }
    2195   }
    2196 
    2197   if (sequenceIn_>=0) {
    2198     lowerIn_ = lower_[sequenceIn_];
    2199     upperIn_ = upper_[sequenceIn_];
    2200     valueIn_ = solution_[sequenceIn_];
    2201     dualIn_ = dj_[sequenceIn_];
    2202 
    2203     if (numberTimesOptimal_) {
    2204       // can we adjust cost back closer to original
    2205       //*** add coding
    2206     }
    2207 #if MODIFYCOST>1   
    2208     // modify cost to hit zero exactly
    2209     // so (dualIn_+modification)==theta_*alpha_
    2210     double modification = theta_*alpha_-dualIn_;
    2211     dualIn_ += modification;
    2212     dj_[sequenceIn_]=dualIn_;
    2213     cost_[sequenceIn_] += modification;
    2214     //int costOffset = numberRows_+numberColumns_;
    2215     //cost_[sequenceIn_+costOffset] += modification; // save change
    2216     //assert (fabs(modification)<1.0e-6);
    2217 #ifdef CLP_DEBUG
    2218     if ((handler_->logLevel()&32)&&fabs(modification)>1.0e-15)
    2219       printf("exact %d new cost %g, change %g\n",sequenceIn_,
    2220              cost_[sequenceIn_],modification);
    2221 #endif
    2222 #endif
    2223    
    2224     if (alpha_<0.0) {
    2225       // as if from upper bound
    2226       directionIn_=-1;
    2227       upperIn_=valueIn_;
    2228     } else {
    2229       // as if from lower bound
    2230       directionIn_=1;
    2231       lowerIn_=valueIn_;
    2232     }
    2233   }
    2234   //if (thisIncrease)
    2235   //dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance];
    2236 
    2237   // clear arrays
    2238 
    2239   for (i=0;i<2;i++) {
    2240     memset(array[i],0,marker[i][0]*sizeof(double));
    2241     memset(array[i]+marker[i][1],0,
    2242            (numberColumns_-marker[i][1])*sizeof(double));
    2243   }
    2244 }
    2245 /* Checks if tentative optimal actually means unbounded
    2246    Returns -3 if not, 2 if is unbounded */
    2247 int
    2248 ClpPdco::checkUnbounded(CoinIndexedVector * ray,
    2249                                    CoinIndexedVector * spare,
    2250                                    double changeCost)
    2251 {
    2252   int status=2; // say unbounded
    2253   factorization_->updateColumn(spare,ray);
    2254   // get reduced cost
    2255   int i;
    2256   int number=ray->getNumElements();
    2257   int * index = ray->getIndices();
    2258   double * array = ray->denseVector();
    2259   for (i=0;i<number;i++) {
    2260     int iRow=index[i];
    2261     int iPivot=pivotVariable_[iRow];
    2262     changeCost -= cost(iPivot)*array[iRow];
    2263   }
    2264   double way;
    2265   if (changeCost>0.0) {
    2266     //try going down
    2267     way=1.0;
    2268   } else if (changeCost<0.0) {
    2269     //try going up
    2270     way=-1.0;
    2271   } else {
    2272 #ifdef CLP_DEBUG
    2273     printf("can't decide on up or down\n");
    2274 #endif
    2275     way=0.0;
    2276     status=-3;
    2277   }
    2278   double movement=1.0e10*way; // some largish number
    2279   double zeroTolerance = 1.0e-14*dualBound_;
    2280   for (i=0;i<number;i++) {
    2281     int iRow=index[i];
    2282     int iPivot=pivotVariable_[iRow];
    2283     double arrayValue = array[iRow];
    2284     if (fabs(arrayValue)<zeroTolerance)
    2285       arrayValue=0.0;
    2286     double newValue=solution(iPivot)+movement*arrayValue;
    2287     if (newValue>upper(iPivot)+primalTolerance_||
    2288         newValue<lower(iPivot)-primalTolerance_)
    2289       status=-3; // not unbounded
    2290   }
    2291   if (status==2) {
    2292     // create ray
    2293     delete [] ray_;
    2294     ray_ = new double [numberColumns_];
    2295     ClpFillN(ray_,numberColumns_,0.0);
    2296     for (i=0;i<number;i++) {
    2297       int iRow=index[i];
    2298       int iPivot=pivotVariable_[iRow];
    2299       double arrayValue = array[iRow];
    2300       if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
    2301         ray_[iPivot] = way* array[iRow];
    2302     }
    2303   }
    2304   ray->clear();
    2305   return status;
    2306 }
    2307 /* Checks if finished.  Updates status */
    2308 void
    2309 ClpPdco::statusOfProblemInDual(int & lastCleaned,int type,
    2310                                       ClpSimplexProgress * progress,
    2311                                       double * givenDuals)
    2312 {
    2313   bool normalType=true;
    2314   if (type==2) {
    2315     // trouble - restore solution
    2316     memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
    2317     memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
    2318            numberRows_*sizeof(double));
    2319     memcpy(columnActivityWork_,savedSolution_ ,
    2320            numberColumns_*sizeof(double));
    2321     forceFactorization_=1; // a bit drastic but ..
    2322     changeMade_++; // say something changed
    2323   }
    2324   int tentativeStatus = problemStatus_;
    2325   double changeCost;
    2326   bool unflagVariables = true;
    2327   if (problemStatus_>-3||factorization_->pivots()) {
    2328     // factorize
    2329     // later on we will need to recover from singularities
    2330     // also we could skip if first time
    2331     // save dual weights
    2332     dualRowPivot_->saveWeights(this,1);
    2333     if (type) {
    2334       // is factorization okay?
    2335       if (internalFactorize(1)) {
    2336         // no - restore previous basis
    2337         unflagVariables = false;
    2338         assert (type==1);
    2339         changeMade_++; // say something changed
    2340         memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
    2341         memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
    2342                numberRows_*sizeof(double));
    2343         memcpy(columnActivityWork_,savedSolution_ ,
    2344                numberColumns_*sizeof(double));
    2345         // get correct bounds on all variables
    2346         double dummyChangeCost=0.0;
    2347         changeBounds(true,rowArray_[2],dummyChangeCost);
    2348         // throw away change
    2349         int i;
    2350         for (i=0;i<4;i++)
    2351           rowArray_[i]->clear();
    2352         // need to reject something
    2353         char x = isColumn(sequenceOut_) ? 'C' :'R';
    2354         handler_->message(CLP_SIMPLEX_FLAG,messages_)
    2355           <<x<<sequenceWithin(sequenceOut_)
    2356           <<CoinMessageEol;
    2357         setFlagged(sequenceOut_);
    2358        
    2359         forceFactorization_=1; // a bit drastic but ..
    2360         type = 2;
    2361         //assert (internalFactorize(1)==0);
    2362         if (internalFactorize(1)) {
    2363           memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
    2364           memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
    2365                  numberRows_*sizeof(double));
    2366           memcpy(columnActivityWork_,savedSolution_ ,
    2367                  numberColumns_*sizeof(double));
    2368           // debug
    2369           assert (internalFactorize(1)==0);
    2370         }
    2371       }
    2372     }
    2373     if (problemStatus_!=-4||factorization_->pivots()>10)
    2374       problemStatus_=-3;
    2375   }
    2376   // at this stage status is -3 or -4 if looks infeasible
    2377   // get primal and dual solutions
    2378   gutsOfSolution(givenDuals,NULL);
    2379   // Double check infeasibility if no action
    2380   if (progress->lastIterationNumber(0)==numberIterations_) {
    2381     if (dualRowPivot_->looksOptimal()) {
    2382       numberPrimalInfeasibilities_ = 0;
    2383       sumPrimalInfeasibilities_ = 0.0;
    2384     }
    2385   }
    2386   // Check if looping
    2387   int loop;
    2388   if (!givenDuals&&type!=2)
    2389     loop = progress->looping();
    2390   else
    2391     loop=-1;
    2392   int situationChanged=0;
    2393   if (loop>=0) {
    2394     problemStatus_ = loop; //exit if in loop
    2395     if (!problemStatus_) {
    2396       // declaring victory
    2397       numberPrimalInfeasibilities_ = 0;
    2398       sumPrimalInfeasibilities_ = 0.0;
    2399     }
    2400     return;
    2401   } else if (loop<-1) {
    2402     // something may have changed
    2403     gutsOfSolution(NULL,NULL);
    2404     situationChanged=1;
    2405   }
    2406   // really for free variables in
    2407   if((progressFlag_&2)!=0) {
    2408     situationChanged=2;
    2409   }
    2410   progressFlag_ = 0; //reset progress flag
    2411   handler_->message(CLP_SIMPLEX_STATUS,messages_)
    2412     <<numberIterations_<<objectiveValue();
    2413   handler_->printing(sumPrimalInfeasibilities_>0.0)
    2414     <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
    2415   handler_->printing(sumDualInfeasibilities_>0.0)
    2416     <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
    2417   handler_->printing(numberDualInfeasibilitiesWithoutFree_
    2418                      <numberDualInfeasibilities_)
    2419                        <<numberDualInfeasibilitiesWithoutFree_;
    2420   handler_->message()<<CoinMessageEol;
    2421 
    2422   // dual bound coming in
    2423   //double saveDualBound = dualBound_;
    2424   while (problemStatus_<=-3) {
    2425     int cleanDuals=0;
    2426     if (situationChanged!=0)
    2427       cleanDuals=1;
    2428     int numberChangedBounds=0;
    2429     int doOriginalTolerance=0;
    2430     if ( lastCleaned==numberIterations_)
    2431       doOriginalTolerance=1;
    2432     // check optimal
    2433     // give code benefit of doubt
    2434     if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
    2435         sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
    2436       // say optimal (with these bounds etc)
    2437       numberDualInfeasibilities_ = 0;
    2438       sumDualInfeasibilities_ = 0.0;
    2439       numberPrimalInfeasibilities_ = 0;
    2440       sumPrimalInfeasibilities_ = 0.0;
    2441     }
    2442     //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
    2443     if (dualFeasible()||problemStatus_==-4) {
    2444       progress->modifyObjective(objectiveValue_
    2445                                -sumDualInfeasibilities_*dualBound_);
    2446       if (primalFeasible()) {
    2447         normalType=false;
    2448         // may be optimal - or may be bounds are wrong
    2449         handler_->message(CLP_DUAL_BOUNDS,messages_)
    2450           <<dualBound_
    2451           <<CoinMessageEol;
    2452         // save solution in case unbounded
    2453         ClpDisjointCopyN(columnActivityWork_,numberColumns_,
    2454                           columnArray_[0]->denseVector());
    2455         ClpDisjointCopyN(rowActivityWork_,numberRows_,
    2456                           rowArray_[2]->denseVector());
    2457         numberChangedBounds=changeBounds(false,rowArray_[3],changeCost);
    2458         if (numberChangedBounds<=0&&!numberDualInfeasibilities_) {
    2459           //looks optimal - do we need to reset tolerance
    2460           if (perturbation_==101) {
    2461             perturbation_=102; // stop any perturbations
    2462             cleanDuals=1;
    2463             createRim(4);
    2464           }
    2465           if (lastCleaned<numberIterations_&&numberTimesOptimal_<4) {
    2466             doOriginalTolerance=2;
    2467             numberTimesOptimal_++;
    2468             changeMade_++; // say something changed
    2469             if (numberTimesOptimal_==1) {
    2470               dualTolerance_ = dblParam_[ClpDualTolerance];
    2471               // better to have small tolerance even if slower
    2472               factorization_->zeroTolerance(1.0e-15);
    2473             } else {
    2474               dualTolerance_ = dblParam_[ClpDualTolerance];
    2475               dualTolerance_ *= pow(2.0,numberTimesOptimal_-1);
    2476             }
    2477             cleanDuals=2; // If nothing changed optimal else primal
    2478           } else {
    2479             problemStatus_=0; // optimal
    2480             if (lastCleaned<numberIterations_) {
    2481               handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
    2482                 <<CoinMessageEol;
    2483             }
    2484           }
    2485         } else {
    2486           cleanDuals=1;
    2487           if (doOriginalTolerance==1) {
    2488             // check unbounded
    2489             // find a variable with bad dj
    2490             int iSequence;
    2491             int iChosen=-1;
    2492             double largest = 100.0*primalTolerance_;
    2493             for (iSequence=0;iSequence<numberRows_+numberColumns_;
    2494                  iSequence++) {
    2495               double djValue = dj_[iSequence];
    2496               double originalLo = originalLower(iSequence);
    2497               double originalUp = originalUpper(iSequence);
    2498               if (fabs(djValue)>fabs(largest)) {
    2499                 if (getStatus(iSequence)!=basic) {
    2500                   if (djValue>0&&originalLo<-1.0e20) {
    2501                     if (djValue>fabs(largest)) {
    2502                       largest=djValue;
    2503                       iChosen=iSequence;
    2504                     }
    2505                   } else if (djValue<0&&originalUp>1.0e20) {
    2506                     if (-djValue>fabs(largest)) {
    2507                       largest=djValue;
    2508                       iChosen=iSequence;
    2509                     }
    2510                   }
    2511                 }
    2512               }
    2513             }
    2514             if (iChosen>=0) {
    2515               int iSave=sequenceIn_;
    2516               sequenceIn_=iChosen;
    2517               unpack(rowArray_[1]);
    2518               sequenceIn_ = iSave;
    2519               // if dual infeasibilities then must be free vector so add in dual
    2520               if (numberDualInfeasibilities_) {
    2521                 if (fabs(changeCost)>1.0e-5)
    2522                   printf("Odd free/unbounded combo\n");
    2523                 changeCost += cost_[iChosen];
    2524               }
    2525               problemStatus_ = checkUnbounded(rowArray_[1],rowArray_[0],
    2526                                               changeCost);
    2527               rowArray_[1]->clear();
    2528             } else {
    2529               problemStatus_=-3;
    2530             }
    2531             if (problemStatus_==2&&perturbation_==101) {
    2532               perturbation_=102; // stop any perturbations
    2533               cleanDuals=1;
    2534               createRim(4);
    2535               problemStatus_=-1;
    2536             }
    2537             if (problemStatus_==2) {
    2538               // it is unbounded - restore solution
    2539               // but first add in changes to non-basic
    2540               int iColumn;
    2541               double * original = columnArray_[0]->denseVector();
    2542               for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    2543                 if(getColumnStatus(iColumn)!= basic)
    2544                   ray_[iColumn] +=
    2545                     columnActivityWork_[iColumn]-original[iColumn];
    2546                 columnActivityWork_[iColumn] = original[iColumn];
    2547               }
    2548               ClpDisjointCopyN(rowArray_[2]->denseVector(),numberRows_,
    2549                                 rowActivityWork_);
    2550             }
    2551           } else {
    2552             doOriginalTolerance=2;
    2553             rowArray_[0]->clear();
    2554           }
    2555         }
    2556         ClpFillN(columnArray_[0]->denseVector(),numberColumns_,0.0);
    2557         ClpFillN(rowArray_[2]->denseVector(),numberRows_,0.0);
    2558       }
    2559       if (problemStatus_==-4) {
    2560         // may be infeasible - or may be bounds are wrong
    2561         numberChangedBounds=changeBounds(false,NULL,changeCost);
    2562         if (perturbation_==101) {
    2563           perturbation_=102; // stop any perturbations
    2564           cleanDuals=1;
    2565           numberChangedBounds=1;
    2566           createRim(4);
    2567         }
    2568         if (numberChangedBounds<=0||dualBound_>1.0e20||
    2569             (largestPrimalError_>1.0&&dualBound_>1.0e17)) {
    2570           problemStatus_=1; // infeasible
    2571         } else {
    2572           normalType=false;
    2573           problemStatus_=-1; //iterate
    2574           cleanDuals=1;
    2575           if (numberChangedBounds<=0)
    2576             doOriginalTolerance=2;
    2577           // and delete ray which has been created
    2578           delete [] ray_;
    2579           ray_ = NULL;
    2580         }
    2581 
    2582       }
    2583     } else {
    2584       cleanDuals=1;
    2585     }
    2586     if (problemStatus_<0) {
    2587       if (doOriginalTolerance==2) {
    2588         // put back original tolerance
    2589         lastCleaned=numberIterations_;
    2590         handler_->message(CLP_DUAL_ORIGINAL,messages_)
    2591           <<CoinMessageEol;
    2592 
    2593         perturbation_=102; // stop any perturbations
    2594 #if 0
    2595         double * xcost = new double[numberRows_+numberColumns_];
    2596         double * xlower = new double[numberRows_+numberColumns_];
    2597         double * xupper = new double[numberRows_+numberColumns_];
    2598         double * xdj = new double[numberRows_+numberColumns_];
    2599         double * xsolution = new double[numberRows_+numberColumns_];
    2600         memcpy(xcost,cost_,(numberRows_+numberColumns_)*sizeof(double));
    2601         memcpy(xlower,lower_,(numberRows_+numberColumns_)*sizeof(double));
    2602         memcpy(xupper,upper_,(numberRows_+numberColumns_)*sizeof(double));
    2603         memcpy(xdj,dj_,(numberRows_+numberColumns_)*sizeof(double));
    2604         memcpy(xsolution,solution_,(numberRows_+numberColumns_)*sizeof(double));
    2605 #endif
    2606         createRim(4);
    2607         // make sure duals are current
    2608         computeDuals(givenDuals);
    2609         checkDualSolution();
    2610 #if 0
    2611         int i;
    2612         for (i=0;i<numberRows_+numberColumns_;i++) {
    2613           if (cost_[i]!=xcost[i])
    2614             printf("** %d old cost %g new %g sol %g\n",
    2615                    i,xcost[i],cost_[i],solution_[i]);
    2616           if (lower_[i]!=xlower[i])
    2617             printf("** %d old lower %g new %g sol %g\n",
    2618                    i,xlower[i],lower_[i],solution_[i]);
    2619           if (upper_[i]!=xupper[i])
    2620             printf("** %d old upper %g new %g sol %g\n",
    2621                    i,xupper[i],upper_[i],solution_[i]);
    2622           if (dj_[i]!=xdj[i])
    2623             printf("** %d old dj %g new %g sol %g\n",
    2624                    i,xdj[i],dj_[i],solution_[i]);
    2625           if (solution_[i]!=xsolution[i])
    2626             printf("** %d old solution %g new %g sol %g\n",
    2627                    i,xsolution[i],solution_[i],solution_[i]);
    2628         }
    2629         //delete [] xcost;
    2630         //delete [] xupper;
    2631         //delete [] xlower;
    2632         //delete [] xdj;
    2633         //delete [] xsolution;
    2634 #endif
    2635         // put back bounds as they were if was optimal
    2636         if (doOriginalTolerance==2) {
    2637           changeMade_++; // say something changed
    2638           changeBounds(true,NULL,changeCost);
    2639           cleanDuals=2;
    2640           //cleanDuals=1;
    2641         }
    2642 #if 0
    2643         //int i;
    2644         for (i=0;i<numberRows_+numberColumns_;i++) {
    2645           if (cost_[i]!=xcost[i])
    2646             printf("** %d old cost %g new %g sol %g\n",
    2647                    i,xcost[i],cost_[i],solution_[i]);
    2648           if (lower_[i]!=xlower[i])
    2649             printf("** %d old lower %g new %g sol %g\n",
    2650                    i,xlower[i],lower_[i],solution_[i]);
    2651           if (upper_[i]!=xupper[i])
    2652             printf("** %d old upper %g new %g sol %g\n",
    2653                    i,xupper[i],upper_[i],solution_[i]);
    2654           if (dj_[i]!=xdj[i])
    2655             printf("** %d old dj %g new %g sol %g\n",
    2656                    i,xdj[i],dj_[i],solution_[i]);
    2657           if (solution_[i]!=xsolution[i])
    2658             printf("** %d old solution %g new %g sol %g\n",
    2659                    i,xsolution[i],solution_[i],solution_[i]);
    2660         }
    2661         delete [] xcost;
    2662         delete [] xupper;
    2663         delete [] xlower;
    2664         delete [] xdj;
    2665         delete [] xsolution;
    2666 #endif
    2667       }
    2668       if (cleanDuals==1||(cleanDuals==2&&!numberDualInfeasibilities_)) {
    2669         // make sure dual feasible
    2670         // look at all rows and columns
    2671         rowArray_[0]->clear();
    2672         columnArray_[0]->clear();
    2673         double objectiveChange=0.0;
    2674 #if 0
    2675         double * xcost = new double[numberRows_+numberColumns_];
    2676         double * xlower = new double[numberRows_+numberColumns_];
    2677         double * xupper = new double[numberRows_+numberColumns_];
    2678         double * xdj = new double[numberRows_+numberColumns_];
    2679         double * xsolution = new double[numberRows_+numberColumns_];
    2680         memcpy(xcost,cost_,(numberRows_+numberColumns_)*sizeof(double));
    2681         memcpy(xlower,lower_,(numberRows_+numberColumns_)*sizeof(double));
    2682         memcpy(xupper,upper_,(numberRows_+numberColumns_)*sizeof(double));
    2683         memcpy(xdj,dj_,(numberRows_+numberColumns_)*sizeof(double));
    2684         memcpy(xsolution,solution_,(numberRows_+numberColumns_)*sizeof(double));
    2685 #endif
    2686         updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
    2687                           0.0,objectiveChange,true);
    2688 #if 0
    2689         int i;
    2690         for (i=0;i<numberRows_+numberColumns_;i++) {
    2691           if (cost_[i]!=xcost[i])
    2692             printf("** %d old cost %g new %g sol %g\n",
    2693                    i,xcost[i],cost_[i],solution_[i]);
    2694           if (lower_[i]!=xlower[i])
    2695             printf("** %d old lower %g new %g sol %g\n",
    2696                    i,xlower[i],lower_[i],solution_[i]);
    2697           if (upper_[i]!=xupper[i])
    2698             printf("** %d old upper %g new %g sol %g\n",
    2699                    i,xupper[i],upper_[i],solution_[i]);
    2700           if (dj_[i]!=xdj[i])
    2701             printf("** %d old dj %g new %g sol %g\n",
    2702                    i,xdj[i],dj_[i],solution_[i]);
    2703           if (solution_[i]!=xsolution[i])
    2704             printf("** %d old solution %g new %g sol %g\n",
    2705                    i,xsolution[i],solution_[i],solution_[i]);
    2706         }
    2707         delete [] xcost;
    2708         delete [] xupper;
    2709         delete [] xlower;
    2710         delete [] xdj;
    2711         delete [] xsolution;
    2712 #endif
    2713         // for now - recompute all
    2714         gutsOfSolution(NULL,NULL);
    2715         //assert(numberDualInfeasibilitiesWithoutFree_==0);
    2716 
    2717         if (numberDualInfeasibilities_||situationChanged==2)
    2718           problemStatus_=-1; // carry on as normal
    2719         situationChanged=0;
    2720       } else {
    2721         // iterate
    2722         if (cleanDuals!=2)
    2723           problemStatus_=-1;
    2724         else
    2725           problemStatus_ = 10; // try primal
    2726       }
    2727     }
    2728   }
    2729   if (type==0||type==1) {
    2730     if (!type) {
    2731       // create save arrays
    2732       delete [] saveStatus_;
    2733       delete [] savedSolution_;
    2734       saveStatus_ = new unsigned char [numberRows_+numberColumns_];
    2735       savedSolution_ = new double [numberRows_+numberColumns_];
    2736     }
    2737     // save arrays
    2738     memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
    2739     memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
    2740            numberRows_*sizeof(double));
    2741     memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
    2742   }
    2743 
    2744   // restore weights (if saved) - also recompute infeasibility list
    2745   if (tentativeStatus>-3)
    2746     dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
    2747   else
    2748     dualRowPivot_->saveWeights(this,3);
    2749   // unflag all variables (we may want to wait a bit?)
    2750   if (tentativeStatus!=-2&&unflagVariables) {
    2751     int iRow;
    2752     for (iRow=0;iRow<numberRows_;iRow++) {
    2753       int iPivot=pivotVariable_[iRow];
    2754       clearFlagged(iPivot);
    2755     }
    2756   }
    2757   // see if cutoff reached
    2758   double limit = 0.0;
    2759   getDblParam(ClpDualObjectiveLimit, limit);
    2760   if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
    2761            optimizationDirection_*limit&&
    2762            !numberAtFakeBound()&&!numberDualInfeasibilities_) {
    2763     problemStatus_=1;
    2764     secondaryStatus_ = 1; // and say was on cutoff
    2765   }
    2766   if (problemStatus_<0&&!changeMade_) {
    2767     problemStatus_=4; // unknown
    2768   }
    2769   lastGoodIteration_ = numberIterations_;
    2770 #if 0
    2771   double thisObj = progress->lastObjective(0);
    2772   double lastObj = progress->lastObjective(1);
    2773   if (lastObj>thisObj+1.0e-6*max(fabs(thisObj),fabs(lastObj))+1.0e-8
    2774       &&givenDuals==NULL) {
    2775     int maxFactor = factorization_->maximumPivots();
    2776     if (maxFactor>10) {
    2777       if (forceFactorization_<0)
    2778         forceFactorization_= maxFactor;
    2779       forceFactorization_ = max (1,(forceFactorization_>>1));
    2780       printf("Reducing factorization frequency\n");
    2781     }
    2782   }
    2783 #endif
    2784 }
    2785 /* While updateDualsInDual sees what effect is of flip
    2786    this does actual flipping.
    2787    If change >0.0 then value in array >0.0 => from lower to upper
    2788 */
    2789 void
    2790 ClpPdco::flipBounds(CoinIndexedVector * rowArray,
    2791                   CoinIndexedVector * columnArray,
    2792                   double change)
    2793 {
    2794   int number;
    2795   int * which;
    2796  
    2797   int iSection;
    2798 
    2799   for (iSection=0;iSection<2;iSection++) {
    2800     int i;
    2801     double * solution = solutionRegion(iSection);
    2802     double * lower = lowerRegion(iSection);
    2803     double * upper = upperRegion(iSection);
    2804     int addSequence;
    2805     if (!iSection) {
    2806       number = rowArray->getNumElements();
    2807       which = rowArray->getIndices();
    2808       addSequence = numberColumns_;
    2809     } else {
    2810       number = columnArray->getNumElements();
    2811       which = columnArray->getIndices();
    2812       addSequence = 0;
    2813     }
    2814    
    2815     for (i=0;i<number;i++) {
    2816       int iSequence = which[i];
    2817       Status status = getStatus(iSequence+addSequence);
    2818 
    2819       switch(status) {
    2820 
    2821       case basic:
    2822       case isFree:
    2823       case superBasic:
    2824       case ClpSimplex::isFixed:
    2825         break;
    2826       case atUpperBound:
    2827         // to lower bound
    2828         setStatus(iSequence+addSequence,atLowerBound);
    2829         solution[iSequence] = lower[iSequence];
    2830         break;
    2831       case atLowerBound:
    2832         // to upper bound
    2833         setStatus(iSequence+addSequence,atUpperBound);
    2834         solution[iSequence] = upper[iSequence];
    2835         break;
    2836       }
    2837     }
    2838   }
    2839   rowArray->setNumElements(0);
    2840   columnArray->setNumElements(0);
    2841 }
    2842 // Restores bound to original bound
    2843 void
    2844 ClpPdco::originalBound( int iSequence)
    2845 {
    2846   if (getFakeBound(iSequence)!=noFake)
    2847     numberFake_--;;
    2848   if (iSequence>=numberColumns_) {
    2849     // rows
    2850     int iRow = iSequence-numberColumns_;
    2851     rowLowerWork_[iRow]=rowLower_[iRow];
    2852     rowUpperWork_[iRow]=rowUpper_[iRow];
    2853     if (rowScale_) {
    2854       if (rowLowerWork_[iRow]>-1.0e50)
    2855         rowLowerWork_[iRow] *= rowScale_[iRow];
    2856       if (rowUpperWork_[iRow]<1.0e50)
    2857         rowUpperWork_[iRow] *= rowScale_[iRow];
    2858     }
    2859   } else {
    2860     // columns
    2861     columnLowerWork_[iSequence]=columnLower_[iSequence];
    2862     columnUpperWork_[iSequence]=columnUpper_[iSequence];
    2863     if (rowScale_) {
    2864       double multiplier = 1.0/columnScale_[iSequence];
    2865       if (columnLowerWork_[iSequence]>-1.0e50)
    2866         columnLowerWork_[iSequence] *= multiplier;
    2867       if (columnUpperWork_[iSequence]<1.0e50)
    2868         columnUpperWork_[iSequence] *= multiplier;
    2869     }
    2870   }
    2871   setFakeBound(iSequence,noFake);
    2872 }
    2873 /* As changeBounds but just changes new bounds for a single variable.
    2874    Returns true if change */
    2875 bool
    2876 ClpPdco::changeBound( int iSequence)
    2877 {
    2878   // old values
    2879   double oldLower=lower_[iSequence];
    2880   double oldUpper=upper_[iSequence];
    2881   double value=solution_[iSequence];
    2882   bool modified=false;
    2883   originalBound(iSequence);
    2884   // original values
    2885   double lowerValue=lower_[iSequence];
    2886   double upperValue=upper_[iSequence];
    2887   // back to altered values
    2888   lower_[iSequence] = oldLower;
    2889   upper_[iSequence] = oldUpper;
    2890   if (getFakeBound(iSequence)!=noFake)
    2891     numberFake_--;;
    2892   if (value==oldLower) {
    2893     if (upperValue > oldLower + dualBound_) {
    2894       upper_[iSequence]=oldLower+dualBound_;
    2895       setFakeBound(iSequence,upperFake);
    2896       modified=true;
    2897       numberFake_++;
    2898     }
    2899   } else if (value==oldUpper) {
    2900     if (lowerValue < oldUpper - dualBound_) {
    2901       lower_[iSequence]=oldUpper-dualBound_;
    2902       setFakeBound(iSequence,lowerFake);
    2903       modified=true;
    2904       numberFake_++;
    2905     }
    2906   } else {
    2907     assert(value==oldLower||value==oldUpper);
    2908   }
    2909   return modified;
    2910 }
    2911 // Perturbs problem
    2912 void
    2913 ClpPdco::perturb()
    2914 {
    2915   if (perturbation_>100)
    2916     return; //perturbed already
    2917   if (perturbation_==100)
    2918     perturbation_=50; // treat as normal
    2919   bool modifyRowCosts=false;
    2920   // dual perturbation
    2921   double perturbation=1.0e-20;
    2922   // maximum fraction of cost to perturb
    2923   double maximumFraction = 1.0e-5;
    2924   double constantPerturbation = 100.0*dualTolerance_;
    2925   const int * lengths = matrix_->getVectorLengths();
    2926   int maxLength=0;
    2927   int minLength=numberRows_;
    2928   if (!numberIterations_&&perturbation_==50) {
    2929     // See if we need to perturb
    2930     double * sort = new double[numberColumns_];
    2931     int i;
    2932     for (i=0;i<numberColumns_;i++) {
    2933       double value = fabs(objectiveWork_[i]);
    2934       sort[i]=value;
    2935     }
    2936     std::sort(sort,sort+numberColumns_);
    2937     int number=1;
    2938     double last = sort[0];
    2939     for (i=1;i<numberColumns_;i++) {
    2940       if (last!=sort[i])
    2941         number++;
    2942       last=sort[i];
    2943     }
    2944     delete [] sort;
    2945     if (number*4>numberColumns_)
    2946       return; // good enough
    2947   }
    2948   int iColumn;
    2949   for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    2950     if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
    2951       int length = lengths[iColumn];
    2952       if (length>2) {
    2953         maxLength = max(maxLength,length);
    2954         minLength = min(minLength,length);
    2955       }
    2956     }
    2957   }
    2958   // If > 70 then do rows
    2959   if (perturbation_>=70) {
    2960     modifyRowCosts=true;
    2961     perturbation_ -= 20;
    2962     printf("Row costs modified, ");
    2963   }
    2964   bool uniformChange=false;
    2965   if (perturbation_>50) {
    2966     // Experiment
    2967     // maximumFraction could be 1.0e-10 to 1.0
    2968     double m[]={1.0e-10,1.0e-9,1.0e-8,1.0e-7,1.0e-6,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,1.0};
    2969     maximumFraction = m[min(perturbation_-51,10)];
    2970   }
    2971   int iRow;
    2972   double smallestNonZero=1.0e100;
    2973   if (perturbation_>=50) {
    2974     perturbation = 1.0e-8;
    2975     bool allSame=true;
    2976     double lastValue=0.0;
    2977     for (iRow=0;iRow<numberRows_;iRow++) {
    2978       double lo = rowLowerWork_[iRow];
    2979       double up = rowUpperWork_[iRow];
    2980       if (lo<up) {
    2981         double value = fabs(rowObjectiveWork_[iRow]);
    2982         perturbation = max(perturbation,value);
    2983         if (value) {
    2984           modifyRowCosts=true;
    2985           smallestNonZero = min(smallestNonZero,value);
    2986         }
    2987       }
    2988       if (lo&&lo>-1.0e10) {
    2989         lo=fabs(lo);
    2990         if (!lastValue)
    2991           lastValue=lo;
    2992         else if (fabs(lo-lastValue)>1.0e-7)
    2993           allSame=false;
    2994       }
    2995       if (up&&up<1.0e10) {
    2996         up=fabs(up);
    2997         if (!lastValue)
    2998           lastValue=up;
    2999         else if (fabs(up-lastValue)>1.0e-7)
    3000           allSame=false;
    3001       }
    3002     }
    3003     for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    3004       double lo = columnLowerWork_[iColumn];
    3005       double up = columnUpperWork_[iColumn];
    3006       if (lo<up) {
    3007         double value =
    3008           fabs(objectiveWork_[iColumn]);
    3009         perturbation = max(perturbation,value);
    3010         if (value) {
    3011           smallestNonZero = min(smallestNonZero,value);
    3012         }
    3013       }
    3014       if (lo&&lo>-1.0e10) {
    3015         lo=fabs(lo);
    3016         if (!lastValue)
    3017           lastValue=lo;
    3018         else if (fabs(lo-lastValue)>1.0e-7)
    3019           allSame=false;
    3020       }
    3021       if (up&&up<1.0e10) {
    3022         up=fabs(up);
    3023         if (!lastValue)
    3024           lastValue=up;
    3025         else if (fabs(up-lastValue)>1.0e-7)
    3026           allSame=false;
    3027       }
    3028     }
    3029     if (allSame) {
    3030       // Really hit perturbation
    3031       maximumFraction = max(1.0e-2*lastValue,maximumFraction);
    3032     }
    3033     perturbation = min(perturbation,smallestNonZero/maximumFraction);
    3034   } else {
    3035     // user is in charge
    3036     maximumFraction = 1.0e-1;
    3037     // but some experiments
    3038     if (perturbation_<=-900) {
    3039       modifyRowCosts=true;
    3040       perturbation_ += 1000;
    3041       printf("Row costs modified, ");
    3042     }
    3043     if (perturbation_<=-10) {
    3044       perturbation_ += 10;
    3045       maximumFraction = 1.0;
    3046       if ((-perturbation_)%100>=10) {
    3047         uniformChange=true;
    3048         perturbation_+=20;
    3049       }
    3050       while (perturbation_<-10) {
    3051         perturbation_ += 100;
    3052         maximumFraction *= 1.0e-1;
    3053       }
    3054     }
    3055     perturbation = pow(10.0,perturbation_);
    3056   }
    3057   double largestZero=0.0;
    3058   double largest=0.0;
    3059   double largestPerCent=0.0;
    3060   // modify costs
    3061   bool printOut=(handler_->logLevel()==63);
    3062   printOut=false;
    3063   if (modifyRowCosts) {
    3064     for (iRow=0;iRow<numberRows_;iRow++) {
    3065       if (rowLowerWork_[iRow]<rowUpperWork_[iRow]) {
    3066         double value = perturbation;
    3067         double currentValue = rowObjectiveWork_[iRow];
    3068         value = min(value,maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-3));
    3069         if (rowLowerWork_[iRow]>-largeValue_) {
    3070           if (fabs(rowLowerWork_[iRow])<fabs(rowUpperWork_[iRow]))
    3071             value *= CoinDrand48();
    3072           else
    3073             value *= -CoinDrand48();
    3074         } else if (rowUpperWork_[iRow]<largeValue_) {
    3075           value *= -CoinDrand48();
    3076         } else {
    3077           value=0.0;
    3078         }
    3079         if (currentValue) {
    3080           largest = max(largest,fabs(value));
    3081           if (fabs(value)>fabs(currentValue)*largestPerCent)
    3082             largestPerCent=fabs(value/currentValue);
    3083         } else {
    3084           largestZero=max(largestZero,fabs(value));
    3085         }
    3086         if (printOut)
    3087           printf("row %d cost %g change %g\n",iRow,rowObjectiveWork_[iRow],value);
    3088         rowObjectiveWork_[iRow] += value;
    3089       }
    3090     }
    3091   }
    3092   // more its but faster double weight[]={1.0e-4,1.0e-2,1.0e-1,1.0,2.0,10.0,100.0,200.0,400.0,600.0,1000.0};
    3093   // good its double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,5.0,10.0,20.0,30.0,40.0,100.0};
    3094   double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,5.0,10.0,20.0,30.0,40.0,100.0};
    3095   //double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,20.0,50.0,100.0,120.0,130.0,140.0,200.0};
    3096   double extraWeight=10.0;
    3097   // Scale back if wanted
    3098   double weight2[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0};
    3099   if (constantPerturbation<99.0*dualTolerance_) {
    3100     perturbation *= 0.1;
    3101     extraWeight=0.5;
    3102     memcpy(weight,weight2,sizeof(weight2));
    3103   }
    3104   // adjust weights if all columns long
    3105   double factor=1.0;
    3106   if (maxLength) {
    3107     factor = 3.0/(double) minLength;
    3108   }
    3109   // Make variables with more elements more expensive
    3110   const double m1 = 0.5;
    3111   for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    3112     if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]&&getStatus(iColumn)!=basic) {
    3113       double value = perturbation;
    3114       double currentValue = objectiveWork_[iColumn];
    3115       value = min(value,constantPerturbation+maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-8));
    3116       //value = min(value,constantPerturbation;+maximumFraction*fabs(currentValue));
    3117       double value2 = constantPerturbation+1.0e-1*smallestNonZero;
    3118       if (uniformChange) {
    3119         value = maximumFraction;
    3120         value2=maximumFraction;
    3121       }
    3122       if (columnLowerWork_[iColumn]>-largeValue_) {
    3123         if (fabs(columnLowerWork_[iColumn])<
    3124             fabs(columnUpperWork_[iColumn])) {
    3125           value *= (1.0-m1 +m1*CoinDrand48());
    3126           value2 *= (1.0-m1 +m1*CoinDrand48());
    3127         } else {
    3128           value *= -(1.0-m1+m1*CoinDrand48());
    3129           value2 *= -(1.0-m1+m1*CoinDrand48());
    3130         }
    3131       } else if (columnUpperWork_[iColumn]<largeValue_) {
    3132         value *= -(1.0-m1+m1*CoinDrand48());
    3133         value2 *= -(1.0-m1+m1*CoinDrand48());
    3134       } else {
    3135         value=0.0;
    3136       }
    3137       if (value) {
    3138         int length = lengths[iColumn];
    3139         if (length>3) {
    3140           length = (int) ((double) length * factor);
    3141           length = max(3,length);
    3142         }
    3143         double multiplier;
    3144 #if 1
    3145         if (length<10)
    3146           multiplier=weight[length];
    3147         else
    3148           multiplier = weight[10];
    3149 #else
    3150         if (length<10)
    3151           multiplier=weight[length];
    3152         else
    3153           multiplier = weight[10]+extraWeight*(length-10);
    3154         multiplier *= 0.5;
    3155 #endif
    3156         value *= multiplier;
    3157         value = min (value,value2);
    3158         if (fabs(value)<=dualTolerance_)
    3159           value=0.0;
    3160         if (currentValue) {
    3161           largest = max(largest,fabs(value));
    3162           if (fabs(value)>fabs(currentValue)*largestPerCent)
    3163             largestPerCent=fabs(value/currentValue);
    3164         } else {
    3165           largestZero=max(largestZero,fabs(value));
    3166         }
    3167         if (printOut)
    3168           printf("col %d cost %g change %g\n",iColumn,objectiveWork_[iColumn],value);
    3169         objectiveWork_[iColumn] += value;
    3170       }
    3171     }
    3172   }
    3173   handler_->message(CLP_SIMPLEX_PERTURB,messages_)
    3174     <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
    3175     <<CoinMessageEol;
    3176   // and zero changes
    3177   //int nTotal = numberRows_+numberColumns_;
    3178   //memset(cost_+nTotal,0,nTotal*sizeof(double));
    3179   // say perturbed
    3180   perturbation_=101;
    3181 
    3182 }
    3183 /* For strong branching.  On input lower and upper are new bounds
    3184    while on output they are change in objective function values
    3185    (>1.0e50 infeasible).
    3186    Return code is 0 if nothing interesting, -1 if infeasible both
    3187    ways and +1 if infeasible one way (check values to see which one(s))
    3188 */
    3189 int ClpPdco::strongBranching(int numberVariables,const int * variables,
    3190                                     double * newLower, double * newUpper,
    3191                                     double ** outputSolution,
    3192                                     int * outputStatus, int * outputIterations,
    3193                                     bool stopOnFirstInfeasible,
    3194                                     bool alwaysFinish)
    3195 {
    3196   int i;
    3197   int returnCode=0;
    3198   double saveObjectiveValue = objectiveValue_;
    3199 #if 1
    3200   algorithm_ = -1;
    3201 
    3202   //scaling(false);
    3203 
    3204   // put in standard form (and make row copy)
    3205   // create modifiable copies of model rim and do optional scaling
    3206   createRim(7+8+16,true);
    3207 
    3208   // change newLower and newUpper if scaled
    3209 
    3210   // Do initial factorization
    3211   // and set certain stuff
    3212   // We can either set increasing rows so ...IsBasic gives pivot row
    3213   // or we can just increment iBasic one by one
    3214   // for now let ...iBasic give pivot row
    3215   factorization_->increasingRows(2);
    3216   // row activities have negative sign
    3217   factorization_->slackValue(-1.0);
    3218   factorization_->zeroTolerance(1.0e-13);
    3219 
    3220   int factorizationStatus = internalFactorize(0);
    3221   if (factorizationStatus<0)
    3222     return 1; // some error
    3223   else if (factorizationStatus)
    3224     handler_->message(CLP_SINGULARITIES,messages_)
    3225       <<factorizationStatus
    3226       <<CoinMessageEol;
    3227  
    3228   // save stuff
    3229   ClpFactorization saveFactorization(*factorization_);
    3230   // save basis and solution
    3231   double * saveSolution = new double[numberRows_+numberColumns_];
    3232   memcpy(saveSolution,solution_,
    3233          (numberRows_+numberColumns_)*sizeof(double));
    3234   unsigned char * saveStatus =
    3235     new unsigned char [numberRows_+numberColumns_];
    3236   memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
    3237   // save bounds as createRim makes clean copies
    3238   double * saveLower = new double[numberRows_+numberColumns_];
    3239   memcpy(saveLower,lower_,
    3240          (numberRows_+numberColumns_)*sizeof(double));
    3241   double * saveUpper = new double[numberRows_+numberColumns_];
    3242   memcpy(saveUpper,upper_,
    3243          (numberRows_+numberColumns_)*sizeof(double));
    3244   double * saveObjective = new double[numberRows_+numberColumns_];
    3245   memcpy(saveObjective,cost_,
    3246          (numberRows_+numberColumns_)*sizeof(double));
    3247   int * savePivot = new int [numberRows_];
    3248   memcpy(savePivot, pivotVariable_, numberRows_*sizeof(int));
    3249   // need to save/restore weights.
    3250 
    3251   int iSolution = 0;
    3252   for (i=0;i<numberVariables;i++) {
    3253     int iColumn = variables[i];
    3254     double objectiveChange;
    3255     double saveBound;
    3256    
    3257     // try down
    3258    
    3259     saveBound = columnUpper_[iColumn];
    3260     // external view - in case really getting optimal
    3261     columnUpper_[iColumn] = newUpper[i];
    3262     if (scalingFlag_<=0)
    3263       upper_[iColumn] = newUpper[i];
    3264     else
    3265       upper_[iColumn] = newUpper[i]/columnScale_[iColumn]; // scale
    3266     // Start of fast iterations
    3267     int status = fastDual(alwaysFinish);
    3268     if (status) {
    3269       // not finished - might be optimal
    3270       checkPrimalSolution(rowActivityWork_,columnActivityWork_);
    3271       double limit = 0.0;
    3272       getDblParam(ClpDualObjectiveLimit, limit);
    3273       if (!numberPrimalInfeasibilities_&&objectiveValue()*optimizationDirection_<
    3274           optimizationDirection_*limit) {
    3275         problemStatus_=0;
    3276       }
    3277       status=problemStatus_;
    3278     }
    3279 
    3280     if (scalingFlag_<=0) {
    3281       memcpy(outputSolution[iSolution],solution_,numberColumns_*sizeof(double));
    3282     } else {
    3283       int j;
    3284       double * sol = outputSolution[iSolution];
    3285       for (j=0;j<numberColumns_;j++)
    3286         sol[j] = solution_[j]*columnScale_[j];
    3287     }
    3288     outputStatus[iSolution]=status;
    3289     outputIterations[iSolution]=numberIterations_;
    3290     iSolution++;
    3291     // restore
    3292     memcpy(solution_,saveSolution,
    3293            (numberRows_+numberColumns_)*sizeof(double));
    3294     memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
    3295     memcpy(lower_,saveLower,
    3296            (numberRows_+numberColumns_)*sizeof(double));
    3297     memcpy(upper_,saveUpper,
    3298            (numberRows_+numberColumns_)*sizeof(double));
    3299     memcpy(cost_,saveObjective,
    3300          (numberRows_+numberColumns_)*sizeof(double));
    3301     columnUpper_[iColumn] = saveBound;
    3302     memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int));
    3303     delete factorization_;
    3304     factorization_ = new ClpFactorization(saveFactorization);
    3305 
    3306     if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
    3307       objectiveChange = objectiveValue_-saveObjectiveValue;
    3308     } else {
    3309       objectiveChange = 1.0e100;
    3310     }
    3311     newUpper[i]=objectiveChange;
    3312 #ifdef CLP_DEBUG
    3313     printf("down on %d costs %g\n",iColumn,objectiveChange);
    3314 #endif
    3315          
    3316     // try up
    3317    
    3318     saveBound = columnLower_[iColumn];
    3319     // external view - in case really getting optimal
    3320     columnLower_[iColumn] = newLower[i];
    3321     if (scalingFlag_<=0)
    3322       lower_[iColumn] = newLower[i];
    3323     else
    3324       lower_[iColumn] = newLower[i]/columnScale_[iColumn]; // scale
    3325     // Start of fast iterations
    3326     status = fastDual(alwaysFinish);
    3327     if (status) {
    3328       // not finished - might be optimal
    3329       checkPrimalSolution(rowActivityWork_,columnActivityWork_);
    3330       double limit = 0.0;
    3331       getDblParam(ClpDualObjectiveLimit, limit);
    3332       if (!numberPrimalInfeasibilities_&&objectiveValue()*optimizationDirection_<
    3333           optimizationDirection_*limit) {
    3334         problemStatus_=0;
    3335       }
    3336       status=problemStatus_;
    3337     }
    3338     if (scalingFlag_<=0) {
    3339       memcpy(outputSolution[iSolution],solution_,numberColumns_*sizeof(double));
    3340     } else {
    3341       int j;
    3342       double * sol = outputSolution[iSolution];
    3343       for (j=0;j<numberColumns_;j++)
    3344         sol[j] = solution_[j]*columnScale_[j];
    3345     }
    3346     outputStatus[iSolution]=status;
    3347     outputIterations[iSolution]=numberIterations_;
    3348     iSolution++;
    3349 
    3350     // restore
    3351     memcpy(solution_,saveSolution,
    3352            (numberRows_+numberColumns_)*sizeof(double));
    3353     memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
    3354     memcpy(lower_,saveLower,
    3355            (numberRows_+numberColumns_)*sizeof(double));
    3356     memcpy(upper_,saveUpper,
    3357            (numberRows_+numberColumns_)*sizeof(double));
    3358     memcpy(cost_,saveObjective,
    3359          (numberRows_+numberColumns_)*sizeof(double));
    3360     columnLower_[iColumn] = saveBound;
    3361     memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int));
    3362     delete factorization_;
    3363     factorization_ = new ClpFactorization(saveFactorization);
    3364 
    3365     if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
    3366       objectiveChange = objectiveValue_-saveObjectiveValue;
    3367     } else {
    3368       objectiveChange = 1.0e100;
    3369     }
    3370     newLower[i]=objectiveChange;
    3371 #ifdef CLP_DEBUG
    3372     printf("up on %d costs %g\n",iColumn,objectiveChange);
    3373 #endif
    3374          
    3375     /* Possibilities are:
    3376        Both sides feasible - store
    3377        Neither side feasible - set objective high and exit if desired
    3378        One side feasible - change bounds and resolve
    3379     */
    3380     if (newUpper[i]<1.0e100) {
    3381       if(newLower[i]<1.0e100) {
    3382         // feasible - no action
    3383       } else {
    3384         // up feasible, down infeasible
    3385         returnCode=1;
    3386         if (stopOnFirstInfeasible)
    3387           break;
    3388       }
    3389     } else {
    3390       if(newLower[i]<1.0e100) {
    3391         // down feasible, up infeasible
    3392         returnCode=1;
    3393         if (stopOnFirstInfeasible)
    3394           break;
    3395       } else {
    3396         // neither side feasible
    3397         returnCode=-1;
    3398         break;
    3399       }
    3400     }
    3401   }
    3402   delete [] saveSolution;
    3403   delete [] saveLower;
    3404   delete [] saveUpper;
    3405   delete [] saveObjective;
    3406   delete [] saveStatus;
    3407   delete [] savePivot;
    3408 
    3409   // Get rid of some arrays and empty factorization
    3410   deleteRim();
    3411 #else
    3412   // save basis and solution
    3413   double * rowSolution = new double[numberRows_];
    3414   memcpy(rowSolution,rowActivity_,numberRows_*sizeof(double));
    3415   double * columnSolution = new double[numberColumns_];
    3416   memcpy(columnSolution,columnActivity_,numberColumns_*sizeof(double));
    3417   unsigned char * saveStatus =
    3418     new unsigned char [numberRows_+numberColumns_];
    3419   if (!status_) {
    3420     // odd but anyway setup all slack basis
    3421     status_ = new unsigned char [numberColumns_+numberRows_];
    3422     memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
    3423   }
    3424   memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
    3425   int iSolution =0;
    3426   for (i=0;i<numberVariables;i++) {
    3427     int iColumn = variables[i];
    3428     double objectiveChange;
    3429    
    3430     // try down
    3431    
    3432     double saveUpper = columnUpper_[iColumn];
    3433     columnUpper_[iColumn] = newUpper[i];
    3434     int status=dual(0);
    3435     memcpy(outputSolution[iSolution],columnActivity_,numberColumns_*sizeof(double));
    3436     outputStatus[iSolution]=status;
    3437     outputIterations[iSolution]=numberIterations_;
    3438     iSolution++;
    3439 
    3440     // restore
    3441     columnUpper_[iColumn] = saveUpper;
    3442     memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
    3443     memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
    3444     memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
    3445 
    3446     if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
    3447       objectiveChange = objectiveValue_-saveObjectiveValue;
    3448     } else {
    3449       objectiveChange = 1.0e100;
    3450     }
    3451     newUpper[i]=objectiveChange;
    3452 #ifdef CLP_DEBUG
    3453     printf("down on %d costs %g\n",iColumn,objectiveChange);
    3454 #endif
    3455          
    3456     // try up
    3457    
    3458     double saveLower = columnLower_[iColumn];
    3459     columnLower_[iColumn] = newLower[i];
    3460     status=dual(0);
    3461     memcpy(outputSolution[iSolution],columnActivity_,numberColumns_*sizeof(double));
    3462     outputStatus[iSolution]=status;
    3463     outputIterations[iSolution]=numberIterations_;
    3464     iSolution++;
    3465 
    3466     // restore
    3467     columnLower_[iColumn] = saveLower;
    3468     memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
    3469     memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
    3470     memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
    3471 
    3472     if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
    3473       objectiveChange = objectiveValue_-saveObjectiveValue;
    3474     } else {
    3475       objectiveChange = 1.0e100;
    3476     }
    3477     newLower[i]=objectiveChange;
    3478 #ifdef CLP_DEBUG
    3479     printf("up on %d costs %g\n",iColumn,objectiveChange);
    3480 #endif
    3481          
    3482     /* Possibilities are:
    3483        Both sides feasible - store
    3484        Neither side feasible - set objective high and exit
    3485        One side feasible - change bounds and resolve
    3486     */
    3487     if (newUpper[i]<1.0e100) {
    3488       if(newLower[i]<1.0e100) {
    3489         // feasible - no action
    3490       } else {
    3491         // up feasible, down infeasible
    3492         returnCode=1;
    3493         if (stopOnFirstInfeasible)
    3494           break;
    3495       }
    3496     } else {
    3497       if(newLower[i]<1.0e100) {
    3498         // down feasible, up infeasible
    3499         returnCode=1;
    3500         if (stopOnFirstInfeasible)
    3501           break;
    3502       } else {
    3503         // neither side feasible
    3504         returnCode=-1;
    3505         break;
    3506       }
    3507     }
    3508   }
    3509   delete [] rowSolution;
    3510   delete [] columnSolution;
    3511   delete [] saveStatus;
    3512 #endif
    3513   objectiveValue_ = saveObjectiveValue;
    3514   return returnCode;
    3515 }
    3516 // treat no pivot as finished (unless interesting)
    3517 int ClpPdco::fastDual(bool alwaysFinish)
    3518 {
    3519   algorithm_ = -1;
    3520   // save data
    3521   ClpDataSave data = saveData();
    3522   dualTolerance_=dblParam_[ClpDualTolerance];
    3523   primalTolerance_=dblParam_[ClpPrimalTolerance];
    3524 
    3525   // save dual bound
    3526   double saveDualBound = dualBound_;
    3527 
    3528   double objectiveChange;
    3529   // for dual we will change bounds using dualBound_
    3530   // for this we need clean basis so it is after factorize
    3531   gutsOfSolution(NULL,NULL);
    3532   numberFake_ =0; // Number of variables at fake bounds
    3533   changeBounds(true,NULL,objectiveChange);
    3534 
    3535   problemStatus_ = -1;
    3536   numberIterations_=0;
    3537   factorization_->sparseThreshold(0);
    3538   factorization_->goSparse();
    3539 
    3540   int lastCleaned=0; // last time objective or bounds cleaned up
    3541 
    3542   // number of times we have declared optimality
    3543   numberTimesOptimal_=0;
    3544 
    3545   // This says whether to restore things etc
    3546   int factorType=0;
    3547   /*
    3548     Status of problem:
    3549     0 - optimal
    3550     1 - infeasible
    3551     2 - unbounded
    3552     -1 - iterating
    3553     -2 - factorization wanted
    3554     -3 - redo checking without factorization
    3555     -4 - looks infeasible
    3556 
    3557     BUT also from whileIterating return code is:
    3558 
    3559    -1 iterations etc
    3560    -2 inaccuracy
    3561    -3 slight inaccuracy (and done iterations)
    3562    +0 looks optimal (might be unbounded - but we will investigate)
    3563    +1 looks infeasible
    3564    +3 max iterations
    3565 
    3566   */
    3567 
    3568   int returnCode = 0;
    3569 
    3570   int iRow,iColumn;
    3571   while (problemStatus_<0) {
    3572     // clear
    3573     for (iRow=0;iRow<4;iRow++) {
    3574       rowArray_[iRow]->clear();
    3575     }   
    3576    
    3577     for (iColumn=0;iColumn<2;iColumn++) {
    3578       columnArray_[iColumn]->clear();
    3579     }   
    3580 
    3581     // give matrix (and model costs and bounds a chance to be
    3582     // refreshed (normally null)
    3583     matrix_->refresh(this);
    3584     // may factorize, checks if problem finished
    3585     // should be able to speed this up on first time
    3586     statusOfProblemInDual(lastCleaned,factorType,progress_,NULL);
    3587 
    3588     // Say good factorization
    3589     factorType=1;
    3590 
    3591     // Do iterations
    3592     if (problemStatus_<0) {
    3593       double * givenPi=NULL;
    3594       returnCode = whileIterating(givenPi);
    3595       if (!alwaysFinish&&(returnCode<1||returnCode==3)) {
    3596         double limit = 0.0;
    3597         getDblParam(ClpDualObjectiveLimit, limit);
    3598         if(fabs(limit)>1.0e30||objectiveValue()*optimizationDirection_<
    3599            optimizationDirection_*limit||
    3600            numberAtFakeBound()) {
    3601           returnCode=1;
    3602           secondaryStatus_ = 1; // and say was on cutoff
    3603           // can't say anything interesting - might as well return
    3604 #ifdef CLP_DEBUG
    3605           printf("returning from fastDual after %d iterations with code %d\n",
    3606                  numberIterations_,returnCode);
    3607 #endif
    3608           break;
    3609         }
    3610       }
    3611       returnCode=0;
    3612     }
    3613   }
    3614 
    3615   // clear
    3616   for (iRow=0;iRow<4;iRow++) {
    3617     rowArray_[iRow]->clear();
    3618   }   
    3619  
    3620   for (iColumn=0;iColumn<2;iColumn++) {
    3621     columnArray_[iColumn]->clear();
    3622   }   
    3623   assert(!numberFake_||returnCode||problemStatus_); // all bounds should be okay
    3624   // Restore any saved stuff
    3625   restoreData(data);
    3626   dualBound_ = saveDualBound;
    3627   return returnCode;
    3628 }
    3629 /* Checks number of variables at fake bounds.  This is used by fastDual
    3630    so can exit gracefully before end */
    3631 int
    3632 ClpPdco::numberAtFakeBound()
    3633 {
    3634   int iSequence;
    3635   int numberFake=0;
    3636  
    3637   for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
    3638     FakeBound bound = getFakeBound(iSequence);
    3639     switch(getStatus(iSequence)) {
    3640 
    3641     case basic:
    3642       break;
    3643     case isFree:
    3644     case superBasic:
    3645     case ClpSimplex::isFixed:
    3646       assert (bound==noFake);
    3647       break;
    3648     case atUpperBound:
    3649       if (bound==upperFake||bound==bothFake)
    3650         numberFake++;
    3651       break;
    3652     case atLowerBound:
    3653       if (bound==lowerFake||bound==bothFake)
    3654         numberFake++;
    3655       break;
    3656     }
    3657   }
    3658   numberFake_ = numberFake;
    3659   return numberFake;
    3660 }
    3661 /* Pivot out a variable and choose an incoing one.  Assumes dual
    3662    feasible - will not go through a reduced cost. 
    3663    Returns step length in theta
    3664    Returns ray in ray_ (or NULL if no pivot)
    3665    Return codes as before but -1 means no acceptable pivot
    3666 */
    3667 int
    3668 ClpPdco::pivotResult()
    3669 {
    3670   abort();
    3671   return 0;
    3672 }
    3673 /*
    3674    Row array has row part of pivot row
    3675    Column array has column part.
    3676    This is used in dual values pass
    3677 */
    3678 int
    3679 ClpPdco::checkPossibleValuesMove(CoinIndexedVector * rowArray,
    3680                                         CoinIndexedVector * columnArray,
    3681                                         double acceptablePivot,
    3682                                         CoinBigIndex * dubiousWeights)
    3683 {
    3684   double * work;
    3685   int number;
    3686   int * which;
    3687   int iSection;
    3688 
    3689   double tolerance = dualTolerance_*1.001;
    3690 
    3691   double thetaDown = 1.0e31;
    3692   double changeDown ;
    3693   double thetaUp = 1.0e31;
    3694   double bestAlphaDown = acceptablePivot*0.99999;
    3695   double bestAlphaUp = acceptablePivot*0.99999;
    3696   int sequenceDown =-1;
    3697   int sequenceUp = sequenceOut_;
    3698 
    3699   double djBasic = dj_[sequenceOut_];
    3700   if (djBasic>0.0) {
    3701     // basic at lower bound so directionOut_ 1 and -1 in pivot row
    3702     // dj will go to zero on other way
    3703     thetaUp = djBasic;
    3704     changeDown = -lower_[sequenceOut_];
    3705   } else {
    3706     // basic at upper bound so directionOut_ -1 and 1 in pivot row
    3707     // dj will go to zero on other way
    3708     thetaUp = -djBasic;
    3709     changeDown = upper_[sequenceOut_];
    3710   }
    3711   bestAlphaUp = 1.0;
    3712   int addSequence;
    3713 
    3714   double alphaUp=0.0;
    3715   double alphaDown=0.0;
    3716 
    3717   for (iSection=0;iSection<2;iSection++) {
    3718 
    3719     int i;
    3720     if (!iSection) {
    3721       work = rowArray->denseVector();
    3722       number = rowArray->getNumElements();
    3723       which = rowArray->getIndices();
    3724       addSequence = numberColumns_;
    3725     } else {
    3726       work = columnArray->denseVector();
    3727       number = columnArray->getNumElements();
    3728       which = columnArray->getIndices();
    3729       addSequence = 0;
    3730     }
    3731    
    3732     for (i=0;i<number;i++) {
    3733       int iSequence = which[i];
    3734       int iSequence2 = iSequence + addSequence;
    3735       double alpha;
    3736       double oldValue;
    3737       double value;
    3738 
    3739       switch(getStatus(iSequence2)) {
    3740          
    3741       case basic:
    3742         break;
    3743       case ClpSimplex::isFixed:
    3744         alpha = work[i];
    3745         changeDown += alpha*upper_[iSequence2];
    3746         break;
    3747       case isFree:
    3748       case superBasic:
    3749         alpha = work[i];
    3750         // dj must be effectively zero as dual feasible
    3751         if (fabs(alpha)>bestAlphaUp) {
    3752           thetaDown = 0.0;
    3753           thetaUp = 0.0;
    3754           bestAlphaDown = fabs(alpha);
    3755           bestAlphaUp = bestAlphaUp;
    3756           sequenceDown =iSequence2;
    3757           sequenceUp = sequenceDown;
    3758           alphaUp = alpha;
    3759           alphaDown = alpha;
    3760         }
    3761         break;
    3762       case atUpperBound:
    3763         alpha = work[i];
    3764         oldValue = dj_[iSequence2];
    3765         changeDown += alpha*upper_[iSequence2];
    3766         if (alpha>=acceptablePivot) {
    3767           // might do other way
    3768           value = oldValue+thetaUp*alpha;
    3769           if (value>-tolerance) {
    3770             if (value>tolerance||fabs(alpha)>bestAlphaUp) {
    3771               thetaUp = -oldValue/alpha;
    3772               bestAlphaUp = fabs(alpha);
    3773               sequenceUp = iSequence2;
    3774               alphaUp = alpha;
    3775             }
    3776           }
    3777         } else if (alpha<=-acceptablePivot) {
    3778           // might do this way
    3779           value = oldValue-thetaDown*alpha;
    3780           if (value>-tolerance) {
    3781             if (value>tolerance||fabs(alpha)>bestAlphaDown) {
    3782               thetaDown = oldValue/alpha;
    3783               bestAlphaDown = fabs(alpha);
    3784               sequenceDown = iSequence2;
    3785               alphaDown = alpha;
    3786             }
    3787           }
    3788         }
    3789         break;
    3790       case atLowerBound:
    3791         alpha = work[i];
    3792         oldValue = dj_[iSequence2];
    3793         changeDown += alpha*lower_[iSequence2];
    3794         if (alpha<=-acceptablePivot) {
    3795           // might do other way
    3796           value = oldValue+thetaUp*alpha;
    3797           if (value<tolerance) {
    3798             if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
    3799               thetaUp = -oldValue/alpha;
    3800               bestAlphaUp = fabs(alpha);
    3801               sequenceUp = iSequence2;
    3802               alphaUp = alpha;
    3803             }
    3804           }
    3805         } else if (alpha>=acceptablePivot) {
    3806           // might do this way
    3807           value = oldValue-thetaDown*alpha;
    3808           if (value<tolerance) {
    3809             if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
    3810               thetaDown = oldValue/alpha;
    3811               bestAlphaDown = fabs(alpha);
    3812               sequenceDown = iSequence2;
    3813               alphaDown = alpha;
    3814             }
    3815           }
    3816         }
    3817         break;
    3818       }
    3819     }
    3820   }
    3821   int returnCode;
    3822   thetaUp *= -1.0;
    3823   double changeUp = -thetaUp * changeDown;
    3824   changeDown = -thetaDown*changeDown;
    3825   if (max(fabs(thetaDown),fabs(thetaUp))<1.0e-8) {
    3826     // largest
    3827     if (fabs(alphaDown)<fabs(alphaUp)) {
    3828       sequenceDown =-1;
    3829     }
    3830   }
    3831   // choose
    3832   if (changeDown>changeUp&&sequenceDown>=0) {
    3833     theta_ = thetaDown;
    3834     sequenceIn_ = sequenceDown;
    3835     alpha_ = alphaDown;
    3836 #ifdef CLP_DEBUG
    3837     if ((handler_->logLevel()&32))
    3838       printf("predicted way - dirout %d, change %g,%g theta %g\n",
    3839              directionOut_,changeDown,changeUp,theta_);
    3840 #endif
    3841     returnCode = 0;
    3842   } else {
    3843     theta_ = thetaUp;
    3844     sequenceIn_ = sequenceUp;
    3845     alpha_ = alphaUp;
    3846     if (sequenceIn_!=sequenceOut_) {
    3847 #ifdef CLP_DEBUG
    3848       if ((handler_->logLevel()&32))
    3849         printf("opposite way - dirout %d, change %g,%g theta %g\n",
    3850                directionOut_,changeDown,changeUp,theta_);
    3851 #endif
    3852       returnCode = 0;
    3853     } else {
    3854 #ifdef CLP_DEBUG
    3855       if ((handler_->logLevel()&32))
    3856         printf("opposite way to zero dj - dirout %d, change %g,%g theta %g\n",
    3857                directionOut_,changeDown,changeUp,theta_);
    3858 #endif
    3859       returnCode = 1;
    3860     }
    3861   }
    3862   return returnCode;
    3863 }
    3864 /*
    3865    This sees if we can move duals in dual values pass.
    3866    This is done before any pivoting
    3867 */
    3868 void ClpPdco::doEasyOnesInValuesPass(double * dj)
    3869 {
    3870   // Get column copy
    3871   CoinPackedMatrix * columnCopy = matrix();
    3872   // Get a row copy in standard format
    3873   CoinPackedMatrix copy;
    3874   copy.reverseOrderedCopyOf(*columnCopy);
    3875   // get matrix data pointers
    3876   const int * column = copy.getIndices();
    3877   const CoinBigIndex * rowStart = copy.getVectorStarts();
    3878   const int * rowLength = copy.getVectorLengths();
    3879   const double * elementByRow = copy.getElements();
    3880   double tolerance = dualTolerance_*1.001;
    3881 
    3882   int iRow;
    3883 #ifdef CLP_DEBUG
    3884   {
    3885     double value5=0.0;
    3886     int i;
    3887     for (i=0;i<numberRows_+numberColumns_;i++) {
    3888       if (dj[i]<-1.0e-6)
    3889         value5 += dj[i]*upper_[i];
    3890       else if (dj[i] >1.0e-6)
    3891         value5 += dj[i]*lower_[i];
    3892     }
    3893     printf("Values objective Value before %g\n",value5);
    3894   }
    3895 #endif
    3896   // for scaled row
    3897   double * scaled=NULL;
    3898   if (rowScale_)
    3899     scaled = new double[numberColumns_];
    3900   for (iRow=0;iRow<numberRows_;iRow++) {
    3901 
    3902     int iSequence = iRow + numberColumns_;
    3903     double djBasic = dj[iSequence];
    3904     if (getRowStatus(iRow)==basic&&fabs(djBasic)>tolerance) {
    3905 
    3906       double changeUp ;
    3907       // always -1 in pivot row
    3908       if (djBasic>0.0) {
    3909         // basic at lower bound
    3910         changeUp = -lower_[iSequence];
    3911       } else {
    3912         // basic at upper bound
    3913         changeUp = upper_[iSequence];
    3914       }
    3915       bool canMove=true;
    3916       int i;
    3917       const double * thisElements = elementByRow + rowStart[iRow];
    3918       const int * thisIndices = column+rowStart[iRow];
    3919       if (rowScale_) {
    3920         // scale row
    3921         double scale = rowScale_[iRow];
    3922         for (i = 0;i<rowLength[iRow];i++) {
    3923           int iColumn = thisIndices[i];
    3924           double alpha = thisElements[i];
    3925           scaled[i] = scale*alpha*columnScale_[iColumn];
    3926         }
    3927         thisElements = scaled;
    3928       }
    3929       for (i = 0;i<rowLength[iRow];i++) {
    3930         int iColumn = thisIndices[i];
    3931         double alpha = thisElements[i];
    3932         double oldValue = dj[iColumn];;
    3933         double value;
    3934 
    3935         switch(getStatus(iColumn)) {
    3936          
    3937         case basic:
    3938           if (dj[iColumn]<-tolerance&&
    3939               fabs(solution_[iColumn]-upper_[iColumn])<1.0e-8) {
    3940             // at ub
    3941             changeUp += alpha*upper_[iColumn];
    3942             // might do other way
    3943             value = oldValue+djBasic*alpha;
    3944             if (value>tolerance)
    3945               canMove=false;
    3946           } else if (dj[iColumn]>tolerance&&
    3947               fabs(solution_[iColumn]-lower_[iColumn])<1.0e-8) {
    3948             changeUp += alpha*lower_[iColumn];
    3949             // might do other way
    3950             value = oldValue+djBasic*alpha;
    3951             if (value<-tolerance)
    3952               canMove=false;
    3953           } else {
    3954             canMove=false;
    3955           }
    3956           break;
    3957         case ClpSimplex::isFixed:
    3958           changeUp += alpha*upper_[iColumn];
    3959           break;
    3960         case isFree:
    3961         case superBasic:
    3962           canMove=false;
    3963         break;
    3964       case atUpperBound:
    3965         changeUp += alpha*upper_[iColumn];
    3966         // might do other way
    3967         value = oldValue+djBasic*alpha;
    3968         if (value>tolerance)
    3969           canMove=false;
    3970         break;
    3971       case atLowerBound:
    3972         changeUp += alpha*lower_[iColumn];
    3973         // might do other way
    3974         value = oldValue+djBasic*alpha;
    3975         if (value<-tolerance)
    3976           canMove=false;
    3977         break;
    3978         }
    3979       }
    3980       if (canMove) {
    3981         if (changeUp*djBasic>1.0e-12||fabs(changeUp)<1.0e-8) {
    3982           // move
    3983           for (i = 0;i<rowLength[iRow];i++) {
    3984             int iColumn = thisIndices[i];
    3985             double alpha = thisElements[i];
    3986             dj[iColumn] += djBasic * alpha;
    3987           }
    3988           dj[iSequence]=0.0;
    3989 #ifdef CLP_DEBUG
    3990           {
    3991             double value5=0.0;
    3992             int i;
    3993             for (i=0;i<numberRows_+numberColumns_;i++) {
    3994               if (dj[i]<-1.0e-6)
    3995                 value5 += dj[i]*upper_[i];
    3996               else if (dj[i] >1.0e-6)
    3997                 value5 += dj[i]*lower_[i];
    3998             }
    3999             printf("Values objective Value after row %d old dj %g %g\n",
    4000                    iRow,djBasic,value5);
    4001           }
    4002 #endif
    4003         }
    4004       }
    4005     }     
    4006   }
    4007   delete [] scaled;
    4008 }
    4009 int
    4010 ClpPdco::nextSuperBasic()
    4011 {
    4012   if (firstFree_>=0) {
    4013     int returnValue=firstFree_;
    4014     int iColumn=firstFree_+1;
    4015     for (;iColumn<numberRows_+numberColumns_;iColumn++) {
    4016       if (getStatus(iColumn)==isFree)
    4017         if (fabs(dj_[iColumn])>1.0e2*dualTolerance_)
    4018           break;
    4019     }
    4020     firstFree_ = iColumn;
    4021     if (firstFree_==numberRows_+numberColumns_)
    4022       firstFree_=-1;
    4023     return returnValue;
    4024   } else {
    4025     return -1;
    4026   }
    4027 }
  • branches/pre/Makefile.Clp

    r213 r214  
    11# Static or shared libraries should be built (STATIC or SHARED)?
    22LibType := SHARED
    3 LibType := STATIC
     3#LibType := STATIC
    44
    55# Select optimization (-O or -g). -O will be automatically bumped up to the
     
    77# between then specify the exact level you want, e.g., -O1 or -O2
    88OptLevel := -g
    9 #OptLevel := -O2
     9OptLevel := -O1
    1010
    1111LIBNAME := Clp
     
    3535LIBSRC += ClpSolve.cpp
    3636LIBSRC += ClpInterior.cpp
     37LIBSRC += ClpPdco.cpp
    3738# and Presolve stuff
    3839LIBSRC += ClpPresolve.cpp
  • branches/pre/Test/Makefile.test

    r208 r214  
    6363endif
    6464#LDFLAGS += -lefence
    65 LDFLAGS += -static
     65#LDFLAGS += -static
    6666CXXFLAGS += -pg
    6767#if DENSE and using ATLAS
  • branches/pre/include/ClpInterior.hpp

    r213 r214  
    2222class ClpNonLinearCost;
    2323class ClpInteriorProgress;
    24 
     24// ******** DATA to be moved into protected section of ClpInterior
     25typedef struct{
     26  double  atolmin;
     27  double  r3norm;
     28  double  LSdamp;
     29  double* deltay;
     30} Info;
     31
     32typedef struct{
     33  double  atolold;
     34  double  atolnew;
     35  double  r3ratio;
     36  int   istop;
     37  int   itncg;
     38} Outfo;
     39 
     40typedef struct{
     41double  gamma;
     42double  delta;
     43int MaxIter;
     44double  FeaTol;
     45double  OptTol;
     46double  StepTol;
     47double  x0min;
     48double  z0min;
     49double  mu0;
     50int   LSmethod;   // 1=Cholesky    2=QR    3=LSQR
     51int   LSproblem;  // See below
     52int LSQRMaxIter;
     53double  LSQRatol1; // Initial  atol
     54double  LSQRatol2; // Smallest atol (unless atol1 is smaller)
     55double  LSQRconlim;
     56int  wait;
     57} Options;
     58class Lsqr;
     59// ***** END
    2560/** This solves LPs using interior point methods
    2661
     
    107142  /** Pdco algorithm - see ClpPdco.hpp for method */
    108143  int pdco();
    109   //@}
    110 
    111   /**@name Needed for functionality of OsiSimplexInterface */
    112   //@{
    113   /// LSQR
    114   void lsqr();
     144  // ** Temporary version
     145  int  pdco( Lsqr *lsqr, Options &options, Info &info, Outfo &outfo);
    115146  //@}
    116147
  • branches/pre/include/ClpPdco.hpp

    r212 r214  
    1111#define ClpPdco_H
    1212
    13 #include "ClpSimplex.hpp"
     13#include "ClpInterior.hpp"
    1414
    15 /** This solves LPs using the dual simplex method
     15/** This solves problems in Primal Dual Convex Optimization
    1616
    17     It inherits from ClpSimplex.  It has no data of its own and
    18     is never created - only cast from a ClpSimplex object at algorithm time.
     17    It inherits from ClpInterior.  It has no data of its own and
     18    is never created - only cast from a ClpInterior object at algorithm time.
    1919
    2020*/
    2121
    22 class ClpPdco : public ClpSimplex {
     22class ClpPdco : public ClpInterior {
    2323
    2424public:
     
    2626  /**@name Description of algorithm */
    2727  //@{
    28   /** Dual algorithm
     28  /** Pdco algorithm
    2929
    3030      Method
    3131
    32      It tries to be a single phase approach with a weight of 1.0 being
    33      given to getting optimal and a weight of updatedDualBound_ being
    34      given to getting dual feasible.  In this version I have used the
    35      idea that this weight can be thought of as a fake bound.  If the
    36      distance between the lower and upper bounds on a variable is less
    37      than the feasibility weight then we are always better off flipping
    38      to other bound to make dual feasible.  If the distance is greater
    39      then we make up a fake bound updatedDualBound_ away from one bound.
    40      If we end up optimal or primal infeasible, we check to see if
    41      bounds okay.  If so we have finished, if not we increase updatedDualBound_
    42      and continue (after checking if unbounded). I am undecided about
    43      free variables - there is coding but I am not sure about it.  At
    44      present I put them in basis anyway.
    45 
    46      The code is designed to take advantage of sparsity so arrays are
    47      seldom zeroed out from scratch or gone over in their entirety.
    48      The only exception is a full scan to find outgoing variable for
    49      Dantzig row choice.  For steepest edge we keep an updated list
    50      of infeasibilities (actually squares). 
    51      On easy problems we don't need full scan - just
    52      pick first reasonable.
    53 
    54      One problem is how to tackle degeneracy and accuracy.  At present
    55      I am using the modification of costs which I put in OSL and some
    56      of what I think is the dual analog of Gill et al.
    57      I am still not sure of the exact details.
    58 
    59      The flow of dual is three while loops as follows:
    60 
    61      while (not finished) {
    62 
    63        while (not clean solution) {
    64 
    65           Factorize and/or clean up solution by flipping variables so
    66           dual feasible.  If looks finished check fake dual bounds.
    67           Repeat until status is iterating (-1) or finished (0,1,2)
    68 
    69        }
    70 
    71        while (status==-1) {
    72 
    73          Iterate until no pivot in or out or time to re-factorize.
    74 
    75          Flow is:
    76 
    77          choose pivot row (outgoing variable).  if none then
    78          we are primal feasible so looks as if done but we need to
    79          break and check bounds etc.
    80 
    81          Get pivot row in tableau
    82 
    83          Choose incoming column.  If we don't find one then we look
    84          primal infeasible so break and check bounds etc.  (Also the
    85          pivot tolerance is larger after any iterations so that may be
    86          reason)
    87 
    88          If we do find incoming column, we may have to adjust costs to
    89          keep going forwards (anti-degeneracy).  Check pivot will be stable
    90          and if unstable throw away iteration and break to re-factorize.
    91          If minor error re-factorize after iteration.
    92 
    93          Update everything (this may involve flipping variables to stay
    94          dual feasible.
    95 
    96        }
    97 
    98      }
    99 
    100      TODO's (or maybe not)
    101 
    102      At present we never check we are going forwards.  I overdid that in
    103      OSL so will try and make a last resort.
    104 
    105      Needs partial scan pivot out option.
    106 
    107      May need other anti-degeneracy measures, especially if we try and use
    108      loose tolerances as a way to solve in fewer iterations.
    109 
    110      I like idea of dynamic scaling.  This gives opportunity to decouple
    111      different implications of scaling for accuracy, iteration count and
    112      feasibility tolerance.
    11332
    11433  */
    11534
    116   int dual(int ifValuesPass);
    117   /** For strong branching.  On input lower and upper are new bounds
    118       while on output they are change in objective function values
    119       (>1.0e50 infeasible).
    120       Return code is 0 if nothing interesting, -1 if infeasible both
    121       ways and +1 if infeasible one way (check values to see which one(s))
    122       Solutions are filled in as well - even down, odd up - also
    123       status and number of iterations
    124   */
    125   int strongBranching(int numberVariables,const int * variables,
    126                       double * newLower, double * newUpper,
    127                       double ** outputSolution,
    128                       int * outputStatus, int * outputIterations,
    129                       bool stopOnFirstInfeasible=true,
    130                       bool alwaysFinish=false);
     35  int pdco();
     36  // ** Temporary version
     37  int  pdco( Lsqr *lsqr, Options &options, Info &info, Outfo &outfo);
     38
    13139  //@}
    13240
    13341  /**@name Functions used in dual */
    13442  //@{
    135   /** This has the flow between re-factorizations
    136       Broken out for clarity and will be used by strong branching
    137 
    138       Reasons to come out:
    139       -1 iterations etc
    140       -2 inaccuracy
    141       -3 slight inaccuracy (and done iterations)
    142       +0 looks optimal (might be unbounded - but we will investigate)
    143       +1 looks infeasible
    144       +3 max iterations
    145 
    146       If givenPi not NULL then in values pass
    147    */
    148   int whileIterating(double * & givenPi);
    149   /** The duals are updated by the given arrays.
    150       Returns number of infeasibilities.
    151       After rowArray and columnArray will just have those which
    152       have been flipped.
    153       Variables may be flipped between bounds to stay dual feasible.
    154       The output vector has movement of primal
    155       solution (row length array) */
    156   int updateDualsInDual(CoinIndexedVector * rowArray,
    157                   CoinIndexedVector * columnArray,
    158                   CoinIndexedVector * outputArray,
    159                   double theta,
    160                   double & objectiveChange,
    161                         bool fullRecompute);
    162   /** The duals are updated by the given arrays.
    163       This is in values pass - so no changes to primal is made
    164   */
    165   void updateDualsInValuesPass(CoinIndexedVector * rowArray,
    166                   CoinIndexedVector * columnArray,
    167                   double theta);
    168   /** While updateDualsInDual sees what effect is of flip
    169       this does actuall flipping.
    170       If change >0.0 then value in array >0.0 => from lower to upper
    171   */
    172   void flipBounds(CoinIndexedVector * rowArray,
    173                   CoinIndexedVector * columnArray,
    174                   double change);
    175   /**
    176       Row array has row part of pivot row
    177       Column array has column part.
    178       This chooses pivot column.
    179       Spare arrays are used to save pivots which will go infeasible
    180       We will check for basic so spare array will never overflow.
    181       If necessary will modify costs
    182       For speed, we may need to go to a bucket approach when many
    183       variables are being flipped
    184 
    185   */
    186   void dualColumn(CoinIndexedVector * rowArray,
    187                   CoinIndexedVector * columnArray,
    188                   CoinIndexedVector * spareArray,
    189                   CoinIndexedVector * spareArray2,
    190                   double accpetablePivot,
    191                   CoinBigIndex * dubiousWeights);
    192   /**
    193       Row array has row part of pivot row
    194       Column array has column part.
    195       This sees what is best thing to do in dual values pass
    196       Returns 0 if theta_ move will put basic variable out to bound,
    197       1 if can change dual on chosen row and leave variable in basis
    198   */
    199   int checkPossibleValuesMove(CoinIndexedVector * rowArray,
    200                                CoinIndexedVector * columnArray,
    201                               double acceptablePivot,
    202                               CoinBigIndex * dubiousWeights);
    203   /**
    204       This sees if we can move duals in dual values pass.
    205       This is done before any pivoting
    206   */
    207   void doEasyOnesInValuesPass(double * givenReducedCosts);
    208   /**
    209       Chooses dual pivot row
    210       Would be faster with separate region to scan
    211       and will have this (with square of infeasibility) when steepest
    212       For easy problems we can just choose one of the first rows we look at
    213      
    214       If alreadyChosen >=0 then in values pass and that row has been
    215       selected
    216   */
    217   void dualRow(int alreadyChosen);
    218   /** Checks if any fake bounds active - if so returns number and modifies
    219       updatedDualBound_ and everything.
    220       Free variables will be left as free
    221       Returns number of bounds changed if >=0
    222       Returns -1 if not initialize and no effect
    223       Fills in changeVector which can be used to see if unbounded
    224       and cost of change vector
    225   */
    226   int changeBounds(bool initialize,CoinIndexedVector * outputArray,
    227                    double & changeCost);
    228   /** As changeBounds but just changes new bounds for a single variable.
    229       Returns true if change */
    230   bool changeBound( int iSequence);
    231   /// Restores bound to original bound
    232   void originalBound(int iSequence);
    233   /** Checks if tentative optimal actually means unbounded in dual
    234       Returns -3 if not, 2 if is unbounded */
    235   int checkUnbounded(CoinIndexedVector * ray,CoinIndexedVector * spare,
    236                      double changeCost);
    237   /**  Refactorizes if necessary
    238        Checks if finished.  Updates status.
    239        lastCleaned refers to iteration at which some objective/feasibility
    240        cleaning too place.
    241 
    242        type - 0 initial so set up save arrays etc
    243             - 1 normal -if good update save
    244             - 2 restoring from saved
    245   */
    246   void statusOfProblemInDual(int & lastCleaned, int type,
    247                              ClpSimplexProgress * progress,
    248                              double * givenDjs);
    249   /// Perturbs problem (method depends on perturbation())
    250   void perturb();
    251   /** Fast iterations.  Misses out a lot of initialization.
    252       Normally stops on maximum iterations, first re-factorization
    253       or tentative optimum.  If looks interesting then continues as
    254       normal.  Returns 0 if finished properly, 1 otherwise.
    255   */
    256   int fastDual(bool alwaysFinish=false);
    257   /** Checks number of variables at fake bounds.  This is used by fastDual
    258       so can exit gracefully before end */
    259   int numberAtFakeBound();
    260 
    261   /** Pivot in a variable and choose an outgoing one.  Assumes dual
    262       feasible - will not go through a reduced cost.  Returns step length in theta
    263       Returns ray in ray_ (or NULL if no pivot)
    264       Return codes as before but -1 means no acceptable pivot
    265   */
    266   int pivotResult();
    267   /** Get next free , -1 if none */
    268   int nextSuperBasic();
    269  
     43  /// LSQR
     44  void lsqr();
    27045  //@}
    27146};
Note: See TracChangeset for help on using the changeset viewer.