Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (4 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpPdco.cpp

    r1972 r2385  
    1010
    1111*/
    12 
    13 
    1412
    1513#include "CoinPragma.hpp"
     
    3129#include <iostream>
    3230
    33 int
    34 ClpPdco::pdco()
    35 {
    36      printf("Dummy pdco solve\n");
    37      return 0;
     31int ClpPdco::pdco()
     32{
     33  printf("Dummy pdco solve\n");
     34  return 0;
    3835}
    3936// ** Temporary version
    40 int
    41 ClpPdco::pdco( ClpPdcoBase * stuff, Options &options, Info &info, Outfo &outfo)
    42 {
    43 //    D1, D2 are positive-definite diagonal matrices defined from d1, d2.
    44 //           In particular, d2 indicates the accuracy required for
    45 //           satisfying each row of Ax = b.
    46 //
    47 // D1 and D2 (via d1 and d2) provide primal and dual regularization
    48 // respectively.  They ensure that the primal and dual solutions
    49 // (x,r) and (y,z) are unique and bounded.
    50 //
    51 // A scalar d1 is equivalent to d1 = ones(n,1), D1 = diag(d1).
    52 // A scalar d2 is equivalent to d2 = ones(m,1), D2 = diag(d2).
    53 // Typically, d1 = d2 = 1e-4.
    54 // These values perturb phi(x) only slightly  (by about 1e-8) and request
    55 // that A*x = b be satisfied quite accurately (to about 1e-8).
    56 // Set d1 = 1e-4, d2 = 1 for least-squares problems with bound constraints.
    57 // The problem is then
    58 //
    59 //    minimize    phi(x) + 1/2 norm(d1*x)^2 + 1/2 norm(A*x - b)^2
    60 //    subject to  bl <= x <= bu.
    61 //
    62 // More generally, d1 and d2 may be n and m vectors containing any positive
    63 // values (preferably not too small, and typically no larger than 1).
    64 // Bigger elements of d1 and d2 improve the stability of the solver.
    65 //
    66 // At an optimal solution, if x(j) is on its lower or upper bound,
    67 // the corresponding z(j) is positive or negative respectively.
    68 // If x(j) is between its bounds, z(j) = 0.
    69 // If bl(j) = bu(j), x(j) is fixed at that value and z(j) may have
    70 // either sign.
    71 //
    72 // Also, r and y satisfy r = D2 y, so that Ax + D2^2 y = b.
    73 // Thus if d2(i) = 1e-4, the i-th row of Ax = b will be satisfied to
    74 // approximately 1e-8.  This determines how large d2(i) can safely be.
    75 //
    76 //
    77 // EXTERNAL FUNCTIONS:
    78 // options         = pdcoSet;                  provided with pdco.m
    79 // [obj,grad,hess] = pdObj( x );               provided by user
    80 //               y = pdMat( name,mode,m,n,x ); provided by user if pdMat
    81 //                                             is a string, not a matrix
    82 //
    83 // INPUT ARGUMENTS:
    84 // pdObj      is a string containing the name of a function pdObj.m
    85 //            or a function_handle for such a function
    86 //            such that [obj,grad,hess] = pdObj(x) defines
    87 //            obj  = phi(x)              : a scalar,
    88 //            grad = gradient of phi(x)  : an n-vector,
    89 //            hess = diag(Hessian of phi): an n-vector.
    90 //         Examples:
    91 //            If phi(x) is the linear function c"x, pdObj should return
    92 //               [obj,grad,hess] = [c"*x, c, zeros(n,1)].
    93 //            If phi(x) is the entropy function E(x) = sum x(j) log x(j),
    94 //               [obj,grad,hess] = [E(x), log(x)+1, 1./x].
    95 // pdMat      may be an ifexplicit m x n matrix A (preferably sparse!),
    96 //            or a string containing the name of a function pdMat.m
    97 //            or a function_handle for such a function
    98 //            such that y = pdMat( name,mode,m,n,x )
    99 //            returns   y = A*x (mode=1)  or  y = A"*x (mode=2).
    100 //            The input parameter "name" will be the string pdMat.
    101 // b          is an m-vector.
    102 // bl         is an n-vector of lower bounds.  Non-existent bounds
    103 //            may be represented by bl(j) = -Inf or bl(j) <= -1e+20.
    104 // bu         is an n-vector of upper bounds.  Non-existent bounds
    105 //            may be represented by bu(j) =  Inf or bu(j) >=  1e+20.
    106 // d1, d2     may be positive scalars or positive vectors (see above).
    107 // options    is a structure that may be set and altered by pdcoSet
    108 //            (type help pdcoSet).
    109 // x0, y0, z0 provide an initial solution.
    110 // xsize, zsize are estimates of the biggest x and z at the solution.
    111 //            They are used to scale (x,y,z).  Good estimates
    112 //            should improve the performance of the barrier method.
    113 //
    114 //
    115 // OUTPUT ARGUMENTS:
    116 // x          is the primal solution.
    117 // y          is the dual solution associated with Ax + D2 r = b.
    118 // z          is the dual solution associated with bl <= x <= bu.
    119 // inform = 0 if a solution is found;
    120 //        = 1 if too many iterations were required;
    121 //        = 2 if the linesearch failed too often.
    122 // PDitns     is the number of Primal-Dual Barrier iterations required.
    123 // CGitns     is the number of Conjugate-Gradient  iterations required
    124 //            if an iterative solver is used (LSQR).
    125 // time       is the cpu time used.
    126 //----------------------------------------------------------------------
    127 
    128 // PRIVATE FUNCTIONS:
    129 //    pdxxxbounds
    130 //    pdxxxdistrib
    131 //    pdxxxlsqr
    132 //    pdxxxlsqrmat
    133 //    pdxxxmat
    134 //    pdxxxmerit
    135 //    pdxxxresid1
    136 //    pdxxxresid2
    137 //    pdxxxstep
    138 //
    139 // GLOBAL VARIABLES:
    140 //    global pdDDD1 pdDDD2 pdDDD3
    141 //
    142 //
    143 // NOTES:
    144 // The matrix A should be reasonably well scaled: norm(A,inf) =~ 1.
    145 // The vector b and objective phi(x) may be of any size, but ensure that
    146 // xsize and zsize are reasonably close to norm(x,inf) and norm(z,inf)
    147 // at the solution.
    148 //
    149 // The files defining pdObj  and pdMat
    150 // must not be called Fname.m or Aname.m!!
    151 //
    152 //
    153 // AUTHOR:
    154 //    Michael Saunders, Systems Optimization Laboratory (SOL),
    155 //    Stanford University, Stanford, California, USA.
    156 //    saunders@stanford.edu
    157 //
    158 // CONTRIBUTORS:
    159 //    Byunggyoo Kim, SOL, Stanford University.
    160 //    bgkim@stanford.edu
    161 //
    162 // DEVELOPMENT:
    163 // 20 Jun 1997: Original version of pdsco.m derived from pdlp0.m.
    164 // 29 Sep 2002: Original version of pdco.m  derived from pdsco.m.
    165 //              Introduced D1, D2 in place of gamma*I, delta*I
    166 //              and allowed for general bounds bl <= x <= bu.
    167 // 06 Oct 2002: Allowed for fixed variabes: bl(j) = bu(j) for any j.
    168 // 15 Oct 2002: Eliminated some work vectors (since m, n might be LARGE).
    169 //              Modularized residuals, linesearch
    170 // 16 Oct 2002: pdxxx..., pdDDD... names rationalized.
    171 //              pdAAA eliminated (global copy of A).
    172 //              Aname is now used directly as an ifexplicit A or a function.
    173 //              NOTE: If Aname is a function, it now has an extra parameter.
    174 // 23 Oct 2002: Fname and Aname can now be function handles.
    175 // 01 Nov 2002: Bug fixed in feval in pdxxxmat.
    176 //-----------------------------------------------------------------------
    177 
    178 //  global pdDDD1 pdDDD2 pdDDD3
    179      double inf = 1.0e30;
    180      double eps = 1.0e-15;
    181      double atolold = -1.0, r3ratio = -1.0, Pinf, Dinf, Cinf, Cinf0;
    182 
    183      printf("\n   --------------------------------------------------------");
    184      printf("\n   pdco.m                            Version of 01 Nov 2002");
    185      printf("\n   Primal-dual barrier method to minimize a convex function");
    186      printf("\n   subject to linear constraints Ax + r = b,  bl <= x <= bu");
    187      printf("\n   --------------------------------------------------------\n");
    188 
    189      int m = numberRows_;
    190      int n = numberColumns_;
    191      bool ifexplicit = true;
    192 
    193      CoinDenseVector<double> b(m, rhs_);
    194      CoinDenseVector<double> x(n, x_);
    195      CoinDenseVector<double> y(m, y_);
    196      CoinDenseVector<double> z(n, dj_);
    197      //delete old arrays
    198      delete [] rhs_;
    199      delete [] x_;
    200      delete [] y_;
    201      delete [] dj_;
    202      rhs_ = NULL;
    203      x_ = NULL;
    204      y_ = NULL;
    205      dj_ = NULL;
    206 
    207      // Save stuff so available elsewhere
    208      pdcoStuff_ = stuff;
    209 
    210      double normb  = b.infNorm();
    211      double normx0 = x.infNorm();
    212      double normy0 = y.infNorm();
    213      double normz0 = z.infNorm();
    214 
    215      printf("\nmax |b | = %8g     max |x0| = %8g", normb , normx0);
    216      printf(                "      xsize   = %8g", xsize_);
    217      printf("\nmax |y0| = %8g     max |z0| = %8g", normy0, normz0);
    218      printf(                "      zsize   = %8g", zsize_);
    219 
    220      //---------------------------------------------------------------------
    221      // Initialize.
    222      //---------------------------------------------------------------------
    223      //true   = 1;
    224      //false  = 0;
    225      //zn     = zeros(n,1);
    226      //int nb     = n + m;
    227      int CGitns = 0;
    228      int inform = 0;
    229      //---------------------------------------------------------------------
    230      //  Only allow scalar d1, d2 for now
    231      //---------------------------------------------------------------------
    232      /*
     37int ClpPdco::pdco(ClpPdcoBase *stuff, 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  double inf = 1.0e30;
     176  double eps = 1.0e-15;
     177  double atolold = -1.0, r3ratio = -1.0, 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 = numberRows_;
     186  int n = numberColumns_;
     187  bool ifexplicit = true;
     188
     189  CoinDenseVector< double > b(m, rhs_);
     190  CoinDenseVector< double > x(n, x_);
     191  CoinDenseVector< double > y(m, y_);
     192  CoinDenseVector< double > z(n, dj_);
     193  //delete old arrays
     194  delete[] rhs_;
     195  delete[] x_;
     196  delete[] y_;
     197  delete[] dj_;
     198  rhs_ = NULL;
     199  x_ = NULL;
     200  y_ = NULL;
     201  dj_ = NULL;
     202
     203  // Save stuff so available elsewhere
     204  pdcoStuff_ = stuff;
     205
     206  double normb = b.infNorm();
     207  double normx0 = x.infNorm();
     208  double normy0 = y.infNorm();
     209  double normz0 = z.infNorm();
     210
     211  printf("\nmax |b | = %8g     max |x0| = %8g", normb, normx0);
     212  printf("      xsize   = %8g", xsize_);
     213  printf("\nmax |y0| = %8g     max |z0| = %8g", normy0, normz0);
     214  printf("      zsize   = %8g", zsize_);
     215
     216  //---------------------------------------------------------------------
     217  // Initialize.
     218  //---------------------------------------------------------------------
     219  //true   = 1;
     220  //false  = 0;
     221  //zn     = zeros(n,1);
     222  //int nb     = n + m;
     223  int CGitns = 0;
     224  int inform = 0;
     225  //---------------------------------------------------------------------
     226  //  Only allow scalar d1, d2 for now
     227  //---------------------------------------------------------------------
     228  /*
    233229     if (d1_->size()==1)
    234230         d1_->resize(n, d1_->getElements()[0]);  // Allow scalar d1, d2
     
    236232         d2->resize(m, d2->getElements()[0]);  // to mean dk * unit vector
    237233      */
    238      assert (stuff->sizeD1() == 1);
    239      double d1 = stuff->getD1();
    240      double d2 = stuff->getD2();
    241 
    242      //---------------------------------------------------------------------
    243      // Grab input options.
    244      //---------------------------------------------------------------------
    245      int  maxitn    = options.MaxIter;
    246      double featol    = options.FeaTol;
    247      double opttol    = options.OptTol;
    248      double steptol  = options.StepTol;
    249      int  stepSame  = 1; /* options.StepSame;   // 1 means stepx == stepz */
    250      double x0min    = options.x0min;
    251      double z0min    = options.z0min;
    252      double mu0      = options.mu0;
    253      int  LSproblem = options.LSproblem; // See below
    254      int  LSmethod  = options.LSmethod;  // 1=Cholesky    2=QR    3=LSQR
    255      int  itnlim    = options.LSQRMaxIter * CoinMin(m, n);
    256      double atol1     = options.LSQRatol1; // Initial  atol
    257      double atol2     = options.LSQRatol2; // Smallest atol,unless atol1 is smaller
    258      double conlim    = options.LSQRconlim;
    259      //int  wait      = options.wait;
    260 
    261      // LSproblem:
    262      //  1 = dy          2 = dy shifted, DLS
    263      // 11 = s          12 =  s shifted, DLS    (dx = Ds)
    264      // 21 = dx
    265      // 31 = 3x3 system, symmetrized by Z^{1/2}
    266      // 32 = 2x2 system, symmetrized by X^{1/2}
    267 
    268      //---------------------------------------------------------------------
    269      // Set other parameters.
    270      //---------------------------------------------------------------------
    271      int  kminor    = 0;      // 1 stops after each iteration
    272      double eta       = 1e-4;  // Linesearch tolerance for "sufficient descent"
    273      double maxf      = 10;    // Linesearch backtrack limit (function evaluations)
    274      double maxfail   = 1;      // Linesearch failure limit (consecutive iterations)
    275      double bigcenter = 1e+3;  // mu is reduced if center < bigcenter.
    276 
    277      // Parameters for LSQR.
    278      double atolmin   = eps;    // Smallest atol if linesearch back-tracks
    279      double btol      = 0;      // Should be small (zero is ok)
    280      double show      = false; // Controls lsqr iteration log
    281      /*
     234  assert(stuff->sizeD1() == 1);
     235  double d1 = stuff->getD1();
     236  double d2 = stuff->getD2();
     237
     238  //---------------------------------------------------------------------
     239  // Grab input options.
     240  //---------------------------------------------------------------------
     241  int maxitn = options.MaxIter;
     242  double featol = options.FeaTol;
     243  double opttol = options.OptTol;
     244  double steptol = options.StepTol;
     245  int stepSame = 1; /* options.StepSame;   // 1 means stepx == stepz */
     246  double x0min = options.x0min;
     247  double z0min = options.z0min;
     248  double mu0 = options.mu0;
     249  int LSproblem = options.LSproblem; // See below
     250  int LSmethod = options.LSmethod; // 1=Cholesky    2=QR    3=LSQR
     251  int itnlim = options.LSQRMaxIter * CoinMin(m, n);
     252  double atol1 = options.LSQRatol1; // Initial  atol
     253  double atol2 = options.LSQRatol2; // Smallest atol,unless atol1 is smaller
     254  double conlim = options.LSQRconlim;
     255  //int  wait      = options.wait;
     256
     257  // LSproblem:
     258  //  1 = dy          2 = dy shifted, DLS
     259  // 11 = s          12 =  s shifted, DLS    (dx = Ds)
     260  // 21 = dx
     261  // 31 = 3x3 system, symmetrized by Z^{1/2}
     262  // 32 = 2x2 system, symmetrized by X^{1/2}
     263
     264  //---------------------------------------------------------------------
     265  // Set other parameters.
     266  //---------------------------------------------------------------------
     267  int kminor = 0; // 1 stops after each iteration
     268  double eta = 1e-4; // Linesearch tolerance for "sufficient descent"
     269  double maxf = 10; // Linesearch backtrack limit (function evaluations)
     270  double maxfail = 1; // Linesearch failure limit (consecutive iterations)
     271  double bigcenter = 1e+3; // mu is reduced if center < bigcenter.
     272
     273  // Parameters for LSQR.
     274  double atolmin = eps; // Smallest atol if linesearch back-tracks
     275  double btol = 0; // Should be small (zero is ok)
     276  double show = false; // Controls lsqr iteration log
     277  /*
    282278     double gamma     = d1->infNorm();
    283279     double delta     = d2->infNorm();
    284280     */
    285      double gamma = d1;
    286      double delta = d2;
    287 
    288      printf("\n\nx0min    = %8g     featol   = %8.1e", x0min, featol);
    289      printf(                  "      d1max   = %8.1e", gamma);
    290      printf(  "\nz0min    = %8g     opttol   = %8.1e", z0min, opttol);
    291      printf(                  "      d2max   = %8.1e", delta);
    292      printf(  "\nmu0      = %8.1e     steptol  = %8g", mu0  , steptol);
    293      printf(                  "     bigcenter= %8g"  , bigcenter);
    294 
    295      printf("\n\nLSQR:");
    296      printf("\natol1    = %8.1e     atol2    = %8.1e", atol1 , atol2 );
    297      printf(                  "      btol    = %8.1e", btol );
    298      printf("\nconlim   = %8.1e     itnlim   = %8d"  , conlim, itnlim);
    299      printf(                  "      show    = %8g"  , show );
    300 
    301 // LSmethod  = 3;  ////// Hardwire LSQR
    302 // LSproblem = 1;  ////// and LS problem defining "dy".
    303      /*
     281  double gamma = d1;
     282  double delta = d2;
     283
     284  printf("\n\nx0min    = %8g     featol   = %8.1e", x0min, featol);
     285  printf("      d1max   = %8.1e", gamma);
     286  printf("\nz0min    = %8g     opttol   = %8.1e", z0min, opttol);
     287  printf("      d2max   = %8.1e", delta);
     288  printf("\nmu0      = %8.1e     steptol  = %8g", mu0, steptol);
     289  printf("     bigcenter= %8g", bigcenter);
     290
     291  printf("\n\nLSQR:");
     292  printf("\natol1    = %8.1e     atol2    = %8.1e", atol1, atol2);
     293  printf("      btol    = %8.1e", btol);
     294  printf("\nconlim   = %8.1e     itnlim   = %8d", conlim, itnlim);
     295  printf("      show    = %8g", show);
     296
     297  // LSmethod  = 3;  ////// Hardwire LSQR
     298  // LSproblem = 1;  ////// and LS problem defining "dy".
     299  /*
    304300       if wait
    305301         printf("\n\nReview parameters... then type "return"\n")
     
    307303       end
    308304     */
    309      if (eta < 0)
    310           printf("\n\nLinesearch disabled by eta < 0");
    311 
    312      //---------------------------------------------------------------------
    313      // All parameters have now been set.
    314      //---------------------------------------------------------------------
    315      double time    = CoinCpuTime();
    316      //bool useChol = (LSmethod == 1);
    317      //bool useQR   = (LSmethod == 2);
    318      bool direct = (LSmethod <= 2 && ifexplicit);
    319      char solver[7];
    320      strcpy(solver, "  LSQR");
    321 
    322 
    323      //---------------------------------------------------------------------
    324      // Categorize bounds and allow for fixed variables by modifying b.
    325      //---------------------------------------------------------------------
    326 
    327      int nlow, nupp, nfix;
    328      int *bptrs[3] = {0};
    329      getBoundTypes(&nlow, &nupp, &nfix, bptrs );
    330      int *low = bptrs[0];
    331      int *upp = bptrs[1];
    332      int *fix = bptrs[2];
    333 
    334      int nU = n;
    335      if (nupp == 0) nU = 1; //Make dummy vectors if no Upper bounds
    336 
    337      //---------------------------------------------------------------------
    338      //  Get pointers to local copy of model bounds
    339      //---------------------------------------------------------------------
    340 
    341      CoinDenseVector<double> bl(n, columnLower_);
    342      double *bl_elts = bl.getElements();
    343      CoinDenseVector<double> bu(nU, columnUpper_); // this is dummy if no UB
    344      double *bu_elts = bu.getElements();
    345 
    346      CoinDenseVector<double> r1(m, 0.0);
    347      double *r1_elts = r1.getElements();
    348      CoinDenseVector<double> x1(n, 0.0);
    349      double *x1_elts = x1.getElements();
    350 
    351      if (nfix > 0) {
    352           for (int k = 0; k < nfix; k++)
    353                x1_elts[fix[k]] = bl[fix[k]];
    354           matVecMult(1, r1, x1);
    355           b = b - r1;
    356           // At some stage, might want to look at normfix = norm(r1,inf);
    357      }
    358 
    359      //---------------------------------------------------------------------
    360      // Scale the input data.
    361      // The scaled variables are
    362      //    xbar     = x/beta,
    363      //    ybar     = y/zeta,
    364      //    zbar     = z/zeta.
    365      // Define
    366      //    theta    = beta*zeta;
    367      // The scaled function is
    368      //    phibar   = ( 1   /theta) fbar(beta*xbar),
    369      //    gradient = (beta /theta) grad,
    370      //    Hessian  = (beta2/theta) hess.
    371      //---------------------------------------------------------------------
    372      double beta = xsize_;
    373      if (beta == 0) beta = 1; // beta scales b, x.
    374      double zeta = zsize_;
    375      if (zeta == 0) zeta = 1; // zeta scales y, z.
    376      double theta  = beta * zeta;                          // theta scales obj.
    377      // (theta could be anything, but theta = beta*zeta makes
    378      // scaled grad = grad/zeta = 1 approximately if zeta is chosen right.)
    379 
    380      for (int k = 0; k < nlow; k++)
    381           bl_elts[low[k]] = bl_elts[low[k]] / beta;
    382      for (int k = 0; k < nupp; k++)
    383           bu_elts[upp[k]] = bu_elts[upp[k]] / beta;
    384      d1     = d1 * ( beta / sqrt(theta) );
    385      d2     = d2 * ( sqrt(theta) / beta );
    386 
    387      double beta2  = beta * beta;
    388      b.scale( (1.0 / beta) );
    389      y.scale( (1.0 / zeta) );
    390      x.scale( (1.0 / beta) );
    391      z.scale( (1.0 / zeta) );
    392 
    393      //---------------------------------------------------------------------
    394      // Initialize vectors that are not fully used if bounds are missing.
    395      //---------------------------------------------------------------------
    396      CoinDenseVector<double> rL(n, 0.0);
    397      CoinDenseVector<double> cL(n, 0.0);
    398      CoinDenseVector<double> z1(n, 0.0);
    399      CoinDenseVector<double> dx1(n, 0.0);
    400      CoinDenseVector<double> dz1(n, 0.0);
    401      CoinDenseVector<double> r2(n, 0.0);
    402 
    403      // Assign upper bd regions (dummy if no UBs)
    404 
    405      CoinDenseVector<double> rU(nU, 0.0);
    406      CoinDenseVector<double> cU(nU, 0.0);
    407      CoinDenseVector<double> x2(nU, 0.0);
    408      CoinDenseVector<double> z2(nU, 0.0);
    409      CoinDenseVector<double> dx2(nU, 0.0);
    410      CoinDenseVector<double> dz2(nU, 0.0);
    411 
    412      //---------------------------------------------------------------------
    413      // Initialize x, y, z, objective, etc.
    414      //---------------------------------------------------------------------
    415      CoinDenseVector<double> dx(n, 0.0);
    416      CoinDenseVector<double> dy(m, 0.0);
    417      CoinDenseVector<double> Pr(m);
    418      CoinDenseVector<double> D(n);
    419      double *D_elts = D.getElements();
    420      CoinDenseVector<double> w(n);
    421      double *w_elts = w.getElements();
    422      CoinDenseVector<double> rhs(m + n);
    423 
    424 
    425      //---------------------------------------------------------------------
    426      // Pull out the element array pointers for efficiency
    427      //---------------------------------------------------------------------
    428      double *x_elts  = x.getElements();
    429      double *x2_elts = x2.getElements();
    430      double *z_elts  = z.getElements();
    431      double *z1_elts = z1.getElements();
    432      double *z2_elts = z2.getElements();
    433 
    434      for (int k = 0; k < nlow; k++) {
    435           x_elts[low[k]]  = CoinMax( x_elts[low[k]], bl[low[k]]);
    436           x1_elts[low[k]] = CoinMax( x_elts[low[k]] - bl[low[k]], x0min  );
    437           z1_elts[low[k]] = CoinMax( z_elts[low[k]], z0min  );
    438      }
    439      for (int k = 0; k < nupp; k++) {
    440           x_elts[upp[k]]  = CoinMin( x_elts[upp[k]], bu[upp[k]]);
    441           x2_elts[upp[k]] = CoinMax(bu[upp[k]] -  x_elts[upp[k]], x0min  );
    442           z2_elts[upp[k]] = CoinMax(-z_elts[upp[k]], z0min  );
    443      }
    444      //////////////////// Assume hessian is diagonal. //////////////////////
    445 
    446 //  [obj,grad,hess] = feval( Fname, (x*beta) );
    447      x.scale(beta);
    448      double obj = getObj(x);
    449      CoinDenseVector<double> grad(n);
    450      getGrad(x, grad);
    451      CoinDenseVector<double> H(n);
    452      getHessian(x , H);
    453      x.scale((1.0 / beta));
    454 
    455      //double * g_elts = grad.getElements();
    456      double * H_elts = H.getElements();
    457 
    458      obj /= theta;                       // Scaled obj.
    459      grad = grad * (beta / theta) + (d1 * d1) * x; // grad includes x regularization.
    460      H  = H * (beta2 / theta) + (d1 * d1)      ; // H    includes x regularization.
    461 
    462 
    463      /*---------------------------------------------------------------------
     305  if (eta < 0)
     306    printf("\n\nLinesearch disabled by eta < 0");
     307
     308  //---------------------------------------------------------------------
     309  // All parameters have now been set.
     310  //---------------------------------------------------------------------
     311  double time = CoinCpuTime();
     312  //bool useChol = (LSmethod == 1);
     313  //bool useQR   = (LSmethod == 2);
     314  bool direct = (LSmethod <= 2 && ifexplicit);
     315  char solver[7];
     316  strcpy(solver, "  LSQR");
     317
     318  //---------------------------------------------------------------------
     319  // Categorize bounds and allow for fixed variables by modifying b.
     320  //---------------------------------------------------------------------
     321
     322  int nlow, nupp, nfix;
     323  int *bptrs[3] = { 0 };
     324  getBoundTypes(&nlow, &nupp, &nfix, bptrs);
     325  int *low = bptrs[0];
     326  int *upp = bptrs[1];
     327  int *fix = bptrs[2];
     328
     329  int nU = n;
     330  if (nupp == 0)
     331    nU = 1; //Make dummy vectors if no Upper bounds
     332
     333  //---------------------------------------------------------------------
     334  //  Get pointers to local copy of model bounds
     335  //---------------------------------------------------------------------
     336
     337  CoinDenseVector< double > bl(n, columnLower_);
     338  double *bl_elts = bl.getElements();
     339  CoinDenseVector< double > bu(nU, columnUpper_); // this is dummy if no UB
     340  double *bu_elts = bu.getElements();
     341
     342  CoinDenseVector< double > r1(m, 0.0);
     343  double *r1_elts = r1.getElements();
     344  CoinDenseVector< double > x1(n, 0.0);
     345  double *x1_elts = x1.getElements();
     346
     347  if (nfix > 0) {
     348    for (int k = 0; k < nfix; k++)
     349      x1_elts[fix[k]] = bl[fix[k]];
     350    matVecMult(1, r1, x1);
     351    b = b - r1;
     352    // At some stage, might want to look at normfix = norm(r1,inf);
     353  }
     354
     355  //---------------------------------------------------------------------
     356  // Scale the input data.
     357  // The scaled variables are
     358  //    xbar     = x/beta,
     359  //    ybar     = y/zeta,
     360  //    zbar     = z/zeta.
     361  // Define
     362  //    theta    = beta*zeta;
     363  // The scaled function is
     364  //    phibar   = ( 1   /theta) fbar(beta*xbar),
     365  //    gradient = (beta /theta) grad,
     366  //    Hessian  = (beta2/theta) hess.
     367  //---------------------------------------------------------------------
     368  double beta = xsize_;
     369  if (beta == 0)
     370    beta = 1; // beta scales b, x.
     371  double zeta = zsize_;
     372  if (zeta == 0)
     373    zeta = 1; // zeta scales y, z.
     374  double theta = beta * zeta; // theta scales obj.
     375  // (theta could be anything, but theta = beta*zeta makes
     376  // scaled grad = grad/zeta = 1 approximately if zeta is chosen right.)
     377
     378  for (int k = 0; k < nlow; k++)
     379    bl_elts[low[k]] = bl_elts[low[k]] / beta;
     380  for (int k = 0; k < nupp; k++)
     381    bu_elts[upp[k]] = bu_elts[upp[k]] / beta;
     382  d1 = d1 * (beta / sqrt(theta));
     383  d2 = d2 * (sqrt(theta) / beta);
     384
     385  double beta2 = beta * beta;
     386  b.scale((1.0 / beta));
     387  y.scale((1.0 / zeta));
     388  x.scale((1.0 / beta));
     389  z.scale((1.0 / zeta));
     390
     391  //---------------------------------------------------------------------
     392  // Initialize vectors that are not fully used if bounds are missing.
     393  //---------------------------------------------------------------------
     394  CoinDenseVector< double > rL(n, 0.0);
     395  CoinDenseVector< double > cL(n, 0.0);
     396  CoinDenseVector< double > z1(n, 0.0);
     397  CoinDenseVector< double > dx1(n, 0.0);
     398  CoinDenseVector< double > dz1(n, 0.0);
     399  CoinDenseVector< double > r2(n, 0.0);
     400
     401  // Assign upper bd regions (dummy if no UBs)
     402
     403  CoinDenseVector< double > rU(nU, 0.0);
     404  CoinDenseVector< double > cU(nU, 0.0);
     405  CoinDenseVector< double > x2(nU, 0.0);
     406  CoinDenseVector< double > z2(nU, 0.0);
     407  CoinDenseVector< double > dx2(nU, 0.0);
     408  CoinDenseVector< double > dz2(nU, 0.0);
     409
     410  //---------------------------------------------------------------------
     411  // Initialize x, y, z, objective, etc.
     412  //---------------------------------------------------------------------
     413  CoinDenseVector< double > dx(n, 0.0);
     414  CoinDenseVector< double > dy(m, 0.0);
     415  CoinDenseVector< double > Pr(m);
     416  CoinDenseVector< double > D(n);
     417  double *D_elts = D.getElements();
     418  CoinDenseVector< double > w(n);
     419  double *w_elts = w.getElements();
     420  CoinDenseVector< double > rhs(m + n);
     421
     422  //---------------------------------------------------------------------
     423  // Pull out the element array pointers for efficiency
     424  //---------------------------------------------------------------------
     425  double *x_elts = x.getElements();
     426  double *x2_elts = x2.getElements();
     427  double *z_elts = z.getElements();
     428  double *z1_elts = z1.getElements();
     429  double *z2_elts = z2.getElements();
     430
     431  for (int k = 0; k < nlow; k++) {
     432    x_elts[low[k]] = CoinMax(x_elts[low[k]], bl[low[k]]);
     433    x1_elts[low[k]] = CoinMax(x_elts[low[k]] - bl[low[k]], x0min);
     434    z1_elts[low[k]] = CoinMax(z_elts[low[k]], z0min);
     435  }
     436  for (int k = 0; k < nupp; k++) {
     437    x_elts[upp[k]] = CoinMin(x_elts[upp[k]], bu[upp[k]]);
     438    x2_elts[upp[k]] = CoinMax(bu[upp[k]] - x_elts[upp[k]], x0min);
     439    z2_elts[upp[k]] = CoinMax(-z_elts[upp[k]], z0min);
     440  }
     441  //////////////////// Assume hessian is diagonal. //////////////////////
     442
     443  //  [obj,grad,hess] = feval( Fname, (x*beta) );
     444  x.scale(beta);
     445  double obj = getObj(x);
     446  CoinDenseVector< double > grad(n);
     447  getGrad(x, grad);
     448  CoinDenseVector< double > H(n);
     449  getHessian(x, H);
     450  x.scale((1.0 / beta));
     451
     452  //double * g_elts = grad.getElements();
     453  double *H_elts = H.getElements();
     454
     455  obj /= theta; // Scaled obj.
     456  grad = grad * (beta / theta) + (d1 * d1) * x; // grad includes x regularization.
     457  H = H * (beta2 / theta) + (d1 * d1); // H    includes x regularization.
     458
     459  /*---------------------------------------------------------------------
    464460     // Compute primal and dual residuals:
    465461     // r1 =  b - Aprod(x) - d2*d2*y;
     
    467463     //  rL =  bl - x + x1;
    468464     //  rU =  x + x2 - bu; */
    469      //---------------------------------------------------------------------
    470      //  [r1,r2,rL,rU,Pinf,Dinf] = ...
    471      //      pdxxxresid1( Aname,fix,low,upp, ...
    472      //                   b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 );
    473      pdxxxresid1( this, nlow, nupp, nfix, low, upp, fix,
    474                   b, bl_elts, bu_elts, d1, d2, grad, rL, rU, x, x1, x2, y, z1, z2,
    475                   r1, r2, &Pinf, &Dinf);
    476      //---------------------------------------------------------------------
    477      // Initialize mu and complementarity residuals:
    478      //    cL   = mu*e - X1*z1.
    479      //    cU   = mu*e - X2*z2.
    480      //
    481      // 25 Jan 2001: Now that b and obj are scaled (and hence x,y,z),
    482      //              we should be able to use mufirst = mu0 (absolute value).
    483      //              0.1 worked poorly on StarTest1 with x0min = z0min = 0.1.
    484      // 29 Jan 2001: We might as well use mu0 = x0min * z0min;
    485      //              so that most variables are centered after a warm start.
    486      // 29 Sep 2002: Use mufirst = mu0*(x0min * z0min),
    487      //              regarding mu0 as a scaling of the initial center.
    488      //---------------------------------------------------------------------
    489      //  double mufirst = mu0*(x0min * z0min);
    490      double mufirst = mu0;  // revert to absolute value
    491      double mulast = 0.1 * opttol;
    492      mulast  = CoinMin( mulast, mufirst );
    493      double mu      = mufirst;
    494      double center, fmerit;
    495      pdxxxresid2( mu, nlow, nupp, low, upp, cL, cU, x1, x2,
    496                   z1, z2, &center, &Cinf, &Cinf0 );
    497      fmerit = pdxxxmerit(nlow, nupp, low, upp, r1, r2, rL, rU, cL, cU );
    498 
    499      // Initialize other things.
    500 
    501      bool  precon  = true;
    502      double PDitns    = 0;
    503      //bool converged = false;
    504      double atol      = atol1;
    505      atol2     = CoinMax( atol2, atolmin );
    506      atolmin  = atol2;
    507      //  pdDDD2    = d2;    // Global vector for diagonal matrix D2
    508 
    509      //  Iteration log.
    510 
    511      int nf      = 0;
    512      int itncg  = 0;
    513      int nfail  = 0;
    514 
    515      printf("\n\nItn   mu   stepx   stepz  Pinf  Dinf");
    516      printf("  Cinf   Objective    nf  center");
    517      if (direct) {
    518           printf("\n");
    519      } else {
    520           printf("  atol   solver   Inexact\n");
    521      }
    522 
    523      double regx = (d1 * x).twoNorm();
    524      double regy = (d2 * y).twoNorm();
    525      //  regterm = twoNorm(d1.*x)^2  +  norm(d2.*y)^2;
    526      double regterm = regx * regx + regy * regy;
    527      double objreg  = obj  + 0.5 * regterm;
    528      double objtrue = objreg * theta;
    529 
    530      printf("\n%3g                     ", PDitns        );
    531      printf("%6.1f%6.1f" , log10(Pinf ), log10(Dinf));
    532      printf("%6.1f%15.7e", log10(Cinf0), objtrue    );
    533      printf("   %8.1f\n"   , center                   );
    534      /*
     465  //---------------------------------------------------------------------
     466  //  [r1,r2,rL,rU,Pinf,Dinf] = ...
     467  //      pdxxxresid1( Aname,fix,low,upp, ...
     468  //                   b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 );
     469  pdxxxresid1(this, nlow, nupp, nfix, low, upp, fix,
     470    b, bl_elts, bu_elts, d1, d2, grad, rL, rU, x, x1, x2, y, z1, z2,
     471    r1, r2, &Pinf, &Dinf);
     472  //---------------------------------------------------------------------
     473  // Initialize mu and complementarity residuals:
     474  //    cL   = mu*e - X1*z1.
     475  //    cU   = mu*e - X2*z2.
     476  //
     477  // 25 Jan 2001: Now that b and obj are scaled (and hence x,y,z),
     478  //              we should be able to use mufirst = mu0 (absolute value).
     479  //              0.1 worked poorly on StarTest1 with x0min = z0min = 0.1.
     480  // 29 Jan 2001: We might as well use mu0 = x0min * z0min;
     481  //              so that most variables are centered after a warm start.
     482  // 29 Sep 2002: Use mufirst = mu0*(x0min * z0min),
     483  //              regarding mu0 as a scaling of the initial center.
     484  //---------------------------------------------------------------------
     485  //  double mufirst = mu0*(x0min * z0min);
     486  double mufirst = mu0; // revert to absolute value
     487  double mulast = 0.1 * opttol;
     488  mulast = CoinMin(mulast, mufirst);
     489  double mu = mufirst;
     490  double center, fmerit;
     491  pdxxxresid2(mu, nlow, nupp, low, upp, cL, cU, x1, x2,
     492    z1, z2, &center, &Cinf, &Cinf0);
     493  fmerit = pdxxxmerit(nlow, nupp, low, upp, r1, r2, rL, rU, cL, cU);
     494
     495  // Initialize other things.
     496
     497  bool precon = true;
     498  double PDitns = 0;
     499  //bool converged = false;
     500  double atol = atol1;
     501  atol2 = CoinMax(atol2, atolmin);
     502  atolmin = atol2;
     503  //  pdDDD2    = d2;    // Global vector for diagonal matrix D2
     504
     505  //  Iteration log.
     506
     507  int nf = 0;
     508  int itncg = 0;
     509  int nfail = 0;
     510
     511  printf("\n\nItn   mu   stepx   stepz  Pinf  Dinf");
     512  printf("  Cinf   Objective    nf  center");
     513  if (direct) {
     514    printf("\n");
     515  } else {
     516    printf("  atol   solver   Inexact\n");
     517  }
     518
     519  double regx = (d1 * x).twoNorm();
     520  double regy = (d2 * y).twoNorm();
     521  //  regterm = twoNorm(d1.*x)^2  +  norm(d2.*y)^2;
     522  double regterm = regx * regx + regy * regy;
     523  double objreg = obj + 0.5 * regterm;
     524  double objtrue = objreg * theta;
     525
     526  printf("\n%3g                     ", PDitns);
     527  printf("%6.1f%6.1f", log10(Pinf), log10(Dinf));
     528  printf("%6.1f%15.7e", log10(Cinf0), objtrue);
     529  printf("   %8.1f\n", center);
     530  /*
    535531     if kminor
    536532       printf("\n\nStart of first minor itn...\n");
     
    538534     end
    539535     */
    540      //---------------------------------------------------------------------
    541      // Main loop.
    542      //---------------------------------------------------------------------
    543      // Lsqr
    544      ClpLsqr thisLsqr(this);
    545      //  while (converged) {
    546      while(PDitns < maxitn) {
    547           PDitns = PDitns + 1;
    548 
    549           // 31 Jan 2001: Set atol according to progress, a la Inexact Newton.
    550           // 07 Feb 2001: 0.1 not small enough for Satellite problem.  Try 0.01.
    551           // 25 Apr 2001: 0.01 seems wasteful for Star problem.
    552           //              Now that starting conditions are better, go back to 0.1.
    553 
    554           double r3norm = CoinMax(Pinf,   CoinMax(Dinf, Cinf));
    555           atol   = CoinMin(atol, r3norm * 0.1);
    556           atol   = CoinMax(atol,  atolmin   );
    557           info.r3norm = r3norm;
    558 
    559           //-------------------------------------------------------------------
    560           //  Define a damped Newton iteration for solving f = 0,
    561           //  keeping  x1, x2, z1, z2 > 0.  We eliminate dx1, dx2, dz1, dz2
    562           //  to obtain the system
    563           //
    564           //     [-H2  A"  ] [ dx ] = [ w ],   H2 = H + D1^2 + X1inv Z1 + X2inv Z2,
    565           //     [ A   D2^2] [ dy ] = [ r1]    w  = r2 - X1inv(cL + Z1 rL)
    566           //                                           + X2inv(cU + Z2 rU),
    567           //
    568           //  which is equivalent to the least-squares problem
    569           //
    570           //     min || [ D A"]dy  -  [  D w   ] ||,   D = H2^{-1/2}.         (*)
    571           //         || [  D2 ]       [D2inv r1] ||
    572           //-------------------------------------------------------------------
    573           for (int k = 0; k < nlow; k++)
    574                H_elts[low[k]] = H_elts[low[k]] + z1[low[k]] / x1[low[k]];
    575           for (int k = 0; k < nupp; k++)
    576                H[upp[k]] = H[upp[k]] + z2[upp[k]] / x2[upp[k]];
    577           w = r2;
    578           for (int k = 0; k < nlow; k++)
    579                w[low[k]] = w[low[k]] - (cL[low[k]] + z1[low[k]] * rL[low[k]]) / x1[low[k]];
    580           for (int k = 0; k < nupp; k++)
    581                w[upp[k]] = w[upp[k]] + (cU[upp[k]] + z2[upp[k]] * rU[upp[k]]) / x2[upp[k]];
    582 
    583           if (LSproblem == 1) {
    584                //-----------------------------------------------------------------
    585                //  Solve (*) for dy.
    586                //-----------------------------------------------------------------
    587                H      = 1.0 / H; // H is now Hinv (NOTE!)
    588                for (int k = 0; k < nfix; k++)
    589                     H[fix[k]] = 0;
    590                for (int k = 0; k < n; k++)
    591                     D_elts[k] = sqrt(H_elts[k]);
    592                thisLsqr.borrowDiag1(D_elts);
    593                thisLsqr.diag2_ = d2;
    594 
    595                if (direct) {
    596                     // Omit direct option for now
    597                } else {// Iterative solve using LSQR.
    598                     //rhs     = [ D.*w; r1./d2 ];
    599                     for (int k = 0; k < n; k++)
    600                          rhs[k] = D_elts[k] * w_elts[k];
    601                     for (int k = 0; k < m; k++)
    602                          rhs[n+k] = r1_elts[k] * (1.0 / d2);
    603                     double damp    = 0;
    604 
    605                     if (precon) {  // Construct diagonal preconditioner for LSQR
    606                          matPrecon(d2, Pr, D);
    607                     }
    608                     /*
     536  //---------------------------------------------------------------------
     537  // Main loop.
     538  //---------------------------------------------------------------------
     539  // Lsqr
     540  ClpLsqr thisLsqr(this);
     541  //  while (converged) {
     542  while (PDitns < maxitn) {
     543    PDitns = PDitns + 1;
     544
     545    // 31 Jan 2001: Set atol according to progress, a la Inexact Newton.
     546    // 07 Feb 2001: 0.1 not small enough for Satellite problem.  Try 0.01.
     547    // 25 Apr 2001: 0.01 seems wasteful for Star problem.
     548    //              Now that starting conditions are better, go back to 0.1.
     549
     550    double r3norm = CoinMax(Pinf, CoinMax(Dinf, Cinf));
     551    atol = CoinMin(atol, r3norm * 0.1);
     552    atol = CoinMax(atol, atolmin);
     553    info.r3norm = r3norm;
     554
     555    //-------------------------------------------------------------------
     556    //  Define a damped Newton iteration for solving f = 0,
     557    //  keeping  x1, x2, z1, z2 > 0.  We eliminate dx1, dx2, dz1, dz2
     558    //  to obtain the system
     559    //
     560    //     [-H2  A"  ] [ dx ] = [ w ],   H2 = H + D1^2 + X1inv Z1 + X2inv Z2,
     561    //     [ A   D2^2] [ dy ] = [ r1]    w  = r2 - X1inv(cL + Z1 rL)
     562    //                                           + X2inv(cU + Z2 rU),
     563    //
     564    //  which is equivalent to the least-squares problem
     565    //
     566    //     min || [ D A"]dy  -  [  D w   ] ||,   D = H2^{-1/2}.         (*)
     567    //         || [  D2 ]       [D2inv r1] ||
     568    //-------------------------------------------------------------------
     569    for (int k = 0; k < nlow; k++)
     570      H_elts[low[k]] = H_elts[low[k]] + z1[low[k]] / x1[low[k]];
     571    for (int k = 0; k < nupp; k++)
     572      H[upp[k]] = H[upp[k]] + z2[upp[k]] / x2[upp[k]];
     573    w = r2;
     574    for (int k = 0; k < nlow; k++)
     575      w[low[k]] = w[low[k]] - (cL[low[k]] + z1[low[k]] * rL[low[k]]) / x1[low[k]];
     576    for (int k = 0; k < nupp; k++)
     577      w[upp[k]] = w[upp[k]] + (cU[upp[k]] + z2[upp[k]] * rU[upp[k]]) / x2[upp[k]];
     578
     579    if (LSproblem == 1) {
     580      //-----------------------------------------------------------------
     581      //  Solve (*) for dy.
     582      //-----------------------------------------------------------------
     583      H = 1.0 / H; // H is now Hinv (NOTE!)
     584      for (int k = 0; k < nfix; k++)
     585        H[fix[k]] = 0;
     586      for (int k = 0; k < n; k++)
     587        D_elts[k] = sqrt(H_elts[k]);
     588      thisLsqr.borrowDiag1(D_elts);
     589      thisLsqr.diag2_ = d2;
     590
     591      if (direct) {
     592        // Omit direct option for now
     593      } else { // Iterative solve using LSQR.
     594        //rhs     = [ D.*w; r1./d2 ];
     595        for (int k = 0; k < n; k++)
     596          rhs[k] = D_elts[k] * w_elts[k];
     597        for (int k = 0; k < m; k++)
     598          rhs[n + k] = r1_elts[k] * (1.0 / d2);
     599        double damp = 0;
     600
     601        if (precon) { // Construct diagonal preconditioner for LSQR
     602          matPrecon(d2, Pr, D);
     603        }
     604        /*
    609605                        rw(7)        = precon;
    610606                            info.atolmin = atolmin;
     
    621617                        thisLsqr.do_lsqr(this);
    622618                        */
    623                     //  New version of lsqr
    624 
    625                     int istop;
    626                     dy.clear();
    627                     show = false;
    628                     info.atolmin = atolmin;
    629                     info.r3norm  = fmerit;  // Must be the 2-norm here.
    630 
    631                     thisLsqr.do_lsqr( rhs, damp, atol, btol, conlim, itnlim,
    632                                       show, info, dy , &istop, &itncg, &outfo, precon, Pr);
    633                     if (precon)
    634                          dy = dy * Pr;
    635 
    636                     if (!precon && itncg > 999999)
    637                          precon = true;
    638 
    639                     if (istop == 3  ||  istop == 7 )  // conlim or itnlim
    640                          printf("\n    LSQR stopped early:  istop = //%d", istop);
    641 
    642 
    643                     atolold   = outfo.atolold;
    644                     atol      = outfo.atolnew;
    645                     r3ratio   = outfo.r3ratio;
    646                }// LSproblem 1
    647 
    648                //      grad      = pdxxxmat( Aname,2,m,n,dy );   // grad = A"dy
    649                grad.clear();
    650                matVecMult(2, grad, dy);
    651                for (int k = 0; k < nfix; k++)
    652                     grad[fix[k]] = 0;                            // grad is a work vector
    653                dx = H * (grad - w);
    654 
    655           } else {
    656                perror( "This LSproblem not yet implemented\n" );
    657           }
    658           //-------------------------------------------------------------------
    659 
    660           CGitns += itncg;
    661 
    662           //-------------------------------------------------------------------
    663           // dx and dy are now known.  Get dx1, dx2, dz1, dz2.
    664           //-------------------------------------------------------------------
    665           for (int k = 0; k < nlow; k++) {
    666                dx1[low[k]] = - rL[low[k]] + dx[low[k]];
    667                dz1[low[k]] =  (cL[low[k]] - z1[low[k]] * dx1[low[k]]) / x1[low[k]];
    668           }
    669           for (int k = 0; k < nupp; k++) {
    670                dx2[upp[k]] = - rU[upp[k]] - dx[upp[k]];
    671                dz2[upp[k]] =  (cU[upp[k]] - z2[upp[k]] * dx2[upp[k]]) / x2[upp[k]];
    672           }
    673           //-------------------------------------------------------------------
    674           // Find the maximum step.
    675           //--------------------------------------------------------------------
    676           double stepx1 = pdxxxstep(nlow, low, x1, dx1 );
    677           double stepx2 = inf;
    678           if (nupp > 0)
    679                stepx2 = pdxxxstep(nupp, upp, x2, dx2 );
    680           double stepz1 = pdxxxstep( z1     , dz1      );
    681           double stepz2 = inf;
    682           if (nupp > 0)
    683                stepz2 = pdxxxstep( z2     , dz2      );
    684           double stepx  = CoinMin( stepx1, stepx2 );
    685           double stepz  = CoinMin( stepz1, stepz2 );
    686           stepx  = CoinMin( steptol * stepx, 1.0 );
    687           stepz  = CoinMin( steptol * stepz, 1.0 );
    688           if (stepSame) {                  // For NLPs, force same step
    689                stepx = CoinMin( stepx, stepz );   // (true Newton method)
    690                stepz = stepx;
    691           }
    692 
    693           //-------------------------------------------------------------------
    694           // Backtracking linesearch.
    695           //-------------------------------------------------------------------
    696           bool fail     =  true;
    697           nf       =  0;
    698 
    699           while (nf < maxf) {
    700                nf      = nf + 1;
    701                x       = x        +  stepx * dx;
    702                y       = y        +  stepz * dy;
    703                for (int k = 0; k < nlow; k++) {
    704                     x1[low[k]] = x1[low[k]]  +  stepx * dx1[low[k]];
    705                     z1[low[k]] = z1[low[k]]  +  stepz * dz1[low[k]];
    706                }
    707                for (int k = 0; k < nupp; k++) {
    708                     x2[upp[k]] = x2[upp[k]]  +  stepx * dx2[upp[k]];
    709                     z2[upp[k]] = z2[upp[k]]  +  stepz * dz2[upp[k]];
    710                }
    711                //      [obj,grad,hess] = feval( Fname, (x*beta) );
    712                x.scale(beta);
    713                obj = getObj(x);
    714                getGrad(x, grad);
    715                getHessian(x, H);
    716                x.scale((1.0 / beta));
    717 
    718                obj        /= theta;
    719                grad       = grad * (beta / theta)  +  d1 * d1 * x;
    720                H          = H * (beta2 / theta)  +  d1 * d1;
    721 
    722                //      [r1,r2,rL,rU,Pinf,Dinf] = ...
    723                pdxxxresid1( this, nlow, nupp, nfix, low, upp, fix,
    724                             b, bl_elts, bu_elts, d1, d2, grad, rL, rU, x, x1, x2,
    725                             y, z1, z2, r1, r2, &Pinf, &Dinf );
    726                //double center, Cinf, Cinf0;
    727                //      [cL,cU,center,Cinf,Cinf0] = ...
    728                pdxxxresid2( mu, nlow, nupp, low, upp, cL, cU, x1, x2, z1, z2,
    729                             &center, &Cinf, &Cinf0);
    730                double fmeritnew = pdxxxmerit(nlow, nupp, low, upp, r1, r2, rL, rU, cL, cU );
    731                double step      = CoinMin( stepx, stepz );
    732 
    733                if (fmeritnew <= (1 - eta * step)*fmerit) {
    734                     fail = false;
    735                     break;
    736                }
    737 
    738                // Merit function didn"t decrease.
    739                // Restore variables to previous values.
    740                // (This introduces a little error, but save lots of space.)
    741 
    742                x       = x        -  stepx * dx;
    743                y       = y        -  stepz * dy;
    744                for (int k = 0; k < nlow; k++) {
    745                     x1[low[k]] = x1[low[k]]  -  stepx * dx1[low[k]];
    746                     z1[low[k]] = z1[low[k]]  -  stepz * dz1[low[k]];
    747                }
    748                for (int k = 0; k < nupp; k++) {
    749                     x2[upp[k]] = x2[upp[k]]  -  stepx * dx2[upp[k]];
    750                     z2[upp[k]] = z2[upp[k]]  -  stepz * dz2[upp[k]];
    751                }
    752                // Back-track.
    753                // If it"s the first time,
    754                // make stepx and stepz the same.
    755 
    756                if (nf == 1 && stepx != stepz) {
    757                     stepx = step;
    758                } else if (nf < maxf) {
    759                     stepx = stepx / 2;
    760                }
    761                stepz = stepx;
    762           }
    763 
    764           if (fail) {
    765                printf("\n     Linesearch failed (nf too big)");
    766                nfail += 1;
    767           } else {
    768                nfail = 0;
    769           }
    770 
    771           //-------------------------------------------------------------------
    772           // Set convergence measures.
    773           //--------------------------------------------------------------------
    774           regx = (d1 * x).twoNorm();
    775           regy = (d2 * y).twoNorm();
    776           regterm = regx * regx + regy * regy;
    777           objreg  = obj  +  0.5 * regterm;
    778           objtrue = objreg * theta;
    779 
    780           bool primalfeas    = Pinf  <=  featol;
    781           bool dualfeas      = Dinf  <=  featol;
    782           bool complementary = Cinf0 <=  opttol;
    783           bool enough        = PDitns >=       4; // Prevent premature termination.
    784           bool converged     = primalfeas  &  dualfeas  &  complementary  &  enough;
    785 
    786           //-------------------------------------------------------------------
    787           // Iteration log.
    788           //-------------------------------------------------------------------
    789           char str1[100], str2[100], str3[100], str4[100], str5[100];
    790           sprintf(str1, "\n%3g%5.1f" , PDitns      , log10(mu)   );
    791           sprintf(str2, "%8.5f%8.5f" , stepx       , stepz       );
    792           if (stepx < 0.0001 || stepz < 0.0001) {
    793                sprintf(str2, " %6.1e %6.1e" , stepx       , stepz       );
    794           }
    795 
    796           sprintf(str3, "%6.1f%6.1f" , log10(Pinf) , log10(Dinf));
    797           sprintf(str4, "%6.1f%15.7e", log10(Cinf0), objtrue     );
    798           sprintf(str5, "%3d%8.1f"   , nf          , center      );
    799           if (center > 99999) {
    800                sprintf(str5, "%3d%8.1e"   , nf          , center      );
    801           }
    802           printf("%s%s%s%s%s", str1, str2, str3, str4, str5);
    803           if (direct) {
    804                // relax
    805           } else {
    806                printf(" %5.1f%7d%7.3f", log10(atolold), itncg, r3ratio);
    807           }
    808           //-------------------------------------------------------------------
    809           // Test for termination.
    810           //-------------------------------------------------------------------
    811           if (kminor) {
    812                printf( "\nStart of next minor itn...\n");
    813                //      keyboard;
    814           }
    815 
    816           if (converged) {
    817                printf("\n   Converged");
    818                break;
    819           } else if (PDitns >= maxitn) {
    820                printf("\n   Too many iterations");
    821                inform = 1;
    822                break;
    823           } else if (nfail  >= maxfail) {
    824                printf("\n   Too many linesearch failures");
    825                inform = 2;
    826                break;
    827           } else {
    828 
    829                // Reduce mu, and reset certain residuals.
    830 
    831                double stepmu  = CoinMin( stepx , stepz   );
    832                stepmu  = CoinMin( stepmu, steptol );
    833                double muold   = mu;
    834                mu      = mu   -  stepmu * mu;
    835                if (center >= bigcenter)
    836                     mu = muold;
    837 
    838                // mutrad = mu0*(sum(Xz)/n); // 24 May 1998: Traditional value, but
    839                // mu     = CoinMin(mu,mutrad ); // it seemed to decrease mu too much.
    840 
    841                mu      = CoinMax(mu, mulast); // 13 Jun 1998: No need for smaller mu.
    842                //      [cL,cU,center,Cinf,Cinf0] = ...
    843                pdxxxresid2( mu, nlow, nupp, low, upp, cL, cU, x1, x2, z1, z2,
    844                             &center, &Cinf, &Cinf0 );
    845                fmerit = pdxxxmerit( nlow, nupp, low, upp, r1, r2, rL, rU, cL, cU );
    846 
    847                // Reduce atol for LSQR (and SYMMLQ).
    848                // NOW DONE AT TOP OF LOOP.
    849 
    850                atolold = atol;
    851                // if atol > atol2
    852                //   atolfac = (mu/mufirst)^0.25;
    853                //   atol    = CoinMax( atol*atolfac, atol2 );
    854                // end
    855 
    856                // atol = CoinMin( atol, mu );     // 22 Jan 2001: a la Inexact Newton.
    857                // atol = CoinMin( atol, 0.5*mu ); // 30 Jan 2001: A bit tighter
    858 
    859                // If the linesearch took more than one function (nf > 1),
    860                // we assume the search direction needed more accuracy
    861                // (though this may be true only for LPs).
    862                // 12 Jun 1998: Ask for more accuracy if nf > 2.
    863                // 24 Nov 2000: Also if the steps are small.
    864                // 30 Jan 2001: Small steps might be ok with warm start.
    865                // 06 Feb 2001: Not necessarily.  Reinstated tests in next line.
    866 
    867                if (nf > 2  ||  CoinMin( stepx, stepz ) <= 0.01)
    868                     atol = atolold * 0.1;
    869           }
    870           //---------------------------------------------------------------------
    871           // End of main loop.
    872           //---------------------------------------------------------------------
    873      }
    874 
    875 
    876      for (int k = 0; k < nfix; k++)
    877           x[fix[k]] = bl[fix[k]];
    878      z      = z1;
    879      if (nupp > 0)
    880           z = z - z2;
    881      printf("\n\nmax |x| =%10.3f", x.infNorm() );
    882      printf("    max |y| =%10.3f", y.infNorm() );
    883      printf("    max |z| =%10.3f", z.infNorm() );
    884      printf("   scaled");
    885 
    886      x.scale(beta);
    887      y.scale(zeta);
    888      z.scale(zeta);   // Unscale x, y, z.
    889 
    890      printf(  "\nmax |x| =%10.3f", x.infNorm() );
    891      printf("    max |y| =%10.3f", y.infNorm() );
    892      printf("    max |z| =%10.3f", z.infNorm() );
    893      printf(" unscaled\n");
    894 
    895      time   = CoinCpuTime() - time;
    896      char str1[100], str2[100];
    897      sprintf(str1, "\nPDitns  =%10g", PDitns );
    898      sprintf(str2, "itns =%10d", CGitns );
    899      //  printf( [str1 " " solver str2] );
    900      printf("    time    =%10.1f\n", time);
    901      /*
     619        //  New version of lsqr
     620
     621        int istop;
     622        dy.clear();
     623        show = false;
     624        info.atolmin = atolmin;
     625        info.r3norm = fmerit; // Must be the 2-norm here.
     626
     627        thisLsqr.do_lsqr(rhs, damp, atol, btol, conlim, itnlim,
     628          show, info, dy, &istop, &itncg, &outfo, precon, Pr);
     629        if (precon)
     630          dy = dy * Pr;
     631
     632        if (!precon && itncg > 999999)
     633          precon = true;
     634
     635        if (istop == 3 || istop == 7) // conlim or itnlim
     636          printf("\n    LSQR stopped early:  istop = //%d", istop);
     637
     638        atolold = outfo.atolold;
     639        atol = outfo.atolnew;
     640        r3ratio = outfo.r3ratio;
     641      } // LSproblem 1
     642
     643      //      grad      = pdxxxmat( Aname,2,m,n,dy );   // grad = A"dy
     644      grad.clear();
     645      matVecMult(2, grad, dy);
     646      for (int k = 0; k < nfix; k++)
     647        grad[fix[k]] = 0; // grad is a work vector
     648      dx = H * (grad - w);
     649
     650    } else {
     651      perror("This LSproblem not yet implemented\n");
     652    }
     653    //-------------------------------------------------------------------
     654
     655    CGitns += itncg;
     656
     657    //-------------------------------------------------------------------
     658    // dx and dy are now known.  Get dx1, dx2, dz1, dz2.
     659    //-------------------------------------------------------------------
     660    for (int k = 0; k < nlow; k++) {
     661      dx1[low[k]] = -rL[low[k]] + dx[low[k]];
     662      dz1[low[k]] = (cL[low[k]] - z1[low[k]] * dx1[low[k]]) / x1[low[k]];
     663    }
     664    for (int k = 0; k < nupp; k++) {
     665      dx2[upp[k]] = -rU[upp[k]] - dx[upp[k]];
     666      dz2[upp[k]] = (cU[upp[k]] - z2[upp[k]] * dx2[upp[k]]) / x2[upp[k]];
     667    }
     668    //-------------------------------------------------------------------
     669    // Find the maximum step.
     670    //--------------------------------------------------------------------
     671    double stepx1 = pdxxxstep(nlow, low, x1, dx1);
     672    double stepx2 = inf;
     673    if (nupp > 0)
     674      stepx2 = pdxxxstep(nupp, upp, x2, dx2);
     675    double stepz1 = pdxxxstep(z1, dz1);
     676    double stepz2 = inf;
     677    if (nupp > 0)
     678      stepz2 = pdxxxstep(z2, dz2);
     679    double stepx = CoinMin(stepx1, stepx2);
     680    double stepz = CoinMin(stepz1, stepz2);
     681    stepx = CoinMin(steptol * stepx, 1.0);
     682    stepz = CoinMin(steptol * stepz, 1.0);
     683    if (stepSame) { // For NLPs, force same step
     684      stepx = CoinMin(stepx, stepz); // (true Newton method)
     685      stepz = stepx;
     686    }
     687
     688    //-------------------------------------------------------------------
     689    // Backtracking linesearch.
     690    //-------------------------------------------------------------------
     691    bool fail = true;
     692    nf = 0;
     693
     694    while (nf < maxf) {
     695      nf = nf + 1;
     696      x = x + stepx * dx;
     697      y = y + stepz * dy;
     698      for (int k = 0; k < nlow; k++) {
     699        x1[low[k]] = x1[low[k]] + stepx * dx1[low[k]];
     700        z1[low[k]] = z1[low[k]] + stepz * dz1[low[k]];
     701      }
     702      for (int k = 0; k < nupp; k++) {
     703        x2[upp[k]] = x2[upp[k]] + stepx * dx2[upp[k]];
     704        z2[upp[k]] = z2[upp[k]] + stepz * dz2[upp[k]];
     705      }
     706      //      [obj,grad,hess] = feval( Fname, (x*beta) );
     707      x.scale(beta);
     708      obj = getObj(x);
     709      getGrad(x, grad);
     710      getHessian(x, H);
     711      x.scale((1.0 / beta));
     712
     713      obj /= theta;
     714      grad = grad * (beta / theta) + d1 * d1 * x;
     715      H = H * (beta2 / theta) + d1 * d1;
     716
     717      //      [r1,r2,rL,rU,Pinf,Dinf] = ...
     718      pdxxxresid1(this, nlow, nupp, nfix, low, upp, fix,
     719        b, bl_elts, bu_elts, d1, d2, grad, rL, rU, x, x1, x2,
     720        y, z1, z2, r1, r2, &Pinf, &Dinf);
     721      //double center, Cinf, Cinf0;
     722      //      [cL,cU,center,Cinf,Cinf0] = ...
     723      pdxxxresid2(mu, nlow, nupp, low, upp, cL, cU, x1, x2, z1, z2,
     724        &center, &Cinf, &Cinf0);
     725      double fmeritnew = pdxxxmerit(nlow, nupp, low, upp, r1, r2, rL, rU, cL, cU);
     726      double step = CoinMin(stepx, stepz);
     727
     728      if (fmeritnew <= (1 - eta * step) * fmerit) {
     729        fail = false;
     730        break;
     731      }
     732
     733      // Merit function didn"t decrease.
     734      // Restore variables to previous values.
     735      // (This introduces a little error, but save lots of space.)
     736
     737      x = x - stepx * dx;
     738      y = y - stepz * dy;
     739      for (int k = 0; k < nlow; k++) {
     740        x1[low[k]] = x1[low[k]] - stepx * dx1[low[k]];
     741        z1[low[k]] = z1[low[k]] - stepz * dz1[low[k]];
     742      }
     743      for (int k = 0; k < nupp; k++) {
     744        x2[upp[k]] = x2[upp[k]] - stepx * dx2[upp[k]];
     745        z2[upp[k]] = z2[upp[k]] - stepz * dz2[upp[k]];
     746      }
     747      // Back-track.
     748      // If it"s the first time,
     749      // make stepx and stepz the same.
     750
     751      if (nf == 1 && stepx != stepz) {
     752        stepx = step;
     753      } else if (nf < maxf) {
     754        stepx = stepx / 2;
     755      }
     756      stepz = stepx;
     757    }
     758
     759    if (fail) {
     760      printf("\n     Linesearch failed (nf too big)");
     761      nfail += 1;
     762    } else {
     763      nfail = 0;
     764    }
     765
     766    //-------------------------------------------------------------------
     767    // Set convergence measures.
     768    //--------------------------------------------------------------------
     769    regx = (d1 * x).twoNorm();
     770    regy = (d2 * y).twoNorm();
     771    regterm = regx * regx + regy * regy;
     772    objreg = obj + 0.5 * regterm;
     773    objtrue = objreg * theta;
     774
     775    bool primalfeas = Pinf <= featol;
     776    bool dualfeas = Dinf <= featol;
     777    bool complementary = Cinf0 <= opttol;
     778    bool enough = PDitns >= 4; // Prevent premature termination.
     779    bool converged = primalfeas & dualfeas & complementary & enough;
     780
     781    //-------------------------------------------------------------------
     782    // Iteration log.
     783    //-------------------------------------------------------------------
     784    char str1[100], str2[100], str3[100], str4[100], str5[100];
     785    sprintf(str1, "\n%3g%5.1f", PDitns, log10(mu));
     786    sprintf(str2, "%8.5f%8.5f", stepx, stepz);
     787    if (stepx < 0.0001 || stepz < 0.0001) {
     788      sprintf(str2, " %6.1e %6.1e", stepx, stepz);
     789    }
     790
     791    sprintf(str3, "%6.1f%6.1f", log10(Pinf), log10(Dinf));
     792    sprintf(str4, "%6.1f%15.7e", log10(Cinf0), objtrue);
     793    sprintf(str5, "%3d%8.1f", nf, center);
     794    if (center > 99999) {
     795      sprintf(str5, "%3d%8.1e", nf, center);
     796    }
     797    printf("%s%s%s%s%s", str1, str2, str3, str4, str5);
     798    if (direct) {
     799      // relax
     800    } else {
     801      printf(" %5.1f%7d%7.3f", log10(atolold), itncg, r3ratio);
     802    }
     803    //-------------------------------------------------------------------
     804    // Test for termination.
     805    //-------------------------------------------------------------------
     806    if (kminor) {
     807      printf("\nStart of next minor itn...\n");
     808      //      keyboard;
     809    }
     810
     811    if (converged) {
     812      printf("\n   Converged");
     813      break;
     814    } else if (PDitns >= maxitn) {
     815      printf("\n   Too many iterations");
     816      inform = 1;
     817      break;
     818    } else if (nfail >= maxfail) {
     819      printf("\n   Too many linesearch failures");
     820      inform = 2;
     821      break;
     822    } else {
     823
     824      // Reduce mu, and reset certain residuals.
     825
     826      double stepmu = CoinMin(stepx, stepz);
     827      stepmu = CoinMin(stepmu, steptol);
     828      double muold = mu;
     829      mu = mu - stepmu * mu;
     830      if (center >= bigcenter)
     831        mu = muold;
     832
     833      // mutrad = mu0*(sum(Xz)/n); // 24 May 1998: Traditional value, but
     834      // mu     = CoinMin(mu,mutrad ); // it seemed to decrease mu too much.
     835
     836      mu = CoinMax(mu, mulast); // 13 Jun 1998: No need for smaller mu.
     837      //      [cL,cU,center,Cinf,Cinf0] = ...
     838      pdxxxresid2(mu, nlow, nupp, low, upp, cL, cU, x1, x2, z1, z2,
     839        &center, &Cinf, &Cinf0);
     840      fmerit = pdxxxmerit(nlow, nupp, low, upp, r1, r2, rL, rU, cL, cU);
     841
     842      // Reduce atol for LSQR (and SYMMLQ).
     843      // NOW DONE AT TOP OF LOOP.
     844
     845      atolold = atol;
     846      // if atol > atol2
     847      //   atolfac = (mu/mufirst)^0.25;
     848      //   atol    = CoinMax( atol*atolfac, atol2 );
     849      // end
     850
     851      // atol = CoinMin( atol, mu );     // 22 Jan 2001: a la Inexact Newton.
     852      // atol = CoinMin( atol, 0.5*mu ); // 30 Jan 2001: A bit tighter
     853
     854      // If the linesearch took more than one function (nf > 1),
     855      // we assume the search direction needed more accuracy
     856      // (though this may be true only for LPs).
     857      // 12 Jun 1998: Ask for more accuracy if nf > 2.
     858      // 24 Nov 2000: Also if the steps are small.
     859      // 30 Jan 2001: Small steps might be ok with warm start.
     860      // 06 Feb 2001: Not necessarily.  Reinstated tests in next line.
     861
     862      if (nf > 2 || CoinMin(stepx, stepz) <= 0.01)
     863        atol = atolold * 0.1;
     864    }
     865    //---------------------------------------------------------------------
     866    // End of main loop.
     867    //---------------------------------------------------------------------
     868  }
     869
     870  for (int k = 0; k < nfix; k++)
     871    x[fix[k]] = bl[fix[k]];
     872  z = z1;
     873  if (nupp > 0)
     874    z = z - z2;
     875  printf("\n\nmax |x| =%10.3f", x.infNorm());
     876  printf("    max |y| =%10.3f", y.infNorm());
     877  printf("    max |z| =%10.3f", z.infNorm());
     878  printf("   scaled");
     879
     880  x.scale(beta);
     881  y.scale(zeta);
     882  z.scale(zeta); // Unscale x, y, z.
     883
     884  printf("\nmax |x| =%10.3f", x.infNorm());
     885  printf("    max |y| =%10.3f", y.infNorm());
     886  printf("    max |z| =%10.3f", z.infNorm());
     887  printf(" unscaled\n");
     888
     889  time = CoinCpuTime() - time;
     890  char str1[100], str2[100];
     891  sprintf(str1, "\nPDitns  =%10g", PDitns);
     892  sprintf(str2, "itns =%10d", CGitns);
     893  //  printf( [str1 " " solver str2] );
     894  printf("    time    =%10.1f\n", time);
     895  /*
    902896     pdxxxdistrib( abs(x),abs(z) );   // Private function
    903897
     
    905899       keyboard;
    906900     */
    907 //-----------------------------------------------------------------------
    908 // End function pdco.m
    909 //-----------------------------------------------------------------------
    910      /*  printf("Solution x values:\n\n");
     901  //-----------------------------------------------------------------------
     902  // End function pdco.m
     903  //-----------------------------------------------------------------------
     904  /*  printf("Solution x values:\n\n");
    911905       for (int k=0; k<n; k++)
    912906         printf(" %d   %e\n", k, x[k]);
    913907     */
    914 // Print distribution
    915      double thresh[9] = { 0.00000001, 0.0000001, 0.000001, 0.00001, 0.0001, 0.001, 0.01, 0.1, 1.00001};
    916      int counts[9] = {0};
    917      for (int ij = 0; ij < n; ij++) {
    918           for (int j = 0; j < 9; j++) {
    919                if(x[ij] < thresh[j]) {
    920                     counts[j] += 1;
    921                     break;
    922                }
    923           }
    924      }
    925      printf ("Distribution of Solution Values\n");
    926      for (int j = 8; j > 1; j--)
    927           printf(" %g  to  %g %d\n", thresh[j-1], thresh[j], counts[j]);
    928      printf("   Less than   %g %d\n", thresh[2], counts[0]);
    929 
    930      return inform;
     908  // Print distribution
     909  double thresh[9] = { 0.00000001, 0.0000001, 0.000001, 0.00001, 0.0001, 0.001, 0.01, 0.1, 1.00001 };
     910  int counts[9] = { 0 };
     911  for (int ij = 0; ij < n; ij++) {
     912    for (int j = 0; j < 9; j++) {
     913      if (x[ij] < thresh[j]) {
     914        counts[j] += 1;
     915        break;
     916      }
     917    }
     918  }
     919  printf("Distribution of Solution Values\n");
     920  for (int j = 8; j > 1; j--)
     921    printf(" %g  to  %g %d\n", thresh[j - 1], thresh[j], counts[j]);
     922  printf("   Less than   %g %d\n", thresh[2], counts[0]);
     923
     924  return inform;
    931925}
    932926// LSQR
    933 void
    934 ClpPdco::lsqr()
    935 {
    936 }
    937 
    938 void ClpPdco::matVecMult( int mode, double* x_elts, double* y_elts)
    939 {
    940      pdcoStuff_->matVecMult(this, mode, x_elts, y_elts);
    941 }
    942 void ClpPdco::matVecMult( int mode, CoinDenseVector<double> &x, double *y_elts)
    943 {
    944      double *x_elts = x.getElements();
    945      matVecMult( mode, x_elts, y_elts);
    946      return;
    947 }
    948 
    949 void ClpPdco::matVecMult( int mode, CoinDenseVector<double> &x, CoinDenseVector<double> &y)
    950 {
    951      double *x_elts = x.getElements();
    952      double *y_elts = y.getElements();
    953      matVecMult( mode, x_elts, y_elts);
    954      return;
    955 }
    956 
    957 void ClpPdco::matVecMult( int mode, CoinDenseVector<double> *x, CoinDenseVector<double> *y)
    958 {
    959      double *x_elts = x->getElements();
    960      double *y_elts = y->getElements();
    961      matVecMult( mode, x_elts, y_elts);
    962      return;
    963 }
    964 void ClpPdco::matPrecon(double delta, double* x_elts, double* y_elts)
    965 {
    966      pdcoStuff_->matPrecon(this, delta, x_elts, y_elts);
    967 }
    968 void ClpPdco::matPrecon(double delta, CoinDenseVector<double> &x, double *y_elts)
    969 {
    970      double *x_elts = x.getElements();
    971      matPrecon(delta, x_elts, y_elts);
    972      return;
    973 }
    974 
    975 void ClpPdco::matPrecon(double delta, CoinDenseVector<double> &x, CoinDenseVector<double> &y)
    976 {
    977      double *x_elts = x.getElements();
    978      double *y_elts = y.getElements();
    979      matPrecon(delta, x_elts, y_elts);
    980      return;
    981 }
    982 
    983 void ClpPdco::matPrecon(double delta, CoinDenseVector<double> *x, CoinDenseVector<double> *y)
    984 {
    985      double *x_elts = x->getElements();
    986      double *y_elts = y->getElements();
    987      matPrecon(delta, x_elts, y_elts);
    988      return;
     927void ClpPdco::lsqr()
     928{
     929}
     930
     931void ClpPdco::matVecMult(int mode, double *x_elts, double *y_elts)
     932{
     933  pdcoStuff_->matVecMult(this, mode, x_elts, y_elts);
     934}
     935void ClpPdco::matVecMult(int mode, CoinDenseVector< double > &x, double *y_elts)
     936{
     937  double *x_elts = x.getElements();
     938  matVecMult(mode, x_elts, y_elts);
     939  return;
     940}
     941
     942void ClpPdco::matVecMult(int mode, CoinDenseVector< double > &x, CoinDenseVector< double > &y)
     943{
     944  double *x_elts = x.getElements();
     945  double *y_elts = y.getElements();
     946  matVecMult(mode, x_elts, y_elts);
     947  return;
     948}
     949
     950void ClpPdco::matVecMult(int mode, CoinDenseVector< double > *x, CoinDenseVector< double > *y)
     951{
     952  double *x_elts = x->getElements();
     953  double *y_elts = y->getElements();
     954  matVecMult(mode, x_elts, y_elts);
     955  return;
     956}
     957void ClpPdco::matPrecon(double delta, double *x_elts, double *y_elts)
     958{
     959  pdcoStuff_->matPrecon(this, delta, x_elts, y_elts);
     960}
     961void ClpPdco::matPrecon(double delta, CoinDenseVector< double > &x, double *y_elts)
     962{
     963  double *x_elts = x.getElements();
     964  matPrecon(delta, x_elts, y_elts);
     965  return;
     966}
     967
     968void ClpPdco::matPrecon(double delta, CoinDenseVector< double > &x, CoinDenseVector< double > &y)
     969{
     970  double *x_elts = x.getElements();
     971  double *y_elts = y.getElements();
     972  matPrecon(delta, x_elts, y_elts);
     973  return;
     974}
     975
     976void ClpPdco::matPrecon(double delta, CoinDenseVector< double > *x, CoinDenseVector< double > *y)
     977{
     978  double *x_elts = x->getElements();
     979  double *y_elts = y->getElements();
     980  matPrecon(delta, x_elts, y_elts);
     981  return;
    989982}
    990983void ClpPdco::getBoundTypes(int *nlow, int *nupp, int *nfix, int **bptrs)
    991984{
    992      *nlow = numberColumns_;
    993      *nupp = *nfix = 0;
    994      int *low_ = (int *)malloc(numberColumns_ * sizeof(int)) ;
    995      for (int k = 0; k < numberColumns_; k++)
    996           low_[k] = k;
    997      bptrs[0] = low_;
    998      return;
    999 }
    1000 
    1001 double ClpPdco::getObj(CoinDenseVector<double> &x)
    1002 {
    1003      return pdcoStuff_->getObj(this, x);
    1004 }
    1005 
    1006 void ClpPdco::getGrad(CoinDenseVector<double> &x, CoinDenseVector<double> &g)
    1007 {
    1008      pdcoStuff_->getGrad(this, x, g);
    1009 }
    1010 
    1011 void ClpPdco::getHessian(CoinDenseVector<double> &x, CoinDenseVector<double> &H)
    1012 {
    1013      pdcoStuff_->getHessian(this, x, H);
    1014 }
     985  *nlow = numberColumns_;
     986  *nupp = *nfix = 0;
     987  int *low_ = (int *)malloc(numberColumns_ * sizeof(int));
     988  for (int k = 0; k < numberColumns_; k++)
     989    low_[k] = k;
     990  bptrs[0] = low_;
     991  return;
     992}
     993
     994double ClpPdco::getObj(CoinDenseVector< double > &x)
     995{
     996  return pdcoStuff_->getObj(this, x);
     997}
     998
     999void ClpPdco::getGrad(CoinDenseVector< double > &x, CoinDenseVector< double > &g)
     1000{
     1001  pdcoStuff_->getGrad(this, x, g);
     1002}
     1003
     1004void ClpPdco::getHessian(CoinDenseVector< double > &x, CoinDenseVector< double > &H)
     1005{
     1006  pdcoStuff_->getHessian(this, x, H);
     1007}
     1008
     1009/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     1010*/
Note: See TracChangeset for help on using the changeset viewer.