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