Clp trunk
|
00001 /* $Id$ */ 00002 // Copyright (C) 2003, International Business Machines 00003 // Corporation and others. All Rights Reserved. 00004 // This code is licensed under the terms of the Eclipse Public License (EPL). 00005 00006 #ifndef ClpHelperFunctions_H 00007 #define ClpHelperFunctions_H 00008 00009 #include "ClpConfig.h" 00010 #ifdef HAVE_CMATH 00011 # include <cmath> 00012 #else 00013 # ifdef HAVE_MATH_H 00014 # include <math.h> 00015 # else 00016 # error "don't have header file for math" 00017 # endif 00018 #endif 00019 00028 double maximumAbsElement(const double * region, int size); 00029 void setElements(double * region, int size, double value); 00030 void multiplyAdd(const double * region1, int size, double multiplier1, 00031 double * region2, double multiplier2); 00032 double innerProduct(const double * region1, int size, const double * region2); 00033 void getNorms(const double * region, int size, double & norm1, double & norm2); 00034 #if COIN_LONG_WORK 00035 // For long double versions 00036 CoinWorkDouble maximumAbsElement(const CoinWorkDouble * region, int size); 00037 void setElements(CoinWorkDouble * region, int size, CoinWorkDouble value); 00038 void multiplyAdd(const CoinWorkDouble * region1, int size, CoinWorkDouble multiplier1, 00039 CoinWorkDouble * region2, CoinWorkDouble multiplier2); 00040 CoinWorkDouble innerProduct(const CoinWorkDouble * region1, int size, const CoinWorkDouble * region2); 00041 void getNorms(const CoinWorkDouble * region, int size, CoinWorkDouble & norm1, CoinWorkDouble & norm2); 00042 inline void 00043 CoinMemcpyN(const double * from, const int size, CoinWorkDouble * to) 00044 { 00045 for (int i = 0; i < size; i++) 00046 to[i] = from[i]; 00047 } 00048 inline void 00049 CoinMemcpyN(const CoinWorkDouble * from, const int size, double * to) 00050 { 00051 for (int i = 0; i < size; i++) 00052 to[i] = static_cast<double>(from[i]); 00053 } 00054 inline CoinWorkDouble 00055 CoinMax(const CoinWorkDouble x1, const double x2) 00056 { 00057 return (x1 > x2) ? x1 : x2; 00058 } 00059 inline CoinWorkDouble 00060 CoinMax(double x1, const CoinWorkDouble x2) 00061 { 00062 return (x1 > x2) ? x1 : x2; 00063 } 00064 inline CoinWorkDouble 00065 CoinMin(const CoinWorkDouble x1, const double x2) 00066 { 00067 return (x1 < x2) ? x1 : x2; 00068 } 00069 inline CoinWorkDouble 00070 CoinMin(double x1, const CoinWorkDouble x2) 00071 { 00072 return (x1 < x2) ? x1 : x2; 00073 } 00074 inline CoinWorkDouble CoinSqrt(CoinWorkDouble x) 00075 { 00076 return sqrtl(x); 00077 } 00078 #else 00079 inline double CoinSqrt(double x) 00080 { 00081 return sqrt(x); 00082 } 00083 #endif 00084 00086 #ifdef ClpPdco_H 00087 00088 00089 inline double pdxxxmerit(int nlow, int nupp, int *low, int *upp, CoinDenseVector <double> &r1, 00090 CoinDenseVector <double> &r2, CoinDenseVector <double> &rL, 00091 CoinDenseVector <double> &rU, CoinDenseVector <double> &cL, 00092 CoinDenseVector <double> &cU ) 00093 { 00094 00095 // Evaluate the merit function for Newton's method. 00096 // It is the 2-norm of the three sets of residuals. 00097 double sum1, sum2; 00098 CoinDenseVector <double> f(6); 00099 f[0] = r1.twoNorm(); 00100 f[1] = r2.twoNorm(); 00101 sum1 = sum2 = 0.0; 00102 for (int k = 0; k < nlow; k++) { 00103 sum1 += rL[low[k]] * rL[low[k]]; 00104 sum2 += cL[low[k]] * cL[low[k]]; 00105 } 00106 f[2] = sqrt(sum1); 00107 f[4] = sqrt(sum2); 00108 sum1 = sum2 = 0.0; 00109 for (int k = 0; k < nupp; k++) { 00110 sum1 += rL[upp[k]] * rL[upp[k]]; 00111 sum2 += cL[upp[k]] * cL[upp[k]]; 00112 } 00113 f[3] = sqrt(sum1); 00114 f[5] = sqrt(sum2); 00115 00116 return f.twoNorm(); 00117 } 00118 00119 //----------------------------------------------------------------------- 00120 // End private function pdxxxmerit 00121 //----------------------------------------------------------------------- 00122 00123 00124 //function [r1,r2,rL,rU,Pinf,Dinf] = ... 00125 // pdxxxresid1( Aname,fix,low,upp, ... 00126 // b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 ) 00127 00128 inline void pdxxxresid1(ClpPdco *model, const int nlow, const int nupp, const int nfix, 00129 int *low, int *upp, int *fix, 00130 CoinDenseVector <double> &b, double *bl, double *bu, double d1, double d2, 00131 CoinDenseVector <double> &grad, CoinDenseVector <double> &rL, 00132 CoinDenseVector <double> &rU, CoinDenseVector <double> &x, 00133 CoinDenseVector <double> &x1, CoinDenseVector <double> &x2, 00134 CoinDenseVector <double> &y, CoinDenseVector <double> &z1, 00135 CoinDenseVector <double> &z2, CoinDenseVector <double> &r1, 00136 CoinDenseVector <double> &r2, double *Pinf, double *Dinf) 00137 { 00138 00139 // Form residuals for the primal and dual equations. 00140 // rL, rU are output, but we input them as full vectors 00141 // initialized (permanently) with any relevant zeros. 00142 00143 // Get some element pointers for efficiency 00144 double *x_elts = x.getElements(); 00145 double *r2_elts = r2.getElements(); 00146 00147 for (int k = 0; k < nfix; k++) 00148 x_elts[fix[k]] = 0; 00149 00150 r1.clear(); 00151 r2.clear(); 00152 model->matVecMult( 1, r1, x ); 00153 model->matVecMult( 2, r2, y ); 00154 for (int k = 0; k < nfix; k++) 00155 r2_elts[fix[k]] = 0; 00156 00157 00158 r1 = b - r1 - d2 * d2 * y; 00159 r2 = grad - r2 - z1; // grad includes d1*d1*x 00160 if (nupp > 0) 00161 r2 = r2 + z2; 00162 00163 for (int k = 0; k < nlow; k++) 00164 rL[low[k]] = bl[low[k]] - x[low[k]] + x1[low[k]]; 00165 for (int k = 0; k < nupp; k++) 00166 rU[upp[k]] = - bu[upp[k]] + x[upp[k]] + x2[upp[k]]; 00167 00168 double normL = 0.0; 00169 double normU = 0.0; 00170 for (int k = 0; k < nlow; k++) 00171 if (rL[low[k]] > normL) normL = rL[low[k]]; 00172 for (int k = 0; k < nupp; k++) 00173 if (rU[upp[k]] > normU) normU = rU[upp[k]]; 00174 00175 *Pinf = CoinMax(normL, normU); 00176 *Pinf = CoinMax( r1.infNorm() , *Pinf ); 00177 *Dinf = r2.infNorm(); 00178 *Pinf = CoinMax( *Pinf, 1e-99 ); 00179 *Dinf = CoinMax( *Dinf, 1e-99 ); 00180 } 00181 00182 //----------------------------------------------------------------------- 00183 // End private function pdxxxresid1 00184 //----------------------------------------------------------------------- 00185 00186 00187 //function [cL,cU,center,Cinf,Cinf0] = ... 00188 // pdxxxresid2( mu,low,upp,cL,cU,x1,x2,z1,z2 ) 00189 00190 inline void pdxxxresid2(double mu, int nlow, int nupp, int *low, int *upp, 00191 CoinDenseVector <double> &cL, CoinDenseVector <double> &cU, 00192 CoinDenseVector <double> &x1, CoinDenseVector <double> &x2, 00193 CoinDenseVector <double> &z1, CoinDenseVector <double> &z2, 00194 double *center, double *Cinf, double *Cinf0) 00195 { 00196 00197 // Form residuals for the complementarity equations. 00198 // cL, cU are output, but we input them as full vectors 00199 // initialized (permanently) with any relevant zeros. 00200 // Cinf is the complementarity residual for X1 z1 = mu e, etc. 00201 // Cinf0 is the same for mu=0 (i.e., for the original problem). 00202 00203 double maxXz = -1e20; 00204 double minXz = 1e20; 00205 00206 double *x1_elts = x1.getElements(); 00207 double *z1_elts = z1.getElements(); 00208 double *cL_elts = cL.getElements(); 00209 for (int k = 0; k < nlow; k++) { 00210 double x1z1 = x1_elts[low[k]] * z1_elts[low[k]]; 00211 cL_elts[low[k]] = mu - x1z1; 00212 if (x1z1 > maxXz) maxXz = x1z1; 00213 if (x1z1 < minXz) minXz = x1z1; 00214 } 00215 00216 double *x2_elts = x2.getElements(); 00217 double *z2_elts = z2.getElements(); 00218 double *cU_elts = cU.getElements(); 00219 for (int k = 0; k < nupp; k++) { 00220 double x2z2 = x2_elts[upp[k]] * z2_elts[upp[k]]; 00221 cU_elts[upp[k]] = mu - x2z2; 00222 if (x2z2 > maxXz) maxXz = x2z2; 00223 if (x2z2 < minXz) minXz = x2z2; 00224 } 00225 00226 maxXz = CoinMax( maxXz, 1e-99 ); 00227 minXz = CoinMax( minXz, 1e-99 ); 00228 *center = maxXz / minXz; 00229 00230 double normL = 0.0; 00231 double normU = 0.0; 00232 for (int k = 0; k < nlow; k++) 00233 if (cL_elts[low[k]] > normL) normL = cL_elts[low[k]]; 00234 for (int k = 0; k < nupp; k++) 00235 if (cU_elts[upp[k]] > normU) normU = cU_elts[upp[k]]; 00236 *Cinf = CoinMax( normL, normU); 00237 *Cinf0 = maxXz; 00238 } 00239 //----------------------------------------------------------------------- 00240 // End private function pdxxxresid2 00241 //----------------------------------------------------------------------- 00242 00243 inline double pdxxxstep( CoinDenseVector <double> &x, CoinDenseVector <double> &dx ) 00244 { 00245 00246 // Assumes x > 0. 00247 // Finds the maximum step such that x + step*dx >= 0. 00248 00249 double step = 1e+20; 00250 00251 int n = x.size(); 00252 double *x_elts = x.getElements(); 00253 double *dx_elts = dx.getElements(); 00254 for (int k = 0; k < n; k++) 00255 if (dx_elts[k] < 0) 00256 if ((x_elts[k] / (-dx_elts[k])) < step) 00257 step = x_elts[k] / (-dx_elts[k]); 00258 return step; 00259 } 00260 //----------------------------------------------------------------------- 00261 // End private function pdxxxstep 00262 //----------------------------------------------------------------------- 00263 00264 inline double pdxxxstep(int nset, int *set, CoinDenseVector <double> &x, CoinDenseVector <double> &dx ) 00265 { 00266 00267 // Assumes x > 0. 00268 // Finds the maximum step such that x + step*dx >= 0. 00269 00270 double step = 1e+20; 00271 00272 int n = x.size(); 00273 double *x_elts = x.getElements(); 00274 double *dx_elts = dx.getElements(); 00275 for (int k = 0; k < n; k++) 00276 if (dx_elts[k] < 0) 00277 if ((x_elts[k] / (-dx_elts[k])) < step) 00278 step = x_elts[k] / (-dx_elts[k]); 00279 return step; 00280 } 00281 //----------------------------------------------------------------------- 00282 // End private function pdxxxstep 00283 //----------------------------------------------------------------------- 00284 #endif 00285 #endif