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

Last change on this file since 754 was 754, checked in by andreasw, 14 years ago

first version

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