Changeset 2385 for trunk/Clp/src/ClpPdco.cpp
 Timestamp:
 Jan 6, 2019 2:43:06 PM (4 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

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