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

- File List
- File Members

# /home/coin/SVN-release/Bcps-0.91.0/Clp/src/ClpHelperFunctions.hpp

Go to the documentation of this file.

    00001 // Copyright (C) 2003, International Business Machines
    00002 // Corporation and others.  All Rights Reserved.
    00003 #ifndef ClpHelperFunctions_H
    00004 #define ClpHelperFunctions_H
    00005
    00014 double maximumAbsElement(const double * region, int size);
    00015 void setElements(double * region, int size, double value);
    00016 void multiplyAdd(const double * region1, int size, double multiplier1,
    00017                  double * region2, double multiplier2);
    00018 double innerProduct(const double * region1, int size, const double * region2);
    00019 void getNorms(const double * region, int size, double & norm1, double & norm2);
    00020
    00022 #ifdef ClpPdco_H
    00023
    00024
    00025 inline double pdxxxmerit(int nlow, int nupp, int *low, int *upp, CoinDenseVector <double> &r1,
    00026                 CoinDenseVector <double> &r2, CoinDenseVector <double> &rL,
    00027                 CoinDenseVector <double> &rU, CoinDenseVector <double> &cL,
    00028                 CoinDenseVector <double> &cU ){
    00029
    00030 // Evaluate the merit function for Newton's method.
    00031 // It is the 2-norm of the three sets of residuals.
    00032   double sum1, sum2;
    00033   CoinDenseVector <double> f(6);
    00034   f[0] = r1.twoNorm();
    00035   f[1] = r2.twoNorm();
    00036   sum1 = sum2 = 0.0;
    00037   for (int k=0; k<nlow; k++){
    00038     sum1 += rL[low[k]]*rL[low[k]];
    00039     sum2 += cL[low[k]]*cL[low[k]];
    00040   }
    00041   f[2] = sqrt(sum1);
    00042   f[4] = sqrt(sum2);
    00043   sum1 = sum2 = 0.0;
    00044   for (int k=0; k<nupp; k++){
    00045     sum1 += rL[upp[k]]*rL[upp[k]];
    00046     sum2 += cL[upp[k]]*cL[upp[k]];
    00047   }
    00048   f[3] = sqrt(sum1);
    00049   f[5] = sqrt(sum2);
    00050
    00051   return f.twoNorm();
    00052 }
    00053
    00054 //-----------------------------------------------------------------------
    00055 // End private function pdxxxmerit
    00056 //-----------------------------------------------------------------------
    00057
    00058
    00059 //function [r1,r2,rL,rU,Pinf,Dinf] =    ...
    00060 //      pdxxxresid1( Aname,fix,low,upp, ...
    00061 //                   b,bl,bu,d1,d2,grad,rL,rU,x,x1,x2,y,z1,z2 )
    00062
    00063 inline void pdxxxresid1(ClpPdco *model, const int nlow, const int nupp, const int nfix,
    00064                  int *low, int *upp, int *fix,
    00065                  CoinDenseVector <double> &b, double *bl, double *bu, double d1, double d2,
    00066                  CoinDenseVector <double> &grad, CoinDenseVector <double> &rL,
    00067                  CoinDenseVector <double> &rU, CoinDenseVector <double> &x,
    00068                  CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
    00069                  CoinDenseVector <double> &y,  CoinDenseVector <double> &z1,
    00070                  CoinDenseVector <double> &z2, CoinDenseVector <double> &r1,
    00071                  CoinDenseVector <double> &r2, double *Pinf, double *Dinf){
    00072
    00073 // Form residuals for the primal and dual equations.
    00074 // rL, rU are output, but we input them as full vectors
    00075 // initialized (permanently) with any relevant zeros.
    00076
    00077 // Get some element pointers for efficiency
    00078   double *x_elts  = x.getElements();
    00079   double *r2_elts = r2.getElements();
    00080
    00081   for (int k=0; k<nfix; k++)
    00082     x_elts[fix[k]]  = 0;
    00083
    00084   r1.clear();
    00085   r2.clear();
    00086   model->matVecMult( 1, r1, x );
    00087   model->matVecMult( 2, r2, y );
    00088   for (int k=0; k<nfix; k++)
    00089     r2_elts[fix[k]]  = 0;
    00090
    00091
    00092   r1      = b    - r1 - d2*d2*y;
    00093   r2      = grad - r2 - z1;              // grad includes d1*d1*x
    00094   if(nupp > 0)
    00095     r2    = r2 + z2;
    00096
    00097   for (int k=0; k<nlow; k++)
    00098     rL[low[k]] = bl[low[k]] - x[low[k]] + x1[low[k]];
    00099   for (int k=0; k<nupp; k++)
    00100     rU[upp[k]] = - bu[upp[k]] + x[upp[k]] + x2[upp[k]];
    00101
    00102   double normL = 0.0;
    00103   double normU = 0.0;
    00104   for (int k=0; k<nlow; k++)
    00105     if (rL[low[k]] > normL) normL = rL[low[k]];
    00106   for (int k=0; k<nupp; k++)
    00107     if (rU[upp[k]] > normU) normU = rU[upp[k]];
    00108
    00109   *Pinf    = CoinMax(normL, normU);
    00110   *Pinf    = CoinMax( r1.infNorm() , *Pinf );
    00111   *Dinf    = r2.infNorm();
    00112   *Pinf    = CoinMax( *Pinf, 1e-99 );
    00113   *Dinf    = CoinMax( *Dinf, 1e-99 );
    00114 }
    00115
    00116 //-----------------------------------------------------------------------
    00117 // End private function pdxxxresid1
    00118 //-----------------------------------------------------------------------
    00119
    00120
    00121 //function [cL,cU,center,Cinf,Cinf0] = ...
    00122 //      pdxxxresid2( mu,low,upp,cL,cU,x1,x2,z1,z2 )
    00123
    00124 inline void pdxxxresid2(double mu, int nlow, int nupp, int *low, int *upp,
    00125                  CoinDenseVector <double> &cL, CoinDenseVector <double> &cU,
    00126                  CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
    00127                  CoinDenseVector <double> &z1, CoinDenseVector <double> &z2,
    00128                  double *center, double *Cinf, double *Cinf0){
    00129
    00130 // Form residuals for the complementarity equations.
    00131 // cL, cU are output, but we input them as full vectors
    00132 // initialized (permanently) with any relevant zeros.
    00133 // Cinf  is the complementarity residual for X1 z1 = mu e, etc.
    00134 // Cinf0 is the same for mu=0 (i.e., for the original problem).
    00135
    00136   double maxXz = -1e20;
    00137   double minXz = 1e20;
    00138
    00139   double *x1_elts = x1.getElements();
    00140   double *z1_elts = z1.getElements();
    00141   double *cL_elts = cL.getElements();
    00142   for (int k=0; k<nlow; k++){
    00143     double x1z1    = x1_elts[low[k]] * z1_elts[low[k]];
    00144     cL_elts[low[k]] = mu - x1z1;
    00145     if(x1z1 > maxXz) maxXz = x1z1;
    00146     if(x1z1 < minXz) minXz = x1z1;
    00147   }
    00148
    00149   double *x2_elts = x2.getElements();
    00150   double *z2_elts = z2.getElements();
    00151   double *cU_elts = cU.getElements();
    00152   for (int k=0; k<nupp; k++){
    00153     double x2z2    = x2_elts[upp[k]] * z2_elts[upp[k]];
    00154     cU_elts[upp[k]] = mu - x2z2;
    00155     if(x2z2 > maxXz) maxXz = x2z2;
    00156     if(x2z2 < minXz) minXz = x2z2;
    00157   }
    00158
    00159   maxXz   = CoinMax( maxXz, 1e-99 );
    00160   minXz   = CoinMax( minXz, 1e-99 );
    00161   *center  = maxXz / minXz;
    00162
    00163   double normL = 0.0;
    00164   double normU = 0.0;
    00165   for (int k=0; k<nlow; k++)
    00166     if (cL_elts[low[k]] > normL) normL = cL_elts[low[k]];
    00167   for (int k=0; k<nupp; k++)
    00168     if (cU_elts[upp[k]] > normU) normU = cU_elts[upp[k]];
    00169   *Cinf    = CoinMax( normL, normU);
    00170   *Cinf0   = maxXz;
    00171 }
    00172 //-----------------------------------------------------------------------
    00173 // End private function pdxxxresid2
    00174 //-----------------------------------------------------------------------
    00175
    00176 inline double  pdxxxstep( CoinDenseVector <double> &x, CoinDenseVector <double> &dx ){
    00177
    00178 // Assumes x > 0.
    00179 // Finds the maximum step such that x + step*dx >= 0.
    00180
    00181   double step     = 1e+20;
    00182
    00183   int n = x.size();
    00184   double *x_elts = x.getElements();
    00185   double *dx_elts = dx.getElements();
    00186   for (int k=0; k<n; k++)
    00187     if (dx_elts[k] < 0)
    00188       if ((x_elts[k]/(-dx_elts[k])) < step)
    00189         step = x_elts[k]/(-dx_elts[k]);
    00190   return step;
    00191 }
    00192 //-----------------------------------------------------------------------
    00193 // End private function pdxxxstep
    00194 //-----------------------------------------------------------------------
    00195
    00196 inline double  pdxxxstep(int nset, int *set, CoinDenseVector <double> &x, CoinDenseVector <double> &dx ){
    00197
    00198 // Assumes x > 0.
    00199 // Finds the maximum step such that x + step*dx >= 0.
    00200
    00201   double step     = 1e+20;
    00202
    00203   int n = x.size();
    00204   double *x_elts = x.getElements();
    00205   double *dx_elts = dx.getElements();
    00206   for (int k=0; k<n; k++)
    00207     if (dx_elts[k] < 0)
    00208       if ((x_elts[k]/(-dx_elts[k])) < step)
    00209         step = x_elts[k]/(-dx_elts[k]);
    00210   return step;
    00211 }
    00212 //-----------------------------------------------------------------------
    00213 // End private function pdxxxstep
    00214 //-----------------------------------------------------------------------
    00215 #endif
    00216 #endif