source: branches/pre/ClpPdco.cpp @ 215

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

Second try at pdco

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