Changeset 2385 for trunk/Clp/src/ClpHelperFunctions.hpp
 Timestamp:
 Jan 6, 2019 2:43:06 PM (4 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/ClpHelperFunctions.hpp
r2034 r2385 9 9 #include "ClpConfig.h" 10 10 #ifdef HAVE_CMATH 11 # 11 #include <cmath> 12 12 #else 13 # 14 # 15 # 16 # 17 # 13 #ifdef HAVE_MATH_H 14 #include <math.h> 15 #else 16 #error "don't have header file for math" 17 #endif 18 18 #endif 19 19 … … 26 26 */ 27 27 28 double maximumAbsElement(const double * 29 void setElements(double * 30 void multiplyAdd(const double * 31 double *region2, double multiplier2);32 double innerProduct(const double * region1, int size, const double *region2);33 void getNorms(const double * region, int size, double & norm1, double &norm2);28 double maximumAbsElement(const double *region, int size); 29 void setElements(double *region, int size, double value); 30 void multiplyAdd(const double *region1, int size, double multiplier1, 31 double *region2, double multiplier2); 32 double innerProduct(const double *region1, int size, const double *region2); 33 void getNorms(const double *region, int size, double &norm1, double &norm2); 34 34 #if COIN_LONG_WORK 35 35 // For long double versions 36 CoinWorkDouble maximumAbsElement(const CoinWorkDouble * 37 void setElements(CoinWorkDouble * 38 void multiplyAdd(const CoinWorkDouble * 39 CoinWorkDouble *region2, CoinWorkDouble multiplier2);40 CoinWorkDouble innerProduct(const CoinWorkDouble * region1, int size, const CoinWorkDouble *region2);41 void getNorms(const CoinWorkDouble * region, int size, CoinWorkDouble & norm1, CoinWorkDouble &norm2);36 CoinWorkDouble maximumAbsElement(const CoinWorkDouble *region, int size); 37 void setElements(CoinWorkDouble *region, int size, CoinWorkDouble value); 38 void multiplyAdd(const CoinWorkDouble *region1, int size, CoinWorkDouble multiplier1, 39 CoinWorkDouble *region2, CoinWorkDouble multiplier2); 40 CoinWorkDouble innerProduct(const CoinWorkDouble *region1, int size, const CoinWorkDouble *region2); 41 void getNorms(const CoinWorkDouble *region, int size, CoinWorkDouble &norm1, CoinWorkDouble &norm2); 42 42 inline void 43 CoinMemcpyN(const double * from, const int size, CoinWorkDouble *to)44 { 45 46 43 CoinMemcpyN(const double *from, const int size, CoinWorkDouble *to) 44 { 45 for (int i = 0; i < size; i++) 46 to[i] = from[i]; 47 47 } 48 48 inline void 49 CoinMemcpyN(const CoinWorkDouble * from, const int size, double *to)50 { 51 52 to[i] = static_cast<double>(from[i]);49 CoinMemcpyN(const CoinWorkDouble *from, const int size, double *to) 50 { 51 for (int i = 0; i < size; i++) 52 to[i] = static_cast< double >(from[i]); 53 53 } 54 54 inline CoinWorkDouble 55 55 CoinMax(const CoinWorkDouble x1, const double x2) 56 56 { 57 57 return (x1 > x2) ? x1 : x2; 58 58 } 59 59 inline CoinWorkDouble 60 60 CoinMax(double x1, const CoinWorkDouble x2) 61 61 { 62 62 return (x1 > x2) ? x1 : x2; 63 63 } 64 64 inline CoinWorkDouble 65 65 CoinMin(const CoinWorkDouble x1, const double x2) 66 66 { 67 67 return (x1 < x2) ? x1 : x2; 68 68 } 69 69 inline CoinWorkDouble 70 70 CoinMin(double x1, const CoinWorkDouble x2) 71 71 { 72 72 return (x1 < x2) ? x1 : x2; 73 73 } 74 74 inline CoinWorkDouble CoinSqrt(CoinWorkDouble x) 75 75 { 76 76 return sqrtl(x); 77 77 } 78 78 #else 79 79 inline double CoinSqrt(double x) 80 80 { 81 81 return sqrt(x); 82 82 } 83 83 #endif 84 84 /// Trace 85 # ifdef NDEBUG 86 # define ClpTraceDebug(expression) {} 87 # else 85 #ifdef NDEBUG 86 #define ClpTraceDebug(expression) \ 87 { \ 88 } 89 #else 88 90 void ClpTracePrint(std::string fileName, std::string message, int line); 89 # define ClpTraceDebug(expression) { \ 90 if (!(expression)) { ClpTracePrint(__FILE__,__STRING(expression),__LINE__); } \ 91 } 92 # endif 91 #define ClpTraceDebug(expression) \ 92 { \ 93 if (!(expression)) { \ 94 ClpTracePrint(__FILE__, __STRING(expression), __LINE__); \ 95 } \ 96 } 97 #endif 93 98 /// Following only included if ClpPdco defined 94 99 #ifdef ClpPdco_H 95 100 96 97 inline double pdxxxmerit(int nlow, int nupp, int *low, int *upp, CoinDenseVector <double> &r1, 98 CoinDenseVector <double> &r2, CoinDenseVector <double> &rL, 99 CoinDenseVector <double> &rU, CoinDenseVector <double> &cL, 100 CoinDenseVector <double> &cU ) 101 { 102 103 // Evaluate the merit function for Newton's method. 104 // It is the 2norm of the three sets of residuals. 105 double sum1, sum2; 106 CoinDenseVector <double> f(6); 107 f[0] = r1.twoNorm(); 108 f[1] = r2.twoNorm(); 109 sum1 = sum2 = 0.0; 110 for (int k = 0; k < nlow; k++) { 111 sum1 += rL[low[k]] * rL[low[k]]; 112 sum2 += cL[low[k]] * cL[low[k]]; 113 } 114 f[2] = sqrt(sum1); 115 f[4] = sqrt(sum2); 116 sum1 = sum2 = 0.0; 117 for (int k = 0; k < nupp; k++) { 118 sum1 += rL[upp[k]] * rL[upp[k]]; 119 sum2 += cL[upp[k]] * cL[upp[k]]; 120 } 121 f[3] = sqrt(sum1); 122 f[5] = sqrt(sum2); 123 124 return f.twoNorm(); 101 inline double pdxxxmerit(int nlow, int nupp, int *low, int *upp, CoinDenseVector< double > &r1, 102 CoinDenseVector< double > &r2, CoinDenseVector< double > &rL, 103 CoinDenseVector< double > &rU, CoinDenseVector< double > &cL, 104 CoinDenseVector< double > &cU) 105 { 106 107 // Evaluate the merit function for Newton's method. 108 // It is the 2norm of the three sets of residuals. 109 double sum1, sum2; 110 CoinDenseVector< double > f(6); 111 f[0] = r1.twoNorm(); 112 f[1] = r2.twoNorm(); 113 sum1 = sum2 = 0.0; 114 for (int k = 0; k < nlow; k++) { 115 sum1 += rL[low[k]] * rL[low[k]]; 116 sum2 += cL[low[k]] * cL[low[k]]; 117 } 118 f[2] = sqrt(sum1); 119 f[4] = sqrt(sum2); 120 sum1 = sum2 = 0.0; 121 for (int k = 0; k < nupp; k++) { 122 sum1 += rL[upp[k]] * rL[upp[k]]; 123 sum2 += cL[upp[k]] * cL[upp[k]]; 124 } 125 f[3] = sqrt(sum1); 126 f[5] = sqrt(sum2); 127 128 return f.twoNorm(); 125 129 } 126 130 … … 128 132 // End private function pdxxxmerit 129 133 // 130 131 134 132 135 //function [r1,r2,rL,rU,Pinf,Dinf] = ... … … 135 138 136 139 inline void pdxxxresid1(ClpPdco *model, const int nlow, const int nupp, const int nfix, 137 int *low, int *upp, int *fix, 138 CoinDenseVector <double> &b, double *bl, double *bu, double d1, double d2, 139 CoinDenseVector <double> &grad, CoinDenseVector <double> &rL, 140 CoinDenseVector <double> &rU, CoinDenseVector <double> &x, 141 CoinDenseVector <double> &x1, CoinDenseVector <double> &x2, 142 CoinDenseVector <double> &y, CoinDenseVector <double> &z1, 143 CoinDenseVector <double> &z2, CoinDenseVector <double> &r1, 144 CoinDenseVector <double> &r2, double *Pinf, double *Dinf) 145 { 146 147 // Form residuals for the primal and dual equations. 148 // rL, rU are output, but we input them as full vectors 149 // initialized (permanently) with any relevant zeros. 150 151 // Get some element pointers for efficiency 152 double *x_elts = x.getElements(); 153 double *r2_elts = r2.getElements(); 154 155 for (int k = 0; k < nfix; k++) 156 x_elts[fix[k]] = 0; 157 158 r1.clear(); 159 r2.clear(); 160 model>matVecMult( 1, r1, x ); 161 model>matVecMult( 2, r2, y ); 162 for (int k = 0; k < nfix; k++) 163 r2_elts[fix[k]] = 0; 164 165 166 r1 = b  r1  d2 * d2 * y; 167 r2 = grad  r2  z1; // grad includes d1*d1*x 168 if (nupp > 0) 169 r2 = r2 + z2; 170 171 for (int k = 0; k < nlow; k++) 172 rL[low[k]] = bl[low[k]]  x[low[k]] + x1[low[k]]; 173 for (int k = 0; k < nupp; k++) 174 rU[upp[k]] =  bu[upp[k]] + x[upp[k]] + x2[upp[k]]; 175 176 double normL = 0.0; 177 double normU = 0.0; 178 for (int k = 0; k < nlow; k++) 179 if (rL[low[k]] > normL) normL = rL[low[k]]; 180 for (int k = 0; k < nupp; k++) 181 if (rU[upp[k]] > normU) normU = rU[upp[k]]; 182 183 *Pinf = CoinMax(normL, normU); 184 *Pinf = CoinMax( r1.infNorm() , *Pinf ); 185 *Dinf = r2.infNorm(); 186 *Pinf = CoinMax( *Pinf, 1e99 ); 187 *Dinf = CoinMax( *Dinf, 1e99 ); 140 int *low, int *upp, int *fix, 141 CoinDenseVector< double > &b, double *bl, double *bu, double d1, double d2, 142 CoinDenseVector< double > &grad, CoinDenseVector< double > &rL, 143 CoinDenseVector< double > &rU, CoinDenseVector< double > &x, 144 CoinDenseVector< double > &x1, CoinDenseVector< double > &x2, 145 CoinDenseVector< double > &y, CoinDenseVector< double > &z1, 146 CoinDenseVector< double > &z2, CoinDenseVector< double > &r1, 147 CoinDenseVector< double > &r2, double *Pinf, double *Dinf) 148 { 149 150 // Form residuals for the primal and dual equations. 151 // rL, rU are output, but we input them as full vectors 152 // initialized (permanently) with any relevant zeros. 153 154 // Get some element pointers for efficiency 155 double *x_elts = x.getElements(); 156 double *r2_elts = r2.getElements(); 157 158 for (int k = 0; k < nfix; k++) 159 x_elts[fix[k]] = 0; 160 161 r1.clear(); 162 r2.clear(); 163 model>matVecMult(1, r1, x); 164 model>matVecMult(2, r2, y); 165 for (int k = 0; k < nfix; k++) 166 r2_elts[fix[k]] = 0; 167 168 r1 = b  r1  d2 * d2 * y; 169 r2 = grad  r2  z1; // grad includes d1*d1*x 170 if (nupp > 0) 171 r2 = r2 + z2; 172 173 for (int k = 0; k < nlow; k++) 174 rL[low[k]] = bl[low[k]]  x[low[k]] + x1[low[k]]; 175 for (int k = 0; k < nupp; k++) 176 rU[upp[k]] = bu[upp[k]] + x[upp[k]] + x2[upp[k]]; 177 178 double normL = 0.0; 179 double normU = 0.0; 180 for (int k = 0; k < nlow; k++) 181 if (rL[low[k]] > normL) 182 normL = rL[low[k]]; 183 for (int k = 0; k < nupp; k++) 184 if (rU[upp[k]] > normU) 185 normU = rU[upp[k]]; 186 187 *Pinf = CoinMax(normL, normU); 188 *Pinf = CoinMax(r1.infNorm(), *Pinf); 189 *Dinf = r2.infNorm(); 190 *Pinf = CoinMax(*Pinf, 1e99); 191 *Dinf = CoinMax(*Dinf, 1e99); 188 192 } 189 193 … … 191 195 // End private function pdxxxresid1 192 196 // 193 194 197 195 198 //function [cL,cU,center,Cinf,Cinf0] = ... … … 197 200 198 201 inline void pdxxxresid2(double mu, int nlow, int nupp, int *low, int *upp, 199 CoinDenseVector <double> &cL, CoinDenseVector <double> &cU, 200 CoinDenseVector <double> &x1, CoinDenseVector <double> &x2, 201 CoinDenseVector <double> &z1, CoinDenseVector <double> &z2, 202 double *center, double *Cinf, double *Cinf0) 203 { 204 205 // Form residuals for the complementarity equations. 206 // cL, cU are output, but we input them as full vectors 207 // initialized (permanently) with any relevant zeros. 208 // Cinf is the complementarity residual for X1 z1 = mu e, etc. 209 // Cinf0 is the same for mu=0 (i.e., for the original problem). 210 211 double maxXz = 1e20; 212 double minXz = 1e20; 213 214 double *x1_elts = x1.getElements(); 215 double *z1_elts = z1.getElements(); 216 double *cL_elts = cL.getElements(); 217 for (int k = 0; k < nlow; k++) { 218 double x1z1 = x1_elts[low[k]] * z1_elts[low[k]]; 219 cL_elts[low[k]] = mu  x1z1; 220 if (x1z1 > maxXz) maxXz = x1z1; 221 if (x1z1 < minXz) minXz = x1z1; 222 } 223 224 double *x2_elts = x2.getElements(); 225 double *z2_elts = z2.getElements(); 226 double *cU_elts = cU.getElements(); 227 for (int k = 0; k < nupp; k++) { 228 double x2z2 = x2_elts[upp[k]] * z2_elts[upp[k]]; 229 cU_elts[upp[k]] = mu  x2z2; 230 if (x2z2 > maxXz) maxXz = x2z2; 231 if (x2z2 < minXz) minXz = x2z2; 232 } 233 234 maxXz = CoinMax( maxXz, 1e99 ); 235 minXz = CoinMax( minXz, 1e99 ); 236 *center = maxXz / minXz; 237 238 double normL = 0.0; 239 double normU = 0.0; 240 for (int k = 0; k < nlow; k++) 241 if (cL_elts[low[k]] > normL) normL = cL_elts[low[k]]; 242 for (int k = 0; k < nupp; k++) 243 if (cU_elts[upp[k]] > normU) normU = cU_elts[upp[k]]; 244 *Cinf = CoinMax( normL, normU); 245 *Cinf0 = maxXz; 202 CoinDenseVector< double > &cL, CoinDenseVector< double > &cU, 203 CoinDenseVector< double > &x1, CoinDenseVector< double > &x2, 204 CoinDenseVector< double > &z1, CoinDenseVector< double > &z2, 205 double *center, double *Cinf, double *Cinf0) 206 { 207 208 // Form residuals for the complementarity equations. 209 // cL, cU are output, but we input them as full vectors 210 // initialized (permanently) with any relevant zeros. 211 // Cinf is the complementarity residual for X1 z1 = mu e, etc. 212 // Cinf0 is the same for mu=0 (i.e., for the original problem). 213 214 double maxXz = 1e20; 215 double minXz = 1e20; 216 217 double *x1_elts = x1.getElements(); 218 double *z1_elts = z1.getElements(); 219 double *cL_elts = cL.getElements(); 220 for (int k = 0; k < nlow; k++) { 221 double x1z1 = x1_elts[low[k]] * z1_elts[low[k]]; 222 cL_elts[low[k]] = mu  x1z1; 223 if (x1z1 > maxXz) 224 maxXz = x1z1; 225 if (x1z1 < minXz) 226 minXz = x1z1; 227 } 228 229 double *x2_elts = x2.getElements(); 230 double *z2_elts = z2.getElements(); 231 double *cU_elts = cU.getElements(); 232 for (int k = 0; k < nupp; k++) { 233 double x2z2 = x2_elts[upp[k]] * z2_elts[upp[k]]; 234 cU_elts[upp[k]] = mu  x2z2; 235 if (x2z2 > maxXz) 236 maxXz = x2z2; 237 if (x2z2 < minXz) 238 minXz = x2z2; 239 } 240 241 maxXz = CoinMax(maxXz, 1e99); 242 minXz = CoinMax(minXz, 1e99); 243 *center = maxXz / minXz; 244 245 double normL = 0.0; 246 double normU = 0.0; 247 for (int k = 0; k < nlow; k++) 248 if (cL_elts[low[k]] > normL) 249 normL = cL_elts[low[k]]; 250 for (int k = 0; k < nupp; k++) 251 if (cU_elts[upp[k]] > normU) 252 normU = cU_elts[upp[k]]; 253 *Cinf = CoinMax(normL, normU); 254 *Cinf0 = maxXz; 246 255 } 247 256 // … … 249 258 // 250 259 251 inline double pdxxxstep( 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 260 261 262 263 264 265 266 260 inline double pdxxxstep(CoinDenseVector< double > &x, CoinDenseVector< double > &dx) 261 { 262 263 // Assumes x > 0. 264 // Finds the maximum step such that x + step*dx >= 0. 265 266 double step = 1e+20; 267 268 int n = x.size(); 269 double *x_elts = x.getElements(); 270 double *dx_elts = dx.getElements(); 271 for (int k = 0; k < n; k++) 272 if (dx_elts[k] < 0) 273 if ((x_elts[k] / (dx_elts[k])) < step) 274 step = x_elts[k] / (dx_elts[k]); 275 return step; 267 276 } 268 277 // … … 270 279 // 271 280 272 inline double pdxxxstep(int nset, int *set, CoinDenseVector <double> &x, CoinDenseVector <double> &dx)273 { 274 275 // Assumes x > 0.276 // Finds the maximum step such that x + step*dx >= 0.277 278 double step= 1e+20;279 280 281 282 283 284 285 286 287 281 inline double pdxxxstep(int nset, int *set, CoinDenseVector< double > &x, CoinDenseVector< double > &dx) 282 { 283 284 // Assumes x > 0. 285 // Finds the maximum step such that x + step*dx >= 0. 286 287 double step = 1e+20; 288 289 int n = x.size(); 290 double *x_elts = x.getElements(); 291 double *dx_elts = dx.getElements(); 292 for (int k = 0; k < n; k++) 293 if (dx_elts[k] < 0) 294 if ((x_elts[k] / (dx_elts[k])) < step) 295 step = x_elts[k] / (dx_elts[k]); 296 return step; 288 297 } 289 298 // … … 292 301 #endif 293 302 #endif 303 304 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 305 */
Note: See TracChangeset
for help on using the changeset viewer.