# source:stable/1.15/Clp/src/ClpPdco.cpp@1941

Last change on this file since 1941 was 1941, checked in by stefan, 7 years ago

sync with trunk rev 1940

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