ClpHelperFunctions.hpp
Go to the documentation of this file.
1 /* $Id: ClpHelperFunctions.hpp 1525 2010-02-26 17:27:59Z mjs $ */
2 // Copyright (C) 2003, International Business Machines
3 // Corporation and others. All Rights Reserved.
4 #ifndef ClpHelperFunctions_H
5 #define ClpHelperFunctions_H
6 
15 double maximumAbsElement(const double * region, int size);
16 void setElements(double * region, int size, double value);
17 void multiplyAdd(const double * region1, int size, double multiplier1,
18  double * region2, double multiplier2);
19 double innerProduct(const double * region1, int size, const double * region2);
20 void getNorms(const double * region, int size, double & norm1, double & norm2);
21 #if COIN_LONG_WORK
22 // For long double versions
23 CoinWorkDouble maximumAbsElement(const CoinWorkDouble * region, int size);
24 void setElements(CoinWorkDouble * region, int size, CoinWorkDouble value);
25 void multiplyAdd(const CoinWorkDouble * region1, int size, CoinWorkDouble multiplier1,
26  CoinWorkDouble * region2, CoinWorkDouble multiplier2);
27 CoinWorkDouble innerProduct(const CoinWorkDouble * region1, int size, const CoinWorkDouble * region2);
28 void getNorms(const CoinWorkDouble * region, int size, CoinWorkDouble & norm1, CoinWorkDouble & norm2);
29 inline void
30 CoinMemcpyN(const double * from, const int size, CoinWorkDouble * to)
31 {
32  for (int i = 0; i < size; i++)
33  to[i] = from[i];
34 }
35 inline void
36 CoinMemcpyN(const CoinWorkDouble * from, const int size, double * to)
37 {
38  for (int i = 0; i < size; i++)
39  to[i] = static_cast<double>(from[i]);
40 }
41 inline CoinWorkDouble
42 CoinMax(const CoinWorkDouble x1, const double x2)
43 {
44  return (x1 > x2) ? x1 : x2;
45 }
46 inline CoinWorkDouble
47 CoinMax(double x1, const CoinWorkDouble x2)
48 {
49  return (x1 > x2) ? x1 : x2;
50 }
51 inline CoinWorkDouble
52 CoinMin(const CoinWorkDouble x1, const double x2)
53 {
54  return (x1 < x2) ? x1 : x2;
55 }
56 inline CoinWorkDouble
57 CoinMin(double x1, const CoinWorkDouble x2)
58 {
59  return (x1 < x2) ? x1 : x2;
60 }
61 inline CoinWorkDouble CoinSqrt(CoinWorkDouble x)
62 {
63  return sqrtl(x);
64 }
65 #else
66 inline double CoinSqrt(double x)
67 {
68  return sqrt(x);
69 }
70 #endif
71 
73 #ifdef ClpPdco_H
74 
75 
76 inline double pdxxxmerit(int nlow, int nupp, int *low, int *upp, CoinDenseVector <double> &r1,
77  CoinDenseVector <double> &r2, CoinDenseVector <double> &rL,
78  CoinDenseVector <double> &rU, CoinDenseVector <double> &cL,
79  CoinDenseVector <double> &cU )
80 {
81 
82 // Evaluate the merit function for Newton's method.
83 // It is the 2-norm of the three sets of residuals.
84  double sum1, sum2;
85  CoinDenseVector <double> f(6);
86  f[0] = r1.twoNorm();
87  f[1] = r2.twoNorm();
88  sum1 = sum2 = 0.0;
89  for (int k = 0; k < nlow; k++) {
90  sum1 += rL[low[k]] * rL[low[k]];
91  sum2 += cL[low[k]] * cL[low[k]];
92  }
93  f[2] = sqrt(sum1);
94  f[4] = sqrt(sum2);
95  sum1 = sum2 = 0.0;
96  for (int k = 0; k < nupp; k++) {
97  sum1 += rL[upp[k]] * rL[upp[k]];
98  sum2 += cL[upp[k]] * cL[upp[k]];
99  }
100  f[3] = sqrt(sum1);
101  f[5] = sqrt(sum2);
102 
103  return f.twoNorm();
104 }
105 
106 //-----------------------------------------------------------------------
107 // End private function pdxxxmerit
108 //-----------------------------------------------------------------------
109 
110 
111 //function [r1,r2,rL,rU,Pinf,Dinf] = ...
112 // pdxxxresid1( Aname,fix,low,upp, ...
113 // b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 )
114 
115 inline void pdxxxresid1(ClpPdco *model, const int nlow, const int nupp, const int nfix,
116  int *low, int *upp, int *fix,
117  CoinDenseVector <double> &b, double *bl, double *bu, double d1, double d2,
118  CoinDenseVector <double> &grad, CoinDenseVector <double> &rL,
119  CoinDenseVector <double> &rU, CoinDenseVector <double> &x,
120  CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
121  CoinDenseVector <double> &y, CoinDenseVector <double> &z1,
122  CoinDenseVector <double> &z2, CoinDenseVector <double> &r1,
123  CoinDenseVector <double> &r2, double *Pinf, double *Dinf)
124 {
125 
126 // Form residuals for the primal and dual equations.
127 // rL, rU are output, but we input them as full vectors
128 // initialized (permanently) with any relevant zeros.
129 
130 // Get some element pointers for efficiency
131  double *x_elts = x.getElements();
132  double *r2_elts = r2.getElements();
133 
134  for (int k = 0; k < nfix; k++)
135  x_elts[fix[k]] = 0;
136 
137  r1.clear();
138  r2.clear();
139  model->matVecMult( 1, r1, x );
140  model->matVecMult( 2, r2, y );
141  for (int k = 0; k < nfix; k++)
142  r2_elts[fix[k]] = 0;
143 
144 
145  r1 = b - r1 - d2 * d2 * y;
146  r2 = grad - r2 - z1; // grad includes d1*d1*x
147  if (nupp > 0)
148  r2 = r2 + z2;
149 
150  for (int k = 0; k < nlow; k++)
151  rL[low[k]] = bl[low[k]] - x[low[k]] + x1[low[k]];
152  for (int k = 0; k < nupp; k++)
153  rU[upp[k]] = - bu[upp[k]] + x[upp[k]] + x2[upp[k]];
154 
155  double normL = 0.0;
156  double normU = 0.0;
157  for (int k = 0; k < nlow; k++)
158  if (rL[low[k]] > normL) normL = rL[low[k]];
159  for (int k = 0; k < nupp; k++)
160  if (rU[upp[k]] > normU) normU = rU[upp[k]];
161 
162  *Pinf = CoinMax(normL, normU);
163  *Pinf = CoinMax( r1.infNorm() , *Pinf );
164  *Dinf = r2.infNorm();
165  *Pinf = CoinMax( *Pinf, 1e-99 );
166  *Dinf = CoinMax( *Dinf, 1e-99 );
167 }
168 
169 //-----------------------------------------------------------------------
170 // End private function pdxxxresid1
171 //-----------------------------------------------------------------------
172 
173 
174 //function [cL,cU,center,Cinf,Cinf0] = ...
175 // pdxxxresid2( mu,low,upp,cL,cU,x1,x2,z1,z2 )
176 
177 inline void pdxxxresid2(double mu, int nlow, int nupp, int *low, int *upp,
178  CoinDenseVector <double> &cL, CoinDenseVector <double> &cU,
179  CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
180  CoinDenseVector <double> &z1, CoinDenseVector <double> &z2,
181  double *center, double *Cinf, double *Cinf0)
182 {
183 
184 // Form residuals for the complementarity equations.
185 // cL, cU are output, but we input them as full vectors
186 // initialized (permanently) with any relevant zeros.
187 // Cinf is the complementarity residual for X1 z1 = mu e, etc.
188 // Cinf0 is the same for mu=0 (i.e., for the original problem).
189 
190  double maxXz = -1e20;
191  double minXz = 1e20;
192 
193  double *x1_elts = x1.getElements();
194  double *z1_elts = z1.getElements();
195  double *cL_elts = cL.getElements();
196  for (int k = 0; k < nlow; k++) {
197  double x1z1 = x1_elts[low[k]] * z1_elts[low[k]];
198  cL_elts[low[k]] = mu - x1z1;
199  if (x1z1 > maxXz) maxXz = x1z1;
200  if (x1z1 < minXz) minXz = x1z1;
201  }
202 
203  double *x2_elts = x2.getElements();
204  double *z2_elts = z2.getElements();
205  double *cU_elts = cU.getElements();
206  for (int k = 0; k < nupp; k++) {
207  double x2z2 = x2_elts[upp[k]] * z2_elts[upp[k]];
208  cU_elts[upp[k]] = mu - x2z2;
209  if (x2z2 > maxXz) maxXz = x2z2;
210  if (x2z2 < minXz) minXz = x2z2;
211  }
212 
213  maxXz = CoinMax( maxXz, 1e-99 );
214  minXz = CoinMax( minXz, 1e-99 );
215  *center = maxXz / minXz;
216 
217  double normL = 0.0;
218  double normU = 0.0;
219  for (int k = 0; k < nlow; k++)
220  if (cL_elts[low[k]] > normL) normL = cL_elts[low[k]];
221  for (int k = 0; k < nupp; k++)
222  if (cU_elts[upp[k]] > normU) normU = cU_elts[upp[k]];
223  *Cinf = CoinMax( normL, normU);
224  *Cinf0 = maxXz;
225 }
226 //-----------------------------------------------------------------------
227 // End private function pdxxxresid2
228 //-----------------------------------------------------------------------
229 
230 inline double pdxxxstep( CoinDenseVector <double> &x, CoinDenseVector <double> &dx )
231 {
232 
233 // Assumes x > 0.
234 // Finds the maximum step such that x + step*dx >= 0.
235 
236  double step = 1e+20;
237 
238  int n = x.size();
239  double *x_elts = x.getElements();
240  double *dx_elts = dx.getElements();
241  for (int k = 0; k < n; k++)
242  if (dx_elts[k] < 0)
243  if ((x_elts[k] / (-dx_elts[k])) < step)
244  step = x_elts[k] / (-dx_elts[k]);
245  return step;
246 }
247 //-----------------------------------------------------------------------
248 // End private function pdxxxstep
249 //-----------------------------------------------------------------------
250 
251 inline double pdxxxstep(int nset, int *set, CoinDenseVector <double> &x, CoinDenseVector <double> &dx )
252 {
253 
254 // Assumes x > 0.
255 // Finds the maximum step such that x + step*dx >= 0.
256 
257  double step = 1e+20;
258 
259  int n = x.size();
260  double *x_elts = x.getElements();
261  double *dx_elts = dx.getElements();
262  for (int k = 0; k < n; k++)
263  if (dx_elts[k] < 0)
264  if ((x_elts[k] / (-dx_elts[k])) < step)
265  step = x_elts[k] / (-dx_elts[k]);
266  return step;
267 }
268 //-----------------------------------------------------------------------
269 // End private function pdxxxstep
270 //-----------------------------------------------------------------------
271 #endif
272 #endif
void matVecMult(int, double *, double *)
LSQR.
void getNorms(const double *region, int size, double &norm1, double &norm2)
void setElements(double *region, int size, double value)
double maximumAbsElement(const double *region, int size)
Note (JJF) I have added some operations on arrays even though they may duplicate CoinDenseVector.
double innerProduct(const double *region1, int size, const double *region2)
This solves problems in Primal Dual Convex Optimization.
Definition: ClpPdco.hpp:22
double CoinSqrt(double x)
void multiplyAdd(const double *region1, int size, double multiplier1, double *region2, double multiplier2)