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

Last change on this file since 1665 was 1665, checked in by lou, 10 years ago

Add EPL license notice in src.

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