• Main Page
  • Classes
  • Files
  • File List
  • File Members

/build/buildd/clp-1.11.1/Clp/src/ClpHelperFunctions.hpp

Go to the documentation of this file.
00001 /* $Id: ClpHelperFunctions.hpp 1458 2009-11-05 12:34:07Z forrest $ */
00002 // Copyright (C) 2003, International Business Machines
00003 // Corporation and others.  All Rights Reserved.
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   // For long double versions 
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 // Evaluate the merit function for Newton's method.
00082 // It is the 2-norm of the three sets of residuals.
00083   double sum1, sum2;
00084   CoinDenseVector <double> f(6);
00085   f[0] = r1.twoNorm();
00086   f[1] = r2.twoNorm();
00087   sum1 = sum2 = 0.0;
00088   for (int k=0; k<nlow; k++){
00089     sum1 += rL[low[k]]*rL[low[k]];
00090     sum2 += cL[low[k]]*cL[low[k]];
00091   }
00092   f[2] = sqrt(sum1);
00093   f[4] = sqrt(sum2);
00094   sum1 = sum2 = 0.0;
00095   for (int k=0; k<nupp; k++){
00096     sum1 += rL[upp[k]]*rL[upp[k]];
00097     sum2 += cL[upp[k]]*cL[upp[k]];
00098   }
00099   f[3] = sqrt(sum1);
00100   f[5] = sqrt(sum2);
00101 
00102   return f.twoNorm();
00103 }
00104 
00105 //-----------------------------------------------------------------------
00106 // End private function pdxxxmerit
00107 //-----------------------------------------------------------------------
00108 
00109 
00110 //function [r1,r2,rL,rU,Pinf,Dinf] =    ...
00111 //      pdxxxresid1( Aname,fix,low,upp, ...
00112 //                   b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 )
00113 
00114 inline void pdxxxresid1(ClpPdco *model, const int nlow, const int nupp, const int nfix,
00115                  int *low, int *upp, int *fix,
00116                  CoinDenseVector <double> &b, double *bl, double *bu, double d1, double d2,
00117                  CoinDenseVector <double> &grad, CoinDenseVector <double> &rL,
00118                  CoinDenseVector <double> &rU, CoinDenseVector <double> &x,
00119                  CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
00120                  CoinDenseVector <double> &y,  CoinDenseVector <double> &z1,
00121                  CoinDenseVector <double> &z2, CoinDenseVector <double> &r1,
00122                  CoinDenseVector <double> &r2, double *Pinf, double *Dinf){
00123 
00124 // Form residuals for the primal and dual equations.
00125 // rL, rU are output, but we input them as full vectors
00126 // initialized (permanently) with any relevant zeros.
00127 
00128 // Get some element pointers for efficiency
00129   double *x_elts  = x.getElements();
00130   double *r2_elts = r2.getElements();
00131   
00132   for (int k=0; k<nfix; k++)
00133     x_elts[fix[k]]  = 0;
00134 
00135   r1.clear();
00136   r2.clear();   
00137   model->matVecMult( 1, r1, x );
00138   model->matVecMult( 2, r2, y );
00139   for (int k=0; k<nfix; k++)
00140     r2_elts[fix[k]]  = 0;
00141 
00142 
00143   r1      = b    - r1 - d2*d2*y;
00144   r2      = grad - r2 - z1;              // grad includes d1*d1*x
00145   if(nupp > 0)       
00146     r2    = r2 + z2;   
00147 
00148   for (int k=0; k<nlow; k++)
00149     rL[low[k]] = bl[low[k]] - x[low[k]] + x1[low[k]];
00150   for (int k=0; k<nupp; k++)
00151     rU[upp[k]] = - bu[upp[k]] + x[upp[k]] + x2[upp[k]];
00152 
00153   double normL = 0.0;
00154   double normU = 0.0;
00155   for (int k=0; k<nlow; k++)
00156     if (rL[low[k]] > normL) normL = rL[low[k]];
00157   for (int k=0; k<nupp; k++)
00158     if (rU[upp[k]] > normU) normU = rU[upp[k]];
00159 
00160   *Pinf    = CoinMax(normL, normU);  
00161   *Pinf    = CoinMax( r1.infNorm() , *Pinf );
00162   *Dinf    = r2.infNorm();
00163   *Pinf    = CoinMax( *Pinf, 1e-99 );
00164   *Dinf    = CoinMax( *Dinf, 1e-99 );
00165 }
00166 
00167 //-----------------------------------------------------------------------
00168 // End private function pdxxxresid1
00169 //-----------------------------------------------------------------------
00170 
00171 
00172 //function [cL,cU,center,Cinf,Cinf0] = ...
00173 //      pdxxxresid2( mu,low,upp,cL,cU,x1,x2,z1,z2 )
00174 
00175 inline void pdxxxresid2(double mu, int nlow, int nupp, int *low, int *upp,
00176                  CoinDenseVector <double> &cL, CoinDenseVector <double> &cU,
00177                  CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
00178                  CoinDenseVector <double> &z1, CoinDenseVector <double> &z2,
00179                  double *center, double *Cinf, double *Cinf0){
00180 
00181 // Form residuals for the complementarity equations.
00182 // cL, cU are output, but we input them as full vectors
00183 // initialized (permanently) with any relevant zeros.
00184 // Cinf  is the complementarity residual for X1 z1 = mu e, etc.
00185 // Cinf0 is the same for mu=0 (i.e., for the original problem).
00186 
00187   double maxXz = -1e20;
00188   double minXz = 1e20;
00189 
00190   double *x1_elts = x1.getElements();
00191   double *z1_elts = z1.getElements();
00192   double *cL_elts = cL.getElements();
00193   for (int k=0; k<nlow; k++){
00194     double x1z1    = x1_elts[low[k]] * z1_elts[low[k]];
00195     cL_elts[low[k]] = mu - x1z1;
00196     if(x1z1 > maxXz) maxXz = x1z1;
00197     if(x1z1 < minXz) minXz = x1z1;
00198   }
00199 
00200   double *x2_elts = x2.getElements();
00201   double *z2_elts = z2.getElements();
00202   double *cU_elts = cU.getElements();
00203   for (int k=0; k<nupp; k++){
00204     double x2z2    = x2_elts[upp[k]] * z2_elts[upp[k]];
00205     cU_elts[upp[k]] = mu - x2z2;
00206     if(x2z2 > maxXz) maxXz = x2z2;
00207     if(x2z2 < minXz) minXz = x2z2;
00208   }
00209 
00210   maxXz   = CoinMax( maxXz, 1e-99 );
00211   minXz   = CoinMax( minXz, 1e-99 );
00212   *center  = maxXz / minXz;
00213 
00214   double normL = 0.0;
00215   double normU = 0.0;
00216   for (int k=0; k<nlow; k++)
00217     if (cL_elts[low[k]] > normL) normL = cL_elts[low[k]];
00218   for (int k=0; k<nupp; k++)
00219     if (cU_elts[upp[k]] > normU) normU = cU_elts[upp[k]];
00220   *Cinf    = CoinMax( normL, normU);
00221   *Cinf0   = maxXz;
00222 }
00223 //-----------------------------------------------------------------------
00224 // End private function pdxxxresid2
00225 //-----------------------------------------------------------------------
00226 
00227 inline double  pdxxxstep( CoinDenseVector <double> &x, CoinDenseVector <double> &dx ){
00228 
00229 // Assumes x > 0.
00230 // Finds the maximum step such that x + step*dx >= 0.
00231 
00232   double step     = 1e+20;
00233 
00234   int n = x.size();
00235   double *x_elts = x.getElements();
00236   double *dx_elts = dx.getElements();
00237   for (int k=0; k<n; k++)
00238     if (dx_elts[k] < 0)
00239       if ((x_elts[k]/(-dx_elts[k])) < step)
00240         step = x_elts[k]/(-dx_elts[k]);
00241   return step;
00242 }
00243 //-----------------------------------------------------------------------
00244 // End private function pdxxxstep
00245 //-----------------------------------------------------------------------
00246 
00247 inline double  pdxxxstep(int nset, int *set, CoinDenseVector <double> &x, CoinDenseVector <double> &dx ){
00248 
00249 // Assumes x > 0.
00250 // Finds the maximum step such that x + step*dx >= 0.
00251 
00252   double step     = 1e+20;
00253 
00254   int n = x.size();
00255   double *x_elts = x.getElements();
00256   double *dx_elts = dx.getElements();
00257   for (int k=0; k<n; k++)
00258     if (dx_elts[k] < 0)
00259       if ((x_elts[k]/(-dx_elts[k])) < step)
00260         step = x_elts[k]/(-dx_elts[k]);
00261   return step;
00262 }
00263 //-----------------------------------------------------------------------
00264 // End private function pdxxxstep
00265 //-----------------------------------------------------------------------
00266 #endif
00267 #endif

Generated on Fri Aug 20 2010 06:29:28 by  doxygen 1.7.1