# source:branches/pre/ClpPdco.cpp@1926

Last change on this file since 1926 was 219, checked in by forrest, 17 years ago

pdco seems to work

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