Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (4 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpHelperFunctions.hpp

    r2034 r2385  
    99#include "ClpConfig.h"
    1010#ifdef HAVE_CMATH
    11 # include <cmath>
     11#include <cmath>
    1212#else
    13 # ifdef HAVE_MATH_H
    14 #  include <math.h>
    15 # else
    16 #  error "don't have header file for math"
    17 # endif
     13#ifdef HAVE_MATH_H
     14#include <math.h>
     15#else
     16#error "don't have header file for math"
     17#endif
    1818#endif
    1919
     
    2626*/
    2727
    28 double maximumAbsElement(const double * region, int size);
    29 void setElements(double * region, int size, double value);
    30 void multiplyAdd(const double * region1, int size, double multiplier1,
    31                  double * region2, double multiplier2);
    32 double innerProduct(const double * region1, int size, const double * region2);
    33 void getNorms(const double * region, int size, double & norm1, double & norm2);
     28double maximumAbsElement(const double *region, int size);
     29void setElements(double *region, int size, double value);
     30void multiplyAdd(const double *region1, int size, double multiplier1,
     31  double *region2, double multiplier2);
     32double innerProduct(const double *region1, int size, const double *region2);
     33void getNorms(const double *region, int size, double &norm1, double &norm2);
    3434#if COIN_LONG_WORK
    3535// For long double versions
    36 CoinWorkDouble maximumAbsElement(const CoinWorkDouble * region, int size);
    37 void setElements(CoinWorkDouble * region, int size, CoinWorkDouble value);
    38 void multiplyAdd(const CoinWorkDouble * region1, int size, CoinWorkDouble multiplier1,
    39                  CoinWorkDouble * region2, CoinWorkDouble multiplier2);
    40 CoinWorkDouble innerProduct(const CoinWorkDouble * region1, int size, const CoinWorkDouble * region2);
    41 void getNorms(const CoinWorkDouble * region, int size, CoinWorkDouble & norm1, CoinWorkDouble & norm2);
     36CoinWorkDouble maximumAbsElement(const CoinWorkDouble *region, int size);
     37void setElements(CoinWorkDouble *region, int size, CoinWorkDouble value);
     38void multiplyAdd(const CoinWorkDouble *region1, int size, CoinWorkDouble multiplier1,
     39  CoinWorkDouble *region2, CoinWorkDouble multiplier2);
     40CoinWorkDouble innerProduct(const CoinWorkDouble *region1, int size, const CoinWorkDouble *region2);
     41void getNorms(const CoinWorkDouble *region, int size, CoinWorkDouble &norm1, CoinWorkDouble &norm2);
    4242inline void
    43 CoinMemcpyN(const double * from, const int size, CoinWorkDouble * to)
    44 {
    45      for (int i = 0; i < size; i++)
    46           to[i] = from[i];
     43CoinMemcpyN(const double *from, const int size, CoinWorkDouble *to)
     44{
     45  for (int i = 0; i < size; i++)
     46    to[i] = from[i];
    4747}
    4848inline void
    49 CoinMemcpyN(const CoinWorkDouble * from, const int size, double * to)
    50 {
    51      for (int i = 0; i < size; i++)
    52           to[i] = static_cast<double>(from[i]);
     49CoinMemcpyN(const CoinWorkDouble *from, const int size, double *to)
     50{
     51  for (int i = 0; i < size; i++)
     52    to[i] = static_cast< double >(from[i]);
    5353}
    5454inline CoinWorkDouble
    5555CoinMax(const CoinWorkDouble x1, const double x2)
    5656{
    57      return (x1 > x2) ? x1 : x2;
     57  return (x1 > x2) ? x1 : x2;
    5858}
    5959inline CoinWorkDouble
    6060CoinMax(double x1, const CoinWorkDouble x2)
    6161{
    62      return (x1 > x2) ? x1 : x2;
     62  return (x1 > x2) ? x1 : x2;
    6363}
    6464inline CoinWorkDouble
    6565CoinMin(const CoinWorkDouble x1, const double x2)
    6666{
    67      return (x1 < x2) ? x1 : x2;
     67  return (x1 < x2) ? x1 : x2;
    6868}
    6969inline CoinWorkDouble
    7070CoinMin(double x1, const CoinWorkDouble x2)
    7171{
    72      return (x1 < x2) ? x1 : x2;
     72  return (x1 < x2) ? x1 : x2;
    7373}
    7474inline CoinWorkDouble CoinSqrt(CoinWorkDouble x)
    7575{
    76      return sqrtl(x);
     76  return sqrtl(x);
    7777}
    7878#else
    7979inline double CoinSqrt(double x)
    8080{
    81      return sqrt(x);
     81  return sqrt(x);
    8282}
    8383#endif
    8484/// Trace
    85 #   ifdef NDEBUG
    86 #      define ClpTraceDebug(expression)         {}
    87 #   else
     85#ifdef NDEBUG
     86#define ClpTraceDebug(expression) \
     87  {                               \
     88  }
     89#else
    8890void ClpTracePrint(std::string fileName, std::string message, int line);
    89 #      define ClpTraceDebug(expression) { \
    90        if (!(expression)) { ClpTracePrint(__FILE__,__STRING(expression),__LINE__); } \
    91   }
    92 #   endif
     91#define ClpTraceDebug(expression)                              \
     92  {                                                            \
     93    if (!(expression)) {                                       \
     94      ClpTracePrint(__FILE__, __STRING(expression), __LINE__); \
     95    }                                                          \
     96  }
     97#endif
    9398/// Following only included if ClpPdco defined
    9499#ifdef ClpPdco_H
    95100
    96 
    97 inline double pdxxxmerit(int nlow, int nupp, int *low, int *upp, CoinDenseVector <double> &r1,
    98                          CoinDenseVector <double> &r2, CoinDenseVector <double> &rL,
    99                          CoinDenseVector <double> &rU, CoinDenseVector <double> &cL,
    100                          CoinDenseVector <double> &cU )
    101 {
    102 
    103 // Evaluate the merit function for Newton's method.
    104 // It is the 2-norm of the three sets of residuals.
    105      double sum1, sum2;
    106      CoinDenseVector <double> f(6);
    107      f[0] = r1.twoNorm();
    108      f[1] = r2.twoNorm();
    109      sum1 = sum2 = 0.0;
    110      for (int k = 0; k < nlow; k++) {
    111           sum1 += rL[low[k]] * rL[low[k]];
    112           sum2 += cL[low[k]] * cL[low[k]];
    113      }
    114      f[2] = sqrt(sum1);
    115      f[4] = sqrt(sum2);
    116      sum1 = sum2 = 0.0;
    117      for (int k = 0; k < nupp; k++) {
    118           sum1 += rL[upp[k]] * rL[upp[k]];
    119           sum2 += cL[upp[k]] * cL[upp[k]];
    120      }
    121      f[3] = sqrt(sum1);
    122      f[5] = sqrt(sum2);
    123 
    124      return f.twoNorm();
     101inline double pdxxxmerit(int nlow, int nupp, int *low, int *upp, CoinDenseVector< double > &r1,
     102  CoinDenseVector< double > &r2, CoinDenseVector< double > &rL,
     103  CoinDenseVector< double > &rU, CoinDenseVector< double > &cL,
     104  CoinDenseVector< double > &cU)
     105{
     106
     107  // Evaluate the merit function for Newton's method.
     108  // It is the 2-norm of the three sets of residuals.
     109  double sum1, sum2;
     110  CoinDenseVector< double > f(6);
     111  f[0] = r1.twoNorm();
     112  f[1] = r2.twoNorm();
     113  sum1 = sum2 = 0.0;
     114  for (int k = 0; k < nlow; k++) {
     115    sum1 += rL[low[k]] * rL[low[k]];
     116    sum2 += cL[low[k]] * cL[low[k]];
     117  }
     118  f[2] = sqrt(sum1);
     119  f[4] = sqrt(sum2);
     120  sum1 = sum2 = 0.0;
     121  for (int k = 0; k < nupp; k++) {
     122    sum1 += rL[upp[k]] * rL[upp[k]];
     123    sum2 += cL[upp[k]] * cL[upp[k]];
     124  }
     125  f[3] = sqrt(sum1);
     126  f[5] = sqrt(sum2);
     127
     128  return f.twoNorm();
    125129}
    126130
     
    128132// End private function pdxxxmerit
    129133//-----------------------------------------------------------------------
    130 
    131134
    132135//function [r1,r2,rL,rU,Pinf,Dinf] =    ...
     
    135138
    136139inline void pdxxxresid1(ClpPdco *model, const int nlow, const int nupp, const int nfix,
    137                         int *low, int *upp, int *fix,
    138                         CoinDenseVector <double> &b, double *bl, double *bu, double d1, double d2,
    139                         CoinDenseVector <double> &grad, CoinDenseVector <double> &rL,
    140                         CoinDenseVector <double> &rU, CoinDenseVector <double> &x,
    141                         CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
    142                         CoinDenseVector <double> &y,  CoinDenseVector <double> &z1,
    143                         CoinDenseVector <double> &z2, CoinDenseVector <double> &r1,
    144                         CoinDenseVector <double> &r2, double *Pinf, double *Dinf)
    145 {
    146 
    147 // Form residuals for the primal and dual equations.
    148 // rL, rU are output, but we input them as full vectors
    149 // initialized (permanently) with any relevant zeros.
    150 
    151 // Get some element pointers for efficiency
    152      double *x_elts  = x.getElements();
    153      double *r2_elts = r2.getElements();
    154 
    155      for (int k = 0; k < nfix; k++)
    156           x_elts[fix[k]]  = 0;
    157 
    158      r1.clear();
    159      r2.clear();
    160      model->matVecMult( 1, r1, x );
    161      model->matVecMult( 2, r2, y );
    162      for (int k = 0; k < nfix; k++)
    163           r2_elts[fix[k]]  = 0;
    164 
    165 
    166      r1      = b    - r1 - d2 * d2 * y;
    167      r2      = grad - r2 - z1;              // grad includes d1*d1*x
    168      if (nupp > 0)
    169           r2    = r2 + z2;
    170 
    171      for (int k = 0; k < nlow; k++)
    172           rL[low[k]] = bl[low[k]] - x[low[k]] + x1[low[k]];
    173      for (int k = 0; k < nupp; k++)
    174           rU[upp[k]] = - bu[upp[k]] + x[upp[k]] + x2[upp[k]];
    175 
    176      double normL = 0.0;
    177      double normU = 0.0;
    178      for (int k = 0; k < nlow; k++)
    179           if (rL[low[k]] > normL) normL = rL[low[k]];
    180      for (int k = 0; k < nupp; k++)
    181           if (rU[upp[k]] > normU) normU = rU[upp[k]];
    182 
    183      *Pinf    = CoinMax(normL, normU);
    184      *Pinf    = CoinMax( r1.infNorm() , *Pinf );
    185      *Dinf    = r2.infNorm();
    186      *Pinf    = CoinMax( *Pinf, 1e-99 );
    187      *Dinf    = CoinMax( *Dinf, 1e-99 );
     140  int *low, int *upp, int *fix,
     141  CoinDenseVector< double > &b, double *bl, double *bu, double d1, double d2,
     142  CoinDenseVector< double > &grad, CoinDenseVector< double > &rL,
     143  CoinDenseVector< double > &rU, CoinDenseVector< double > &x,
     144  CoinDenseVector< double > &x1, CoinDenseVector< double > &x2,
     145  CoinDenseVector< double > &y, CoinDenseVector< double > &z1,
     146  CoinDenseVector< double > &z2, CoinDenseVector< double > &r1,
     147  CoinDenseVector< double > &r2, double *Pinf, double *Dinf)
     148{
     149
     150  // Form residuals for the primal and dual equations.
     151  // rL, rU are output, but we input them as full vectors
     152  // initialized (permanently) with any relevant zeros.
     153
     154  // Get some element pointers for efficiency
     155  double *x_elts = x.getElements();
     156  double *r2_elts = r2.getElements();
     157
     158  for (int k = 0; k < nfix; k++)
     159    x_elts[fix[k]] = 0;
     160
     161  r1.clear();
     162  r2.clear();
     163  model->matVecMult(1, r1, x);
     164  model->matVecMult(2, r2, y);
     165  for (int k = 0; k < nfix; k++)
     166    r2_elts[fix[k]] = 0;
     167
     168  r1 = b - r1 - d2 * d2 * y;
     169  r2 = grad - r2 - z1; // grad includes d1*d1*x
     170  if (nupp > 0)
     171    r2 = r2 + z2;
     172
     173  for (int k = 0; k < nlow; k++)
     174    rL[low[k]] = bl[low[k]] - x[low[k]] + x1[low[k]];
     175  for (int k = 0; k < nupp; k++)
     176    rU[upp[k]] = -bu[upp[k]] + x[upp[k]] + x2[upp[k]];
     177
     178  double normL = 0.0;
     179  double normU = 0.0;
     180  for (int k = 0; k < nlow; k++)
     181    if (rL[low[k]] > normL)
     182      normL = rL[low[k]];
     183  for (int k = 0; k < nupp; k++)
     184    if (rU[upp[k]] > normU)
     185      normU = rU[upp[k]];
     186
     187  *Pinf = CoinMax(normL, normU);
     188  *Pinf = CoinMax(r1.infNorm(), *Pinf);
     189  *Dinf = r2.infNorm();
     190  *Pinf = CoinMax(*Pinf, 1e-99);
     191  *Dinf = CoinMax(*Dinf, 1e-99);
    188192}
    189193
     
    191195// End private function pdxxxresid1
    192196//-----------------------------------------------------------------------
    193 
    194197
    195198//function [cL,cU,center,Cinf,Cinf0] = ...
     
    197200
    198201inline void pdxxxresid2(double mu, int nlow, int nupp, int *low, int *upp,
    199                         CoinDenseVector <double> &cL, CoinDenseVector <double> &cU,
    200                         CoinDenseVector <double> &x1, CoinDenseVector <double> &x2,
    201                         CoinDenseVector <double> &z1, CoinDenseVector <double> &z2,
    202                         double *center, double *Cinf, double *Cinf0)
    203 {
    204 
    205 // Form residuals for the complementarity equations.
    206 // cL, cU are output, but we input them as full vectors
    207 // initialized (permanently) with any relevant zeros.
    208 // Cinf  is the complementarity residual for X1 z1 = mu e, etc.
    209 // Cinf0 is the same for mu=0 (i.e., for the original problem).
    210 
    211      double maxXz = -1e20;
    212      double minXz = 1e20;
    213 
    214      double *x1_elts = x1.getElements();
    215      double *z1_elts = z1.getElements();
    216      double *cL_elts = cL.getElements();
    217      for (int k = 0; k < nlow; k++) {
    218           double x1z1    = x1_elts[low[k]] * z1_elts[low[k]];
    219           cL_elts[low[k]] = mu - x1z1;
    220           if (x1z1 > maxXz) maxXz = x1z1;
    221           if (x1z1 < minXz) minXz = x1z1;
    222      }
    223 
    224      double *x2_elts = x2.getElements();
    225      double *z2_elts = z2.getElements();
    226      double *cU_elts = cU.getElements();
    227      for (int k = 0; k < nupp; k++) {
    228           double x2z2    = x2_elts[upp[k]] * z2_elts[upp[k]];
    229           cU_elts[upp[k]] = mu - x2z2;
    230           if (x2z2 > maxXz) maxXz = x2z2;
    231           if (x2z2 < minXz) minXz = x2z2;
    232      }
    233 
    234      maxXz   = CoinMax( maxXz, 1e-99 );
    235      minXz   = CoinMax( minXz, 1e-99 );
    236      *center  = maxXz / minXz;
    237 
    238      double normL = 0.0;
    239      double normU = 0.0;
    240      for (int k = 0; k < nlow; k++)
    241           if (cL_elts[low[k]] > normL) normL = cL_elts[low[k]];
    242      for (int k = 0; k < nupp; k++)
    243           if (cU_elts[upp[k]] > normU) normU = cU_elts[upp[k]];
    244      *Cinf    = CoinMax( normL, normU);
    245      *Cinf0   = maxXz;
     202  CoinDenseVector< double > &cL, CoinDenseVector< double > &cU,
     203  CoinDenseVector< double > &x1, CoinDenseVector< double > &x2,
     204  CoinDenseVector< double > &z1, CoinDenseVector< double > &z2,
     205  double *center, double *Cinf, double *Cinf0)
     206{
     207
     208  // Form residuals for the complementarity equations.
     209  // cL, cU are output, but we input them as full vectors
     210  // initialized (permanently) with any relevant zeros.
     211  // Cinf  is the complementarity residual for X1 z1 = mu e, etc.
     212  // Cinf0 is the same for mu=0 (i.e., for the original problem).
     213
     214  double maxXz = -1e20;
     215  double minXz = 1e20;
     216
     217  double *x1_elts = x1.getElements();
     218  double *z1_elts = z1.getElements();
     219  double *cL_elts = cL.getElements();
     220  for (int k = 0; k < nlow; k++) {
     221    double x1z1 = x1_elts[low[k]] * z1_elts[low[k]];
     222    cL_elts[low[k]] = mu - x1z1;
     223    if (x1z1 > maxXz)
     224      maxXz = x1z1;
     225    if (x1z1 < minXz)
     226      minXz = x1z1;
     227  }
     228
     229  double *x2_elts = x2.getElements();
     230  double *z2_elts = z2.getElements();
     231  double *cU_elts = cU.getElements();
     232  for (int k = 0; k < nupp; k++) {
     233    double x2z2 = x2_elts[upp[k]] * z2_elts[upp[k]];
     234    cU_elts[upp[k]] = mu - x2z2;
     235    if (x2z2 > maxXz)
     236      maxXz = x2z2;
     237    if (x2z2 < minXz)
     238      minXz = x2z2;
     239  }
     240
     241  maxXz = CoinMax(maxXz, 1e-99);
     242  minXz = CoinMax(minXz, 1e-99);
     243  *center = maxXz / minXz;
     244
     245  double normL = 0.0;
     246  double normU = 0.0;
     247  for (int k = 0; k < nlow; k++)
     248    if (cL_elts[low[k]] > normL)
     249      normL = cL_elts[low[k]];
     250  for (int k = 0; k < nupp; k++)
     251    if (cU_elts[upp[k]] > normU)
     252      normU = cU_elts[upp[k]];
     253  *Cinf = CoinMax(normL, normU);
     254  *Cinf0 = maxXz;
    246255}
    247256//-----------------------------------------------------------------------
     
    249258//-----------------------------------------------------------------------
    250259
    251 inline double  pdxxxstep( CoinDenseVector <double> &x, CoinDenseVector <double> &dx )
    252 {
    253 
    254 // Assumes x > 0.
    255 // Finds the maximum step such that x + step*dx >= 0.
    256 
    257      double step    = 1e+20;
    258 
    259      int n = x.size();
    260      double *x_elts = x.getElements();
    261      double *dx_elts = dx.getElements();
    262      for (int k = 0; k < n; k++)
    263           if (dx_elts[k] < 0)
    264                if ((x_elts[k] / (-dx_elts[k])) < step)
    265                     step = x_elts[k] / (-dx_elts[k]);
    266      return step;
     260inline double pdxxxstep(CoinDenseVector< double > &x, CoinDenseVector< double > &dx)
     261{
     262
     263  // Assumes x > 0.
     264  // Finds the maximum step such that x + step*dx >= 0.
     265
     266  double step = 1e+20;
     267
     268  int n = x.size();
     269  double *x_elts = x.getElements();
     270  double *dx_elts = dx.getElements();
     271  for (int k = 0; k < n; k++)
     272    if (dx_elts[k] < 0)
     273      if ((x_elts[k] / (-dx_elts[k])) < step)
     274        step = x_elts[k] / (-dx_elts[k]);
     275  return step;
    267276}
    268277//-----------------------------------------------------------------------
     
    270279//-----------------------------------------------------------------------
    271280
    272 inline double  pdxxxstep(int nset, int *set, CoinDenseVector <double> &x, CoinDenseVector <double> &dx )
    273 {
    274 
    275 // Assumes x > 0.
    276 // Finds the maximum step such that x + step*dx >= 0.
    277 
    278      double step    = 1e+20;
    279 
    280      int n = x.size();
    281      double *x_elts = x.getElements();
    282      double *dx_elts = dx.getElements();
    283      for (int k = 0; k < n; k++)
    284           if (dx_elts[k] < 0)
    285                if ((x_elts[k] / (-dx_elts[k])) < step)
    286                     step = x_elts[k] / (-dx_elts[k]);
    287      return step;
     281inline double pdxxxstep(int nset, int *set, CoinDenseVector< double > &x, CoinDenseVector< double > &dx)
     282{
     283
     284  // Assumes x > 0.
     285  // Finds the maximum step such that x + step*dx >= 0.
     286
     287  double step = 1e+20;
     288
     289  int n = x.size();
     290  double *x_elts = x.getElements();
     291  double *dx_elts = dx.getElements();
     292  for (int k = 0; k < n; k++)
     293    if (dx_elts[k] < 0)
     294      if ((x_elts[k] / (-dx_elts[k])) < step)
     295        step = x_elts[k] / (-dx_elts[k]);
     296  return step;
    288297}
    289298//-----------------------------------------------------------------------
     
    292301#endif
    293302#endif
     303
     304/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     305*/
Note: See TracChangeset for help on using the changeset viewer.