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

Last change on this file since 1879 was 1665, checked in by lou, 9 years ago

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