Clp 1.17.10
Loading...
Searching...
No Matches
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
29double maximumAbsElement(const double *region, int size);
30void setElements(double *region, int size, double value);
31void multiplyAdd(const double *region1, int size, double multiplier1,
32 double *region2, double multiplier2);
33double innerProduct(const double *region1, int size, const double *region2);
34void getNorms(const double *region, int size, double &norm1, double &norm2);
35#if COIN_LONG_WORK
36// For long double versions
37CoinWorkDouble maximumAbsElement(const CoinWorkDouble *region, int size);
38void setElements(CoinWorkDouble *region, int size, CoinWorkDouble value);
39void multiplyAdd(const CoinWorkDouble *region1, int size, CoinWorkDouble multiplier1,
40 CoinWorkDouble *region2, CoinWorkDouble multiplier2);
41CoinWorkDouble innerProduct(const CoinWorkDouble *region1, int size, const CoinWorkDouble *region2);
42void getNorms(const CoinWorkDouble *region, int size, CoinWorkDouble &norm1, CoinWorkDouble &norm2);
43inline void
44CoinMemcpyN(const double *from, const int size, CoinWorkDouble *to)
45{
46 for (int i = 0; i < size; i++)
47 to[i] = from[i];
48}
49inline void
50CoinMemcpyN(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}
55inline CoinWorkDouble
56CoinMax(const CoinWorkDouble x1, const double x2)
57{
58 return (x1 > x2) ? x1 : x2;
59}
60inline CoinWorkDouble
61CoinMax(double x1, const CoinWorkDouble x2)
62{
63 return (x1 > x2) ? x1 : x2;
64}
65inline CoinWorkDouble
66CoinMin(const CoinWorkDouble x1, const double x2)
67{
68 return (x1 < x2) ? x1 : x2;
69}
70inline CoinWorkDouble
71CoinMin(double x1, const CoinWorkDouble x2)
72{
73 return (x1 < x2) ? x1 : x2;
74}
75inline CoinWorkDouble CoinSqrt(CoinWorkDouble x)
76{
77 return sqrtl(x);
78}
79#else
80inline double CoinSqrt(double x)
81{
82 return sqrt(x);
83}
84#endif
86#ifdef NDEBUG
87#define ClpTraceDebug(expression) \
88 { \
89 }
90#else
91void 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
100#ifdef ClpPdco_H
101
102inline 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
140inline 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
202inline 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
261inline 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
282inline 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*/
double CoinSqrt(double x)
double innerProduct(const double *region1, int size, const double *region2)
double maximumAbsElement(const double *region, int size)
Note (JJF) I have added some operations on arrays even though they may duplicate CoinDenseVector.
void getNorms(const double *region, int size, double &norm1, double &norm2)
void ClpTracePrint(std::string fileName, std::string message, int line)
Trace.
void multiplyAdd(const double *region1, int size, double multiplier1, double *region2, double multiplier2)
void setElements(double *region, int size, double value)
This solves problems in Primal Dual Convex Optimization.
Definition ClpPdco.hpp:22
void matVecMult(int, double *, double *)