 Timestamp:
 Oct 3, 2003 11:41:53 AM (17 years ago)
 Location:
 branches/pre
 Files:

 8 added
 6 edited
Legend:
 Unmodified
 Added
 Removed

branches/pre/ClpInterior.cpp
r214 r215 11 11 #include "CoinHelperFunctions.hpp" 12 12 #include "ClpInterior.hpp" 13 #include "ClpPackedMatrix.hpp" 14 #include "CoinIndexedVector.hpp" 13 #include "ClpMatrixBase.hpp" 14 #include "ClpLsqr.hpp" 15 #include "ClpPdcoBase.hpp" 16 #include "CoinDenseVector.hpp" 15 17 #include "ClpMessage.hpp" 16 18 #include "ClpLinearObjective.hpp" … … 29 31 sumDualInfeasibilities_(0.0), 30 32 sumPrimalInfeasibilities_(0.0), 33 xsize_(0.0), 34 zsize_(0.0), 31 35 lower_(NULL), 32 36 rowLowerWork_(NULL), … … 35 39 rowUpperWork_(NULL), 36 40 columnUpperWork_(NULL), 37 cost_(NULL) 41 cost_(NULL), 42 rhs_(NULL), 43 x_(NULL), 44 y_(NULL), 45 dj_(NULL), 46 lsqrObject_(NULL), 47 pdcoStuff_(NULL) 38 48 { 39 49 solveType_=2; // say interior based life form … … 42 52 // Subproblem constructor 43 53 ClpInterior::ClpInterior ( const ClpModel * rhs, 44 45 46 54 int numberRows, const int * whichRow, 55 int numberColumns, const int * whichColumn, 56 bool dropNames, bool dropIntegers) 47 57 : ClpModel(rhs, numberRows, whichRow, 48 58 numberColumns,whichColumn,dropNames,dropIntegers), 49 largestPrimalError_(0.0), 50 largestDualError_(0.0), 51 sumDualInfeasibilities_(0.0), 52 sumPrimalInfeasibilities_(0.0), 53 lower_(NULL), 54 rowLowerWork_(NULL), 55 columnLowerWork_(NULL), 56 upper_(NULL), 57 rowUpperWork_(NULL), 58 columnUpperWork_(NULL), 59 cost_(NULL) 59 largestPrimalError_(0.0), 60 largestDualError_(0.0), 61 sumDualInfeasibilities_(0.0), 62 sumPrimalInfeasibilities_(0.0), 63 xsize_(0.0), 64 zsize_(0.0), 65 lower_(NULL), 66 rowLowerWork_(NULL), 67 columnLowerWork_(NULL), 68 upper_(NULL), 69 rowUpperWork_(NULL), 70 columnUpperWork_(NULL), 71 cost_(NULL), 72 rhs_(NULL), 73 x_(NULL), 74 y_(NULL), 75 dj_(NULL), 76 lsqrObject_(NULL), 77 pdcoStuff_(NULL) 60 78 { 61 79 solveType_=2; // say interior based life form … … 86 104 sumDualInfeasibilities_(0.0), 87 105 sumPrimalInfeasibilities_(0.0), 106 xsize_(0.0), 107 zsize_(0.0), 88 108 lower_(NULL), 89 109 rowLowerWork_(NULL), … … 92 112 rowUpperWork_(NULL), 93 113 columnUpperWork_(NULL), 94 cost_(NULL) 114 cost_(NULL), 115 rhs_(NULL), 116 x_(NULL), 117 y_(NULL), 118 dj_(NULL), 119 lsqrObject_(NULL), 120 pdcoStuff_(NULL) 95 121 { 96 122 gutsOfDelete(); … … 105 131 sumDualInfeasibilities_(0.0), 106 132 sumPrimalInfeasibilities_(0.0), 133 xsize_(0.0), 134 zsize_(0.0), 107 135 lower_(NULL), 108 136 rowLowerWork_(NULL), … … 111 139 rowUpperWork_(NULL), 112 140 columnUpperWork_(NULL), 113 cost_(NULL) 141 cost_(NULL), 142 rhs_(NULL), 143 x_(NULL), 144 y_(NULL), 145 dj_(NULL), 146 lsqrObject_(NULL), 147 pdcoStuff_(NULL) 114 148 { 115 149 solveType_=2; // say interior based life form … … 137 171 //cost_ = ClpCopyOfArray(rhs.cost_,2*(numberColumns_+numberRows_)); 138 172 cost_ = ClpCopyOfArray(rhs.cost_,numberColumns_); 173 rhs_ = ClpCopyOfArray(rhs.rhs_,numberRows_); 174 x_ = ClpCopyOfArray(rhs.x_,numberColumns_); 175 y_ = ClpCopyOfArray(rhs.y_,numberRows_); 176 dj_ = ClpCopyOfArray(rhs.dj_,numberColumns_); 177 lsqrObject_= new ClpLsqr(*rhs.lsqrObject_); 178 pdcoStuff_ = rhs.pdcoStuff_>clone(); 139 179 largestPrimalError_ = rhs.largestPrimalError_; 140 180 largestDualError_ = rhs.largestDualError_; 141 181 sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_; 142 182 sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_; 183 xsize_ = rhs.xsize_; 184 zsize_ = rhs.zsize_; 143 185 solveType_=rhs.solveType_; 144 186 } … … 157 199 delete [] cost_; 158 200 cost_=NULL; 201 delete [] rhs_; 202 rhs_ = NULL; 203 delete [] x_; 204 x_ = NULL; 205 delete [] y_; 206 y_ = NULL; 207 delete [] dj_; 208 dj_ = NULL; 209 delete lsqrObject_; 210 lsqrObject_ = NULL; 211 delete pdcoStuff_; 212 pdcoStuff_=NULL; 159 213 } 160 214 bool … … 402 456 // ** Temporary version 403 457 int 404 ClpInterior::pdco( Lsqr *lsqr, Options &options, Info &info, Outfo &outfo)405 { 406 return ((ClpPdco *) this)>pdco( lsqr,options,info,outfo);407 } 458 ClpInterior::pdco( ClpPdcoBase * stuff, Options &options, Info &info, Outfo &outfo) 459 { 460 return ((ClpPdco *) this)>pdco(stuff,options,info,outfo); 461 } 
branches/pre/ClpPdco.cpp
r214 r215 1 // Copyright (C) 200 2, International Business Machines1 // Copyright (C) 2003, International Business Machines 2 2 // Corporation and others. All Rights Reserved. 3 3 … … 15 15 #include <math.h> 16 16 17 #include "CoinDenseVector.hpp" 18 #include "ClpPdco.hpp" 19 #include "ClpPdcoBase.hpp" 17 20 #include "CoinHelperFunctions.hpp" 18 #include "ClpPdco.hpp"19 #include "CoinDenseVector.hpp"20 21 #include "ClpHelperFunctions.hpp" 21 #include "CoinPackedMatrix.hpp" 22 #include "ClpLsqr.hpp" 23 #include "CoinTime.hpp" 22 24 #include "ClpMessage.hpp" 23 25 #include <cfloat> … … 35 37 // ** Temporary version 36 38 int 37 ClpPdco::pdco( Lsqr *lsqr, Options &options, Info &info, Outfo &outfo)39 ClpPdco::pdco( ClpPdcoBase * stuff, Options &options, Info &info, Outfo &outfo) 38 40 { 39 41 // D1, D2 are positivedefinite diagonal matrices defined from d1, d2. … … 173 175 174 176 // global pdDDD1 pdDDD2 pdDDD3 175 realinf = 1.0e30;176 realeps = 1.0e15;177 realatolold, r3ratio, Pinf, Dinf, Cinf, Cinf0;177 double inf = 1.0e30; 178 double eps = 1.0e15; 179 double atolold, r3ratio, Pinf, Dinf, Cinf, Cinf0; 178 180 179 181 printf("\n "); … … 183 185 printf("\n \n"); 184 186 185 int m = model>rowSize();186 int n = model>colSize();187 int m = numberRows_; 188 int n = numberColumns_; 187 189 bool ifexplicit = true; 188 190 189 CoinDenseVector b(m, model>rhs);190 CoinDenseVector x(n, model>x);191 CoinDenseVector y(m, model>y);192 CoinDenseVector z(n, model>dj);191 CoinDenseVector b(m, rhs_); 192 CoinDenseVector x(n, x_); 193 CoinDenseVector y(m, y_); 194 CoinDenseVector z(n, dj_); 193 195 //delete old arrays 194 delete model>rhs; 195 delete model>x; 196 delete model>y; 197 delete model>dj; 198 199 real normb = b.infNorm(); 200 real normx0 = x.infNorm(); 201 real normy0 = y.infNorm(); 202 real normz0 = z.infNorm(); 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(); 203 209 204 210 printf("\nmax b  = %8g max x0 = %8g", normb , normx0); 205 printf( " xsize = %8g", model>xsize);211 printf( " xsize = %8g", xsize_); 206 212 printf("\nmax y0 = %8g max z0 = %8g", normy0, normz0); 207 printf( " zsize = %8g", model>zsize);213 printf( " zsize = %8g", zsize_); 208 214 209 215 // … … 221 227 // 222 228 /* 223 if ( model>d1>size()==1)224 d1 >resize(n, d1>getElements()[0]); // Allow scalar d1, d2225 if ( model>d2>size()==1)229 if (d1_>size()==1) 230 d1_>resize(n, d1_>getElements()[0]); // Allow scalar d1, d2 231 if (d2_>size()==1) 226 232 d2>resize(m, d2>getElements()[0]); // to mean dk * unit vector 227 233 */ 228 real d1 = model>d1; 229 real d2 = model>d2; 234 assert (stuff>sizeD1()==1); 235 double d1 = stuff>getD1(); 236 double d2 = stuff>getD2(); 230 237 231 238 // … … 233 240 // 234 241 int maxitn = options.MaxIter; 235 realfeatol = options.FeaTol;236 realopttol = options.OptTol;237 realsteptol = options.StepTol;242 double featol = options.FeaTol; 243 double opttol = options.OptTol; 244 double steptol = options.StepTol; 238 245 int stepSame = 1; /* options.StepSame; // 1 means stepx == stepz */ 239 realx0min = options.x0min;240 realz0min = options.z0min;241 realmu0 = options.mu0;246 double x0min = options.x0min; 247 double z0min = options.z0min; 248 double mu0 = options.mu0; 242 249 int LSproblem = options.LSproblem; // See below 243 250 int LSmethod = options.LSmethod; // 1=Cholesky 2=QR 3=LSQR 244 251 int itnlim = options.LSQRMaxIter * min(m,n); 245 realatol1 = options.LSQRatol1; // Initial atol246 realatol2 = options.LSQRatol2; // Smallest atol,unless atol1 is smaller247 realconlim = options.LSQRconlim;252 double atol1 = options.LSQRatol1; // Initial atol 253 double atol2 = options.LSQRatol2; // Smallest atol,unless atol1 is smaller 254 double conlim = options.LSQRconlim; 248 255 int wait = options.wait; 249 256 … … 259 266 // 260 267 int kminor = 0; // 1 stops after each iteration 261 realeta = 1e4; // Linesearch tolerance for "sufficient descent"262 realmaxf = 10; // Linesearch backtrack limit (function evaluations)263 realmaxfail = 1; // Linesearch failure limit (consecutive iterations)264 realbigcenter = 1e+3; // mu is reduced if center < bigcenter.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. 265 272 266 273 // Parameters for LSQR. 267 realatolmin = eps; // Smallest atol if linesearch backtracks268 realbtol = 0; // Should be small (zero is ok)269 realshow = false; // Controls lsqr iteration log274 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 270 277 /* 271 realgamma = d1>infNorm();272 realdelta = d2>infNorm();278 double gamma = d1>infNorm(); 279 double delta = d2>infNorm(); 273 280 */ 274 realgamma = d1;275 realdelta = d2;281 double gamma = d1; 282 double delta = d2; 276 283 277 284 printf("\n\nx0min = %8g featol = %8.1e", x0min, featol); … … 302 309 // All parameters have now been set. 303 310 // 304 //real time = cputime(); 305 double time = now(); 311 double time = CoinCpuTime(); 306 312 bool useChol = (LSmethod == 1); 307 313 bool useQR = (LSmethod == 2); … … 317 323 int nlow, nupp, nfix; 318 324 int *bptrs[3] = {0}; 319 model>getBoundTypes(&nlow, &nupp, &nfix, bptrs );325 getBoundTypes(&nlow, &nupp, &nfix, bptrs ); 320 326 int *low = bptrs[0]; 321 327 int *upp = bptrs[1]; … … 329 335 // 330 336 331 CoinDenseVector bl(n, model>colLower); 332 real *bl_elts = bl.getElements(); 333 CoinDenseVector bu(nU, model>colUpper); // this is dummy if no UB 334 real *bu_elts = bu.getElements(); 335 delete model>colLower; 336 delete model>colUpper; 337 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(); 337 341 338 342 CoinDenseVector r1(m, 0.0); 339 real*r1_elts = r1.getElements();343 double *r1_elts = r1.getElements(); 340 344 CoinDenseVector x1(n, 0.0); 341 real*x1_elts = x1.getElements();345 double *x1_elts = x1.getElements(); 342 346 343 347 if (nfix > 0){ 344 348 for (int k=0; k<nfix; k++) 345 349 x1_elts[fix[k]] = bl[fix[k]]; 346 m odel>matVecMult(1, r1, x1);350 matVecMult(1, r1, x1); 347 351 b = b  r1; 348 352 // At some stage, might want to look at normfix = norm(r1,inf); … … 362 366 // Hessian = (beta2/theta) hess. 363 367 // 364 real beta = model>xsize; if (beta==0) beta = 1; // beta scales b, x.365 real zeta = model>zsize; if (zeta==0) zeta = 1; // zeta scales y, z.366 realtheta = beta*zeta; // theta scales obj.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. 367 371 // (theta could be anything, but theta = beta*zeta makes 368 372 // scaled grad = grad/zeta = 1 approximately if zeta is chosen right.) … … 375 379 d2 = d2 * ( sqrt(theta)/beta ); 376 380 377 realbeta2 = beta*beta;381 double beta2 = beta*beta; 378 382 b.scale( (1.0/beta) ); 379 383 y.scale( (1.0/zeta) ); … … 407 411 CoinDenseVector Pr(m); 408 412 CoinDenseVector D(n); 409 real*D_elts = D.getElements();413 double *D_elts = D.getElements(); 410 414 CoinDenseVector w(n); 411 real*w_elts = w.getElements();415 double *w_elts = w.getElements(); 412 416 CoinDenseVector rhs(m+n); 413 417 … … 416 420 // Pull out the element array pointers for efficiency 417 421 // 418 real*x_elts = x.getElements();419 real*x2_elts = x2.getElements();420 real*z_elts = z.getElements();421 real*z1_elts = z1.getElements();422 real*z2_elts = z2.getElements();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(); 423 427 424 428 for (int k=0; k<nlow; k++){ … … 436 440 // [obj,grad,hess] = feval( Fname, (x*beta) ); 437 441 x.scale(beta); 438 real obj = model>getObj(x);442 double obj = getObj(x); 439 443 CoinDenseVector grad(n); 440 model>getGrad(x, grad);444 getGrad(x, grad); 441 445 CoinDenseVector H(n); 442 model>getHessian(x ,H);446 getHessian(x ,H); 443 447 x.scale((1.0/beta)); 444 448 445 real* g_elts=grad.getElements();446 real* H_elts=H.getElements();449 double * g_elts=grad.getElements(); 450 double * H_elts=H.getElements(); 447 451 448 452 obj /= theta; // Scaled obj. … … 461 465 // pdxxxresid1( Aname,fix,low,upp, ... 462 466 // b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 ); 463 pdxxxresid1( model, nlow, nupp, nfix, low, upp, fix,467 pdxxxresid1( this, nlow, nupp, nfix, low, upp, fix, 464 468 b,bl_elts,bu_elts,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2, 465 469 r1, r2, &Pinf, &Dinf); … … 477 481 // regarding mu0 as a scaling of the initial center. 478 482 // 479 // realmufirst = mu0*(x0min * z0min);480 realmufirst = mu0; // revert to absolute value481 realmulast = 0.1 * opttol;483 // double mufirst = mu0*(x0min * z0min); 484 double mufirst = mu0; // revert to absolute value 485 double mulast = 0.1 * opttol; 482 486 mulast = min( mulast, mufirst ); 483 realmu = mufirst;484 realcenter, fmerit;487 double mu = mufirst; 488 double center, fmerit; 485 489 pdxxxresid2( mu, nlow, nupp, low,upp, cL, cU, x1, x2, 486 490 z1, z2, ¢er, &Cinf, &Cinf0 ); … … 490 494 491 495 bool precon = true; 492 realPDitns = 0;496 double PDitns = 0; 493 497 bool converged = false; 494 realatol = atol1;498 double atol = atol1; 495 499 atol2 = max( atol2, atolmin ); 496 500 atolmin = atol2; … … 499 503 // Iteration log. 500 504 501 realstepx = 0;502 realstepz = 0;505 double stepx = 0; 506 double stepz = 0; 503 507 int nf = 0; 504 508 int itncg = 0; … … 513 517 } 514 518 515 realregx = (d1*x).twoNorm();516 realregy = (d2*y).twoNorm();519 double regx = (d1*x).twoNorm(); 520 double regy = (d2*y).twoNorm(); 517 521 // regterm = twoNorm(d1.*x)^2 + norm(d2.*y)^2; 518 realregterm = regx*regx + regy*regy;519 realobjreg = obj + 0.5*regterm;520 realobjtrue = objreg * theta;522 double regterm = regx*regx + regy*regy; 523 double objreg = obj + 0.5*regterm; 524 double objtrue = objreg * theta; 521 525 522 526 printf("\n%3g ", PDitns ); … … 533 537 // Main loop. 534 538 // 539 // Lsqr 540 ClpLsqr thisLsqr(this); 535 541 536 542 // while (converged) { … … 543 549 // Now that starting conditions are better, go back to 0.1. 544 550 545 realr3norm = max(Pinf, max(Dinf, Cinf));551 double r3norm = max(Pinf, max(Dinf, Cinf)); 546 552 atol = min(atol, r3norm*0.1); 547 553 atol = max(atol, atolmin ); … … 581 587 for (int k=0; k<n; k++) 582 588 D_elts[k]= sqrt(H_elts[k]); 583 584 lsqr>diag1 = D_elts; 585 lsqr>diag2 = d2; 589 thisLsqr.setDiag1(D_elts); 590 thisLsqr.diag2_ = d2; 586 591 587 592 if (direct){ … … 593 598 for (int k=0; k<m; k++) 594 599 rhs[n+k] = r1_elts[k]*(1.0/d2); 595 realdamp = 0;600 double damp = 0; 596 601 597 602 if (precon){ // Construct diagonal preconditioner for LSQR 598 m odel>matPrecon(d2, Pr, D);603 matPrecon(d2, Pr, D); 599 604 } 600 605 /* … … 608 613 609 614 610 lsqr>input>rhs_vec = &rhs;611 lsqr>input>sol_vec = &dy;612 lsqr>input>rel_mat_err = atol;613 lsqr>do_lsqr(model);615 thisLsqr.input>rhs_vec = &rhs; 616 thisLsqr.input>sol_vec = &dy; 617 thisLsqr.input>rel_mat_err = atol; 618 thisLsqr.do_lsqr(this); 614 619 */ 615 620 // New version of lsqr … … 621 626 info.r3norm = fmerit; // Must be the 2norm here. 622 627 623 lsqr>do_lsqr(model,rhs, damp, atol, btol, conlim, itnlim,628 thisLsqr.do_lsqr( rhs, damp, atol, btol, conlim, itnlim, 624 629 show, info, dy , &istop, &itncg, &outfo, precon, Pr); 625 630 if (precon) … … 630 635 631 636 if (istop == 3  istop == 7 ) // conlim or itnlim 632 printf("\n LSQR stopped early: istop = // 3d", istop);637 printf("\n LSQR stopped early: istop = //%d", istop); 633 638 634 639 … … 640 645 // grad = pdxxxmat( Aname,2,m,n,dy ); // grad = A"dy 641 646 grad.clear(); 642 m odel>matVecMult(2, grad, dy);647 matVecMult(2, grad, dy); 643 648 for (int k=0; k<nfix; k++) 644 649 grad[fix[k]] = 0; // grad is a work vector … … 666 671 // Find the maximum step. 667 672 // 668 realstepx1 = pdxxxstep(nlow, low, x1, dx1 );669 realstepx2 = inf;673 double stepx1 = pdxxxstep(nlow, low, x1, dx1 ); 674 double stepx2 = inf; 670 675 if (nupp > 0) 671 676 stepx2 = pdxxxstep(nupp, upp, x2, dx2 ); 672 realstepz1 = pdxxxstep( z1 , dz1 );673 realstepz2 = inf;677 double stepz1 = pdxxxstep( z1 , dz1 ); 678 double stepz2 = inf; 674 679 if (nupp > 0) 675 680 stepz2 = pdxxxstep( z2 , dz2 ); 676 realstepx = min( stepx1, stepx2 );677 realstepz = min( stepz1, stepz2 );678 stepx = min( steptol*stepx, 1 );679 stepz = min( steptol*stepz, 1 );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 ); 680 685 if (stepSame){ // For NLPs, force same step 681 686 stepx = min( stepx, stepz ); // (true Newton method) … … 703 708 // [obj,grad,hess] = feval( Fname, (x*beta) ); 704 709 x.scale(beta); 705 obj = model>getObj(x);706 model>getGrad(x, grad);707 model>getHessian(x, H);710 obj = getObj(x); 711 getGrad(x, grad); 712 getHessian(x, H); 708 713 x.scale((1.0/beta)); 709 714 … … 713 718 714 719 // [r1,r2,rL,rU,Pinf,Dinf] = ... 715 pdxxxresid1( model, nlow, nupp, nfix, low, upp, fix,720 pdxxxresid1( this, nlow, nupp, nfix, low, upp, fix, 716 721 b,bl_elts,bu_elts,d1,d2,grad,rL,rU,x,x1,x2, 717 722 y,z1,z2, r1, r2, &Pinf, &Dinf ); 718 // realcenter, Cinf, Cinf0;723 //double center, Cinf, Cinf0; 719 724 // [cL,cU,center,Cinf,Cinf0] = ... 720 725 pdxxxresid2( mu, nlow, nupp, low,upp,cL,cU,x1,x2,z1,z2, 721 726 ¢er, &Cinf, &Cinf0); 722 realfmeritnew = pdxxxmerit(nlow, nupp, low,upp,r1,r2,rL,rU,cL,cU );723 realstep = min( stepx, stepz );727 double fmeritnew = pdxxxmerit(nlow, nupp, low,upp,r1,r2,rL,rU,cL,cU ); 728 double step = min( stepx, stepz ); 724 729 725 730 if (fmeritnew <= (1  eta*step)*fmerit){ … … 822 827 // Reduce mu, and reset certain residuals. 823 828 824 realstepmu = min( stepx , stepz );829 double stepmu = min( stepx , stepz ); 825 830 stepmu = min( stepmu, steptol ); 826 realmuold = mu;831 double muold = mu; 827 832 mu = mu  stepmu * mu; 828 833 if (center >= bigcenter) … … 884 889 printf(" unscaled\n"); 885 890 886 time = now()  time;891 time = CoinCpuTime()  time; 887 892 char str1[100],str2[100]; 888 893 sprintf(str1, "\nPDitns =%10g", PDitns ); 889 sprintf(str2, "itns =%10 g", CGitns );894 sprintf(str2, "itns =%10d", CGitns ); 890 895 // printf( [str1 " " solver str2] ); 891 896 printf(" time =%10.1f\n", time); … … 919 924 printf(" Less than %f %d\n",thresh[2],counts[0]); 920 925 921 922 return 0; 926 return 0; 923 927 } 924 928 // LSQR … … 927 931 { 928 932 } 933 934 void ClpPdco::matVecMult( int mode, double* x_elts, double* y_elts){ 935 pdcoStuff_>matVecMult(this,mode,x_elts,y_elts); 936 } 937 void 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 943 void 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 950 void 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 } 956 void ClpPdco::matPrecon(double delta, double* x_elts, double* y_elts){ 957 pdcoStuff_>matPrecon(this,delta,x_elts,y_elts); 958 } 959 void 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 965 void 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 972 void 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 } 978 void 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 988 double ClpPdco::getObj(CoinDenseVector &x){ 989 return pdcoStuff_>getObj(this, x); 990 } 991 992 void ClpPdco::getGrad(CoinDenseVector &x, CoinDenseVector &g){ 993 pdcoStuff_>getGrad(this, x,g); 994 } 995 996 void ClpPdco::getHessian(CoinDenseVector &x, CoinDenseVector &H){ 997 pdcoStuff_>getHessian(this, x,H); 998 } 
branches/pre/Makefile.Clp
r214 r215 7 7 # between then specify the exact level you want, e.g., O1 or O2 8 8 OptLevel := g 9 OptLevel := O19 #OptLevel := O1 10 10 11 11 LIBNAME := Clp … … 36 36 LIBSRC += ClpInterior.cpp 37 37 LIBSRC += ClpPdco.cpp 38 LIBSRC += ClpPdcoBase.cpp 39 LIBSRC += ClpLsqr.cpp 38 40 # and Presolve stuff 39 41 LIBSRC += ClpPresolve.cpp 
branches/pre/include/ClpHelperFunctions.hpp
r214 r215 38 38 // b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 ) 39 39 40 inline void pdxxxresid1( Model*model, const int nlow, const int nupp, const int nfix,40 inline void pdxxxresid1(ClpPdco *model, const int nlow, const int nupp, const int nfix, 41 41 int *low, int *upp, int *fix, 42 42 CoinDenseVector &b, double *bl, double *bu, double d1, double d2, 
branches/pre/include/ClpInterior.hpp
r214 r215 16 16 #include "ClpMatrixBase.hpp" 17 17 #include "ClpSolve.hpp" 18 class ClpDualRowPivot; 19 class ClpPrimalColumnPivot; 20 class ClpFactorization; 21 class CoinIndexedVector; 22 class ClpNonLinearCost; 23 class ClpInteriorProgress; 18 class ClpLsqr; 19 class ClpPdcoBase; 24 20 // ******** DATA to be moved into protected section of ClpInterior 25 21 typedef struct{ … … 143 139 int pdco(); 144 140 // ** Temporary version 145 int pdco( Lsqr *lsqr, Options &options, Info &info, Outfo &outfo);141 int pdco( ClpPdcoBase * stuff, Options &options, Info &info, Outfo &outfo); 146 142 //@} 147 143 … … 245 241 /// Sum of primal infeasibilities 246 242 double sumPrimalInfeasibilities_; 243 public: 244 double xsize_; 245 double zsize_; 246 protected: 247 247 /// Working copy of lower bounds (Owner of arrays below) 248 248 double * lower_; … … 259 259 /// Working copy of objective 260 260 double * cost_; 261 public: 262 /// Rhs 263 double * rhs_; 264 double * x_; 265 double * y_; 266 double * dj_; 267 protected: 268 /// Pointer to Lsqr object 269 ClpLsqr * lsqrObject_; 270 /// Pointer to stuff 271 ClpPdcoBase * pdcoStuff_; 261 272 /// Which algorithm being used 262 273 int algorithm_; 
branches/pre/include/ClpPdco.hpp
r214 r215 19 19 20 20 */ 21 22 21 class ClpPdco : public ClpInterior { 23 22 … … 35 34 int pdco(); 36 35 // ** Temporary version 37 int pdco( Lsqr *lsqr, Options &options, Info &info, Outfo &outfo);36 int pdco( ClpPdcoBase * stuff, Options &options, Info &info, Outfo &outfo); 38 37 39 38 //@} 40 39 41 /**@name Functions used in dual*/40 /**@name Functions used in pdco */ 42 41 //@{ 43 42 /// LSQR 44 43 void lsqr(); 44 45 void matVecMult( int, double *, double *); 46 47 void matVecMult( int, CoinDenseVector &, double *); 48 49 void matVecMult( int, CoinDenseVector &, CoinDenseVector &); 50 51 void matVecMult( int, CoinDenseVector *, CoinDenseVector *); 52 53 void getBoundTypes( int *, int *, int *, int**); 54 55 void getGrad(CoinDenseVector &x, CoinDenseVector &grad); 56 57 void getHessian(CoinDenseVector &x, CoinDenseVector &H); 58 59 double getObj(CoinDenseVector &x); 60 61 void matPrecon( double, double *, double *); 62 63 void matPrecon( double, CoinDenseVector &, double *); 64 65 void matPrecon( double, CoinDenseVector &, CoinDenseVector &); 66 67 void matPrecon( double, CoinDenseVector *, CoinDenseVector *); 45 68 //@} 69 46 70 }; 47 71 #endif
Note: See TracChangeset
for help on using the changeset viewer.