source: trunk/Clp/src/ClpHelperFunctions.hpp @ 2470

Last change on this file since 2470 was 2385, checked in by unxusr, 10 months ago

formatting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 9.3 KB
Line 
1/* $Id: ClpHelperFunctions.hpp 2385 2019-01-06 19:43:06Z 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
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
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
84/// Trace
85#ifdef NDEBUG
86#define ClpTraceDebug(expression) \
87  {                               \
88  }
89#else
90void 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
101inline 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
139inline 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
201inline 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
260inline 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
281inline 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*/
Note: See TracBrowser for help on using the repository browser.