- Main Page
- Modules
- Namespaces
- Classes
- Files
- Related Pages

- File List
- File Members

# /home/coin/SVN-release/Blis-0.91.2/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