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