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

Last change on this file since 1722 was 1722, checked in by stefan, 10 years ago

adjust to changes in CoinUtils? header files

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 9.5 KB
Line 
1/* $Id: ClpHelperFunctions.hpp 1722 2011-04-17 09:58:37Z 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
85/// Following only included if ClpPdco defined
86#ifdef ClpPdco_H
87
88
89inline double pdxxxmerit(int nlow, int nupp, int *low, int *upp, CoinDenseVector <double> &r1,
90                         CoinDenseVector <double> &r2, CoinDenseVector <double> &rL,
91                         CoinDenseVector <double> &rU, CoinDenseVector <double> &cL,
92                         CoinDenseVector <double> &cU )
93{
94
95// Evaluate the merit function for Newton's method.
96// It is the 2-norm of the three sets of residuals.
97     double sum1, sum2;
98     CoinDenseVector <double> f(6);
99     f[0] = r1.twoNorm();
100     f[1] = r2.twoNorm();
101     sum1 = sum2 = 0.0;
102     for (int k = 0; k < nlow; k++) {
103          sum1 += rL[low[k]] * rL[low[k]];
104          sum2 += cL[low[k]] * cL[low[k]];
105     }
106     f[2] = sqrt(sum1);
107     f[4] = sqrt(sum2);
108     sum1 = sum2 = 0.0;
109     for (int k = 0; k < nupp; k++) {
110          sum1 += rL[upp[k]] * rL[upp[k]];
111          sum2 += cL[upp[k]] * cL[upp[k]];
112     }
113     f[3] = sqrt(sum1);
114     f[5] = sqrt(sum2);
115
116     return f.twoNorm();
117}
118
119//-----------------------------------------------------------------------
120// End private function pdxxxmerit
121//-----------------------------------------------------------------------
122
123
124//function [r1,r2,rL,rU,Pinf,Dinf] =    ...
125//      pdxxxresid1( Aname,fix,low,upp, ...
126//                   b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 )
127
128inline void pdxxxresid1(ClpPdco *model, const int nlow, const int nupp, const int nfix,
129                        int *low, int *upp, int *fix,
130                        CoinDenseVector <double> &b, double *bl, double *bu, double d1, double d2,
131                        CoinDenseVector <double> &grad, CoinDenseVector <double> &rL,
132                        CoinDenseVector <double> &rU, CoinDenseVector <double> &x,
133                        CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
134                        CoinDenseVector <double> &y,  CoinDenseVector <double> &z1,
135                        CoinDenseVector <double> &z2, CoinDenseVector <double> &r1,
136                        CoinDenseVector <double> &r2, double *Pinf, double *Dinf)
137{
138
139// Form residuals for the primal and dual equations.
140// rL, rU are output, but we input them as full vectors
141// initialized (permanently) with any relevant zeros.
142
143// Get some element pointers for efficiency
144     double *x_elts  = x.getElements();
145     double *r2_elts = r2.getElements();
146
147     for (int k = 0; k < nfix; k++)
148          x_elts[fix[k]]  = 0;
149
150     r1.clear();
151     r2.clear();
152     model->matVecMult( 1, r1, x );
153     model->matVecMult( 2, r2, y );
154     for (int k = 0; k < nfix; k++)
155          r2_elts[fix[k]]  = 0;
156
157
158     r1      = b    - r1 - d2 * d2 * y;
159     r2      = grad - r2 - z1;              // grad includes d1*d1*x
160     if (nupp > 0)
161          r2    = r2 + z2;
162
163     for (int k = 0; k < nlow; k++)
164          rL[low[k]] = bl[low[k]] - x[low[k]] + x1[low[k]];
165     for (int k = 0; k < nupp; k++)
166          rU[upp[k]] = - bu[upp[k]] + x[upp[k]] + x2[upp[k]];
167
168     double normL = 0.0;
169     double normU = 0.0;
170     for (int k = 0; k < nlow; k++)
171          if (rL[low[k]] > normL) normL = rL[low[k]];
172     for (int k = 0; k < nupp; k++)
173          if (rU[upp[k]] > normU) normU = rU[upp[k]];
174
175     *Pinf    = CoinMax(normL, normU);
176     *Pinf    = CoinMax( r1.infNorm() , *Pinf );
177     *Dinf    = r2.infNorm();
178     *Pinf    = CoinMax( *Pinf, 1e-99 );
179     *Dinf    = CoinMax( *Dinf, 1e-99 );
180}
181
182//-----------------------------------------------------------------------
183// End private function pdxxxresid1
184//-----------------------------------------------------------------------
185
186
187//function [cL,cU,center,Cinf,Cinf0] = ...
188//      pdxxxresid2( mu,low,upp,cL,cU,x1,x2,z1,z2 )
189
190inline void pdxxxresid2(double mu, int nlow, int nupp, int *low, int *upp,
191                        CoinDenseVector <double> &cL, CoinDenseVector <double> &cU,
192                        CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
193                        CoinDenseVector <double> &z1, CoinDenseVector <double> &z2,
194                        double *center, double *Cinf, double *Cinf0)
195{
196
197// Form residuals for the complementarity equations.
198// cL, cU are output, but we input them as full vectors
199// initialized (permanently) with any relevant zeros.
200// Cinf  is the complementarity residual for X1 z1 = mu e, etc.
201// Cinf0 is the same for mu=0 (i.e., for the original problem).
202
203     double maxXz = -1e20;
204     double minXz = 1e20;
205
206     double *x1_elts = x1.getElements();
207     double *z1_elts = z1.getElements();
208     double *cL_elts = cL.getElements();
209     for (int k = 0; k < nlow; k++) {
210          double x1z1    = x1_elts[low[k]] * z1_elts[low[k]];
211          cL_elts[low[k]] = mu - x1z1;
212          if (x1z1 > maxXz) maxXz = x1z1;
213          if (x1z1 < minXz) minXz = x1z1;
214     }
215
216     double *x2_elts = x2.getElements();
217     double *z2_elts = z2.getElements();
218     double *cU_elts = cU.getElements();
219     for (int k = 0; k < nupp; k++) {
220          double x2z2    = x2_elts[upp[k]] * z2_elts[upp[k]];
221          cU_elts[upp[k]] = mu - x2z2;
222          if (x2z2 > maxXz) maxXz = x2z2;
223          if (x2z2 < minXz) minXz = x2z2;
224     }
225
226     maxXz   = CoinMax( maxXz, 1e-99 );
227     minXz   = CoinMax( minXz, 1e-99 );
228     *center  = maxXz / minXz;
229
230     double normL = 0.0;
231     double normU = 0.0;
232     for (int k = 0; k < nlow; k++)
233          if (cL_elts[low[k]] > normL) normL = cL_elts[low[k]];
234     for (int k = 0; k < nupp; k++)
235          if (cU_elts[upp[k]] > normU) normU = cU_elts[upp[k]];
236     *Cinf    = CoinMax( normL, normU);
237     *Cinf0   = maxXz;
238}
239//-----------------------------------------------------------------------
240// End private function pdxxxresid2
241//-----------------------------------------------------------------------
242
243inline double  pdxxxstep( CoinDenseVector <double> &x, CoinDenseVector <double> &dx )
244{
245
246// Assumes x > 0.
247// Finds the maximum step such that x + step*dx >= 0.
248
249     double step     = 1e+20;
250
251     int n = x.size();
252     double *x_elts = x.getElements();
253     double *dx_elts = dx.getElements();
254     for (int k = 0; k < n; k++)
255          if (dx_elts[k] < 0)
256               if ((x_elts[k] / (-dx_elts[k])) < step)
257                    step = x_elts[k] / (-dx_elts[k]);
258     return step;
259}
260//-----------------------------------------------------------------------
261// End private function pdxxxstep
262//-----------------------------------------------------------------------
263
264inline double  pdxxxstep(int nset, int *set, CoinDenseVector <double> &x, CoinDenseVector <double> &dx )
265{
266
267// Assumes x > 0.
268// Finds the maximum step such that x + step*dx >= 0.
269
270     double step     = 1e+20;
271
272     int n = x.size();
273     double *x_elts = x.getElements();
274     double *dx_elts = dx.getElements();
275     for (int k = 0; k < n; k++)
276          if (dx_elts[k] < 0)
277               if ((x_elts[k] / (-dx_elts[k])) < step)
278                    step = x_elts[k] / (-dx_elts[k]);
279     return step;
280}
281//-----------------------------------------------------------------------
282// End private function pdxxxstep
283//-----------------------------------------------------------------------
284#endif
285#endif
Note: See TracBrowser for help on using the repository browser.