1 /* $Id: ClpHelperFunctions.hpp 2385 2019-01-06 19:43:06Z unxusr $ */
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
20 /**
21 Note (JJF) I have added some operations on arrays even though they may
22 duplicate CoinDenseVector. I think the use of templates was a mistake
23 as I don't think inline generic code can take as much advantage of
24 parallelism or machine architectures or memory hierarchies.
25
26 */
27
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
CoinMemcpyN(const double * from,const int size,CoinWorkDouble * to)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
CoinMemcpyN(const CoinWorkDouble * from,const int size,double * to)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
CoinMax(const CoinWorkDouble x1,const double x2)55 CoinMax(const CoinWorkDouble x1, const double x2)
56 {
57 return (x1 > x2) ? x1 : x2;
58 }
59 inline CoinWorkDouble
CoinMax(double x1,const CoinWorkDouble x2)60 CoinMax(double x1, const CoinWorkDouble x2)
61 {
62 return (x1 > x2) ? x1 : x2;
63 }
64 inline CoinWorkDouble
CoinMin(const CoinWorkDouble x1,const double x2)65 CoinMin(const CoinWorkDouble x1, const double x2)
66 {
67 return (x1 < x2) ? x1 : x2;
68 }
69 inline CoinWorkDouble
CoinMin(double x1,const CoinWorkDouble x2)70 CoinMin(double x1, const CoinWorkDouble x2)
71 {
72 return (x1 < x2) ? x1 : x2;
73 }
CoinSqrt(CoinWorkDouble x)74 inline CoinWorkDouble CoinSqrt(CoinWorkDouble x)
75 {
76 return sqrtl(x);
77 }
78 #else
CoinSqrt(double x)79 inline double CoinSqrt(double x)
80 {
81 return sqrt(x);
82 }
83 #endif
84 /// Trace
85 #ifdef NDEBUG
86 #define ClpTraceDebug(expression) \
87 { \
88 }
89 #else
90 void ClpTracePrint(std::string fileName, std::string message, int line);
91 #define ClpTraceDebug(expression) \
92 { \
93 if (!(expression)) { \
94 ClpTracePrint(__FILE__, __STRING(expression), __LINE__); \
95 } \
96 }
97 #endif
98 /// Following only included if ClpPdco defined
99 #ifdef ClpPdco_H
100
pdxxxmerit(int nlow,int nupp,int * low,int * upp,CoinDenseVector<double> & r1,CoinDenseVector<double> & r2,CoinDenseVector<double> & rL,CoinDenseVector<double> & rU,CoinDenseVector<double> & cL,CoinDenseVector<double> & cU)101 inline double pdxxxmerit(int nlow, int nupp, int *low, int *upp, CoinDenseVector< double > &r1,
102 CoinDenseVector< double > &r2, CoinDenseVector< double > &rL,
103 CoinDenseVector< double > &rU, CoinDenseVector< double > &cL,
104 CoinDenseVector< double > &cU)
105 {
106
107 // Evaluate the merit function for Newton's method.
108 // It is the 2-norm of the three sets of residuals.
109 double sum1, sum2;
110 CoinDenseVector< double > f(6);
111 f[0] = r1.twoNorm();
112 f[1] = r2.twoNorm();
113 sum1 = sum2 = 0.0;
114 for (int k = 0; k < nlow; k++) {
115 sum1 += rL[low[k]] * rL[low[k]];
116 sum2 += cL[low[k]] * cL[low[k]];
117 }
118 f[2] = sqrt(sum1);
119 f[4] = sqrt(sum2);
120 sum1 = sum2 = 0.0;
121 for (int k = 0; k < nupp; k++) {
122 sum1 += rL[upp[k]] * rL[upp[k]];
123 sum2 += cL[upp[k]] * cL[upp[k]];
124 }
125 f[3] = sqrt(sum1);
126 f[5] = sqrt(sum2);
127
128 return f.twoNorm();
129 }
130
131 //-----------------------------------------------------------------------
132 // End private function pdxxxmerit
133 //-----------------------------------------------------------------------
134
135 //function [r1,r2,rL,rU,Pinf,Dinf] = ...
136 // pdxxxresid1( Aname,fix,low,upp, ...
137 // b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 )
138
pdxxxresid1(ClpPdco * model,const int nlow,const int nupp,const int nfix,int * low,int * upp,int * fix,CoinDenseVector<double> & b,double * bl,double * bu,double d1,double d2,CoinDenseVector<double> & grad,CoinDenseVector<double> & rL,CoinDenseVector<double> & rU,CoinDenseVector<double> & x,CoinDenseVector<double> & x1,CoinDenseVector<double> & x2,CoinDenseVector<double> & y,CoinDenseVector<double> & z1,CoinDenseVector<double> & z2,CoinDenseVector<double> & r1,CoinDenseVector<double> & r2,double * Pinf,double * Dinf)139 inline void pdxxxresid1(ClpPdco *model, const int nlow, const int nupp, const int nfix,
140 int *low, int *upp, int *fix,
141 CoinDenseVector< double > &b, double *bl, double *bu, double d1, double d2,
142 CoinDenseVector< double > &grad, CoinDenseVector< double > &rL,
143 CoinDenseVector< double > &rU, CoinDenseVector< double > &x,
144 CoinDenseVector< double > &x1, CoinDenseVector< double > &x2,
145 CoinDenseVector< double > &y, CoinDenseVector< double > &z1,
146 CoinDenseVector< double > &z2, CoinDenseVector< double > &r1,
147 CoinDenseVector< double > &r2, double *Pinf, double *Dinf)
148 {
149
150 // Form residuals for the primal and dual equations.
151 // rL, rU are output, but we input them as full vectors
152 // initialized (permanently) with any relevant zeros.
153
154 // Get some element pointers for efficiency
155 double *x_elts = x.getElements();
156 double *r2_elts = r2.getElements();
157
158 for (int k = 0; k < nfix; k++)
159 x_elts[fix[k]] = 0;
160
161 r1.clear();
162 r2.clear();
163 model->matVecMult(1, r1, x);
164 model->matVecMult(2, r2, y);
165 for (int k = 0; k < nfix; k++)
166 r2_elts[fix[k]] = 0;
167
168 r1 = b - r1 - d2 * d2 * y;
169 r2 = grad - r2 - z1; // grad includes d1*d1*x
170 if (nupp > 0)
171 r2 = r2 + z2;
172
173 for (int k = 0; k < nlow; k++)
174 rL[low[k]] = bl[low[k]] - x[low[k]] + x1[low[k]];
175 for (int k = 0; k < nupp; k++)
176 rU[upp[k]] = -bu[upp[k]] + x[upp[k]] + x2[upp[k]];
177
178 double normL = 0.0;
179 double normU = 0.0;
180 for (int k = 0; k < nlow; k++)
181 if (rL[low[k]] > normL)
182 normL = rL[low[k]];
183 for (int k = 0; k < nupp; k++)
184 if (rU[upp[k]] > normU)
185 normU = rU[upp[k]];
186
187 *Pinf = CoinMax(normL, normU);
188 *Pinf = CoinMax(r1.infNorm(), *Pinf);
189 *Dinf = r2.infNorm();
190 *Pinf = CoinMax(*Pinf, 1e-99);
191 *Dinf = CoinMax(*Dinf, 1e-99);
192 }
193
194 //-----------------------------------------------------------------------
195 // End private function pdxxxresid1
196 //-----------------------------------------------------------------------
197
198 //function [cL,cU,center,Cinf,Cinf0] = ...
199 // pdxxxresid2( mu,low,upp,cL,cU,x1,x2,z1,z2 )
200
pdxxxresid2(double mu,int nlow,int nupp,int * low,int * upp,CoinDenseVector<double> & cL,CoinDenseVector<double> & cU,CoinDenseVector<double> & x1,CoinDenseVector<double> & x2,CoinDenseVector<double> & z1,CoinDenseVector<double> & z2,double * center,double * Cinf,double * Cinf0)201 inline void pdxxxresid2(double mu, int nlow, int nupp, int *low, int *upp,
202 CoinDenseVector< double > &cL, CoinDenseVector< double > &cU,
203 CoinDenseVector< double > &x1, CoinDenseVector< double > &x2,
204 CoinDenseVector< double > &z1, CoinDenseVector< double > &z2,
205 double *center, double *Cinf, double *Cinf0)
206 {
207
208 // Form residuals for the complementarity equations.
209 // cL, cU are output, but we input them as full vectors
210 // initialized (permanently) with any relevant zeros.
211 // Cinf is the complementarity residual for X1 z1 = mu e, etc.
212 // Cinf0 is the same for mu=0 (i.e., for the original problem).
213
214 double maxXz = -1e20;
215 double minXz = 1e20;
216
217 double *x1_elts = x1.getElements();
218 double *z1_elts = z1.getElements();
219 double *cL_elts = cL.getElements();
220 for (int k = 0; k < nlow; k++) {
221 double x1z1 = x1_elts[low[k]] * z1_elts[low[k]];
222 cL_elts[low[k]] = mu - x1z1;
223 if (x1z1 > maxXz)
224 maxXz = x1z1;
225 if (x1z1 < minXz)
226 minXz = x1z1;
227 }
228
229 double *x2_elts = x2.getElements();
230 double *z2_elts = z2.getElements();
231 double *cU_elts = cU.getElements();
232 for (int k = 0; k < nupp; k++) {
233 double x2z2 = x2_elts[upp[k]] * z2_elts[upp[k]];
234 cU_elts[upp[k]] = mu - x2z2;
235 if (x2z2 > maxXz)
236 maxXz = x2z2;
237 if (x2z2 < minXz)
238 minXz = x2z2;
239 }
240
241 maxXz = CoinMax(maxXz, 1e-99);
242 minXz = CoinMax(minXz, 1e-99);
243 *center = maxXz / minXz;
244
245 double normL = 0.0;
246 double normU = 0.0;
247 for (int k = 0; k < nlow; k++)
248 if (cL_elts[low[k]] > normL)
249 normL = cL_elts[low[k]];
250 for (int k = 0; k < nupp; k++)
251 if (cU_elts[upp[k]] > normU)
252 normU = cU_elts[upp[k]];
253 *Cinf = CoinMax(normL, normU);
254 *Cinf0 = maxXz;
255 }
256 //-----------------------------------------------------------------------
257 // End private function pdxxxresid2
258 //-----------------------------------------------------------------------
259
pdxxxstep(CoinDenseVector<double> & x,CoinDenseVector<double> & dx)260 inline double pdxxxstep(CoinDenseVector< double > &x, CoinDenseVector< double > &dx)
261 {
262
263 // Assumes x > 0.
264 // Finds the maximum step such that x + step*dx >= 0.
265
266 double step = 1e+20;
267
268 int n = x.size();
269 double *x_elts = x.getElements();
270 double *dx_elts = dx.getElements();
271 for (int k = 0; k < n; k++)
272 if (dx_elts[k] < 0)
273 if ((x_elts[k] / (-dx_elts[k])) < step)
274 step = x_elts[k] / (-dx_elts[k]);
275 return step;
276 }
277 //-----------------------------------------------------------------------
278 // End private function pdxxxstep
279 //-----------------------------------------------------------------------
280
pdxxxstep(int nset,int * set,CoinDenseVector<double> & x,CoinDenseVector<double> & dx)281 inline double pdxxxstep(int nset, int *set, CoinDenseVector< double > &x, CoinDenseVector< double > &dx)
282 {
283
284 // Assumes x > 0.
285 // Finds the maximum step such that x + step*dx >= 0.
286
287 double step = 1e+20;
288
289 int n = x.size();
290 double *x_elts = x.getElements();
291 double *dx_elts = dx.getElements();
292 for (int k = 0; k < n; k++)
293 if (dx_elts[k] < 0)
294 if ((x_elts[k] / (-dx_elts[k])) < step)
295 step = x_elts[k] / (-dx_elts[k]);
296 return step;
297 }
298 //-----------------------------------------------------------------------
299 // End private function pdxxxstep
300 //-----------------------------------------------------------------------
301 #endif
302 #endif
303
304 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
305 */
306