Changeset 215 for branches


Ignore:
Timestamp:
Oct 3, 2003 11:41:53 AM (16 years ago)
Author:
forrest
Message:

Second try at pdco

Location:
branches/pre
Files:
8 added
6 edited

Legend:

Unmodified
Added
Removed
  • branches/pre/ClpInterior.cpp

    r214 r215  
    1111#include "CoinHelperFunctions.hpp"
    1212#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"
    1517#include "ClpMessage.hpp"
    1618#include "ClpLinearObjective.hpp"
     
    2931  sumDualInfeasibilities_(0.0),
    3032  sumPrimalInfeasibilities_(0.0),
     33  xsize_(0.0),
     34  zsize_(0.0),
    3135  lower_(NULL),
    3236  rowLowerWork_(NULL),
     
    3539  rowUpperWork_(NULL),
    3640  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)
    3848{
    3949  solveType_=2; // say interior based life form
     
    4252// Subproblem constructor
    4353ClpInterior::ClpInterior ( const ClpModel * rhs,
    44                      int numberRows, const int * whichRow,
    45                      int numberColumns, const int * whichColumn,
    46                      bool dropNames, bool dropIntegers)
     54                           int numberRows, const int * whichRow,
     55                           int numberColumns, const int * whichColumn,
     56                           bool dropNames, bool dropIntegers)
    4757  : ClpModel(rhs, numberRows, whichRow,
    4858             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)
    6078{
    6179  solveType_=2; // say interior based life form
     
    86104  sumDualInfeasibilities_(0.0),
    87105  sumPrimalInfeasibilities_(0.0),
     106  xsize_(0.0),
     107  zsize_(0.0),
    88108  lower_(NULL),
    89109  rowLowerWork_(NULL),
     
    92112  rowUpperWork_(NULL),
    93113  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)
    95121{
    96122  gutsOfDelete();
     
    105131  sumDualInfeasibilities_(0.0),
    106132  sumPrimalInfeasibilities_(0.0),
     133  xsize_(0.0),
     134  zsize_(0.0),
    107135  lower_(NULL),
    108136  rowLowerWork_(NULL),
     
    111139  rowUpperWork_(NULL),
    112140  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)
    114148{
    115149  solveType_=2; // say interior based life form
     
    137171  //cost_ = ClpCopyOfArray(rhs.cost_,2*(numberColumns_+numberRows_));
    138172  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();
    139179  largestPrimalError_ = rhs.largestPrimalError_;
    140180  largestDualError_ = rhs.largestDualError_;
    141181  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
    142182  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
     183  xsize_ = rhs.xsize_;
     184  zsize_ = rhs.zsize_;
    143185  solveType_=rhs.solveType_;
    144186}
     
    157199  delete [] cost_;
    158200  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;
    159213}
    160214bool
     
    402456// ** Temporary version
    403457int 
    404 ClpInterior::pdco( Lsqr *lsqr, Options &options, Info &info, Outfo &outfo)
    405 {
    406   return ((ClpPdco *) this)->pdco(lsqr,options,info,outfo);
    407 }
     458ClpInterior::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) 2002, International Business Machines
     1// Copyright (C) 2003, International Business Machines
    22// Corporation and others.  All Rights Reserved.
    33
     
    1515#include <math.h>
    1616
     17#include "CoinDenseVector.hpp"
     18#include "ClpPdco.hpp"
     19#include "ClpPdcoBase.hpp"
    1720#include "CoinHelperFunctions.hpp"
    18 #include "ClpPdco.hpp"
    19 #include "CoinDenseVector.hpp"
    2021#include "ClpHelperFunctions.hpp"
    21 #include "CoinPackedMatrix.hpp"
     22#include "ClpLsqr.hpp"
     23#include "CoinTime.hpp"
    2224#include "ClpMessage.hpp"
    2325#include <cfloat>
     
    3537// ** Temporary version
    3638int 
    37 ClpPdco::pdco( Lsqr *lsqr, Options &options, Info &info, Outfo &outfo)
     39ClpPdco::pdco( ClpPdcoBase * stuff, Options &options, Info &info, Outfo &outfo)
    3840{
    3941//    D1, D2 are positive-definite diagonal matrices defined from d1, d2.
     
    173175
    174176//  global pdDDD1 pdDDD2 pdDDD3
    175   real inf = 1.0e30;
    176   real eps = 1.0e-15;
    177   real atolold, r3ratio, Pinf, Dinf, Cinf, Cinf0;
     177  double inf = 1.0e30;
     178  double eps = 1.0e-15;
     179  double atolold, r3ratio, Pinf, Dinf, Cinf, Cinf0;
    178180 
    179181  printf("\n   --------------------------------------------------------");
     
    183185  printf("\n   --------------------------------------------------------\n");
    184186
    185   int m = model->rowSize();
    186   int n = model->colSize();
     187  int m = numberRows_;
     188  int n = numberColumns_;
    187189  bool ifexplicit = true;
    188190
    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_);
    193195  //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();
    203209
    204210  printf("\nmax |b | = %8g     max |x0| = %8g", normb , normx0);
    205   printf(                "      xsize   = %8g", model->xsize);
     211  printf(                "      xsize   = %8g", xsize_);
    206212  printf("\nmax |y0| = %8g     max |z0| = %8g", normy0, normz0);
    207   printf(                "      zsize   = %8g", model->zsize);
     213  printf(                "      zsize   = %8g", zsize_);
    208214
    209215  //---------------------------------------------------------------------
     
    221227  //---------------------------------------------------------------------
    222228  /*
    223   if (model->d1->size()==1)
    224       d1->resize(n, d1->getElements()[0]);  // Allow scalar d1, d2
    225   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)
    226232      d2->resize(m, d2->getElements()[0]);  // to mean dk * unit vector
    227233 */
    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();
    230237
    231238  //---------------------------------------------------------------------
     
    233240  //---------------------------------------------------------------------
    234241  int  maxitn    = options.MaxIter;
    235   real featol    = options.FeaTol;
    236   real opttol    = options.OptTol;
    237   real steptol   = options.StepTol;
     242  double featol    = options.FeaTol;
     243  double opttol    = options.OptTol;
     244  double steptol   = options.StepTol;
    238245  int  stepSame  = 1;  /* options.StepSame;   // 1 means stepx == stepz */
    239   real x0min     = options.x0min;
    240   real z0min     = options.z0min;
    241   real mu0       = options.mu0;
     246  double x0min     = options.x0min;
     247  double z0min     = options.z0min;
     248  double mu0       = options.mu0;
    242249  int  LSproblem = options.LSproblem;  // See below
    243250  int  LSmethod  = options.LSmethod;   // 1=Cholesky    2=QR    3=LSQR
    244251  int  itnlim    = options.LSQRMaxIter * min(m,n);
    245   real atol1     = options.LSQRatol1;  // Initial  atol
    246   real atol2     = options.LSQRatol2;  // Smallest atol,unless atol1 is smaller
    247   real conlim    = options.LSQRconlim;
     252  double atol1     = options.LSQRatol1;  // Initial  atol
     253  double atol2     = options.LSQRatol2;  // Smallest atol,unless atol1 is smaller
     254  double conlim    = options.LSQRconlim;
    248255  int  wait      = options.wait;
    249256
     
    259266  //---------------------------------------------------------------------
    260267  int  kminor    = 0;      // 1 stops after each iteration
    261   real eta       = 1e-4;   // Linesearch tolerance for "sufficient descent"
    262   real maxf      = 10;     // Linesearch backtrack limit (function evaluations)
    263   real maxfail   = 1;      // Linesearch failure limit (consecutive iterations)
    264   real bigcenter = 1e+3;   // mu is reduced if center < bigcenter.
     268  double eta       = 1e-4;   // Linesearch tolerance for "sufficient descent"
     269  double maxf      = 10;     // Linesearch backtrack limit (function evaluations)
     270  double maxfail   = 1;      // Linesearch failure limit (consecutive iterations)
     271  double bigcenter = 1e+3;   // mu is reduced if center < bigcenter.
    265272
    266273  // Parameters for LSQR.
    267   real atolmin   = eps;    // Smallest atol if linesearch back-tracks
    268   real btol      = 0;      // Should be small (zero is ok)
    269   real show      = false;  // Controls lsqr iteration log
     274  double atolmin   = eps;    // Smallest atol if linesearch back-tracks
     275  double btol      = 0;      // Should be small (zero is ok)
     276  double show      = false;  // Controls lsqr iteration log
    270277  /*
    271   real gamma     = d1->infNorm();
    272   real delta     = d2->infNorm();
     278  double gamma     = d1->infNorm();
     279  double delta     = d2->infNorm();
    273280  */
    274   real gamma = d1;
    275   real delta = d2;
     281  double gamma = d1;
     282  double delta = d2;
    276283
    277284  printf("\n\nx0min    = %8g     featol   = %8.1e", x0min, featol);
     
    302309  // All parameters have now been set.
    303310  //---------------------------------------------------------------------
    304   //real time    = cputime();
    305   double time = now();
     311  double time    = CoinCpuTime();
    306312  bool useChol = (LSmethod == 1);
    307313  bool useQR   = (LSmethod == 2);
     
    317323  int nlow, nupp, nfix;
    318324  int *bptrs[3] = {0};
    319   model->getBoundTypes(&nlow, &nupp, &nfix, bptrs );
     325  getBoundTypes(&nlow, &nupp, &nfix, bptrs );
    320326  int *low = bptrs[0];
    321327  int *upp = bptrs[1];
     
    329335  //---------------------------------------------------------------------
    330336
    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();
    337341
    338342  CoinDenseVector r1(m, 0.0);
    339   real *r1_elts = r1.getElements();
     343  double *r1_elts = r1.getElements();
    340344  CoinDenseVector x1(n, 0.0);
    341   real *x1_elts = x1.getElements();
     345  double *x1_elts = x1.getElements();
    342346
    343347  if (nfix > 0){
    344348    for (int k=0; k<nfix; k++)
    345349      x1_elts[fix[k]] = bl[fix[k]];
    346     model->matVecMult(1, r1, x1);
     350    matVecMult(1, r1, x1);
    347351    b = b - r1;
    348352      // At some stage, might want to look at normfix = norm(r1,inf);
     
    362366  //    Hessian  = (beta2/theta) hess.
    363367  //---------------------------------------------------------------------
    364   real beta = model->xsize;   if (beta==0) beta = 1;   // beta scales b, x.
    365   real zeta = model->zsize;   if (zeta==0) zeta = 1;   // zeta scales y, z.
    366   real theta  = beta*zeta;                            // theta scales obj.
     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.
    367371  // (theta could be anything, but theta = beta*zeta makes
    368372  // scaled grad = grad/zeta = 1 approximately if zeta is chosen right.)
     
    375379  d2     = d2 * ( sqrt(theta)/beta );
    376380
    377   real beta2  = beta*beta;
     381  double beta2  = beta*beta;
    378382  b.scale( (1.0/beta) );
    379383  y.scale( (1.0/zeta) );
     
    407411  CoinDenseVector Pr(m);
    408412  CoinDenseVector D(n);
    409   real *D_elts = D.getElements();
     413  double *D_elts = D.getElements();
    410414  CoinDenseVector w(n);
    411   real *w_elts = w.getElements();
     415  double *w_elts = w.getElements();
    412416  CoinDenseVector rhs(m+n);
    413417   
     
    416420  // Pull out the element array pointers for efficiency
    417421  //---------------------------------------------------------------------
    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();
    423427
    424428  for (int k=0; k<nlow; k++){
     
    436440//  [obj,grad,hess] = feval( Fname, (x*beta) );
    437441  x.scale(beta);
    438   real obj = model->getObj(x);
     442  double obj = getObj(x);
    439443  CoinDenseVector grad(n);
    440   model->getGrad(x, grad);
     444  getGrad(x, grad);
    441445  CoinDenseVector H(n);
    442   model->getHessian(x ,H);
     446  getHessian(x ,H);
    443447  x.scale((1.0/beta));
    444448 
    445   real * g_elts=grad.getElements();
    446   real * H_elts=H.getElements();
     449  double * g_elts=grad.getElements();
     450  double * H_elts=H.getElements();
    447451
    448452  obj /= theta;                       // Scaled obj.
     
    461465  //      pdxxxresid1( Aname,fix,low,upp, ...
    462466  //                   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,
    464468               b,bl_elts,bu_elts,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2,
    465469               r1, r2, &Pinf, &Dinf);
     
    477481  //              regarding mu0 as a scaling of the initial center.
    478482  //---------------------------------------------------------------------
    479   //  real mufirst = mu0*(x0min * z0min);
    480   real mufirst = mu0;   // revert to absolute value
    481   real mulast  = 0.1 * opttol;
     483  //  double mufirst = mu0*(x0min * z0min);
     484  double mufirst = mu0;   // revert to absolute value
     485  double mulast  = 0.1 * opttol;
    482486  mulast  = min( mulast, mufirst );
    483   real mu      = mufirst;
    484   real center,  fmerit;
     487  double mu      = mufirst;
     488  double center,  fmerit;
    485489  pdxxxresid2( mu, nlow, nupp, low,upp, cL, cU, x1, x2,
    486490                        z1, z2, &center, &Cinf, &Cinf0 );
     
    490494
    491495  bool  precon   = true;       
    492   real PDitns    = 0;
     496  double PDitns    = 0;
    493497  bool converged = false;
    494   real atol      = atol1;
     498  double atol      = atol1;
    495499  atol2     = max( atol2, atolmin );
    496500  atolmin   = atol2;
     
    499503  //  Iteration log.
    500504
    501   real stepx   = 0;
    502   real stepz   = 0;
     505  double stepx   = 0;
     506  double stepz   = 0;
    503507  int nf      = 0;
    504508  int itncg   = 0;
     
    513517  }
    514518
    515   real regx = (d1*x).twoNorm();
    516   real regy = (d2*y).twoNorm();
     519  double regx = (d1*x).twoNorm();
     520  double regy = (d2*y).twoNorm();
    517521  //  regterm = twoNorm(d1.*x)^2  +  norm(d2.*y)^2;
    518   real regterm = regx*regx + regy*regy;
    519   real objreg  = obj  +  0.5*regterm;
    520   real objtrue = objreg * theta;
     522  double regterm = regx*regx + regy*regy;
     523  double objreg  = obj  +  0.5*regterm;
     524  double objtrue = objreg * theta;
    521525
    522526  printf("\n%3g                 ", PDitns        );
     
    533537  // Main loop.
    534538  //---------------------------------------------------------------------
     539  // Lsqr
     540  ClpLsqr  thisLsqr(this);
    535541
    536542  //  while (converged) {
     
    543549    //              Now that starting conditions are better, go back to 0.1.
    544550
    545     real r3norm = max(Pinf,   max(Dinf,  Cinf));
     551    double r3norm = max(Pinf,   max(Dinf,  Cinf));
    546552    atol   = min(atol,  r3norm*0.1);
    547553    atol   = max(atol,  atolmin   );
     
    581587      for (int k=0; k<n; k++)
    582588        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;
    586591
    587592      if (direct){
     
    593598        for (int k=0; k<m; k++)
    594599          rhs[n+k] = r1_elts[k]*(1.0/d2);
    595         real damp    = 0;
     600        double damp    = 0;
    596601       
    597602        if (precon){    // Construct diagonal preconditioner for LSQR
    598           model->matPrecon(d2, Pr, D);
     603          matPrecon(d2, Pr, D);
    599604        } 
    600605            /*
     
    608613       
    609614
    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);
    614619        */
    615620        //  New version of lsqr
     
    621626        info.r3norm  = fmerit;  // Must be the 2-norm here.
    622627
    623         lsqr->do_lsqr(model, rhs, damp, atol, btol, conlim, itnlim,
     628        thisLsqr.do_lsqr( rhs, damp, atol, btol, conlim, itnlim,
    624629                      show, info, dy , &istop, &itncg, &outfo, precon, Pr);
    625630        if (precon)
     
    630635
    631636        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);
    633638       
    634639
     
    640645      //      grad      = pdxxxmat( Aname,2,m,n,dy );   // grad = A"dy
    641646      grad.clear();
    642       model->matVecMult(2, grad, dy);
     647      matVecMult(2, grad, dy);
    643648      for (int k=0; k<nfix; k++)
    644649        grad[fix[k]] = 0;                            // grad is a work vector
     
    666671    // Find the maximum step.
    667672    //--------------------------------------------------------------------
    668     real stepx1 = pdxxxstep(nlow, low, x1, dx1 );
    669     real stepx2 = inf;
     673    double stepx1 = pdxxxstep(nlow, low, x1, dx1 );
     674    double stepx2 = inf;
    670675    if (nupp > 0)
    671676      stepx2 = pdxxxstep(nupp, upp, x2, dx2 );
    672     real stepz1 = pdxxxstep( z1     , dz1      );
    673     real stepz2 = inf;
     677    double stepz1 = pdxxxstep( z1     , dz1      );
     678    double stepz2 = inf;
    674679    if (nupp > 0)
    675680      stepz2 = pdxxxstep( z2     , dz2      );
    676     real stepx  = min( stepx1, stepx2 );
    677     real stepz  = min( stepz1, stepz2 );
    678     stepx  = min( steptol*stepx, 1 );
    679     stepz  = min( steptol*stepz, 1 );
     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 );
    680685    if (stepSame){                   // For NLPs, force same step
    681686      stepx = min( stepx, stepz );   // (true Newton method)
     
    703708      //      [obj,grad,hess] = feval( Fname, (x*beta) );
    704709      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);
    708713      x.scale((1.0/beta));
    709714     
     
    713718     
    714719      //      [r1,r2,rL,rU,Pinf,Dinf] = ...
    715       pdxxxresid1( model, nlow, nupp, nfix, low, upp, fix,
     720      pdxxxresid1( this, nlow, nupp, nfix, low, upp, fix,
    716721                   b,bl_elts,bu_elts,d1,d2,grad,rL,rU,x,x1,x2,
    717722                   y,z1,z2, r1, r2, &Pinf, &Dinf );
    718       //real center, Cinf, Cinf0;
     723      //double center, Cinf, Cinf0;
    719724      //      [cL,cU,center,Cinf,Cinf0] = ...
    720725                  pdxxxresid2( mu, nlow, nupp, low,upp,cL,cU,x1,x2,z1,z2,
    721726                               &center, &Cinf, &Cinf0);
    722       real fmeritnew = pdxxxmerit(nlow, nupp, low,upp,r1,r2,rL,rU,cL,cU );
    723       real step      = min( stepx, stepz );
     727      double fmeritnew = pdxxxmerit(nlow, nupp, low,upp,r1,r2,rL,rU,cL,cU );
     728      double step      = min( stepx, stepz );
    724729
    725730      if (fmeritnew <= (1 - eta*step)*fmerit){
     
    822827      // Reduce mu, and reset certain residuals.
    823828
    824       real stepmu  = min( stepx , stepz   );
     829      double stepmu  = min( stepx , stepz   );
    825830      stepmu  = min( stepmu, steptol );
    826       real muold   = mu;
     831      double muold   = mu;
    827832      mu      = mu   -  stepmu * mu;
    828833      if (center >= bigcenter)
     
    884889  printf(" unscaled\n");
    885890
    886   time   = now() - time;
     891  time   = CoinCpuTime() - time;
    887892  char str1[100],str2[100];
    888893  sprintf(str1, "\nPDitns  =%10g", PDitns );
    889   sprintf(str2, "itns =%10g", CGitns );
     894  sprintf(str2, "itns =%10d", CGitns );
    890895  //  printf( [str1 " " solver str2] );
    891896  printf("    time    =%10.1f\n", time);
     
    919924  printf("   Less than   %f %d\n",thresh[2],counts[0]);
    920925
    921 
    922 return 0;
     926  return 0;
    923927}
    924928// LSQR
     
    927931{
    928932}
     933
     934void ClpPdco::matVecMult( int mode, double* x_elts, double* y_elts){
     935  pdcoStuff_->matVecMult(this,mode,x_elts,y_elts);
     936}
     937void ClpPdco::matVecMult( int mode, CoinDenseVector &x, double *y_elts){
     938  double *x_elts = x.getElements();
     939  matVecMult( mode, x_elts, y_elts);
     940  return;
     941}
     942
     943void ClpPdco::matVecMult( int mode, CoinDenseVector &x, CoinDenseVector &y){
     944  double *x_elts = x.getElements();
     945  double *y_elts = y.getElements();
     946  matVecMult( mode, x_elts, y_elts);
     947  return;
     948}
     949
     950void ClpPdco::matVecMult( int mode, CoinDenseVector *x, CoinDenseVector *y){
     951  double *x_elts = x->getElements();
     952  double *y_elts = y->getElements();
     953  matVecMult( mode, x_elts, y_elts);
     954  return;
     955}
     956void ClpPdco::matPrecon(double delta, double* x_elts, double* y_elts){
     957  pdcoStuff_->matPrecon(this,delta,x_elts,y_elts);
     958}
     959void ClpPdco::matPrecon(double delta, CoinDenseVector &x, double *y_elts){
     960  double *x_elts = x.getElements();
     961  matPrecon(delta, x_elts, y_elts);
     962  return;
     963}
     964   
     965void ClpPdco::matPrecon(double delta, CoinDenseVector &x, CoinDenseVector &y){
     966  double *x_elts = x.getElements();
     967  double *y_elts = y.getElements();
     968  matPrecon(delta, x_elts, y_elts);
     969  return;
     970}
     971
     972void ClpPdco::matPrecon(double delta, CoinDenseVector *x, CoinDenseVector *y){
     973  double *x_elts = x->getElements();
     974  double *y_elts = y->getElements();
     975  matPrecon(delta, x_elts, y_elts);
     976  return;
     977}
     978void ClpPdco::getBoundTypes(int *nlow, int *nupp, int *nfix, int **bptrs){
     979  *nlow = numberColumns_;
     980  *nupp = *nfix = 0;
     981  int *low_ = (int *)malloc(numberColumns_*sizeof(int)) ;
     982  for (int k=0; k <numberColumns_; k++)
     983    low_[k] = k;
     984  bptrs[0]= low_;
     985  return;
     986}
     987
     988double ClpPdco::getObj(CoinDenseVector &x){
     989  return pdcoStuff_->getObj(this, x);
     990}
     991
     992void ClpPdco::getGrad(CoinDenseVector &x, CoinDenseVector &g){
     993  pdcoStuff_->getGrad(this, x,g);
     994}
     995
     996void ClpPdco::getHessian(CoinDenseVector &x, CoinDenseVector &H){
     997  pdcoStuff_->getHessian(this, x,H);
     998}
  • branches/pre/Makefile.Clp

    r214 r215  
    77# between then specify the exact level you want, e.g., -O1 or -O2
    88OptLevel := -g
    9 OptLevel := -O1
     9#OptLevel := -O1
    1010
    1111LIBNAME := Clp
     
    3636LIBSRC += ClpInterior.cpp
    3737LIBSRC += ClpPdco.cpp
     38LIBSRC += ClpPdcoBase.cpp
     39LIBSRC += ClpLsqr.cpp
    3840# and Presolve stuff
    3941LIBSRC += ClpPresolve.cpp
  • branches/pre/include/ClpHelperFunctions.hpp

    r214 r215  
    3838//                   b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 )
    3939
    40 inline void pdxxxresid1(Model *model, const int nlow, const int nupp, const int nfix,
     40inline void pdxxxresid1(ClpPdco *model, const int nlow, const int nupp, const int nfix,
    4141                 int *low, int *upp, int *fix,
    4242                 CoinDenseVector &b, double *bl, double *bu, double d1, double d2,
  • branches/pre/include/ClpInterior.hpp

    r214 r215  
    1616#include "ClpMatrixBase.hpp"
    1717#include "ClpSolve.hpp"
    18 class ClpDualRowPivot;
    19 class ClpPrimalColumnPivot;
    20 class ClpFactorization;
    21 class CoinIndexedVector;
    22 class ClpNonLinearCost;
    23 class ClpInteriorProgress;
     18class ClpLsqr;
     19class ClpPdcoBase;
    2420// ******** DATA to be moved into protected section of ClpInterior
    2521typedef struct{
     
    143139  int pdco();
    144140  // ** Temporary version
    145   int  pdco( Lsqr *lsqr, Options &options, Info &info, Outfo &outfo);
     141  int  pdco( ClpPdcoBase * stuff, Options &options, Info &info, Outfo &outfo);
    146142  //@}
    147143
     
    245241  /// Sum of primal infeasibilities
    246242  double sumPrimalInfeasibilities_;
     243public:
     244  double xsize_;
     245  double zsize_;
     246protected:
    247247  /// Working copy of lower bounds (Owner of arrays below)
    248248  double * lower_;
     
    259259  /// Working copy of objective
    260260  double * cost_;
     261public:
     262  /// Rhs
     263  double * rhs_;
     264  double * x_;
     265  double * y_;
     266  double * dj_;
     267protected:
     268  /// Pointer to Lsqr object
     269  ClpLsqr * lsqrObject_;
     270  /// Pointer to stuff
     271  ClpPdcoBase * pdcoStuff_;
    261272  /// Which algorithm being used
    262273  int algorithm_;
  • branches/pre/include/ClpPdco.hpp

    r214 r215  
    1919
    2020*/
    21 
    2221class ClpPdco : public ClpInterior {
    2322
     
    3534  int pdco();
    3635  // ** Temporary version
    37   int  pdco( Lsqr *lsqr, Options &options, Info &info, Outfo &outfo);
     36  int  pdco( ClpPdcoBase * stuff, Options &options, Info &info, Outfo &outfo);
    3837
    3938  //@}
    4039
    41   /**@name Functions used in dual */
     40  /**@name Functions used in pdco */
    4241  //@{
    4342  /// LSQR
    4443  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 *);
    4568  //@}
     69
    4670};
    4771#endif
Note: See TracChangeset for help on using the changeset viewer.