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

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 9.3 KB
Line 
1/* $Id: ClpHelperFunctions.hpp 1525 2010-02-26 17:27:59Z mjs $ */
2// Copyright (C) 2003, International Business Machines
3// Corporation and others.  All Rights Reserved.
4#ifndef ClpHelperFunctions_H
5#define ClpHelperFunctions_H
6
7/**
8    Note (JJF) I have added some operations on arrays even though they may
9    duplicate CoinDenseVector.  I think the use of templates was a mistake
10    as I don't think inline generic code can take as much advantage of
11    parallelism or machine architectures or memory hierarchies.
12
13*/
14
15double maximumAbsElement(const double * region, int size);
16void setElements(double * region, int size, double value);
17void multiplyAdd(const double * region1, int size, double multiplier1,
18                 double * region2, double multiplier2);
19double innerProduct(const double * region1, int size, const double * region2);
20void getNorms(const double * region, int size, double & norm1, double & norm2);
21#if COIN_LONG_WORK
22// For long double versions
23CoinWorkDouble maximumAbsElement(const CoinWorkDouble * region, int size);
24void setElements(CoinWorkDouble * region, int size, CoinWorkDouble value);
25void multiplyAdd(const CoinWorkDouble * region1, int size, CoinWorkDouble multiplier1,
26                 CoinWorkDouble * region2, CoinWorkDouble multiplier2);
27CoinWorkDouble innerProduct(const CoinWorkDouble * region1, int size, const CoinWorkDouble * region2);
28void getNorms(const CoinWorkDouble * region, int size, CoinWorkDouble & norm1, CoinWorkDouble & norm2);
29inline void
30CoinMemcpyN(const double * from, const int size, CoinWorkDouble * to)
31{
32     for (int i = 0; i < size; i++)
33          to[i] = from[i];
34}
35inline void
36CoinMemcpyN(const CoinWorkDouble * from, const int size, double * to)
37{
38     for (int i = 0; i < size; i++)
39          to[i] = static_cast<double>(from[i]);
40}
41inline CoinWorkDouble
42CoinMax(const CoinWorkDouble x1, const double x2)
43{
44     return (x1 > x2) ? x1 : x2;
45}
46inline CoinWorkDouble
47CoinMax(double x1, const CoinWorkDouble x2)
48{
49     return (x1 > x2) ? x1 : x2;
50}
51inline CoinWorkDouble
52CoinMin(const CoinWorkDouble x1, const double x2)
53{
54     return (x1 < x2) ? x1 : x2;
55}
56inline CoinWorkDouble
57CoinMin(double x1, const CoinWorkDouble x2)
58{
59     return (x1 < x2) ? x1 : x2;
60}
61inline CoinWorkDouble CoinSqrt(CoinWorkDouble x)
62{
63     return sqrtl(x);
64}
65#else
66inline double CoinSqrt(double x)
67{
68     return sqrt(x);
69}
70#endif
71
72/// Following only included if ClpPdco defined
73#ifdef ClpPdco_H
74
75
76inline double pdxxxmerit(int nlow, int nupp, int *low, int *upp, CoinDenseVector <double> &r1,
77                         CoinDenseVector <double> &r2, CoinDenseVector <double> &rL,
78                         CoinDenseVector <double> &rU, CoinDenseVector <double> &cL,
79                         CoinDenseVector <double> &cU )
80{
81
82// Evaluate the merit function for Newton's method.
83// It is the 2-norm of the three sets of residuals.
84     double sum1, sum2;
85     CoinDenseVector <double> f(6);
86     f[0] = r1.twoNorm();
87     f[1] = r2.twoNorm();
88     sum1 = sum2 = 0.0;
89     for (int k = 0; k < nlow; k++) {
90          sum1 += rL[low[k]] * rL[low[k]];
91          sum2 += cL[low[k]] * cL[low[k]];
92     }
93     f[2] = sqrt(sum1);
94     f[4] = sqrt(sum2);
95     sum1 = sum2 = 0.0;
96     for (int k = 0; k < nupp; k++) {
97          sum1 += rL[upp[k]] * rL[upp[k]];
98          sum2 += cL[upp[k]] * cL[upp[k]];
99     }
100     f[3] = sqrt(sum1);
101     f[5] = sqrt(sum2);
102
103     return f.twoNorm();
104}
105
106//-----------------------------------------------------------------------
107// End private function pdxxxmerit
108//-----------------------------------------------------------------------
109
110
111//function [r1,r2,rL,rU,Pinf,Dinf] =    ...
112//      pdxxxresid1( Aname,fix,low,upp, ...
113//                   b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 )
114
115inline void pdxxxresid1(ClpPdco *model, const int nlow, const int nupp, const int nfix,
116                        int *low, int *upp, int *fix,
117                        CoinDenseVector <double> &b, double *bl, double *bu, double d1, double d2,
118                        CoinDenseVector <double> &grad, CoinDenseVector <double> &rL,
119                        CoinDenseVector <double> &rU, CoinDenseVector <double> &x,
120                        CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
121                        CoinDenseVector <double> &y,  CoinDenseVector <double> &z1,
122                        CoinDenseVector <double> &z2, CoinDenseVector <double> &r1,
123                        CoinDenseVector <double> &r2, double *Pinf, double *Dinf)
124{
125
126// Form residuals for the primal and dual equations.
127// rL, rU are output, but we input them as full vectors
128// initialized (permanently) with any relevant zeros.
129
130// Get some element pointers for efficiency
131     double *x_elts  = x.getElements();
132     double *r2_elts = r2.getElements();
133
134     for (int k = 0; k < nfix; k++)
135          x_elts[fix[k]]  = 0;
136
137     r1.clear();
138     r2.clear();
139     model->matVecMult( 1, r1, x );
140     model->matVecMult( 2, r2, y );
141     for (int k = 0; k < nfix; k++)
142          r2_elts[fix[k]]  = 0;
143
144
145     r1      = b    - r1 - d2 * d2 * y;
146     r2      = grad - r2 - z1;              // grad includes d1*d1*x
147     if (nupp > 0)
148          r2    = r2 + z2;
149
150     for (int k = 0; k < nlow; k++)
151          rL[low[k]] = bl[low[k]] - x[low[k]] + x1[low[k]];
152     for (int k = 0; k < nupp; k++)
153          rU[upp[k]] = - bu[upp[k]] + x[upp[k]] + x2[upp[k]];
154
155     double normL = 0.0;
156     double normU = 0.0;
157     for (int k = 0; k < nlow; k++)
158          if (rL[low[k]] > normL) normL = rL[low[k]];
159     for (int k = 0; k < nupp; k++)
160          if (rU[upp[k]] > normU) normU = rU[upp[k]];
161
162     *Pinf    = CoinMax(normL, normU);
163     *Pinf    = CoinMax( r1.infNorm() , *Pinf );
164     *Dinf    = r2.infNorm();
165     *Pinf    = CoinMax( *Pinf, 1e-99 );
166     *Dinf    = CoinMax( *Dinf, 1e-99 );
167}
168
169//-----------------------------------------------------------------------
170// End private function pdxxxresid1
171//-----------------------------------------------------------------------
172
173
174//function [cL,cU,center,Cinf,Cinf0] = ...
175//      pdxxxresid2( mu,low,upp,cL,cU,x1,x2,z1,z2 )
176
177inline void pdxxxresid2(double mu, int nlow, int nupp, int *low, int *upp,
178                        CoinDenseVector <double> &cL, CoinDenseVector <double> &cU,
179                        CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
180                        CoinDenseVector <double> &z1, CoinDenseVector <double> &z2,
181                        double *center, double *Cinf, double *Cinf0)
182{
183
184// Form residuals for the complementarity equations.
185// cL, cU are output, but we input them as full vectors
186// initialized (permanently) with any relevant zeros.
187// Cinf  is the complementarity residual for X1 z1 = mu e, etc.
188// Cinf0 is the same for mu=0 (i.e., for the original problem).
189
190     double maxXz = -1e20;
191     double minXz = 1e20;
192
193     double *x1_elts = x1.getElements();
194     double *z1_elts = z1.getElements();
195     double *cL_elts = cL.getElements();
196     for (int k = 0; k < nlow; k++) {
197          double x1z1    = x1_elts[low[k]] * z1_elts[low[k]];
198          cL_elts[low[k]] = mu - x1z1;
199          if (x1z1 > maxXz) maxXz = x1z1;
200          if (x1z1 < minXz) minXz = x1z1;
201     }
202
203     double *x2_elts = x2.getElements();
204     double *z2_elts = z2.getElements();
205     double *cU_elts = cU.getElements();
206     for (int k = 0; k < nupp; k++) {
207          double x2z2    = x2_elts[upp[k]] * z2_elts[upp[k]];
208          cU_elts[upp[k]] = mu - x2z2;
209          if (x2z2 > maxXz) maxXz = x2z2;
210          if (x2z2 < minXz) minXz = x2z2;
211     }
212
213     maxXz   = CoinMax( maxXz, 1e-99 );
214     minXz   = CoinMax( minXz, 1e-99 );
215     *center  = maxXz / minXz;
216
217     double normL = 0.0;
218     double normU = 0.0;
219     for (int k = 0; k < nlow; k++)
220          if (cL_elts[low[k]] > normL) normL = cL_elts[low[k]];
221     for (int k = 0; k < nupp; k++)
222          if (cU_elts[upp[k]] > normU) normU = cU_elts[upp[k]];
223     *Cinf    = CoinMax( normL, normU);
224     *Cinf0   = maxXz;
225}
226//-----------------------------------------------------------------------
227// End private function pdxxxresid2
228//-----------------------------------------------------------------------
229
230inline double  pdxxxstep( CoinDenseVector <double> &x, CoinDenseVector <double> &dx )
231{
232
233// Assumes x > 0.
234// Finds the maximum step such that x + step*dx >= 0.
235
236     double step     = 1e+20;
237
238     int n = x.size();
239     double *x_elts = x.getElements();
240     double *dx_elts = dx.getElements();
241     for (int k = 0; k < n; k++)
242          if (dx_elts[k] < 0)
243               if ((x_elts[k] / (-dx_elts[k])) < step)
244                    step = x_elts[k] / (-dx_elts[k]);
245     return step;
246}
247//-----------------------------------------------------------------------
248// End private function pdxxxstep
249//-----------------------------------------------------------------------
250
251inline double  pdxxxstep(int nset, int *set, 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#endif
272#endif
Note: See TracBrowser for help on using the repository browser.