ClpHelperFunctions.hpp

Go to the documentation of this file.
00001 /* $Id: ClpHelperFunctions.hpp 2034 2014-06-12 16:10:26Z stefan $ */
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
00085 #   ifdef NDEBUG
00086 #      define ClpTraceDebug(expression)         {}
00087 #   else
00088 void ClpTracePrint(std::string fileName, std::string message, int line);
00089 #      define ClpTraceDebug(expression) { \
00090        if (!(expression)) { ClpTracePrint(__FILE__,__STRING(expression),__LINE__); } \
00091   }
00092 #   endif
00094 #ifdef ClpPdco_H
00095 
00096 
00097 inline double pdxxxmerit(int nlow, int nupp, int *low, int *upp, CoinDenseVector <double> &r1,
00098                          CoinDenseVector <double> &r2, CoinDenseVector <double> &rL,
00099                          CoinDenseVector <double> &rU, CoinDenseVector <double> &cL,
00100                          CoinDenseVector <double> &cU )
00101 {
00102 
00103 // Evaluate the merit function for Newton's method.
00104 // It is the 2-norm of the three sets of residuals.
00105      double sum1, sum2;
00106      CoinDenseVector <double> f(6);
00107      f[0] = r1.twoNorm();
00108      f[1] = r2.twoNorm();
00109      sum1 = sum2 = 0.0;
00110      for (int k = 0; k < nlow; k++) {
00111           sum1 += rL[low[k]] * rL[low[k]];
00112           sum2 += cL[low[k]] * cL[low[k]];
00113      }
00114      f[2] = sqrt(sum1);
00115      f[4] = sqrt(sum2);
00116      sum1 = sum2 = 0.0;
00117      for (int k = 0; k < nupp; k++) {
00118           sum1 += rL[upp[k]] * rL[upp[k]];
00119           sum2 += cL[upp[k]] * cL[upp[k]];
00120      }
00121      f[3] = sqrt(sum1);
00122      f[5] = sqrt(sum2);
00123 
00124      return f.twoNorm();
00125 }
00126 
00127 //-----------------------------------------------------------------------
00128 // End private function pdxxxmerit
00129 //-----------------------------------------------------------------------
00130 
00131 
00132 //function [r1,r2,rL,rU,Pinf,Dinf] =    ...
00133 //      pdxxxresid1( Aname,fix,low,upp, ...
00134 //                   b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 )
00135 
00136 inline void pdxxxresid1(ClpPdco *model, const int nlow, const int nupp, const int nfix,
00137                         int *low, int *upp, int *fix,
00138                         CoinDenseVector <double> &b, double *bl, double *bu, double d1, double d2,
00139                         CoinDenseVector <double> &grad, CoinDenseVector <double> &rL,
00140                         CoinDenseVector <double> &rU, CoinDenseVector <double> &x,
00141                         CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
00142                         CoinDenseVector <double> &y,  CoinDenseVector <double> &z1,
00143                         CoinDenseVector <double> &z2, CoinDenseVector <double> &r1,
00144                         CoinDenseVector <double> &r2, double *Pinf, double *Dinf)
00145 {
00146 
00147 // Form residuals for the primal and dual equations.
00148 // rL, rU are output, but we input them as full vectors
00149 // initialized (permanently) with any relevant zeros.
00150 
00151 // Get some element pointers for efficiency
00152      double *x_elts  = x.getElements();
00153      double *r2_elts = r2.getElements();
00154 
00155      for (int k = 0; k < nfix; k++)
00156           x_elts[fix[k]]  = 0;
00157 
00158      r1.clear();
00159      r2.clear();
00160      model->matVecMult( 1, r1, x );
00161      model->matVecMult( 2, r2, y );
00162      for (int k = 0; k < nfix; k++)
00163           r2_elts[fix[k]]  = 0;
00164 
00165 
00166      r1      = b    - r1 - d2 * d2 * y;
00167      r2      = grad - r2 - z1;              // grad includes d1*d1*x
00168      if (nupp > 0)
00169           r2    = r2 + z2;
00170 
00171      for (int k = 0; k < nlow; k++)
00172           rL[low[k]] = bl[low[k]] - x[low[k]] + x1[low[k]];
00173      for (int k = 0; k < nupp; k++)
00174           rU[upp[k]] = - bu[upp[k]] + x[upp[k]] + x2[upp[k]];
00175 
00176      double normL = 0.0;
00177      double normU = 0.0;
00178      for (int k = 0; k < nlow; k++)
00179           if (rL[low[k]] > normL) normL = rL[low[k]];
00180      for (int k = 0; k < nupp; k++)
00181           if (rU[upp[k]] > normU) normU = rU[upp[k]];
00182 
00183      *Pinf    = CoinMax(normL, normU);
00184      *Pinf    = CoinMax( r1.infNorm() , *Pinf );
00185      *Dinf    = r2.infNorm();
00186      *Pinf    = CoinMax( *Pinf, 1e-99 );
00187      *Dinf    = CoinMax( *Dinf, 1e-99 );
00188 }
00189 
00190 //-----------------------------------------------------------------------
00191 // End private function pdxxxresid1
00192 //-----------------------------------------------------------------------
00193 
00194 
00195 //function [cL,cU,center,Cinf,Cinf0] = ...
00196 //      pdxxxresid2( mu,low,upp,cL,cU,x1,x2,z1,z2 )
00197 
00198 inline void pdxxxresid2(double mu, int nlow, int nupp, int *low, int *upp,
00199                         CoinDenseVector <double> &cL, CoinDenseVector <double> &cU,
00200                         CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
00201                         CoinDenseVector <double> &z1, CoinDenseVector <double> &z2,
00202                         double *center, double *Cinf, double *Cinf0)
00203 {
00204 
00205 // Form residuals for the complementarity equations.
00206 // cL, cU are output, but we input them as full vectors
00207 // initialized (permanently) with any relevant zeros.
00208 // Cinf  is the complementarity residual for X1 z1 = mu e, etc.
00209 // Cinf0 is the same for mu=0 (i.e., for the original problem).
00210 
00211      double maxXz = -1e20;
00212      double minXz = 1e20;
00213 
00214      double *x1_elts = x1.getElements();
00215      double *z1_elts = z1.getElements();
00216      double *cL_elts = cL.getElements();
00217      for (int k = 0; k < nlow; k++) {
00218           double x1z1    = x1_elts[low[k]] * z1_elts[low[k]];
00219           cL_elts[low[k]] = mu - x1z1;
00220           if (x1z1 > maxXz) maxXz = x1z1;
00221           if (x1z1 < minXz) minXz = x1z1;
00222      }
00223 
00224      double *x2_elts = x2.getElements();
00225      double *z2_elts = z2.getElements();
00226      double *cU_elts = cU.getElements();
00227      for (int k = 0; k < nupp; k++) {
00228           double x2z2    = x2_elts[upp[k]] * z2_elts[upp[k]];
00229           cU_elts[upp[k]] = mu - x2z2;
00230           if (x2z2 > maxXz) maxXz = x2z2;
00231           if (x2z2 < minXz) minXz = x2z2;
00232      }
00233 
00234      maxXz   = CoinMax( maxXz, 1e-99 );
00235      minXz   = CoinMax( minXz, 1e-99 );
00236      *center  = maxXz / minXz;
00237 
00238      double normL = 0.0;
00239      double normU = 0.0;
00240      for (int k = 0; k < nlow; k++)
00241           if (cL_elts[low[k]] > normL) normL = cL_elts[low[k]];
00242      for (int k = 0; k < nupp; k++)
00243           if (cU_elts[upp[k]] > normU) normU = cU_elts[upp[k]];
00244      *Cinf    = CoinMax( normL, normU);
00245      *Cinf0   = maxXz;
00246 }
00247 //-----------------------------------------------------------------------
00248 // End private function pdxxxresid2
00249 //-----------------------------------------------------------------------
00250 
00251 inline double  pdxxxstep( 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 
00272 inline double  pdxxxstep(int nset, int *set, CoinDenseVector <double> &x, CoinDenseVector <double> &dx )
00273 {
00274 
00275 // Assumes x > 0.
00276 // Finds the maximum step such that x + step*dx >= 0.
00277 
00278      double step     = 1e+20;
00279 
00280      int n = x.size();
00281      double *x_elts = x.getElements();
00282      double *dx_elts = dx.getElements();
00283      for (int k = 0; k < n; k++)
00284           if (dx_elts[k] < 0)
00285                if ((x_elts[k] / (-dx_elts[k])) < step)
00286                     step = x_elts[k] / (-dx_elts[k]);
00287      return step;
00288 }
00289 //-----------------------------------------------------------------------
00290 // End private function pdxxxstep
00291 //-----------------------------------------------------------------------
00292 #endif
00293 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 10 Mar 2015 for Clp by  doxygen 1.6.1