ClpHelperFunctions.hpp
Go to the documentation of this file.
1 /* $Id: ClpHelperFunctions.hpp 2034 2014-06-12 16:10:26Z stefan $ */
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 
28 double maximumAbsElement(const double * region, int size);
29 void setElements(double * region, int size, double value);
30 void multiplyAdd(const double * region1, int size, double multiplier1,
31  double * region2, double multiplier2);
32 double innerProduct(const double * region1, int size, const double * region2);
33 void getNorms(const double * region, int size, double & norm1, double & norm2);
34 #if COIN_LONG_WORK
35 // For long double versions
36 CoinWorkDouble maximumAbsElement(const CoinWorkDouble * region, int size);
37 void setElements(CoinWorkDouble * region, int size, CoinWorkDouble value);
38 void multiplyAdd(const CoinWorkDouble * region1, int size, CoinWorkDouble multiplier1,
39  CoinWorkDouble * region2, CoinWorkDouble multiplier2);
40 CoinWorkDouble innerProduct(const CoinWorkDouble * region1, int size, const CoinWorkDouble * region2);
41 void getNorms(const CoinWorkDouble * region, int size, CoinWorkDouble & norm1, CoinWorkDouble & norm2);
42 inline void
43 CoinMemcpyN(const double * from, const int size, CoinWorkDouble * to)
44 {
45  for (int i = 0; i < size; i++)
46  to[i] = from[i];
47 }
48 inline void
49 CoinMemcpyN(const CoinWorkDouble * from, const int size, double * to)
50 {
51  for (int i = 0; i < size; i++)
52  to[i] = static_cast<double>(from[i]);
53 }
54 inline CoinWorkDouble
55 CoinMax(const CoinWorkDouble x1, const double x2)
56 {
57  return (x1 > x2) ? x1 : x2;
58 }
59 inline CoinWorkDouble
60 CoinMax(double x1, const CoinWorkDouble x2)
61 {
62  return (x1 > x2) ? x1 : x2;
63 }
64 inline CoinWorkDouble
65 CoinMin(const CoinWorkDouble x1, const double x2)
66 {
67  return (x1 < x2) ? x1 : x2;
68 }
69 inline CoinWorkDouble
70 CoinMin(double x1, const CoinWorkDouble x2)
71 {
72  return (x1 < x2) ? x1 : x2;
73 }
74 inline CoinWorkDouble CoinSqrt(CoinWorkDouble x)
75 {
76  return sqrtl(x);
77 }
78 #else
79 inline double CoinSqrt(double x)
80 {
81  return sqrt(x);
82 }
83 #endif
84 # ifdef NDEBUG
86 # define ClpTraceDebug(expression) {}
87 # else
88 void ClpTracePrint(std::string fileName, std::string message, int line);
89 # define ClpTraceDebug(expression) { \
90  if (!(expression)) { ClpTracePrint(__FILE__,__STRING(expression),__LINE__); } \
91  }
92 # endif
93 #ifdef ClpPdco_H
95 
96 
97 inline double pdxxxmerit(int nlow, int nupp, int *low, int *upp, CoinDenseVector <double> &r1,
98  CoinDenseVector <double> &r2, CoinDenseVector <double> &rL,
99  CoinDenseVector <double> &rU, CoinDenseVector <double> &cL,
100  CoinDenseVector <double> &cU )
101 {
102 
103 // Evaluate the merit function for Newton's method.
104 // It is the 2-norm of the three sets of residuals.
105  double sum1, sum2;
106  CoinDenseVector <double> f(6);
107  f[0] = r1.twoNorm();
108  f[1] = r2.twoNorm();
109  sum1 = sum2 = 0.0;
110  for (int k = 0; k < nlow; k++) {
111  sum1 += rL[low[k]] * rL[low[k]];
112  sum2 += cL[low[k]] * cL[low[k]];
113  }
114  f[2] = sqrt(sum1);
115  f[4] = sqrt(sum2);
116  sum1 = sum2 = 0.0;
117  for (int k = 0; k < nupp; k++) {
118  sum1 += rL[upp[k]] * rL[upp[k]];
119  sum2 += cL[upp[k]] * cL[upp[k]];
120  }
121  f[3] = sqrt(sum1);
122  f[5] = sqrt(sum2);
123 
124  return f.twoNorm();
125 }
126 
127 //-----------------------------------------------------------------------
128 // End private function pdxxxmerit
129 //-----------------------------------------------------------------------
130 
131 
132 //function [r1,r2,rL,rU,Pinf,Dinf] = ...
133 // pdxxxresid1( Aname,fix,low,upp, ...
134 // b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 )
135 
136 inline void pdxxxresid1(ClpPdco *model, const int nlow, const int nupp, const int nfix,
137  int *low, int *upp, int *fix,
138  CoinDenseVector <double> &b, double *bl, double *bu, double d1, double d2,
139  CoinDenseVector <double> &grad, CoinDenseVector <double> &rL,
140  CoinDenseVector <double> &rU, CoinDenseVector <double> &x,
141  CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
142  CoinDenseVector <double> &y, CoinDenseVector <double> &z1,
143  CoinDenseVector <double> &z2, CoinDenseVector <double> &r1,
144  CoinDenseVector <double> &r2, double *Pinf, double *Dinf)
145 {
146 
147 // Form residuals for the primal and dual equations.
148 // rL, rU are output, but we input them as full vectors
149 // initialized (permanently) with any relevant zeros.
150 
151 // Get some element pointers for efficiency
152  double *x_elts = x.getElements();
153  double *r2_elts = r2.getElements();
154 
155  for (int k = 0; k < nfix; k++)
156  x_elts[fix[k]] = 0;
157 
158  r1.clear();
159  r2.clear();
160  model->matVecMult( 1, r1, x );
161  model->matVecMult( 2, r2, y );
162  for (int k = 0; k < nfix; k++)
163  r2_elts[fix[k]] = 0;
164 
165 
166  r1 = b - r1 - d2 * d2 * y;
167  r2 = grad - r2 - z1; // grad includes d1*d1*x
168  if (nupp > 0)
169  r2 = r2 + z2;
170 
171  for (int k = 0; k < nlow; k++)
172  rL[low[k]] = bl[low[k]] - x[low[k]] + x1[low[k]];
173  for (int k = 0; k < nupp; k++)
174  rU[upp[k]] = - bu[upp[k]] + x[upp[k]] + x2[upp[k]];
175 
176  double normL = 0.0;
177  double normU = 0.0;
178  for (int k = 0; k < nlow; k++)
179  if (rL[low[k]] > normL) normL = rL[low[k]];
180  for (int k = 0; k < nupp; k++)
181  if (rU[upp[k]] > normU) normU = rU[upp[k]];
182 
183  *Pinf = CoinMax(normL, normU);
184  *Pinf = CoinMax( r1.infNorm() , *Pinf );
185  *Dinf = r2.infNorm();
186  *Pinf = CoinMax( *Pinf, 1e-99 );
187  *Dinf = CoinMax( *Dinf, 1e-99 );
188 }
189 
190 //-----------------------------------------------------------------------
191 // End private function pdxxxresid1
192 //-----------------------------------------------------------------------
193 
194 
195 //function [cL,cU,center,Cinf,Cinf0] = ...
196 // pdxxxresid2( mu,low,upp,cL,cU,x1,x2,z1,z2 )
197 
198 inline void pdxxxresid2(double mu, int nlow, int nupp, int *low, int *upp,
199  CoinDenseVector <double> &cL, CoinDenseVector <double> &cU,
200  CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
201  CoinDenseVector <double> &z1, CoinDenseVector <double> &z2,
202  double *center, double *Cinf, double *Cinf0)
203 {
204 
205 // Form residuals for the complementarity equations.
206 // cL, cU are output, but we input them as full vectors
207 // initialized (permanently) with any relevant zeros.
208 // Cinf is the complementarity residual for X1 z1 = mu e, etc.
209 // Cinf0 is the same for mu=0 (i.e., for the original problem).
210 
211  double maxXz = -1e20;
212  double minXz = 1e20;
213 
214  double *x1_elts = x1.getElements();
215  double *z1_elts = z1.getElements();
216  double *cL_elts = cL.getElements();
217  for (int k = 0; k < nlow; k++) {
218  double x1z1 = x1_elts[low[k]] * z1_elts[low[k]];
219  cL_elts[low[k]] = mu - x1z1;
220  if (x1z1 > maxXz) maxXz = x1z1;
221  if (x1z1 < minXz) minXz = x1z1;
222  }
223 
224  double *x2_elts = x2.getElements();
225  double *z2_elts = z2.getElements();
226  double *cU_elts = cU.getElements();
227  for (int k = 0; k < nupp; k++) {
228  double x2z2 = x2_elts[upp[k]] * z2_elts[upp[k]];
229  cU_elts[upp[k]] = mu - x2z2;
230  if (x2z2 > maxXz) maxXz = x2z2;
231  if (x2z2 < minXz) minXz = x2z2;
232  }
233 
234  maxXz = CoinMax( maxXz, 1e-99 );
235  minXz = CoinMax( minXz, 1e-99 );
236  *center = maxXz / minXz;
237 
238  double normL = 0.0;
239  double normU = 0.0;
240  for (int k = 0; k < nlow; k++)
241  if (cL_elts[low[k]] > normL) normL = cL_elts[low[k]];
242  for (int k = 0; k < nupp; k++)
243  if (cU_elts[upp[k]] > normU) normU = cU_elts[upp[k]];
244  *Cinf = CoinMax( normL, normU);
245  *Cinf0 = maxXz;
246 }
247 //-----------------------------------------------------------------------
248 // End private function pdxxxresid2
249 //-----------------------------------------------------------------------
250 
251 inline double pdxxxstep( 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 
272 inline double pdxxxstep(int nset, int *set, CoinDenseVector <double> &x, CoinDenseVector <double> &dx )
273 {
274 
275 // Assumes x > 0.
276 // Finds the maximum step such that x + step*dx >= 0.
277 
278  double step = 1e+20;
279 
280  int n = x.size();
281  double *x_elts = x.getElements();
282  double *dx_elts = dx.getElements();
283  for (int k = 0; k < n; k++)
284  if (dx_elts[k] < 0)
285  if ((x_elts[k] / (-dx_elts[k])) < step)
286  step = x_elts[k] / (-dx_elts[k]);
287  return step;
288 }
289 //-----------------------------------------------------------------------
290 // End private function pdxxxstep
291 //-----------------------------------------------------------------------
292 #endif
293 #endif
CoinMax
T CoinMax(register const T x1, register const T x2)
Return the larger (according to operator&lt;() of the arguments.
Definition: CoinHelperFunctions.hpp:548
ClpPdco::matVecMult
void matVecMult(int, double *, double *)
LSQR.
getNorms
void getNorms(const double *region, int size, double &norm1, double &norm2)
CoinDenseVector::getElements
const T * getElements() const
Get element values.
Definition: CoinDenseVector.hpp:84
CoinDenseVector
Dense Vector.
Definition: CoinDenseVector.hpp:67
setElements
void setElements(double *region, int size, double value)
CoinDenseVector::infNorm
T infNorm() const
infinity-norm of vector
Definition: CoinDenseVector.hpp:144
ClpConfig.h
ClpTracePrint
void ClpTracePrint(std::string fileName, std::string message, int line)
Trace.
CoinWorkDouble
double CoinWorkDouble
Definition: CoinTypes.hpp:53
maximumAbsElement
double maximumAbsElement(const double *region, int size)
Note (JJF) I have added some operations on arrays even though they may duplicate CoinDenseVector.
CoinMemcpyN
void CoinMemcpyN(register const T *from, const int size, register T *to)
This helper function copies an array to another location.
Definition: CoinHelperFunctions.hpp:259
CoinDenseVector::clear
void clear()
Reset the vector (i.e. set all elemenets to zero)
innerProduct
double innerProduct(const double *region1, int size, const double *region2)
CoinDenseVector::size
int size() const
Get the size.
Definition: CoinDenseVector.hpp:82
ClpPdco
This solves problems in Primal Dual Convex Optimization.
Definition: ClpPdco.hpp:22
CoinMin
T CoinMin(register const T x1, register const T x2)
Return the smaller (according to operator&lt;() of the arguments.
Definition: CoinHelperFunctions.hpp:559
CoinSqrt
double CoinSqrt(double x)
Definition: ClpHelperFunctions.hpp:79
CoinDenseVector::twoNorm
double twoNorm() const
2-norm of vector
Definition: CoinDenseVector.hpp:135
multiplyAdd
void multiplyAdd(const double *region1, int size, double multiplier1, double *region2, double multiplier2)
Generated by   doxygen 1.8.5