source: releases/1.8.0/Clp/src/ClpPdco.cpp @ 2257

Last change on this file since 2257 was 1197, checked in by forrest, 12 years ago

many changes to try and improve performance

File size: 34.6 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#ifdef COIN_DO_PDCO
4
5/* Pdco algorithm
6   
7Method
8
9
10*/
11
12
13
14#include "CoinPragma.hpp"
15
16#include <math.h>
17
18#include "CoinDenseVector.hpp"
19#include "ClpPdco.hpp"
20#include "ClpPdcoBase.hpp"
21#include "CoinHelperFunctions.hpp"
22#include "ClpHelperFunctions.hpp"
23#include "ClpLsqr.hpp"
24#include "CoinTime.hpp"
25#include "ClpMessage.hpp"
26#include <cfloat>
27#include <cassert>
28#include <string>
29#include <stdio.h>
30#include <iostream>
31
32int 
33ClpPdco::pdco()
34{
35  printf("Dummy pdco solve\n");
36  return 0;
37}
38// ** Temporary version
39int 
40ClpPdco::pdco( ClpPdcoBase * stuff, Options &options, Info &info, Outfo &outfo)
41{
42//    D1, D2 are positive-definite diagonal matrices defined from d1, d2.
43//           In particular, d2 indicates the accuracy required for
44//           satisfying each row of Ax = b.
45//
46// D1 and D2 (via d1 and d2) provide primal and dual regularization
47// respectively.  They ensure that the primal and dual solutions
48// (x,r) and (y,z) are unique and bounded.
49//
50// A scalar d1 is equivalent to d1 = ones(n,1), D1 = diag(d1).
51// A scalar d2 is equivalent to d2 = ones(m,1), D2 = diag(d2).
52// Typically, d1 = d2 = 1e-4.
53// These values perturb phi(x) only slightly  (by about 1e-8) and request
54// that A*x = b be satisfied quite accurately (to about 1e-8).
55// Set d1 = 1e-4, d2 = 1 for least-squares problems with bound constraints.
56// The problem is then
57//
58//    minimize    phi(x) + 1/2 norm(d1*x)^2 + 1/2 norm(A*x - b)^2
59//    subject to  bl <= x <= bu.
60//
61// More generally, d1 and d2 may be n and m vectors containing any positive
62// values (preferably not too small, and typically no larger than 1).
63// Bigger elements of d1 and d2 improve the stability of the solver.
64//
65// At an optimal solution, if x(j) is on its lower or upper bound,
66// the corresponding z(j) is positive or negative respectively.
67// If x(j) is between its bounds, z(j) = 0.
68// If bl(j) = bu(j), x(j) is fixed at that value and z(j) may have
69// either sign.
70//
71// Also, r and y satisfy r = D2 y, so that Ax + D2^2 y = b.
72// Thus if d2(i) = 1e-4, the i-th row of Ax = b will be satisfied to
73// approximately 1e-8.  This determines how large d2(i) can safely be.
74//
75//
76// EXTERNAL FUNCTIONS:
77// options         = pdcoSet;                  provided with pdco.m
78// [obj,grad,hess] = pdObj( x );               provided by user
79//               y = pdMat( name,mode,m,n,x ); provided by user if pdMat
80//                                             is a string, not a matrix
81//
82// INPUT ARGUMENTS:
83// pdObj      is a string containing the name of a function pdObj.m
84//            or a function_handle for such a function
85//            such that [obj,grad,hess] = pdObj(x) defines
86//            obj  = phi(x)              : a scalar,
87//            grad = gradient of phi(x)  : an n-vector,
88//            hess = diag(Hessian of phi): an n-vector.
89//         Examples:
90//            If phi(x) is the linear function c"x, pdObj should return
91//               [obj,grad,hess] = [c"*x, c, zeros(n,1)].
92//            If phi(x) is the entropy function E(x) = sum x(j) log x(j),
93//               [obj,grad,hess] = [E(x), log(x)+1, 1./x].
94// pdMat      may be an ifexplicit m x n matrix A (preferably sparse!),
95//            or a string containing the name of a function pdMat.m
96//            or a function_handle for such a function
97//            such that y = pdMat( name,mode,m,n,x )
98//            returns   y = A*x (mode=1)  or  y = A"*x (mode=2).
99//            The input parameter "name" will be the string pdMat.
100// b          is an m-vector.
101// bl         is an n-vector of lower bounds.  Non-existent bounds
102//            may be represented by bl(j) = -Inf or bl(j) <= -1e+20.
103// bu         is an n-vector of upper bounds.  Non-existent bounds
104//            may be represented by bu(j) =  Inf or bu(j) >=  1e+20.
105// d1, d2     may be positive scalars or positive vectors (see above).
106// options    is a structure that may be set and altered by pdcoSet
107//            (type help pdcoSet).
108// x0, y0, z0 provide an initial solution.
109// xsize, zsize are estimates of the biggest x and z at the solution.
110//            They are used to scale (x,y,z).  Good estimates
111//            should improve the performance of the barrier method.
112//
113//
114// OUTPUT ARGUMENTS:
115// x          is the primal solution.
116// y          is the dual solution associated with Ax + D2 r = b.
117// z          is the dual solution associated with bl <= x <= bu.
118// inform = 0 if a solution is found;
119//        = 1 if too many iterations were required;
120//        = 2 if the linesearch failed too often.
121// PDitns     is the number of Primal-Dual Barrier iterations required.
122// CGitns     is the number of Conjugate-Gradient  iterations required
123//            if an iterative solver is used (LSQR).
124// time       is the cpu time used.
125//----------------------------------------------------------------------
126
127// PRIVATE FUNCTIONS:
128//    pdxxxbounds
129//    pdxxxdistrib
130//    pdxxxlsqr
131//    pdxxxlsqrmat
132//    pdxxxmat
133//    pdxxxmerit
134//    pdxxxresid1
135//    pdxxxresid2
136//    pdxxxstep
137//
138// GLOBAL VARIABLES:
139//    global pdDDD1 pdDDD2 pdDDD3
140//
141//
142// NOTES:
143// The matrix A should be reasonably well scaled: norm(A,inf) =~ 1.
144// The vector b and objective phi(x) may be of any size, but ensure that
145// xsize and zsize are reasonably close to norm(x,inf) and norm(z,inf)
146// at the solution.
147//
148// The files defining pdObj  and pdMat
149// must not be called Fname.m or Aname.m!!
150//
151//
152// AUTHOR:
153//    Michael Saunders, Systems Optimization Laboratory (SOL),
154//    Stanford University, Stanford, California, USA.
155//    saunders@stanford.edu
156//
157// CONTRIBUTORS:
158//    Byunggyoo Kim, SOL, Stanford University.
159//    bgkim@stanford.edu
160//
161// DEVELOPMENT:
162// 20 Jun 1997: Original version of pdsco.m derived from pdlp0.m.
163// 29 Sep 2002: Original version of pdco.m  derived from pdsco.m.
164//              Introduced D1, D2 in place of gamma*I, delta*I
165//              and allowed for general bounds bl <= x <= bu.
166// 06 Oct 2002: Allowed for fixed variabes: bl(j) = bu(j) for any j.
167// 15 Oct 2002: Eliminated some work vectors (since m, n might be LARGE).
168//              Modularized residuals, linesearch
169// 16 Oct 2002: pdxxx..., pdDDD... names rationalized.
170//              pdAAA eliminated (global copy of A).
171//              Aname is now used directly as an ifexplicit A or a function.
172//              NOTE: If Aname is a function, it now has an extra parameter.
173// 23 Oct 2002: Fname and Aname can now be function handles.
174// 01 Nov 2002: Bug fixed in feval in pdxxxmat.
175//-----------------------------------------------------------------------
176
177//  global pdDDD1 pdDDD2 pdDDD3
178  double inf = 1.0e30;
179  double eps = 1.0e-15;
180  double atolold, r3ratio, Pinf, Dinf, Cinf, Cinf0;
181 
182  printf("\n   --------------------------------------------------------");
183  printf("\n   pdco.m                            Version of 01 Nov 2002");
184  printf("\n   Primal-dual barrier method to minimize a convex function");
185  printf("\n   subject to linear constraints Ax + r = b,  bl <= x <= bu");
186  printf("\n   --------------------------------------------------------\n");
187
188  int m = numberRows_;
189  int n = numberColumns_;
190  bool ifexplicit = true;
191
192  CoinDenseVector<double> b(m, rhs_); 
193  CoinDenseVector<double> x(n, x_);
194  CoinDenseVector<double> y(m, y_);
195  CoinDenseVector<double> z(n, dj_);
196  //delete old arrays
197  delete [] rhs_;
198  delete [] x_;
199  delete [] y_;
200  delete [] dj_;
201  rhs_=NULL;
202  x_=NULL;
203  y_=NULL;
204  dj_=NULL;
205
206  // Save stuff so available elsewhere
207  pdcoStuff_ = stuff;
208
209  double normb  = b.infNorm();
210  double normx0 = x.infNorm();
211  double normy0 = y.infNorm();
212  double normz0 = z.infNorm();
213
214  printf("\nmax |b | = %8g     max |x0| = %8g", normb , normx0);
215  printf(                "      xsize   = %8g", xsize_);
216  printf("\nmax |y0| = %8g     max |z0| = %8g", normy0, normz0);
217  printf(                "      zsize   = %8g", zsize_);
218
219  //---------------------------------------------------------------------
220  // Initialize.
221  //---------------------------------------------------------------------
222  //true   = 1;
223  //false  = 0;
224  //zn     = zeros(n,1);
225  int nb     = n + m;
226  int nkkt   = nb;
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 * min(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),
369  //    gradient = (beta /theta) grad,
370  //    Hessian  = (beta2/theta) hess.
371  //---------------------------------------------------------------------
372  double beta = xsize_;   if (beta==0) beta = 1;   // beta scales b, x.
373  double zeta = zsize_;   if (zeta==0) zeta = 1;   // zeta scales y, z.
374  double theta  = beta*zeta;                            // theta scales obj.
375  // (theta could be anything, but theta = beta*zeta makes
376  // scaled grad = grad/zeta = 1 approximately if zeta is chosen right.)
377
378  for (int k=0; k<nlow; k++)
379    bl_elts[low[k]] = bl_elts[low[k]]/beta;
380  for (int k=0; k<nupp; k++)
381    bu_elts[upp[k]] = bu_elts[upp[k]]/beta;
382  d1     = d1 * ( beta/sqrt(theta) );
383  d2     = d2 * ( sqrt(theta)/beta );
384
385  double beta2  = beta*beta;
386  b.scale( (1.0/beta) );
387  y.scale( (1.0/zeta) );
388  x.scale( (1.0/beta) ); 
389  z.scale( (1.0/zeta) );
390 
391  //---------------------------------------------------------------------
392  // Initialize vectors that are not fully used if bounds are missing.
393  //---------------------------------------------------------------------
394  CoinDenseVector<double> rL(n, 0.0);
395  CoinDenseVector<double> cL(n, 0.0);
396  CoinDenseVector<double> z1(n, 0.0);
397  CoinDenseVector<double> dx1(n, 0.0);
398  CoinDenseVector<double> dz1(n, 0.0);
399  CoinDenseVector<double> r2(n, 0.0);
400
401  // Assign upper bd regions (dummy if no UBs)
402
403  CoinDenseVector<double> rU(nU, 0.0);
404  CoinDenseVector<double> cU(nU, 0.0);
405  CoinDenseVector<double> x2(nU, 0.0);
406  CoinDenseVector<double> z2(nU, 0.0);
407  CoinDenseVector<double> dx2(nU, 0.0);   
408  CoinDenseVector<double> dz2(nU, 0.0);
409
410  //---------------------------------------------------------------------
411  // Initialize x, y, z, objective, etc.
412  //---------------------------------------------------------------------
413  CoinDenseVector<double> dx(n, 0.0);
414  CoinDenseVector<double> dy(m, 0.0);
415  CoinDenseVector<double> Pr(m);
416  CoinDenseVector<double> D(n);
417  double *D_elts = D.getElements();
418  CoinDenseVector<double> w(n);
419  double *w_elts = w.getElements();
420  CoinDenseVector<double> rhs(m+n);
421   
422
423  //---------------------------------------------------------------------
424  // Pull out the element array pointers for efficiency
425  //---------------------------------------------------------------------
426  double *x_elts  = x.getElements();
427  double *x2_elts = x2.getElements();
428  double *z_elts  = z.getElements();
429  double *z1_elts = z1.getElements();
430  double *z2_elts = z2.getElements();
431
432  for (int k=0; k<nlow; k++){
433    x_elts[low[k]]  = max( x_elts[low[k]], bl[low[k]]);
434    x1_elts[low[k]] = max( x_elts[low[k]] - bl[low[k]], x0min  );
435    z1_elts[low[k]] = max( z_elts[low[k]], z0min  );
436  }
437  for (int k=0; k<nupp; k++){
438    x_elts[upp[k]]  = min( x_elts[upp[k]], bu[upp[k]]);
439    x2_elts[upp[k]] = max(bu[upp[k]] -  x_elts[upp[k]], x0min  );
440    z2_elts[upp[k]] = max(-z_elts[upp[k]], z0min  );
441  }
442  //////////////////// Assume hessian is diagonal. //////////////////////
443
444//  [obj,grad,hess] = feval( Fname, (x*beta) );
445  x.scale(beta);
446  double obj = getObj(x);
447  CoinDenseVector<double> grad(n);
448  getGrad(x, grad);
449  CoinDenseVector<double> H(n);
450  getHessian(x ,H);
451  x.scale((1.0/beta));
452 
453  double * g_elts=grad.getElements();
454  double * H_elts=H.getElements();
455
456  obj /= theta;                       // Scaled obj.
457  grad = grad*(beta/theta) + (d1*d1)*x; // grad includes x regularization.
458  H  = H*(beta2/theta) + (d1*d1)      ;     // H    includes x regularization.
459 
460
461  /*---------------------------------------------------------------------
462  // Compute primal and dual residuals:
463  // r1 =  b - Aprod(x) - d2*d2*y;
464  // r2 =  grad - Atprod(y) + z2 - z1; 
465  //  rL =  bl - x + x1;
466  //  rU =  x + x2 - bu; */
467  //---------------------------------------------------------------------
468  //  [r1,r2,rL,rU,Pinf,Dinf] = ...
469  //      pdxxxresid1( Aname,fix,low,upp, ...
470  //                   b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 );
471  pdxxxresid1( this, nlow, nupp, nfix, low, upp, fix, 
472               b,bl_elts,bu_elts,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2,
473               r1, r2, &Pinf, &Dinf);
474  //---------------------------------------------------------------------
475  // Initialize mu and complementarity residuals:
476  //    cL   = mu*e - X1*z1.
477  //    cU   = mu*e - X2*z2.
478  //
479  // 25 Jan 2001: Now that b and obj are scaled (and hence x,y,z),
480  //              we should be able to use mufirst = mu0 (absolute value).
481  //              0.1 worked poorly on StarTest1 with x0min = z0min = 0.1.
482  // 29 Jan 2001: We might as well use mu0 = x0min * z0min;
483  //              so that most variables are centered after a warm start.
484  // 29 Sep 2002: Use mufirst = mu0*(x0min * z0min),
485  //              regarding mu0 as a scaling of the initial center.
486  //---------------------------------------------------------------------
487  //  double mufirst = mu0*(x0min * z0min);
488  double mufirst = mu0;   // revert to absolute value
489  double mulast  = 0.1 * opttol;
490  mulast  = min( mulast, mufirst );
491  double mu      = mufirst;
492  double center,  fmerit;
493  pdxxxresid2( mu, nlow, nupp, low,upp, cL, cU, x1, x2,
494                        z1, z2, &center, &Cinf, &Cinf0 );
495  fmerit = pdxxxmerit(nlow, nupp, low, upp, r1, r2, rL, rU, cL, cU );
496
497  // Initialize other things.
498
499  bool  precon   = true; 
500  double PDitns    = 0;
501  bool converged = false;
502  double atol      = atol1;
503  atol2     = max( atol2, atolmin );
504  atolmin   = atol2;
505  //  pdDDD2    = d2;    // Global vector for diagonal matrix D2
506
507  //  Iteration log.
508
509  double stepx   = 0;
510  double stepz   = 0;
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 = max(Pinf,   max(Dinf,  Cinf));
555    atol   = min(atol,  r3norm*0.1);
556    atol   = max(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, itn;
626        dy.clear();
627        show = false;
628        info.atolmin = atolmin;
629        info.r3norm  = fmerit;  // Must be the 2-norm here.
630
631        thisLsqr.do_lsqr( rhs, damp, atol, btol, conlim, itnlim, 
632                      show, info, dy , &istop, &itncg, &outfo, precon, Pr);
633        if (precon)
634                  dy = dy*Pr;
635
636        if (!precon && itncg>999999)
637          precon=true;
638
639        if (istop == 3  ||  istop == 7 )  // conlim or itnlim
640          printf("\n    LSQR stopped early:  istop = //%d", istop);
641       
642
643        atolold   = outfo.atolold;
644        atol      = outfo.atolnew;
645        r3ratio   = outfo.r3ratio;
646      }// LSproblem 1
647
648      //      grad      = pdxxxmat( Aname,2,m,n,dy );   // grad = A"dy
649      grad.clear();
650      matVecMult(2, grad, dy);
651      for (int k=0; k<nfix; k++)
652        grad[fix[k]] = 0;                            // grad is a work vector
653      dx = H * (grad - w);
654   
655    }else{
656      perror( "This LSproblem not yet implemented\n" );
657    }
658    //-------------------------------------------------------------------
659
660    CGitns += itncg;
661
662    //-------------------------------------------------------------------
663    // dx and dy are now known.  Get dx1, dx2, dz1, dz2.
664    //-------------------------------------------------------------------
665    for (int k=0; k<nlow; k++){   
666      dx1[low[k]] = - rL[low[k]] + dx[low[k]];
667      dz1[low[k]] =  (cL[low[k]] - z1[low[k]]*dx1[low[k]]) / x1[low[k]];
668    }
669    for (int k=0; k<nupp; k++){   
670      dx2[upp[k]] = - rU[upp[k]] - dx[upp[k]];
671      dz2[upp[k]] =  (cU[upp[k]] - z2[upp[k]]*dx2[upp[k]]) / x2[upp[k]];
672    }
673    //-------------------------------------------------------------------
674    // Find the maximum step.
675    //--------------------------------------------------------------------
676    double stepx1 = pdxxxstep(nlow, low, x1, dx1 );
677    double stepx2 = inf;
678    if (nupp > 0)
679      stepx2 = pdxxxstep(nupp, upp, x2, dx2 );
680    double stepz1 = pdxxxstep( z1     , dz1      );
681    double stepz2 = inf;
682    if (nupp > 0)
683      stepz2 = pdxxxstep( z2     , dz2      );
684    double stepx  = min( stepx1, stepx2 );
685    double stepz  = min( stepz1, stepz2 );
686    stepx  = min( steptol*stepx, 1.0 );
687    stepz  = min( steptol*stepz, 1.0 );
688    if (stepSame){                   // For NLPs, force same step
689      stepx = min( stepx, stepz );   // (true Newton method)
690      stepz = stepx;
691    }
692   
693    //-------------------------------------------------------------------
694    // Backtracking linesearch.
695    //-------------------------------------------------------------------
696    bool fail     =  true;
697    nf       =  0;
698
699    while (nf < maxf){
700      nf      = nf + 1;
701      x       = x        +  stepx * dx;
702      y       = y        +  stepz * dy;
703      for (int k=0; k<nlow; k++){   
704        x1[low[k]] = x1[low[k]]  +  stepx * dx1[low[k]];
705        z1[low[k]] = z1[low[k]]  +  stepz * dz1[low[k]];
706      }
707      for (int k=0; k<nupp; k++){   
708        x2[upp[k]] = x2[upp[k]]  +  stepx * dx2[upp[k]];
709        z2[upp[k]] = z2[upp[k]]  +  stepz * dz2[upp[k]];
710      }
711      //      [obj,grad,hess] = feval( Fname, (x*beta) );
712      x.scale(beta);
713      obj = getObj(x);
714      getGrad(x, grad);
715      getHessian(x, H);
716      x.scale((1.0/beta));
717     
718      obj        /= theta;
719      grad       = grad*(beta /theta)  +  d1*d1*x;
720      H          = H*(beta2/theta)  +  d1*d1;
721     
722      //      [r1,r2,rL,rU,Pinf,Dinf] = ...
723      pdxxxresid1( this, nlow, nupp, nfix, low, upp, fix,
724                   b,bl_elts,bu_elts,d1,d2,grad,rL,rU,x,x1,x2,
725                   y,z1,z2, r1, r2, &Pinf, &Dinf );
726      //double center, Cinf, Cinf0;
727      //      [cL,cU,center,Cinf,Cinf0] = ...
728                  pdxxxresid2( mu, nlow, nupp, low,upp,cL,cU,x1,x2,z1,z2,
729                               &center, &Cinf, &Cinf0);
730      double fmeritnew = pdxxxmerit(nlow, nupp, low,upp,r1,r2,rL,rU,cL,cU );
731      double step      = min( 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  = min( stepx , stepz   );
832      stepmu  = min( 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     = min(mu,mutrad ); // it seemed to decrease mu too much.
840
841      mu      = max(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    = max( atol*atolfac, atol2 );
854      // end
855
856      // atol = min( atol, mu );     // 22 Jan 2001: a la Inexact Newton.
857      // atol = min( 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  ||  min( 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);   y.scale(zeta);   z.scale(zeta);   // Unscale x, y, z.
887
888  printf(  "\nmax |x| =%10.3f", x.infNorm() );
889  printf("    max |y| =%10.3f", y.infNorm() );
890  printf("    max |z| =%10.3f", z.infNorm() );
891  printf(" unscaled\n");
892
893  time   = CoinCpuTime() - time;
894  char str1[100],str2[100];
895  sprintf(str1, "\nPDitns  =%10g", PDitns );
896  sprintf(str2, "itns =%10d", CGitns );
897  //  printf( [str1 " " solver str2] );
898  printf("    time    =%10.1f\n", time);
899  /*
900  pdxxxdistrib( abs(x),abs(z) );   // Private function
901 
902  if (wait)
903    keyboard;
904  */
905//-----------------------------------------------------------------------
906// End function pdco.m
907//-----------------------------------------------------------------------
908/*  printf("Solution x values:\n\n");
909  for (int k=0; k<n; k++)
910    printf(" %d   %e\n", k, x[k]);
911*/
912// Print distribution
913  float thresh[9]={ 0.00000001, 0.0000001, 0.000001,0.00001, 0.0001, 0.001, 0.01, 0.1, 1.00001};
914  int counts[9]={0};
915  for (int ij=0; ij<n; ij++){
916    for (int j=0; j<9; j++){
917      if(x[ij] < thresh[j]){
918        counts[j] += 1;
919        break;
920      }
921    }
922  }
923  printf ("Distribution of Solution Values\n");
924  for (int j=8; j>1; j--)
925    printf(" %f  to  %f %d\n",thresh[j-1],thresh[j],counts[j]);
926  printf("   Less than   %f %d\n",thresh[2],counts[0]);
927
928  return 0;
929}
930// LSQR
931void 
932ClpPdco::lsqr()
933{
934}
935
936void ClpPdco::matVecMult( int mode, double* x_elts, double* y_elts){
937  pdcoStuff_->matVecMult(this,mode,x_elts,y_elts);
938}
939void ClpPdco::matVecMult( int mode, CoinDenseVector<double> &x, double *y_elts){
940  double *x_elts = x.getElements();
941  matVecMult( mode, x_elts, y_elts);
942  return;
943}
944
945void ClpPdco::matVecMult( int mode, CoinDenseVector<double> &x, CoinDenseVector<double> &y){
946  double *x_elts = x.getElements();
947  double *y_elts = y.getElements();
948  matVecMult( mode, x_elts, y_elts);
949  return;
950}
951
952void ClpPdco::matVecMult( int mode, CoinDenseVector<double> *x, CoinDenseVector<double> *y){
953  double *x_elts = x->getElements();
954  double *y_elts = y->getElements();
955  matVecMult( mode, x_elts, y_elts);
956  return;
957}
958void ClpPdco::matPrecon(double delta, double* x_elts, double* y_elts){
959  pdcoStuff_->matPrecon(this,delta,x_elts,y_elts);
960}
961void ClpPdco::matPrecon(double delta, CoinDenseVector<double> &x, double *y_elts){
962  double *x_elts = x.getElements();
963  matPrecon(delta, x_elts, y_elts);
964  return;
965}
966   
967void ClpPdco::matPrecon(double delta, CoinDenseVector<double> &x, CoinDenseVector<double> &y){
968  double *x_elts = x.getElements();
969  double *y_elts = y.getElements();
970  matPrecon(delta, x_elts, y_elts);
971  return;
972}
973
974void ClpPdco::matPrecon(double delta, CoinDenseVector<double> *x, CoinDenseVector<double> *y){
975  double *x_elts = x->getElements();
976  double *y_elts = y->getElements();
977  matPrecon(delta, x_elts, y_elts);
978  return;
979}
980void ClpPdco::getBoundTypes(int *nlow, int *nupp, int *nfix, int **bptrs){
981  *nlow = numberColumns_;
982  *nupp = *nfix = 0;
983  int *low_ = (int *)malloc(numberColumns_*sizeof(int)) ;
984  for (int k=0; k <numberColumns_; k++)
985    low_[k] = k;
986  bptrs[0]= low_;
987  return;
988}
989
990double ClpPdco::getObj(CoinDenseVector<double> &x){
991  return pdcoStuff_->getObj(this, x);
992}
993
994void ClpPdco::getGrad(CoinDenseVector<double> &x, CoinDenseVector<double> &g){
995  pdcoStuff_->getGrad(this, x,g);
996}
997
998void ClpPdco::getHessian(CoinDenseVector<double> &x, CoinDenseVector<double> &H){
999  pdcoStuff_->getHessian(this, x,H);
1000}
1001#endif
Note: See TracBrowser for help on using the repository browser.