source: branches/pre/ClpPdco.cpp @ 214

Last change on this file since 214 was 214, checked in by forrest, 16 years ago

More pdco

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