- Main Page
- Namespaces
- Classes
- Files
- Directories

- Search for |

- File List
- File Members

# ClpHelperFunctions.hpp

Go to the documentation of this file.

    00001 /* $Id: ClpHelperFunctions.hpp 1753 2011-06-19 16:27: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
    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