00001
00002
00003
00004 #ifndef ClpHelperFunctions_H
00005 #define ClpHelperFunctions_H
00006
00015 double maximumAbsElement(const double * region, int size);
00016 void setElements(double * region, int size, double value);
00017 void multiplyAdd(const double * region1, int size, double multiplier1,
00018 double * region2, double multiplier2);
00019 double innerProduct(const double * region1, int size, const double * region2);
00020 void getNorms(const double * region, int size, double & norm1, double & norm2);
00021 #if COIN_LONG_WORK
00022
00023 CoinWorkDouble maximumAbsElement(const CoinWorkDouble * region, int size);
00024 void setElements(CoinWorkDouble * region, int size, CoinWorkDouble value);
00025 void multiplyAdd(const CoinWorkDouble * region1, int size, CoinWorkDouble multiplier1,
00026 CoinWorkDouble * region2, CoinWorkDouble multiplier2);
00027 CoinWorkDouble innerProduct(const CoinWorkDouble * region1, int size, const CoinWorkDouble * region2);
00028 void getNorms(const CoinWorkDouble * region, int size, CoinWorkDouble & norm1, CoinWorkDouble & norm2);
00029 inline void
00030 CoinMemcpyN(const double * from, const int size, CoinWorkDouble * to)
00031 {
00032 for (int i = 0; i < size; i++)
00033 to[i] = from[i];
00034 }
00035 inline void
00036 CoinMemcpyN(const CoinWorkDouble * from, const int size, double * to)
00037 {
00038 for (int i = 0; i < size; i++)
00039 to[i] = static_cast<double>(from[i]);
00040 }
00041 inline CoinWorkDouble
00042 CoinMax(const CoinWorkDouble x1, const double x2)
00043 {
00044 return (x1 > x2) ? x1 : x2;
00045 }
00046 inline CoinWorkDouble
00047 CoinMax(double x1, const CoinWorkDouble x2)
00048 {
00049 return (x1 > x2) ? x1 : x2;
00050 }
00051 inline CoinWorkDouble
00052 CoinMin(const CoinWorkDouble x1, const double x2)
00053 {
00054 return (x1 < x2) ? x1 : x2;
00055 }
00056 inline CoinWorkDouble
00057 CoinMin(double x1, const CoinWorkDouble x2)
00058 {
00059 return (x1 < x2) ? x1 : x2;
00060 }
00061 inline CoinWorkDouble CoinSqrt(CoinWorkDouble x)
00062 {
00063 return sqrtl(x);
00064 }
00065 #else
00066 inline double CoinSqrt(double x)
00067 {
00068 return sqrt(x);
00069 }
00070 #endif
00071
00073 #ifdef ClpPdco_H
00074
00075
00076 inline double pdxxxmerit(int nlow, int nupp, int *low, int *upp, CoinDenseVector <double> &r1,
00077 CoinDenseVector <double> &r2, CoinDenseVector <double> &rL,
00078 CoinDenseVector <double> &rU, CoinDenseVector <double> &cL,
00079 CoinDenseVector <double> &cU )
00080 {
00081
00082
00083
00084 double sum1, sum2;
00085 CoinDenseVector <double> f(6);
00086 f[0] = r1.twoNorm();
00087 f[1] = r2.twoNorm();
00088 sum1 = sum2 = 0.0;
00089 for (int k = 0; k < nlow; k++) {
00090 sum1 += rL[low[k]] * rL[low[k]];
00091 sum2 += cL[low[k]] * cL[low[k]];
00092 }
00093 f[2] = sqrt(sum1);
00094 f[4] = sqrt(sum2);
00095 sum1 = sum2 = 0.0;
00096 for (int k = 0; k < nupp; k++) {
00097 sum1 += rL[upp[k]] * rL[upp[k]];
00098 sum2 += cL[upp[k]] * cL[upp[k]];
00099 }
00100 f[3] = sqrt(sum1);
00101 f[5] = sqrt(sum2);
00102
00103 return f.twoNorm();
00104 }
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115 inline void pdxxxresid1(ClpPdco *model, const int nlow, const int nupp, const int nfix,
00116 int *low, int *upp, int *fix,
00117 CoinDenseVector <double> &b, double *bl, double *bu, double d1, double d2,
00118 CoinDenseVector <double> &grad, CoinDenseVector <double> &rL,
00119 CoinDenseVector <double> &rU, CoinDenseVector <double> &x,
00120 CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
00121 CoinDenseVector <double> &y, CoinDenseVector <double> &z1,
00122 CoinDenseVector <double> &z2, CoinDenseVector <double> &r1,
00123 CoinDenseVector <double> &r2, double *Pinf, double *Dinf)
00124 {
00125
00126
00127
00128
00129
00130
00131 double *x_elts = x.getElements();
00132 double *r2_elts = r2.getElements();
00133
00134 for (int k = 0; k < nfix; k++)
00135 x_elts[fix[k]] = 0;
00136
00137 r1.clear();
00138 r2.clear();
00139 model->matVecMult( 1, r1, x );
00140 model->matVecMult( 2, r2, y );
00141 for (int k = 0; k < nfix; k++)
00142 r2_elts[fix[k]] = 0;
00143
00144
00145 r1 = b - r1 - d2 * d2 * y;
00146 r2 = grad - r2 - z1;
00147 if (nupp > 0)
00148 r2 = r2 + z2;
00149
00150 for (int k = 0; k < nlow; k++)
00151 rL[low[k]] = bl[low[k]] - x[low[k]] + x1[low[k]];
00152 for (int k = 0; k < nupp; k++)
00153 rU[upp[k]] = - bu[upp[k]] + x[upp[k]] + x2[upp[k]];
00154
00155 double normL = 0.0;
00156 double normU = 0.0;
00157 for (int k = 0; k < nlow; k++)
00158 if (rL[low[k]] > normL) normL = rL[low[k]];
00159 for (int k = 0; k < nupp; k++)
00160 if (rU[upp[k]] > normU) normU = rU[upp[k]];
00161
00162 *Pinf = CoinMax(normL, normU);
00163 *Pinf = CoinMax( r1.infNorm() , *Pinf );
00164 *Dinf = r2.infNorm();
00165 *Pinf = CoinMax( *Pinf, 1e-99 );
00166 *Dinf = CoinMax( *Dinf, 1e-99 );
00167 }
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177 inline void pdxxxresid2(double mu, int nlow, int nupp, int *low, int *upp,
00178 CoinDenseVector <double> &cL, CoinDenseVector <double> &cU,
00179 CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
00180 CoinDenseVector <double> &z1, CoinDenseVector <double> &z2,
00181 double *center, double *Cinf, double *Cinf0)
00182 {
00183
00184
00185
00186
00187
00188
00189
00190 double maxXz = -1e20;
00191 double minXz = 1e20;
00192
00193 double *x1_elts = x1.getElements();
00194 double *z1_elts = z1.getElements();
00195 double *cL_elts = cL.getElements();
00196 for (int k = 0; k < nlow; k++) {
00197 double x1z1 = x1_elts[low[k]] * z1_elts[low[k]];
00198 cL_elts[low[k]] = mu - x1z1;
00199 if (x1z1 > maxXz) maxXz = x1z1;
00200 if (x1z1 < minXz) minXz = x1z1;
00201 }
00202
00203 double *x2_elts = x2.getElements();
00204 double *z2_elts = z2.getElements();
00205 double *cU_elts = cU.getElements();
00206 for (int k = 0; k < nupp; k++) {
00207 double x2z2 = x2_elts[upp[k]] * z2_elts[upp[k]];
00208 cU_elts[upp[k]] = mu - x2z2;
00209 if (x2z2 > maxXz) maxXz = x2z2;
00210 if (x2z2 < minXz) minXz = x2z2;
00211 }
00212
00213 maxXz = CoinMax( maxXz, 1e-99 );
00214 minXz = CoinMax( minXz, 1e-99 );
00215 *center = maxXz / minXz;
00216
00217 double normL = 0.0;
00218 double normU = 0.0;
00219 for (int k = 0; k < nlow; k++)
00220 if (cL_elts[low[k]] > normL) normL = cL_elts[low[k]];
00221 for (int k = 0; k < nupp; k++)
00222 if (cU_elts[upp[k]] > normU) normU = cU_elts[upp[k]];
00223 *Cinf = CoinMax( normL, normU);
00224 *Cinf0 = maxXz;
00225 }
00226
00227
00228
00229
00230 inline double pdxxxstep( CoinDenseVector <double> &x, CoinDenseVector <double> &dx )
00231 {
00232
00233
00234
00235
00236 double step = 1e+20;
00237
00238 int n = x.size();
00239 double *x_elts = x.getElements();
00240 double *dx_elts = dx.getElements();
00241 for (int k = 0; k < n; k++)
00242 if (dx_elts[k] < 0)
00243 if ((x_elts[k] / (-dx_elts[k])) < step)
00244 step = x_elts[k] / (-dx_elts[k]);
00245 return step;
00246 }
00247
00248
00249
00250
00251 inline double pdxxxstep(int nset, int *set, CoinDenseVector <double> &x, CoinDenseVector <double> &dx )
00252 {
00253
00254
00255
00256
00257 double step = 1e+20;
00258
00259 int n = x.size();
00260 double *x_elts = x.getElements();
00261 double *dx_elts = dx.getElements();
00262 for (int k = 0; k < n; k++)
00263 if (dx_elts[k] < 0)
00264 if ((x_elts[k] / (-dx_elts[k])) < step)
00265 step = x_elts[k] / (-dx_elts[k]);
00266 return step;
00267 }
00268
00269
00270
00271 #endif
00272 #endif