source: branches/pre/include/ClpHelperFunctions.hpp @ 215

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

Second try at pdco

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.0 KB
Line 
1
2inline double pdxxxmerit(int nlow, int nupp, int *low, int *upp, CoinDenseVector &r1,
3                CoinDenseVector &r2, CoinDenseVector &rL,
4                CoinDenseVector &rU, CoinDenseVector &cL,
5                CoinDenseVector &cU ){
6
7// Evaluate the merit function for Newton's method.
8// It is the 2-norm of the three sets of residuals.
9  double sum1, sum2;
10  CoinDenseVector f(6);
11  f[0] = r1.twoNorm();
12  f[1] = r2.twoNorm();
13  sum1 = sum2 = 0.0;
14  for (int k=0; k<nlow; k++){
15    sum1 += rL[low[k]]*rL[low[k]];
16    sum2 += cL[low[k]]*cL[low[k]];
17  }
18  f[2] = sqrt(sum1);
19  f[4] = sqrt(sum2);
20  sum1 = sum2 = 0.0;
21  for (int k=0; k<nupp; k++){
22    sum1 += rL[upp[k]]*rL[upp[k]];
23    sum2 += cL[upp[k]]*cL[upp[k]];
24  }
25  f[3] = sqrt(sum1);
26  f[5] = sqrt(sum2);
27
28  return f.twoNorm();
29}
30
31//-----------------------------------------------------------------------
32// End private function pdxxxmerit
33//-----------------------------------------------------------------------
34
35
36//function [r1,r2,rL,rU,Pinf,Dinf] =    ...
37//      pdxxxresid1( Aname,fix,low,upp, ...
38//                   b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 )
39
40inline void pdxxxresid1(ClpPdco *model, const int nlow, const int nupp, const int nfix,
41                 int *low, int *upp, int *fix,
42                 CoinDenseVector &b, double *bl, double *bu, double d1, double d2,
43                 CoinDenseVector &grad, CoinDenseVector &rL,
44                 CoinDenseVector &rU, CoinDenseVector &x,
45                 CoinDenseVector &x1, CoinDenseVector &x2,
46                 CoinDenseVector &y,  CoinDenseVector &z1,
47                 CoinDenseVector &z2, CoinDenseVector &r1,
48                 CoinDenseVector &r2, double *Pinf, double *Dinf){
49
50// Form residuals for the primal and dual equations.
51// rL, rU are output, but we input them as full vectors
52// initialized (permanently) with any relevant zeros.
53
54// Get some element pointers for efficiency
55  double *x_elts  = x.getElements();
56  double *x1_elts = x2.getElements();
57  double *x2_elts = x2.getElements();
58  double *z1_elts = z1.getElements();
59  double *z2_elts = z2.getElements();
60  double *r2_elts = r2.getElements();
61 
62  for (int k=0; k<nfix; k++)
63    x_elts[fix[k]]  = 0;
64
65  r1.clear();
66  r2.clear();   
67  model->matVecMult( 1, r1, x );
68  model->matVecMult( 2, r2, y );
69  for (int k=0; k<nfix; k++)
70    r2_elts[fix[k]]  = 0;
71
72
73  r1      = b    - r1 - d2*d2*y;
74  r2      = grad - r2 - z1;              // grad includes d1*d1*x
75  if(nupp > 0)       
76    r2    = r2 + z2;   
77
78  for (int k=0; k<nlow; k++)
79    rL[low[k]] = bl[low[k]] - x[low[k]] + x1[low[k]];
80  for (int k=0; k<nupp; k++)
81    rU[upp[k]] = - bu[upp[k]] + x[upp[k]] + x2[upp[k]];
82
83  double normL = 0.0;
84  double normU = 0.0;
85  for (int k=0; k<nlow; k++)
86    if (rL[low[k]] > normL) normL = rL[low[k]];
87  for (int k=0; k<nupp; k++)
88    if (rU[upp[k]] > normU) normU = rU[upp[k]];
89
90  *Pinf    = max(normL, normU); 
91  *Pinf    = max( r1.infNorm() , *Pinf );
92  *Dinf    = r2.infNorm();
93  *Pinf    = max( *Pinf, 1e-99 );
94  *Dinf    = max( *Dinf, 1e-99 );
95}
96
97//-----------------------------------------------------------------------
98// End private function pdxxxresid1
99//-----------------------------------------------------------------------
100
101
102//function [cL,cU,center,Cinf,Cinf0] = ...
103//      pdxxxresid2( mu,low,upp,cL,cU,x1,x2,z1,z2 )
104
105inline void pdxxxresid2(double mu, int nlow, int nupp, int *low, int *upp,
106                 CoinDenseVector &cL, CoinDenseVector &cU,
107                 CoinDenseVector &x1, CoinDenseVector &x2,
108                 CoinDenseVector &z1, CoinDenseVector &z2,
109                 double *center, double *Cinf, double *Cinf0){
110
111// Form residuals for the complementarity equations.
112// cL, cU are output, but we input them as full vectors
113// initialized (permanently) with any relevant zeros.
114// Cinf  is the complementarity residual for X1 z1 = mu e, etc.
115// Cinf0 is the same for mu=0 (i.e., for the original problem).
116
117  double maxXz = -1e20;
118  double minXz = 1e20;
119
120  double *x1_elts = x1.getElements();
121  double *z1_elts = z1.getElements();
122  double *cL_elts = cL.getElements();
123  for (int k=0; k<nlow; k++){
124    double x1z1    = x1_elts[low[k]] * z1_elts[low[k]];
125    cL_elts[low[k]] = mu - x1z1;
126    if(x1z1 > maxXz) maxXz = x1z1;
127    if(x1z1 < minXz) minXz = x1z1;
128  }
129
130  double *x2_elts = x2.getElements();
131  double *z2_elts = z2.getElements();
132  double *cU_elts = cU.getElements();
133  for (int k=0; k<nupp; k++){
134    double x2z2    = x2_elts[upp[k]] * z2_elts[upp[k]];
135    cU_elts[upp[k]] = mu - x2z2;
136    if(x2z2 > maxXz) maxXz = x2z2;
137    if(x2z2 < minXz) minXz = x2z2;
138  }
139
140  maxXz   = max( maxXz, 1e-99 );
141  minXz   = max( minXz, 1e-99 );
142  *center  = maxXz / minXz;
143
144  double normL = 0.0;
145  double normU = 0.0;
146  for (int k=0; k<nlow; k++)
147    if (cL_elts[low[k]] > normL) normL = cL_elts[low[k]];
148  for (int k=0; k<nupp; k++)
149    if (cU_elts[upp[k]] > normU) normU = cU_elts[upp[k]];
150  *Cinf    = max( normL, normU);
151  *Cinf0   = maxXz;
152}
153//-----------------------------------------------------------------------
154// End private function pdxxxresid2
155//-----------------------------------------------------------------------
156
157inline double  pdxxxstep( CoinDenseVector &x, CoinDenseVector &dx ){
158
159// Assumes x > 0.
160// Finds the maximum step such that x + step*dx >= 0.
161
162  double step     = 1e+20;
163
164  int n = x.size();
165  double *x_elts = x.getElements();
166  double *dx_elts = dx.getElements();
167  for (int k=0; k<n; k++)
168    if (dx_elts[k] < 0)
169      if ((x_elts[k]/(-dx_elts[k])) < step)
170        step = x_elts[k]/(-dx_elts[k]);
171  return step;
172}
173//-----------------------------------------------------------------------
174// End private function pdxxxstep
175//-----------------------------------------------------------------------
176
177inline double  pdxxxstep(int nset, int *set, CoinDenseVector &x, CoinDenseVector &dx ){
178
179// Assumes x > 0.
180// Finds the maximum step such that x + step*dx >= 0.
181
182  double step     = 1e+20;
183
184  int n = x.size();
185  double *x_elts = x.getElements();
186  double *dx_elts = dx.getElements();
187  for (int k=0; k<n; k++)
188    if (dx_elts[k] < 0)
189      if ((x_elts[k]/(-dx_elts[k])) < step)
190        step = x_elts[k]/(-dx_elts[k]);
191  return step;
192}
193//-----------------------------------------------------------------------
194// End private function pdxxxstep
195//-----------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.