source: trunk/Couenne/src/expression/operators/bounds/CouenneExprBMul.hpp @ 591

Last change on this file since 591 was 591, checked in by pbelotti, 9 years ago

probing now applied even in absence of NLP solution

  • Property svn:keywords set to Author Date Id Revision
File size: 3.5 KB
Line 
1/* $Id: CouenneExprBMul.hpp 591 2011-05-31 17:03:03Z pbelotti $
2 *
3 * Name:    exprBMul.hpp
4 * Author:  Pietro Belotti
5 * Purpose: definition of operators to compute lower/upper bounds of multiplications
6 *
7 * (C) Carnegie-Mellon University, 2006.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11#ifndef COUENNE_EXPRBMUL_H
12#define COUENNE_EXPRBMUL_H
13
14#include "CouenneExprOp.hpp"
15#include "CouenneConfig.h"
16#include "CoinHelperFunctions.hpp"
17#include "CoinFinite.hpp"
18
19namespace Couenne {
20
21#define MUL_ZERO 1e-20
22#define MUL_INF  sqrt (COIN_DBL_MAX)
23
24/// product that avoids NaN's
25inline CouNumber safeProd (register CouNumber a, register CouNumber b) {
26
27  if (a >  MUL_INF) return (b < -MUL_ZERO) ? -COIN_DBL_MAX : (b > MUL_ZERO) ?  COIN_DBL_MAX : 0.;
28  if (a < -MUL_INF) return (b < -MUL_ZERO) ?  COIN_DBL_MAX : (b > MUL_ZERO) ? -COIN_DBL_MAX : 0.;
29
30  if (b >  MUL_INF) return (a < -MUL_ZERO) ? -COIN_DBL_MAX : (a > MUL_ZERO) ?  COIN_DBL_MAX : 0.;
31  if (b < -MUL_INF) return (a < -MUL_ZERO) ?  COIN_DBL_MAX : (a > MUL_ZERO) ? -COIN_DBL_MAX : 0.;
32
33  return a*b;
34}
35
36
37/// class to compute lower bound of a product based on the bounds of
38/// both factors
39
40class exprLBMul: public exprOp {
41
42 public:
43
44  /// Constructors, destructor
45  exprLBMul  (expression **al, int n): 
46    exprOp (al, n) {} //< non-leaf expression, with argument list
47
48  /// cloning method
49  expression *clone (Domain *d = NULL) const
50    {return new exprLBMul (clonearglist (d), nargs_);}
51
52  /// function for the evaluation of the expression
53  CouNumber operator () ();
54
55  /// print position (PRE, INSIDE, POST)
56  enum pos printPos () const
57    {return PRE;}
58
59  /// print operator
60  std::string printOp () const
61    {return "LB_Mul";}
62};
63
64
65/// compute sum
66
67inline CouNumber exprLBMul::operator () () {
68
69  register CouNumber n = (*(arglist_ [0])) ();
70  register CouNumber N = (*(arglist_ [1])) ();
71  register CouNumber d = (*(arglist_ [2])) ();
72  register CouNumber D = (*(arglist_ [3])) ();
73
74  if (d>=0)
75    if   (n>=0) return safeProd (n,d);
76    else        return safeProd (n,D);
77  else // d <= 0
78    if (N>0) {
79      CouNumber Nd = safeProd (N,d), nD;
80      if (n<0 && D>0 && 
81          (Nd > (nD = safeProd (n,D)))) return nD;
82      else                              return Nd;
83    }
84    else 
85      if (D>0) return safeProd (n,D);
86      else     return safeProd (N,D);
87}
88
89
90/// class to compute upper bound of a product based on the bounds of
91/// both factors
92
93class exprUBMul: public exprOp {
94
95 public:
96
97  /// Constructors, destructor
98  exprUBMul  (expression **al, int n): 
99    exprOp (al, n) {} //< non-leaf expression, with argument list
100
101  /// cloning method
102  expression *clone (Domain *d = NULL) const
103    {return new exprUBMul (clonearglist (d), nargs_);}
104
105  /// function for the evaluation of the expression
106  CouNumber operator () ();
107
108  /// print position (PRE, INSIDE, POST)
109  enum pos printPos () const
110    {return PRE;}
111
112  /// print operator
113  std::string printOp () const
114    {return "UB_Mul";}
115};
116
117
118/// compute sum
119
120inline CouNumber exprUBMul::operator () () {
121
122  //  exprOp:: operator () ();
123
124  register CouNumber n = (*(arglist_ [0])) ();
125  register CouNumber N = (*(arglist_ [1])) ();
126  register CouNumber d = (*(arglist_ [2])) ();
127  register CouNumber D = (*(arglist_ [3])) ();
128
129  if (d>0)
130    if (N<0) return safeProd (N,d);
131    else     return safeProd (N,D);
132  else // d <= 0
133    if (n<0) {
134      CouNumber nd = safeProd (n,d), ND;
135      if (N>0 && D>0 && 
136          ((ND = safeProd (N,D)) > nd)) return ND;
137      else                              return nd;
138    }
139    else 
140      if (D>0) return safeProd (N,D);
141      else     return safeProd (n,D);
142}
143
144}
145
146#endif
Note: See TracBrowser for help on using the repository browser.