/build/buildd/clp-1.12.0/Clp/src/ClpHelperFunctions.hpp
Go to the documentation of this file.
00001 /* $Id: ClpHelperFunctions.hpp 1525 2010-02-26 17:27:59Z mjs $ */
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 
00082 // Evaluate the merit function for Newton's method.
00083 // It is the 2-norm of the three sets of residuals.
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 // End private function pdxxxmerit
00108 //-----------------------------------------------------------------------
00109 
00110 
00111 //function [r1,r2,rL,rU,Pinf,Dinf] =    ...
00112 //      pdxxxresid1( Aname,fix,low,upp, ...
00113 //                   b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 )
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 // Form residuals for the primal and dual equations.
00127 // rL, rU are output, but we input them as full vectors
00128 // initialized (permanently) with any relevant zeros.
00129 
00130 // Get some element pointers for efficiency
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;              // grad includes d1*d1*x
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 // End private function pdxxxresid1
00171 //-----------------------------------------------------------------------
00172 
00173 
00174 //function [cL,cU,center,Cinf,Cinf0] = ...
00175 //      pdxxxresid2( mu,low,upp,cL,cU,x1,x2,z1,z2 )
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 // Form residuals for the complementarity equations.
00185 // cL, cU are output, but we input them as full vectors
00186 // initialized (permanently) with any relevant zeros.
00187 // Cinf  is the complementarity residual for X1 z1 = mu e, etc.
00188 // Cinf0 is the same for mu=0 (i.e., for the original problem).
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 // End private function pdxxxresid2
00228 //-----------------------------------------------------------------------
00229 
00230 inline double  pdxxxstep( CoinDenseVector <double> &x, CoinDenseVector <double> &dx )
00231 {
00232 
00233 // Assumes x > 0.
00234 // Finds the maximum step such that x + step*dx >= 0.
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 // End private function pdxxxstep
00249 //-----------------------------------------------------------------------
00250 
00251 inline double  pdxxxstep(int nset, int *set, CoinDenseVector <double> &x, CoinDenseVector <double> &dx )
00252 {
00253 
00254 // Assumes x > 0.
00255 // Finds the maximum step such that x + step*dx >= 0.
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 // End private function pdxxxstep
00270 //-----------------------------------------------------------------------
00271 #endif
00272 #endif