# source:trunk/Clp/src/ClpPdco.cpp@2385

Last change on this file since 2385 was 2385, checked in by unxusr, 4 months ago

formatting

• Property svn:keywords set to `Id`
File size: 35.0 KB
Line
1/* \$Id: ClpPdco.cpp 2385 2019-01-06 19:43:06Z unxusr \$ */
5
6/* Pdco algorithm
7
8Method
9
10
11*/
12
13#include "CoinPragma.hpp"
14
15#include <math.h>
16
17#include "CoinDenseVector.hpp"
18#include "ClpPdco.hpp"
19#include "ClpPdcoBase.hpp"
20#include "CoinHelperFunctions.hpp"
21#include "ClpHelperFunctions.hpp"
22#include "ClpLsqr.hpp"
23#include "CoinTime.hpp"
24#include "ClpMessage.hpp"
25#include <cfloat>
26#include <cassert>
27#include <string>
28#include <stdio.h>
29#include <iostream>
30
31int ClpPdco::pdco()
32{
33  printf("Dummy pdco solve\n");
34  return 0;
35}
36// ** Temporary version
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,
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  /*
229     if (d1_->size()==1)
230         d1_->resize(n, d1_->getElements()[0]);  // Allow scalar d1, d2
231     if (d2_->size()==1)
232         d2->resize(m, d2->getElements()[0]);  // to mean dk * unit vector
233      */
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  /*
278     double gamma     = d1->infNorm();
279     double delta     = d2->infNorm();
280     */
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  /*
300       if wait
301         printf("\n\nReview parameters... then type "return"\n")
302         keyboard
303       end
304     */
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),
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
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);
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  /*---------------------------------------------------------------------
460     // Compute primal and dual residuals:
461     // r1 =  b - Aprod(x) - d2*d2*y;
462     // r2 =  grad - Atprod(y) + z2 - z1;
463     //  rL =  bl - x + x1;
464     //  rU =  x + x2 - bu; */
465  //---------------------------------------------------------------------
466  //  [r1,r2,rL,rU,Pinf,Dinf] = ...
467  //      pdxxxresid1( Aname,fix,low,upp, ...
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  /*
531     if kminor
532       printf("\n\nStart of first minor itn...\n");
533       keyboard
534     end
535     */
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        /*
605                        rw(7)        = precon;
606                            info.atolmin = atolmin;
607                            info.r3norm  = fmerit;  // Must be the 2-norm here.
608
609                            [ dy, istop, itncg, outfo ] = ...
610                       pdxxxlsqr( nb,m,"pdxxxlsqrmat",Aname,rw,rhs,damp, ...
611                                  atol,btol,conlim,itnlim,show,info );
612
613
614                        thisLsqr.input->rhs_vec = &rhs;
615                        thisLsqr.input->sol_vec = &dy;
616                        thisLsqr.input->rel_mat_err = atol;
617                        thisLsqr.do_lsqr(this);
618                        */
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
646      for (int k = 0; k < nfix; k++)
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);
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  /*
896     pdxxxdistrib( abs(x),abs(z) );   // Private function
897
898     if (wait)
899       keyboard;
900     */
901  //-----------------------------------------------------------------------
902  // End function pdco.m
903  //-----------------------------------------------------------------------
904  /*  printf("Solution x values:\n\n");
905       for (int k=0; k<n; k++)
906         printf(" %d   %e\n", k, x[k]);
907     */
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;
925}
926// LSQR
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;
982}
983void ClpPdco::getBoundTypes(int *nlow, int *nupp, int *nfix, int **bptrs)
984{
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{