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