/home/coin/SVN-release/OS-2.4.1/Couenne/src/convex/operators/conv-exprTrilinear-gencuts.cpp

Go to the documentation of this file.
00001 /* $Id: conv-exprTrilinear-gencuts.cpp 597 2011-06-02 10:09:33Z pbelotti $
00002  *
00003  * Name:       conv-exprTrilinear-gencuts.cpp.cpp
00004  * Source:     GNU C++
00005  * Author:     Sonia Cafieri
00006  * Purpose:    generate inequalities defining the convex envelope of a 
00007  *             trilinear monomial
00008  * History:    Nov 2010 work started
00009  */
00010 
00011 #include "CouenneCutGenerator.hpp"
00012 
00013 #include "CouenneTypes.hpp"
00014 #include "CouenneExprMul.hpp"
00015 #include "CouenneExprTrilinear.hpp"
00016 #include "CouenneProblem.hpp"
00017 #include "CouenneExprAux.hpp"
00018 
00019 #include <vector>
00020 
00021 //#define DEBUG
00022 using namespace Couenne;
00023 
00024 #define EPSILONT 1.e-6
00025 
00026 //typedef CouNumber double;
00027 
00029 void permutation3(int **ind,int *ibnd)
00030 {
00031   ind[0][0] = ibnd[0]; ind[0][1] = ibnd[1]; ind[0][2] = ibnd[2];
00032   ind[1][0] = ibnd[0]; ind[1][1] = ibnd[2]; ind[1][2] = ibnd[1];
00033   ind[2][0] = ibnd[1]; ind[2][1] = ibnd[0]; ind[2][2] = ibnd[2];
00034   ind[3][0] = ibnd[1]; ind[3][1] = ibnd[2]; ind[3][2] = ibnd[0];
00035   ind[4][0] = ibnd[2]; ind[4][1] = ibnd[0]; ind[4][2] = ibnd[1];
00036   ind[5][0] = ibnd[2]; ind[5][1] = ibnd[1]; ind[5][2] = ibnd[0];
00037 }
00038 
00039 
00041 void TriLinCuts (double *vlb, double *vub, int *varIndices,
00042                  std::vector <std::vector <int> >    &cutIndices,
00043                  std::vector <std::vector <double> > &cutCoeff,
00044                  std::vector <double>                &cutLb,
00045                  std::vector <double>                &cutUb) {
00046 
00047   // var indices
00048   int v1, v2, v3, v4;
00049   // number of cuts
00050   int defcons_size = 20;
00051   // bounds on cuts
00052   double *bnd = new double [12];    
00053 
00054   CouNumber cf = 1.;
00055 
00056   int **ind;
00057   ind = new int*[6];
00058   for(int i=0; i<6; i++) {
00059     ind[i] = new int[6];
00060   }
00061      
00062   int *ibnd; 
00063   ibnd = new int[3];
00064   ibnd[0] = varIndices[0]; ibnd[1] = varIndices[1]; ibnd[2] = varIndices[2];
00065 #ifdef DEBUG
00066 std::cout << "ibnd[0] =" << ibnd[0] << "  ibnd[1] =" << ibnd[1] << "  ibnd[2] =" << ibnd[2] << std::endl;
00067 std::cout << "vlb[ibnd[0]] =" << vlb[ibnd[0]] << "  vub[ibnd[0]] =" << vub[ibnd[0]] << std::endl;
00068 std::cout << "vlb[ibnd[1]] =" << vlb[ibnd[1]] << "  vub[ibnd[1]] =" << vub[ibnd[1]] << std::endl;
00069 std::cout << "vlb[ibnd[2]] =" << vlb[ibnd[2]] << "  vub[ibnd[2]] =" << vub[ibnd[2]] << std::endl;
00070 #endif
00071  
00072   // compute the 6 permutations of the 3 variables 
00073   permutation3(ind,ibnd);
00074 
00075   int i, flag=0, idx=0;
00076   i = 0;
00077   while(i < 6 && flag == 0) {
00078     if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] >=0 && vlb[ind[i][2]] <=0 && vub[ind[i][2]] >=0) {
00079       idx = i;   // store the index of the permutation satisfying the condition
00080       flag = 1;  // this is case 1
00081     }
00082     i++; 
00083   }
00084   i = 0;
00085   while(i < 6 && flag == 0) {
00086     if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] <=0 && vlb[ind[i][2]] <=0 
00087        && vub[ind[i][1]] >=0 && vub[ind[i][2]] >=0) {
00088       idx = i;   // store the index of the permutation satisfying the condition
00089       flag = 2;  // this is case 2
00090     }
00091     i++; 
00092   }
00093   i = 0;
00094   while(i < 6 && flag == 0) {
00095     if(vlb[ind[i][0]] <=0 && vlb[ind[i][1]] <=0 && vlb[ind[i][2]] <=0 
00096        && vub[ind[i][0]] >=0 && vub[ind[i][1]] >=0 && vub[ind[i][2]] >=0) {
00097       idx = i;   // store the index of the permutation satisfying the condition
00098       flag = 3;  // this is case 3
00099     }
00100     i++; 
00101   }
00102   i = 0;
00103   while(i < 6 && flag == 0) {
00104     if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] <=0 && vlb[ind[i][2]] <=0 
00105        && vub[ind[i][1]] >=0 && vub[ind[i][2]] <=0) {
00106       idx = i;   // store the index of the permutation satisfying the condition
00107       flag = 4;  // this is case 4
00108     }
00109     i++; 
00110   }
00111   i = 0;
00112   while(i < 6 && flag == 0) {
00113     if(vlb[ind[i][0]] <=0 && vlb[ind[i][1]] <=0 && vlb[ind[i][2]] <=0 
00114        && vub[ind[i][0]] >=0 && vub[ind[i][1]] >=0 && vub[ind[i][2]] <=0) {
00115       idx = i;   // store the index of the permutation satisfying the condition
00116       flag = 5;  // this is case 5
00117     }
00118     i++; 
00119   }
00120   i = 0;
00121   while(i < 6 && flag == 0) {
00122     if(vlb[ind[i][0]] <=0 && vub[ind[i][0]] >=0 && vub[ind[i][1]] <=0 && vub[ind[i][2]] <=0) {
00123       idx = i;   // store the index of the permutation satisfying the condition
00124       flag = 6;  // this is case 6
00125     }
00126     i++; 
00127   }
00128   i = 0;
00129   while(i < 6 && flag == 0) {
00130     if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] >=0 && vlb[ind[i][2]] >=0) {
00131       idx = i;   // store the index of the permutation satisfying the condition
00132       flag = 7;  // this is case 7
00133     }
00134     i++; 
00135   }
00136   i = 0;
00137   while(i < 6 && flag == 0) {
00138     if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] >=0 && vlb[ind[i][2]] <=0 && vub[ind[i][2]] <=0) {
00139       idx = i;   // store the index of the permutation satisfying the condition
00140       flag = 8;  // this is case 8
00141     }
00142     i++; 
00143   }
00144   i = 0;
00145   while(i < 6 && flag == 0) {
00146     if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] <=0 && vub[ind[i][1]] <=0 
00147        && vlb[ind[i][2]] <=0 && vub[ind[i][2]] <=0) {
00148       idx = i;   // store the index of the permutation satisfying the condition
00149       flag = 9;  // this is case 9
00150     }
00151     i++; 
00152   }
00153   i = 0;
00154   while(i < 6 && flag == 0) {
00155     if(vlb[ind[i][0]] <=0 && vub[ind[i][0]] <=0 && vlb[ind[i][1]] <=0 && vub[ind[i][1]] <=0 
00156        && vlb[ind[i][2]] <=0 && vub[ind[i][2]] <=0) {
00157       idx = i;   // store the index of the permutation satisfying the condition
00158       flag = 10;  // this is case 10
00159     }
00160     i++; 
00161   }
00162 
00163   if (flag==0) {
00164     std::cout << "ERROR: case not implemented" << std::endl;
00165     exit(0);
00166   }
00167 
00168   // var indices
00169   v1 = ind[idx][0]; 
00170   v2 = ind[idx][1]; 
00171   v3 = ind[idx][2];
00172   v4 = varIndices [3];
00173 
00174   // lower and upper bound on variables
00175   double xL1 = cf*vlb[v1];  double  xU1 = cf*vub[v1];
00176   double xL2 = vlb[v2];  double xU2 = vub[v2];
00177   double xL3 = vlb[v3];  double xU3 = vub[v3];
00178 
00179 #define prepareVectors(a) {                             \
00180                                                         \
00181     int size = (int) (cutIndices.size ());              \
00182                                                         \
00183     for (int i=0; i<a; i++) {                           \
00184                                                         \
00185       cutIndices. push_back (std::vector <int>    ());  \
00186       cutCoeff.   push_back (std::vector <double> ());  \
00187                                                         \
00188       cutLb.      push_back (-COUENNE_INFINITY);        \
00189       cutUb.      push_back ( COUENNE_INFINITY);        \
00190                                                         \
00191       for (int j=0; j<4; j++) {                         \
00192                                                         \
00193         cutIndices [size+i].push_back (-1);             \
00194         cutCoeff   [size+i].push_back (0.);             \
00195       }                                                 \
00196     }                                                   \
00197 }
00198 
00199   /*----------------------------------------------------------------------------------------*/
00200 
00201   // case 1
00202   if(flag == 1) {
00203 #ifdef DEBUG
00204     std::cout << " -- case 1 --" << std::endl;
00205 #endif
00206 
00207     double theta  = xL1*xU2*xU3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xU3; 
00208     double theta1 = xU1*xL2*xL3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3; 
00209 
00210     defcons_size = 12;
00211 
00212     prepareVectors(defcons_size);
00213 
00214     for(int ii = 0; ii < defcons_size; ii++) {
00215 
00216       cutIndices [ii][0] = v1; 
00217       cutIndices [ii][1] = v2; 
00218       cutIndices [ii][2] = v3; 
00219       cutIndices [ii][3] = v4;     
00220 
00221       cutCoeff [ii][3] = 1.;
00222     }
00223  
00224     cutCoeff [0][0] = -xU2*xU3; cutCoeff [0][1] = -xU1*xU3; cutCoeff [0][2] = -xU1*xU2;  bnd[0] = - 2.*xU1*xU2*xU3;
00225     cutCoeff [1][0] = -xU2*xL3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xU2;  bnd[1] = - xL1*xU2*xL3 - xL1*xU2*xU3;
00226     cutCoeff [2][0] = -xU2*xL3; cutCoeff [2][1] = -xL1*xL3; cutCoeff [2][2] = -xL1*xL2;  bnd[2] = - xL1*xU2*xL3 - xL1*xL2*xL3;
00227     cutCoeff [3][0] = -xL2*xU3; cutCoeff [3][1] = -xU1*xL3; cutCoeff [3][2] = -xU1*xL2;  bnd[3] = - xU1*xL2*xU3 - xU1*xL2*xL3; 
00228     cutCoeff [4][0] = -xL2*xL3; cutCoeff [4][1] = -xU1*xL3; cutCoeff [4][2] = -xL1*xL2;  bnd[4] = - xU1*xL2*xL3 - xL1*xL2*xL3;
00229     cutCoeff [5][0] = -xL2*xU3; cutCoeff [5][1] = -xL1*xU3; cutCoeff [5][2] = -(theta/(xU3-xL3));
00230     bnd[5] = (-(theta*xL3)/(xU3-xL3)) - xL1*xU2*xU3 - xU1*xL2*xU3 + xU1*xU2*xL3; 
00231 
00232     cutCoeff [6][0] = -xU2*xL3; cutCoeff [6][1] = -xU1*xL3; cutCoeff [6][2] = -xU1*xU2;  bnd[6] = - 2.*xU1*xU2*xL3;
00233     cutCoeff [7][0] = -xL2*xL3; cutCoeff [7][1] = -xU1*xU3; cutCoeff [7][2] = -xU1*xL2;  bnd[7] = - xU1*xL2*xU3 - xU1*xL2*xL3;
00234     cutCoeff [8][0] = -xU2*xU3; cutCoeff [8][1] = -xL1*xU3; cutCoeff [8][2] = -xL1*xL2;  bnd[8] = - xL1*xU2*xU3 - xL1*xL2*xU3;
00235     cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xL1*xL3; cutCoeff [9][2] = -xL1*xU2;  bnd[9] = - xL1*xU2*xU3 - xL1*xU2*xL3;
00236     cutCoeff [10][0] = -xL2*xU3; cutCoeff [10][1] = -xU1*xU3; cutCoeff [10][2] = -xL1*xL2;  bnd[10] = - xU1*xL2*xU3 - xL1*xL2*xU3;
00237     cutCoeff [11][0] = -xL2*xL3; cutCoeff [11][1] = -xL1*xL3; cutCoeff [11][2] = -(theta1/(xL3-xU3)); cutCoeff [11][3] = 1.;
00238     bnd[11] = (-(theta1*xU3)/(xL3-xU3)) - xU1*xL2*xL3 - xL1*xU2*xL3 + xU1*xU2*xU3;
00239   } // end if case 1
00240 
00241   /*----------------------------------------------------------------------------------------*/
00242 
00243   // case 2
00244   if(flag == 2) {
00245 #ifdef DEBUG
00246     std::cout << " -- case 2 --" << std::endl;
00247 #endif
00248 
00249     defcons_size = 12;
00250 
00251     prepareVectors (defcons_size);
00252  
00253     // compute the 6 permutations of the 3 variables 
00254     ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3; 
00255     permutation3(ind,ibnd);
00256     int i, flagg=0, idx=0;
00257     i = 0;
00258     while(i < 6 && flagg == 0) {
00259       if(vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]
00260          <= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]])
00261         {
00262           idx = i;   // store the index of the permutation satisfying the condition
00263           flagg = 1;  // condition is satisfied
00264         }
00265       i++; 
00266     }
00267     if (flagg==0) {
00268       std::cout << "ERROR!!!" << std::endl; exit(0);
00269     }
00270     v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
00271 
00272     double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
00273     double xL2(vlb[v2]); double xU2(vub[v2]);
00274     double xL3(vlb[v3]); double xU3(vub[v3]);
00275 
00276     for(int ii = 0; ii < defcons_size; ii++) {
00277 
00278       cutIndices [ii][0] = v1; 
00279       cutIndices [ii][1] = v2; 
00280       cutIndices [ii][2] = v3; 
00281       cutIndices [ii][3] = v4;     
00282       cutCoeff [ii][3] = 1.;
00283     }
00284 
00285     double theta1 = xL1*xL2*xL3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xU3;
00286     double theta2 = xU1*xL2*xU3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xL1*xU2*xU3;
00287 
00288     cutCoeff [0][0] = -xU2*xU3; cutCoeff [0][1] = -xU1*xU3; cutCoeff [0][2] = -xU1*xU2;  bnd[0] = - 2.*xU1*xU2*xU3; 
00289     cutCoeff [1][0] = -xL2*xL3; cutCoeff [1][1] = -xU1*xL3; cutCoeff [1][2] = -xU1*xL2;  bnd[1] = - 2.*xU1*xL2*xL3;
00290     cutCoeff [2][0] = -xU2*xL3; cutCoeff [2][1] = -xL1*xU3; cutCoeff [2][2] = -xL1*xU2;  bnd[2] = - xL1*xU2*xL3 - xL1*xU2*xU3;
00291     cutCoeff [3][0] = -xU2*xL3; cutCoeff [3][1] = -xL1*xL3; cutCoeff [3][2] = -xL1*xL2;  bnd[3] = - xL1*xU2*xL3 - xL1*xL2*xL3;
00292     cutCoeff [4][0] = -xL2*xU3; cutCoeff [4][1] = -(theta1/(xL2-xU2)); cutCoeff [4][2] = -xL1*xL2;  
00293     bnd[4] = (-(theta1*xU2)/(xL2-xU2)) - xL1*xL2*xL3 - xU1*xL2*xU3 + xU1*xU2*xL3;
00294     cutCoeff [5][0] = -xL2*xU3; cutCoeff [5][1] = -xL1*xU3; cutCoeff [5][2] = -(theta2/(xU3-xL3));
00295     bnd[5] = (-(theta2*xL3)/(xU3-xL3)) - xU1*xL2*xU3 - xL1*xU2*xU3 + xU1*xU2*xL3;
00296       
00297     if (  vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
00298           >= vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3]) {
00299 
00300       double theta1c = xU1*xL2*xL3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3;
00301       double theta2c = xU1*xL2*xL3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xL2*xU3;
00302 
00303       cutCoeff [6][0] = -xL2*xU3; cutCoeff [6][1] = -xU1*xU3; cutCoeff [6][2] = -xU1*xL2;  bnd[6] = - 2.*xU1*xL2*xU3;
00304       cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xU1*xL3; cutCoeff [7][2] = -xU1*xU2;  bnd[7] = - 2.*xU1*xU2*xL3;
00305       cutCoeff [8][0] = -xU2*xU3; cutCoeff [8][1] = -xL1*xU3; cutCoeff [8][2] = -xL1*xL2;  bnd[8] = - xL1*xU2*xU3 - xL1*xL2*xU3;
00306       cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xL1*xL3; cutCoeff [9][2] = -xL1*xU2;  bnd[9] = - xL1*xU2*xU3 - xL1*xU2*xL3;
00307       cutCoeff [10][0] = -xL2*xL3; cutCoeff [10][1] = -xL1*xL3; cutCoeff [10][2] = -(theta1c/(xL3-xU3));  
00308       bnd[10] = (-(theta1c*xU3)/(xL3-xU3)) - xU1*xL2*xL3 - xL1*xU2*xL3 + xU1*xU2*xU3;
00309       cutCoeff [11][0] = -xL2*xL3; cutCoeff [11][1] = -(theta2c/(xL2-xU2)); cutCoeff [11][2] = -xL1*xL2;
00310       bnd[11] = (-(theta2c*xU2)/(xL2-xU2)) - xU1*xL2*xL3 - xL1*xL2*xU3 + xU1*xU2*xU3;
00311 
00312     } else {
00313       double theta1c = xU1*xU2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xL2*xU3;
00314       double theta2c = xU1*xU2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xU2*xL3;
00315 
00316       cutCoeff [6][0] = -xL2*xU3; cutCoeff [6][1] = -xU1*xU3; cutCoeff [6][2] = -xU1*xL2;  bnd[6] = - 2.*xU1*xL2*xU3;
00317       cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xU1*xL3; cutCoeff [7][2] = -xU1*xU2;  bnd[7] = - 2.*xU1*xU2*xL3;
00318       cutCoeff [8][0] = -xL2*xL3; cutCoeff [8][1] = -xL1*xL3; cutCoeff [8][2] = -xL1*xU2;  bnd[8] = - xL1*xL2*xL3 - xL1*xU2*xL3;
00319       cutCoeff [9][0] = -xL2*xL3; cutCoeff [9][1] = -xL1*xU3; cutCoeff [9][2] = -xL1*xL2;  bnd[9] = - xL1*xL2*xL3 - xL1*xL2*xU3; 
00320       cutCoeff [10][0] = -xU2*xU3; cutCoeff [10][1] = -xL1*xU3; cutCoeff [10][2] = -(theta1c/(xU3-xL3));  
00321       bnd[10] = (-(theta1c*xL3)/(xU3-xL3)) - xU1*xU2*xU3 - xL1*xL2*xU3 + xU1*xL2*xL3; 
00322       cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -(theta2c/(xU2-xL2)); cutCoeff [11][2] = -xL1*xU2;
00323       bnd[11] = (-(theta2c*xL2)/(xU2-xL2)) - xU1*xU2*xU3 - xL1*xU2*xL3 + xU1*xL2*xL3; 
00324     }
00325       
00326   } // end if case 2
00327 
00328   /*----------------------------------------------------------------------------------------*/
00329 
00330   // case 3
00331   if(flag == 3) {
00332 #ifdef DEBUG
00333     std::cout << " -- case 3 --" << std::endl;
00334 #endif
00335 
00336     int last;
00337 
00338     if(vub[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3] + vlb[v1]*vub[v2]*vub[v3]
00339        <= vlb[v1]*vlb[v2]*vlb[v3] + 2.*vub[v1]*vub[v2]*vub[v3] &&
00340        vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3]
00341        <= vlb[v1]*vub[v2]*vub[v3] + 2.*vub[v1]*vlb[v2]*vlb[v3] &&
00342        vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3]
00343        <= vub[v1]*vlb[v2]*vub[v3] + 2.*vlb[v1]*vub[v2]*vlb[v3] &&
00344        vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3] + vlb[v1]*vub[v2]*vub[v3]
00345        <= vub[v1]*vub[v2]*vlb[v3] + 2.*vlb[v1]*vlb[v2]*vub[v3] ) {
00346     
00347       double theta3x = 0.5*(xL1*xU2*xU3 + xL1*xL2*xL3 - xU1*xU2*xL3 - xU1*xL2*xU3)/(xL1-xU1);
00348       double theta3y = 0.5*(xU1*xL2*xU3 + xL1*xL2*xL3 - xU1*xU2*xL3 - xL1*xU2*xU3)/(xL2-xU2);
00349       double theta3z = 0.5*(xU1*xU2*xL3 + xL1*xL2*xL3 - xU1*xL2*xU3 - xL1*xU2*xU3)/(xL3-xU3);
00350       double theta3c = xL1*xL2*xL3 - theta3x*xL1 - theta3y*xL2 - theta3z*xL3;
00351 
00352       defcons_size = 5;
00353       prepareVectors (defcons_size);
00354 
00355       for(int ii = 0; ii < 5; ii++) {
00356 
00357         cutIndices [ii][0] = v1; 
00358         cutIndices [ii][1] = v2; 
00359         cutIndices [ii][2] = v3; 
00360         cutIndices [ii][3] = v4;     
00361 
00362         cutCoeff [ii][3] = 1.;
00363       }
00364 
00365       cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xU2;  bnd[0] = - 2.*xL1*xU2*xL3;
00366       cutCoeff [1][0] = -xL2*xU3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xL2;  bnd[1] = - 2.*xL1*xL2*xU3;
00367       cutCoeff [2][0] = -xU2*xU3; cutCoeff [2][1] = -xU1*xU3; cutCoeff [2][2] = -xU1*xU2;  bnd[2] = - 2.*xU1*xU2*xU3;
00368       cutCoeff [3][0] = -xL2*xL3; cutCoeff [3][1] = -xU1*xL3; cutCoeff [3][2] = -xU1*xL2;  bnd[3] = - 2.*xU1*xL2*xL3;
00369       cutCoeff [4][0] = -theta3x; cutCoeff [4][1] = -theta3y; cutCoeff [4][2] = -theta3z;  bnd[4] = theta3c;
00370  
00371       last=4;
00372 
00373     } else if (vub[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3] + vlb[v1]*vub[v2]*vub[v3] 
00374                >= vlb[v1]*vlb[v2]*vlb[v3] + 2.*vub[v1]*vub[v2]*vub[v3]) {
00375 #ifdef DEBUG
00376 std::cout << "else if " << std::endl;
00377 #endif
00378          
00379       double theta1 = xU1*xU2*xL3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xU1*xL2*xU3;
00380       double theta2 = xU1*xL2*xU3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xL1*xU2*xU3;
00381       double theta3 = xU1*xU2*xL3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xL1*xU2*xU3;
00382 
00383       defcons_size = 6;
00384       prepareVectors (defcons_size);
00385 
00386       for(int ii = 0; ii < 6; ii++) {
00387 
00388         cutIndices [ii][0] = v1; 
00389         cutIndices [ii][1] = v2; 
00390         cutIndices [ii][2] = v3; 
00391         cutIndices [ii][3] = v4;     
00392 
00393         cutCoeff [ii][3] = 1.;
00394       }
00395 
00396       cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xU2;  bnd[0] = - 2.*xL1*xU2*xL3;
00397       cutCoeff [1][0] = -xL2*xU3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xL2;  bnd[1] = - 2.*xL1*xL2*xU3;
00398       cutCoeff [2][0] = -xL2*xL3; cutCoeff [2][1] = -xU1*xL3; cutCoeff [2][2] = -xU1*xL2;  bnd[2] = - 2.*xU1*xL2*xL3;
00399       cutCoeff [3][0] = -(theta1/(xU1-xL1)); cutCoeff [3][1] = -xU1*xU3; cutCoeff [3][2] = -xU1*xU2;  
00400       bnd[3] = (-(theta1*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xL2*xL3;
00401       cutCoeff [4][0] = -xU2*xU3; cutCoeff [4][1] = -xU1*xU3; cutCoeff [4][2] = -(theta2/(xU3-xL3));  
00402       bnd[4] = (-(theta2*xL3)/(xU3-xL3)) - xU1*xL2*xU3 - xL1*xU2*xU3 + xL1*xL2*xL3;
00403       cutCoeff [5][0] = -xU2*xU3; cutCoeff [5][1] = -(theta3/(xU2-xL2)); cutCoeff [5][2] = -xU1*xU2;
00404       bnd[5] = (-(theta3*xL2)/(xU2-xL2)) - xU1*xU2*xL3 - xL1*xU2*xU3 + xL1*xL2*xL3;
00405 
00406       last=5;
00407 
00408     } else {
00409 #ifdef DEBUG
00410 std::cout << "else " << std::endl;
00411 #endif
00412       // compute the 6 permutations of the 3 variables 
00413       ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3; 
00414       permutation3(ind,ibnd);
00415       int i, flagg=0, idx=0;
00416       i = 0;
00417       while(i < 6 && flagg == 0) {
00418         if (((vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]] 
00419              >= vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]] + 2.*vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]]) ))
00420           {
00421             idx = i;   // store the index of the permutation satisfying the condition
00422             flagg = 1;  // condition is satisfied
00423           }
00424         i++; 
00425       }
00426       if (flagg==0) {
00427         std::cout << "ERROR!!!" << std::endl; exit(0);
00428       }
00429       v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
00430 
00431       double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
00432       double xL2(vlb[v2]); double xU2(vub[v2]);
00433       double xL3(vlb[v3]); double xU3(vub[v3]);
00434 
00435       //} else if (vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3] 
00436       //         >= vlb[v1]*vub[v2]*vub[v3] + 2.*vub[v1]*vlb[v2]*vlb[v3]) {
00437 
00438       defcons_size = 6;
00439       prepareVectors (defcons_size);
00440 
00441       for(int ii = 0; ii < 6; ii++) {
00442 
00443         cutIndices [ii][0] = v1; 
00444         cutIndices [ii][1] = v2; 
00445         cutIndices [ii][2] = v3; 
00446         cutIndices [ii][3] = v4;     
00447 
00448         cutCoeff [ii][3] = 1.;
00449       }
00450 
00451       double theta1 = xL1*xL2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xL2*xU3;
00452       double theta2 = xU1*xU2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xL2*xU3;
00453       double theta3 = xL1*xL2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xU2*xL3;
00454 
00455       cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xU2;  bnd[0] = - 2.*xL1*xU2*xL3;
00456       cutCoeff [1][0] = -xL2*xU3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xL2;  bnd[1] = - 2.*xL1*xL2*xU3;
00457       cutCoeff [2][0] = -xU2*xU3; cutCoeff [2][1] = -xU1*xU3; cutCoeff [2][2] = -xU1*xU2;  bnd[2] = - 2.*xU1*xU2*xU3;
00458       cutCoeff [3][0] = -xL2*xL3; cutCoeff [3][1] = -(theta1/(xL2-xU2)); cutCoeff [3][2] = -xU1*xL2;  
00459       bnd[3] = (-(theta1*xU2)/(xL2-xU2)) - xL1*xL2*xL3 - xU1*xL2*xU3 + xL1*xU2*xU3;
00460       cutCoeff [4][0] = -(theta2/(xU1-xL1)); cutCoeff [4][1] = -xU1*xL3; cutCoeff [4][2] = -xU1*xL2;  
00461       bnd[4] = (-(theta2*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xU2*xU3;
00462       cutCoeff [5][0] = -xL2*xL3; cutCoeff [5][1] = -xU1*xL3; cutCoeff [5][2] = -(theta3/(xL3-xU3));
00463       bnd[5] = (-(theta3*xU3)/(xL3-xU3)) - xL1*xL2*xL3 - xU1*xU2*xL3 + xL1*xU2*xU3;
00464 
00465       last=5;
00466     }
00467 
00468     if(vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]
00469        >= vub[v1]*vub[v2]*vub[v3] + 2.*vlb[v1]*vlb[v2]*vlb[v3] &&
00470        vlb[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3] + vub[v1]*vub[v2]*vub[v3]
00471        >= vub[v1]*vlb[v2]*vlb[v3] + 2.*vlb[v1]*vub[v2]*vub[v3] &&
00472        vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3] + vub[v1]*vub[v2]*vub[v3]
00473        >= vlb[v1]*vub[v2]*vlb[v3] + 2.*vub[v1]*vlb[v2]*vub[v3] &&
00474        vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
00475        >= vlb[v1]*vlb[v2]*vub[v3] + 2.*vub[v1]*vub[v2]*vlb[v3] ) {
00476    
00477 #ifdef DEBUG
00478 std::cout << "2 - if " << std::endl;
00479 #endif
00480       double theta3x = 0.5*(xU1*xL2*xL3 + xU1*xU2*xU3 - xL1*xL2*xU3 - xL1*xU2*xL3)/(xU1-xL1);
00481       double theta3y = 0.5*(xL1*xU2*xL3 + xU1*xU2*xU3 - xL1*xL2*xU3 - xU1*xL2*xL3)/(xU2-xL2);
00482       double theta3z = 0.5*(xL1*xL2*xU3 + xU1*xU2*xU3 - xL1*xU2*xL3 - xU1*xL2*xL3)/(xU3-xL3);
00483       double theta3c = xU1*xU2*xU3 - theta3x*xU1 - theta3y*xU2 - theta3z*xU3;
00484 
00485       prepareVectors (5);
00486 
00487       for(int ii = last+1; ii <= last+5; ii++) {
00488 
00489         cutIndices [ii][0] = v1; 
00490         cutIndices [ii][1] = v2; 
00491         cutIndices [ii][2] = v3; 
00492         cutIndices [ii][3] = v4;     
00493         cutCoeff [ii][3] = 1.;
00494       }
00495 
00496       cutCoeff [last+1][0] = -xL2*xL3; cutCoeff [last+1][1] = -xL1*xL3; cutCoeff [last+1][2] = -xL1*xL2;  bnd[last+1] = - 2.*xL1*xL2*xL3;
00497       cutCoeff [last+2][0] = -xU2*xL3; cutCoeff [last+2][1] = -xU1*xL3; cutCoeff [last+2][2] = -xU1*xU2;  bnd[last+2] = - 2.*xU1*xU2*xL3;
00498       cutCoeff [last+3][0] = -xL2*xU3; cutCoeff [last+3][1] = -xU1*xU3; cutCoeff [last+3][2] = -xU1*xL2;  bnd[last+3] = - 2.*xU1*xL2*xU3;
00499       cutCoeff [last+4][0] = -xU2*xU3; cutCoeff [last+4][1] = -xL1*xU3; cutCoeff [last+4][2] = -xL1*xU2;  bnd[last+4] = - 2.*xL1*xU2*xU3;
00500       cutCoeff [last+5][0] = -theta3x; cutCoeff [last+5][1] = -theta3y; cutCoeff [last+5][2] = -theta3z;  bnd[last+5] = theta3c;
00501 
00502       defcons_size = last+6;
00503 
00504     } else if (vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3] 
00505                <= vub[v1]*vub[v2]*vub[v3] + 2.*vlb[v1]*vlb[v2]*vlb[v3]) {
00506 #ifdef DEBUG
00507 std::cout << "2 - else if" << std::endl;
00508 #endif
00509 
00510       double theta1 = xU1*xL2*xL3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3;
00511       double theta2 = xL1*xL2*xU3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3;
00512       double theta3 = xL1*xL2*xU3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xU1*xL2*xL3;
00513 
00514       prepareVectors (6);
00515 
00516       for(int ii = last+1; ii <= last+6; ii++) {
00517 
00518         cutIndices [ii][0] = v1; 
00519         cutIndices [ii][1] = v2; 
00520         cutIndices [ii][2] = v3; 
00521         cutIndices [ii][3] = v4;     
00522 
00523         cutCoeff [ii][3] = 1.;
00524       }
00525         
00526       cutCoeff [last+1][0] = -xU2*xL3; cutCoeff [last+1][1] = -xU1*xL3; cutCoeff [last+1][2] = -xU1*xU2;  bnd[last+1] = - 2.*xU1*xU2*xL3;
00527       cutCoeff [last+2][0] = -xL2*xU3; cutCoeff [last+2][1] = -xU1*xU3; cutCoeff [last+2][2] = -xU1*xL2;  bnd[last+2] = - 2.*xU1*xL2*xU3;       
00528       cutCoeff [last+3][0] = -xU2*xU3; cutCoeff [last+3][1] = -xL1*xU3; cutCoeff [last+3][2] = -xL1*xU2;  bnd[last+3] = - 2.*xL1*xU2*xU3;
00529       cutCoeff [last+4][0] = -xL2*xL3; cutCoeff [last+4][1] = -xL1*xL3; cutCoeff [last+4][2] = -(theta1/(xL3-xU3));  
00530       bnd[last+4] = (-(theta1*xU3)/(xL3-xU3)) - xU1*xL2*xL3 - xL1*xU2*xL3 + xU1*xU2*xU3;
00531       cutCoeff [last+5][0] = -(theta2/(xL1-xU1)); cutCoeff [last+5][1] = -xL1*xL3; cutCoeff [last+5][2] = -xL1*xL2;
00532       bnd[last+5] = (-(theta2*xU1)/(xL1-xU1)) - xL1*xL2*xU3 - xL1*xU2*xL3 + xU1*xU2*xU3;
00533       cutCoeff [last+6][0] = - xL2*xL3; cutCoeff [last+6][1] = -(theta3/(xL2-xU2)); cutCoeff [last+6][2] = -xL1*xL2;
00534       bnd[last+6] = (-(theta3*xU2)/(xL2-xU2)) - xL1*xL2*xU3 - xU1*xL2*xL3 + xU1*xU2*xU3;
00535 
00536       defcons_size = last+7;
00537 
00538     } else //if (vlb[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3] + vub[v1]*vub[v2]*vub[v3]
00539            //    <= vub[v1]*vlb[v2]*vlb[v3] + 2.*vlb[v1]*vub[v2]*vub[v3]) 
00540       {
00541 #ifdef DEBUG
00542 std::cout << "2 - another else if" << std::endl;
00543 std::cout << "v1 = " << v1 << " v2 =" << v2 << "  v3 =" << v3 << std::endl;
00544 #endif
00545       // compute the 6 permutations of the 3 variables 
00546       //ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3; 
00547       permutation3(ind,ibnd);
00548       int i, flagg=0, idx=0;
00549       i = 0;
00550       while(i < 6 && flagg == 0) {
00551         if (vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00552                <= vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + 2.*vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]])
00553           {
00554             idx = i;   // store the index of the permutation satisfying the condition
00555             flagg = 1;  // condition is satisfied
00556           }
00557         i++; 
00558       }
00559       if (flagg==0) {
00560         std::cout << "ERROR!!!" << std::endl; exit(0);
00561       }
00562       v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
00563 
00564       double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
00565       double xL2(vlb[v2]); double xU2(vub[v2]);
00566       double xL3(vlb[v3]); double xU3(vub[v3]);
00567 
00568       double theta1 = xL1*xL2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xU2*xL3;
00569       double theta2 = xU1*xU2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xU2*xL3;
00570       double theta3 = xL1*xL2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xU1*xU2*xU3;
00571 
00572       prepareVectors (6);
00573 
00574       for(int ii = last+1; ii <= last+6; ii++) {
00575 
00576         cutIndices [ii][0] = v1; 
00577         cutIndices [ii][1] = v2; 
00578         cutIndices [ii][2] = v3; 
00579         cutIndices [ii][3] = v4;     
00580 
00581         cutCoeff [ii][3] = 1.;
00582       }
00583        
00584       cutCoeff [last+1][0] = -xL2*xL3; cutCoeff [last+1][1] = -xL1*xL3; cutCoeff [last+1][2] = -xL1*xL2;  bnd[last+1] = - 2.*xL1*xL2*xL3;
00585       cutCoeff [last+2][0] = -xU2*xL3; cutCoeff [last+2][1] = -xU1*xL3; cutCoeff [last+2][2] = -xU1*xU2;  bnd[last+2] = - 2.*xU1*xU2*xL3;
00586       cutCoeff [last+3][0] = -xL2*xU3; cutCoeff [last+3][1] = -xU1*xU3; cutCoeff [last+3][2] = -xU1*xL2;  bnd[last+3] = - 2.*xU1*xL2*xU3;
00587       cutCoeff [last+4][0] = -(theta1/(xL1-xU1)); cutCoeff [last+4][1] = -xL1*xU3; cutCoeff [last+4][2] = -xL1*xU2;  
00588       bnd[last+4] = (-(theta1*xU1)/(xL1-xU1)) - xL1*xL2*xU3 - xL1*xU2*xL3 + xU1*xL2*xL3; 
00589       cutCoeff [last+5][0] = -xU2*xU3; cutCoeff [last+5][1] = -(theta2/(xU2-xL2)); cutCoeff [last+5][2] = -xL1*xU2;  
00590       bnd[last+5] = (-(theta2*xL2)/(xU2-xL2)) - xU1*xU2*xU3 - xL1*xU2*xL3 + xU1*xL2*xL3; 
00591       cutCoeff [last+6][0] = -xU2*xU3; cutCoeff [last+6][1] = -xL1*xU3; cutCoeff [last+6][2] = -(theta3/(xU3-xL3));  
00592       bnd[last+6] = (-(theta3*xL3)/(xU3-xL3)) - xL1*xL2*xU3 - xU1*xU2*xU3 + xU1*xL2*xL3;
00593 
00594       defcons_size = last+7;
00595     } 
00596     
00597   } // end if case 3
00598 
00599   /*----------------------------------------------------------------------------------------*/
00600 
00601   // case 4
00602   if(flag == 4) {
00603 #ifdef DEBUG
00604     std::cout << " -- case 4 --" << std::endl;
00605 #endif
00606 
00607     double theta  = xU1*xL2*xU3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xL1*xL2*xL3;
00608     double theta1 = xL1*xU2*xL3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xU1*xU2*xU3;
00609 
00610     defcons_size = 12;
00611 
00612     prepareVectors (defcons_size);
00613 
00614     for(int ii = 0; ii < defcons_size; ii++) {
00615       cutIndices [ii][0] = v1; 
00616       cutIndices [ii][1] = v2; 
00617       cutIndices [ii][2] = v3; 
00618       cutIndices [ii][3] = v4;     
00619       cutCoeff [ii][3] = 1.;
00620     }
00621 
00622     cutCoeff [0][0] = -xL2*xL3; cutCoeff [0][1] = -xU1*xL3; cutCoeff [0][2] = -xU1*xL2;  bnd[0] = - 2.*xU1*xL2*xL3;
00623     cutCoeff [1][0] = -xU2*xL3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xU2;  bnd[1] = - xL1*xU2*xL3 - xL1*xU2*xU3;
00624     cutCoeff [2][0] = -xU2*xL3; cutCoeff [2][1] = -xL1*xL3; cutCoeff [2][2] = -xL1*xL2;  bnd[2] = - xL1*xU2*xL3 - xL1*xL2*xL3;
00625     cutCoeff [3][0] = -xL2*xU3; cutCoeff [3][1] = -xU1*xU3; cutCoeff [3][2] = -xU1*xU2;  bnd[3] = - xU1*xL2*xU3 - xU1*xU2*xU3;
00626     cutCoeff [4][0] = -xU2*xU3; cutCoeff [4][1] = -xL1*xU3; cutCoeff [4][2] = -xU1*xU2;  bnd[4] = - xL1*xU2*xU3 - xU1*xU2*xU3; 
00627     cutCoeff [5][0] = -xL2*xU3; cutCoeff [5][1] = -(theta/(xL2-xU2)); cutCoeff [5][2] = -xL1*xL2;
00628     bnd[5] = (-(theta*xU2)/(xL2-xU2)) - xU1*xL2*xU3 - xL1*xL2*xL3 + xU1*xU2*xL3;
00629  
00630     cutCoeff [6][0] = -xU2*xL3; cutCoeff [6][1] = -xU1*xL3; cutCoeff [6][2] = -xU1*xU2;  bnd[6] = - 2.*xU1*xU2*xL3;
00631     cutCoeff [7][0] = -xU2*xU3; cutCoeff [7][1] = -xU1*xU3; cutCoeff [7][2] = -xU1*xL2;  bnd[7] = - xU1*xU2*xU3 - xU1*xL2*xU3;
00632     cutCoeff [8][0] = -xL2*xL3; cutCoeff [8][1] = -xL1*xL3; cutCoeff [8][2] = -xL1*xU2;  bnd[8] = - xL1*xU2*xL3 - xL1*xL2*xL3;
00633     cutCoeff [9][0] = -xL2*xL3; cutCoeff [9][1] = -xL1*xU3; cutCoeff [9][2] = -xL1*xL2;  bnd[9] = - xL1*xL2*xU3 - xL1*xL2*xL3;
00634     cutCoeff [10][0] = -xL2*xU3; cutCoeff [10][1] = -xL1*xU3; cutCoeff [10][2] = -xU1*xL2;  bnd[10] = - xL1*xL2*xU3 - xU1*xL2*xU3; 
00635     cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -(theta1/(xU2-xL2)); cutCoeff [11][2] = -xL1*xU2;
00636     bnd[11] = (-(theta1*xL2)/(xU2-xL2)) - xL1*xU2*xL3 - xU1*xU2*xU3 + xU1*xL2*xL3;
00637       
00638   } // end if case 4
00639 
00640   /*----------------------------------------------------------------------------------------*/
00641 
00642   // case 5
00643   if(flag == 5) {
00644 #ifdef DEBUG
00645     std::cout << " -- case 5 --" << std::endl;
00646  std::cout << "v1 = " << v1 << " v2 =" << v2 << "  v3 =" << v3 << std::endl;
00647 #endif
00648 
00649     defcons_size = 12;
00650     prepareVectors (defcons_size);
00651 
00652     // compute the permutations of the 3 variables 
00653     ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3; 
00654    ind[0][0] = ibnd[0]; ind[0][1] = ibnd[1]; ind[0][2] = ibnd[2];
00655    ind[1][0] = ibnd[1]; ind[1][1] = ibnd[0]; ind[1][2] = ibnd[2];
00656     int i, flagg=0, idx=0;
00657     i = 0;
00658     while(i < 2 && flagg == 0) {
00659       if(vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00660          >= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]) 
00661         {
00662           idx = i;   // store the index of the permutation satisfying the condition
00663           flagg = 1;  // condition is satisfied
00664         }
00665       i++; 
00666     }
00667     if (flagg==0) {
00668       std::cout << "ERROR!!!" << std::endl; exit(0);
00669     }
00670     v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
00671 #ifdef DEBUG
00672  std::cout << "v1 = " << v1 << " v2 =" << v2 << "  v3 =" << v3 << std::endl;
00673 #endif
00674 
00675     double xL1 = cf*vlb[v1]; double xU1 = cf*vub[v1];
00676     double xL2 = vlb[v2]; double xU2 = vub[v2];
00677     double xL3 = vlb[v3]; double xU3 = vub[v3];
00678 
00679     for(int ii = 0; ii < defcons_size; ii++) {
00680 
00681       cutIndices [ii][0] = v1; 
00682       cutIndices [ii][1] = v2; 
00683       cutIndices [ii][2] = v3; 
00684       cutIndices [ii][3] = v4;     
00685       cutCoeff [ii][3] = 1.;
00686     }
00687 
00688     if(vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
00689        <= vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]) {
00690 #ifdef DEBUG
00691     std::cout << " -- 5 if --" << std::endl;
00692 #endif
00693 
00694       double theta1 = xU1*xU2*xL3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xU1*xL2*xU3;
00695       double theta2 = xU1*xU2*xL3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xL1*xU2*xU3;
00696 
00697       cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xU2;  bnd[0] = - 2.*xL1*xU2*xL3; 
00698       cutCoeff [1][0] = -xL2*xL3; cutCoeff [1][1] = -xU1*xL3; cutCoeff [1][2] = -xU1*xL2;  bnd[1] = - 2.*xU1*xL2*xL3;
00699       cutCoeff [2][0] = -xU2*xU3; cutCoeff [2][1] = -xL1*xU3; cutCoeff [2][2] = -xL1*xL2;  bnd[2] = - xL1*xU2*xU3 - xL1*xL2*xU3;
00700       cutCoeff [3][0] = -xL2*xU3; cutCoeff [3][1] = -xU1*xU3; cutCoeff [3][2] = -xL1*xL2;  bnd[3] = - xU1*xL2*xU3 - xL1*xL2*xU3;
00701       cutCoeff [4][0] = -(theta1/(xU1-xL1)); cutCoeff [4][1] = -xU1*xU3; cutCoeff [4][2] = -xU1*xU2;  
00702       bnd[4] = (-(theta1*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xL2*xL3;
00703       cutCoeff [5][0] = -xU2*xU3; cutCoeff [5][1] = -(theta2/(xU2-xL2)); cutCoeff [5][2] = -xU1*xU2;
00704       bnd[5] = (-(theta2*xL2)/(xU2-xL2)) - xU1*xU2*xL3 - xL1*xU2*xU3 + xL1*xL2*xL3;
00705 
00706     } else {
00707 #ifdef DEBUG
00708     std::cout << " -- 5 else --" << std::endl;
00709 #endif
00710       double theta1 = xL1*xL2*xL3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xL1*xU2*xU3;
00711       double theta2 = xL1*xL2*xL3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xU3;
00712 
00713       cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xU2;  bnd[0] = - 2.*xL1*xU2*xL3;
00714       cutCoeff [1][0] = -xL2*xL3; cutCoeff [1][1] = -xU1*xL3; cutCoeff [1][2] = -xU1*xL2;  bnd[1] = - 2.*xU1*xL2*xL3;
00715       cutCoeff [2][0] = -xL2*xU3; cutCoeff [2][1] = -xU1*xU3; cutCoeff [2][2] = -xU1*xU2;  bnd[2] = - xU1*xL2*xU3 - xU1*xU2*xU3; 
00716       cutCoeff [3][0] = -xU2*xU3; cutCoeff [3][1] = -xL1*xU3; cutCoeff [3][2] = -xU1*xU2;  bnd[3] = - xL1*xU2*xU3 - xU1*xU2*xU3; 
00717       cutCoeff [4][0] = -(theta1/(xL1-xU1)); cutCoeff [4][1] = -xL1*xU3; cutCoeff [4][2] = -xL1*xL2;  
00718       bnd[4] = (-(theta1*xU1)/(xL1-xU1)) - xL1*xL2*xL3 - xL1*xU2*xU3 + xU1*xU2*xL3;
00719       cutCoeff [5][0] = -xL2*xU3; cutCoeff [5][1] = -(theta2/(xL2-xU2)); cutCoeff [5][2] = -xL1*xL2;
00720       bnd[5] = (-(theta2*xU2)/(xL2-xU2)) - xL1*xL2*xL3 - xU1*xL2*xU3 + xU1*xU2*xL3;
00721     }
00722 
00723     double theta1c = xL1*xU2*xL3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xL2*xU3;
00724     double theta2c = xU1*xU2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xU2*xL3;
00725 
00726     cutCoeff [6][0] = -xL2*xL3; cutCoeff [6][1] = -xL1*xL3; cutCoeff [6][2] = -xL1*xL2;  bnd[6] = - 2.*xL1*xL2*xL3;
00727     cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xU1*xL3; cutCoeff [7][2] = -xU1*xU2;  bnd[7] = - 2.*xU1*xU2*xL3;
00728     cutCoeff [8][0] = -xL2*xU3; cutCoeff [8][1] = -xL1*xU3; cutCoeff [8][2] = -xU1*xL2;  bnd[8] = - xL1*xL2*xU3 - xU1*xL2*xU3;
00729     cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xU1*xU3; cutCoeff [9][2] = -xU1*xL2;  bnd[9] = - xU1*xU2*xU3 - xU1*xL2*xU3;
00730     cutCoeff [10][0] = -(theta1c/(xL1-xU1)); cutCoeff [10][1] = -xL1*xU3; cutCoeff [10][2] = -xL1*xU2;  
00731     bnd[10] = (-(theta1c*xU1)/(xL1-xU1)) - xL1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xL3;
00732     cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -(theta2c/(xU2-xL2)); cutCoeff [11][2] = -xL1*xU2;
00733     bnd[11] = (-(theta2c*xL2)/(xU2-xL2)) - xU1*xU2*xU3 - xL1*xU2*xL3 + xU1*xL2*xL3;
00734 
00735   } // end if case 5
00736 
00737   /*----------------------------------------------------------------------------------------*/
00738 
00739   // case 6
00740   if(flag == 6) {
00741 #ifdef DEBUG
00742     std::cout << " -- case 6 --" << std::endl;
00743 #endif
00744 
00745     double theta = xU1*xU2*xL3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xU1*xL2*xU3;
00746     double theta1 = xL1*xU2*xL3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xL2*xU3;
00747 
00748     defcons_size = 12;
00749     prepareVectors (defcons_size);
00750 
00751     for (int ii = 0; ii < defcons_size; ii++) {
00752 
00753       cutIndices [ii][0] = v1; 
00754       cutIndices [ii][1] = v2; 
00755       cutIndices [ii][2] = v3; 
00756       cutIndices [ii][3] = v4;     
00757 
00758       cutCoeff [ii][3] = 1.;
00759     }
00760     
00761     cutCoeff [0][0] = -xL2*xL3; cutCoeff [0][1] = -xU1*xL3; cutCoeff [0][2] = -xU1*xL2;  bnd[0] = - 2.*xU1*xL2*xL3;
00762     cutCoeff [1][0] = -xU2*xU3; cutCoeff [1][1] = -xL1*xL3; cutCoeff [1][2] = -xL1*xU2;  bnd[1] = - xL1*xU2*xU3 - xL1*xU2*xL3;
00763     cutCoeff [2][0] = -xU2*xU3; cutCoeff [2][1] = -xL1*xU3; cutCoeff [2][2] = -xL1*xL2;  bnd[2] = - xL1*xU2*xU3 - xL1*xL2*xU3;
00764     cutCoeff [3][0] = -xL2*xU3; cutCoeff [3][1] = -xU1*xU3; cutCoeff [3][2] = -xL1*xL2;  bnd[3] = - xU1*xL2*xU3 - xL1*xL2*xU3;
00765     cutCoeff [4][0] = -xU2*xL3; cutCoeff [4][1] = -xL1*xL3; cutCoeff [4][2] = -xU1*xU2;  bnd[4] = - xL1*xU2*xL3 - xU1*xU2*xL3;
00766     cutCoeff [5][0] = -(theta/(xU1-xL1)); cutCoeff [5][1] = -xU1*xU3; cutCoeff [5][2] = -xU1*xU2;
00767     bnd[5] = (-(theta*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xL2*xL3;
00768 
00769     cutCoeff [6][0] = -xL2*xL3; cutCoeff [6][1] = -xL1*xL3; cutCoeff [6][2] = -xL1*xL2;  bnd[6] = - 2.*xL1*xL2*xL3;
00770     cutCoeff [7][0] = -xL2*xU3; cutCoeff [7][1] = -xL1*xU3; cutCoeff [7][2] = -xU1*xL2;  bnd[7] = - xL1*xL2*xU3 - xU1*xL2*xU3; 
00771     cutCoeff [8][0] = -xU2*xU3; cutCoeff [8][1] = -xU1*xU3; cutCoeff [8][2] = -xU1*xL2;  bnd[8] = - xU1*xU2*xU3 - xU1*xL2*xU3; 
00772     cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xU1*xL3; cutCoeff [9][2] = -xU1*xU2;  bnd[9] = - xU1*xU2*xU3 - xU1*xU2*xL3;
00773     cutCoeff [10][0] = -xU2*xL3; cutCoeff [10][1] = -xU1*xL3; cutCoeff [10][2] = -xL1*xU2;  bnd[10] = - xL1*xU2*xL3 - xU1*xU2*xL3; 
00774     cutCoeff [11][0] = -(theta1/(xL1-xU1)); cutCoeff [11][1] = -xL1*xU3; cutCoeff [11][2] = -xL1*xU2;
00775     bnd[11] = (-(theta1*xU1)/(xL1-xU1)) - xL1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xL3;
00776   } // end if case 6
00777     
00778   /*----------------------------------------------------------------------------------------*/
00779 
00780   // case 7
00781   if(flag == 7) {
00782 #ifdef DEBUG
00783     std::cout << " -- case 7 --" << std::endl;
00784 #endif
00785 
00786     defcons_size = 12;
00787     prepareVectors (defcons_size);
00788 
00789     if ((vlb[v1]<=EPSILONT && vlb[v2]<=EPSILONT && vlb[v3]<=EPSILONT) ||
00790         (vlb[v1]==vlb[v2] && vlb[v1]==vlb[v3] && vub[v1]==vub[v2] && vub[v1]==vub[v3])) {
00791 #ifdef DEBUG
00792       std::cout << " -- epsilonT --" << std::endl; 
00793 #endif
00794     } else {
00795     // compute the 6 permutations of the 3 variables 
00796     ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3; 
00797     permutation3(ind,ibnd);
00798     int i, flagg=0, idx=0;
00799     i = 0;
00800     while(i < 6 && flagg == 0) {
00801       if(vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00802          <= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]] &&
00803          vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00804          <= vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]) 
00805         {
00806           idx = i;   // store the index of the permutation satisfying the condition
00807           flagg = 1;  // condition is satisfied
00808         }
00809       i++; 
00810     }
00811     if (flagg==0) {
00812       std::cout << "ERROR!!!" << std::endl; exit(0);
00813     }
00814     }
00815     v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
00816 
00817     double xL1 = cf*vlb[v1]; double xU1 = cf*vub[v1];
00818     double xL2 = vlb[v2]; double xU2 = vub[v2];
00819     double xL3 = vlb[v3]; double xU3 = vub[v3];
00820 
00821     //if(vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3]
00822     // <= vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3] &&
00823     // vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3]
00824     // <= vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]) {
00825 
00826     double theta1 = xU1*xU2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xL2*xU3;
00827     double theta2 = xL1*xL2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xU2*xL3;
00828 
00829     for(int ii = 0; ii < defcons_size; ii++) {
00830       cutIndices [ii][0] = v1; 
00831       cutIndices [ii][1] = v2; 
00832       cutIndices [ii][2] = v3; 
00833       cutIndices [ii][3] = v4;     
00834 
00835       cutCoeff [ii][3] = 1.;
00836     }
00837     
00838     cutCoeff [0][0] = -xL2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xL2;  bnd[0] = - 2.*xL1*xL2*xL3; 
00839     cutCoeff [1][0] = -xU2*xU3; cutCoeff [1][1] = -xU1*xU3; cutCoeff [1][2] = -xU1*xU2;  bnd[1] = - 2.*xU1*xU2*xU3;
00840     cutCoeff [2][0] = -xL2*xU3; cutCoeff [2][1] = -xL1*xU3; cutCoeff [2][2] = -xU1*xL2;  bnd[2] = - xL1*xL2*xU3 - xU1*xL2*xU3;
00841     cutCoeff [3][0] = -xU2*xL3; cutCoeff [3][1] = -xU1*xL3; cutCoeff [3][2] = -xL1*xU2;  bnd[3] = - xU1*xU2*xL3 - xL1*xU2*xL3;
00842     cutCoeff [4][0] = -(theta1/(xU1-xL1)); cutCoeff [4][1] = -xU1*xL3; cutCoeff [4][2] = -xU1*xL2;  
00843     bnd[4] = (-(theta1*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xU2*xU3;
00844     cutCoeff [5][0] = -(theta2/(xL1-xU1)); cutCoeff [5][1] = -xL1*xU3; cutCoeff [5][2] = -xL1*xU2;
00845     bnd[5] = (-(theta2*xU1)/(xL1-xU1)) - xL1*xL2*xU3 - xL1*xU2*xL3 + xU1*xL2*xL3; 
00846     //}
00847     cutCoeff [6][0] = -xL2*xL3; cutCoeff [6][1] = -xU1*xL3; cutCoeff [6][2] = -xU1*xU2;  bnd[6] = - xU1*xU2*xL3 - xU1*xL2*xL3; 
00848     cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xL1*xL3; cutCoeff [7][2] = -xU1*xU2;  bnd[7] = - xU1*xU2*xL3 - xL1*xU2*xL3;
00849     cutCoeff [8][0] = -xL2*xL3; cutCoeff [8][1] = -xU1*xU3; cutCoeff [8][2] = -xU1*xL2;  bnd[8] = - xU1*xL2*xU3 - xU1*xL2*xL3; 
00850     cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xL1*xL3; cutCoeff [9][2] = -xL1*xU2;  bnd[9] = - xL1*xU2*xU3 - xL1*xU2*xL3;
00851     cutCoeff [10][0] = -xL2*xU3; cutCoeff [10][1] = -xU1*xU3; cutCoeff [10][2] = -xL1*xL2;  bnd[10] = - xU1*xL2*xU3 - xL1*xL2*xU3;
00852     cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -xL1*xU3; cutCoeff [11][2] = -xL1*xL2;  bnd[11] = - xL1*xU2*xU3 - xL1*xL2*xU3; 
00853   } // end if case 7
00854 
00855   /*----------------------------------------------------------------------------------------*/
00856       
00857   // case 8
00858   if(flag == 8) {
00859 #ifdef DEBUG
00860     std::cout << " -- case 8 --" << std::endl;
00861 #endif
00862 
00863     defcons_size = 12;
00864     prepareVectors (defcons_size);
00865 
00866     // compute the 6 permutations of the 3 variables 
00867     ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3; 
00868     permutation3(ind,ibnd);
00869     int i, flagg=0, idx=0;
00870     i = 0;
00871     while(i < 6 && flagg == 0) {
00872       if(vub[ind[i][2]] <=0  &&
00873          ((vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00874            >= vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]] &&
00875            vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00876            >= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]) ||
00877           (vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00878            >= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]])) ) 
00879         {
00880           idx = i;   // store the index of the permutation satisfying the condition
00881           flagg = 1;  // condition is satisfied
00882         }
00883       i++; 
00884     }
00885     if (flagg==0) {
00886       std::cout << "ERROR!!!" << std::endl; exit(0);
00887     }
00888     v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
00889 
00890     double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
00891     double xL2(vlb[v2]); double xU2(vub[v2]);
00892     double xL3(vlb[v3]); double xU3(vub[v3]);
00893 
00894     for(int ii = 0; ii < defcons_size; ii++) {
00895 
00896       cutIndices [ii][0] = v1; 
00897       cutIndices [ii][1] = v2; 
00898       cutIndices [ii][2] = v3; 
00899       cutIndices [ii][3] = v4;     
00900       cutCoeff [ii][3] = 1.;
00901     }
00902     
00903     // compute the 6 permutations of the 3 variables 
00904     //if(vub[v3]<=0) {
00905     cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xL2;  bnd[0] = - xL1*xU2*xL3 - xL1*xL2*xL3; 
00906     cutCoeff [1][0] = -xU2*xL3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xU2;  bnd[1] = - xL1*xU2*xL3 - xL1*xU2*xU3; 
00907     cutCoeff [2][0] = -xL2*xU3; cutCoeff [2][1] = -xU1*xL3; cutCoeff [2][2] = -xU1*xL2;  bnd[2] = - xU1*xL2*xU3 - xU1*xL2*xL3;
00908     cutCoeff [3][0] = -xL2*xU3; cutCoeff [3][1] = -xU1*xU3; cutCoeff [3][2] = -xU1*xU2;  bnd[3] = - xU1*xL2*xU3 - xU1*xU2*xU3;
00909     cutCoeff [4][0] = -xL2*xL3; cutCoeff [4][1] = -xU1*xL3; cutCoeff [4][2] = -xL1*xL2;  bnd[4] = - xU1*xL2*xL3 - xL1*xL2*xL3;
00910     cutCoeff [5][0] = -xU2*xU3; cutCoeff [5][1] = -xL1*xU3; cutCoeff [5][2] = -xU1*xU2;  bnd[5] = - xU1*xU2*xU3 - xL1*xU2*xU3;
00911 
00912     if(vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
00913        >= vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3] &&
00914        vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
00915        >= vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3]) {
00916 
00917       double theta1c = xU1*xL2*xL3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3;
00918       double theta2c = xU1*xL2*xU3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xL1*xU2*xU3;
00919 
00920       cutCoeff [6][0] = -xL2*xU3; cutCoeff [6][1] = -xL1*xU3; cutCoeff [6][2] = -xL1*xL2;  bnd[6] = - 2.*xL1*xL2*xU3;
00921       cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xU1*xL3; cutCoeff [7][2] = -xU1*xU2;  bnd[7] = - 2.*xU1*xU2*xL3;
00922       cutCoeff [8][0] = -xL2*xL3; cutCoeff [8][1] = -xU1*xU3; cutCoeff [8][2] = -xU1*xL2;  bnd[8] = - xU1*xL2*xU3 - xU1*xL2*xL3;
00923       cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xL1*xL3; cutCoeff [9][2] = -xL1*xU2;  bnd[9] = - xL1*xU2*xU3 - xL1*xU2*xL3;
00924       cutCoeff [10][0] = -xL2*xL3; cutCoeff [10][1] = -xL1*xL3; cutCoeff [10][2] = -(theta1c/(xL3-xU3));  
00925       bnd[10] = (-(theta1c*xU3)/(xL3-xU3)) - xU1*xL2*xL3 - xL1*xU2*xL3 + xU1*xU2*xU3;
00926       cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -xU1*xU3; cutCoeff [11][2] = -(theta2c/(xU3-xL3));
00927       bnd[11] = (-(theta2c*xL3)/(xU3-xL3)) - xU1*xL2*xU3 - xL1*xU2*xU3 + xL1*xL2*xL3;
00928 
00929     } else if(vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3]
00930               >= vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3]) {
00931 
00932       double theta1c = xL1*xL2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xL2*xU3;
00933       double theta2c = xL1*xU2*xL3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xU1*xU2*xU3;
00934 
00935       cutCoeff [6][0] = -xL2*xU3; cutCoeff [6][1] = -xL1*xU3; cutCoeff [6][2] = -xL1*xL2;  bnd[6] = - 2.*xL1*xL2*xU3;
00936       cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xU1*xL3; cutCoeff [7][2] = -xU1*xU2;  bnd[7] = - 2.*xU1*xU2*xL3;
00937       cutCoeff [8][0] = -xL2*xL3; cutCoeff [8][1] = -xL1*xL3; cutCoeff [8][2] = -xL1*xU2;  bnd[8] = - xL1*xL2*xL3 - xL1*xU2*xL3;
00938       cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xU1*xU3; cutCoeff [9][2] = -xU1*xL2;  bnd[9] = - xU1*xL2*xU3 - xU1*xU2*xU3;
00939       cutCoeff [10][0] = -xL2*xL3; cutCoeff [10][1] = -(theta1c/(xL2-xU2)); cutCoeff [10][2] = -xU1*xL2;  
00940       bnd[10] = (-(theta1c*xU2)/(xL2-xU2)) - xL1*xL2*xL3 - xU1*xL2*xU3 + xL1*xU2*xU3; 
00941       cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -(theta2c/(xU2-xL2)); cutCoeff [11][2] = -xL1*xU2;
00942       bnd[11] = (-(theta2c*xL2)/(xU2-xL2)) - xL1*xU2*xL3 - xU1*xU2*xU3 + xU1*xL2*xL3;
00943     }
00944     //} 
00945   } // end if case 8
00946 
00947   /*----------------------------------------------------------------------------------------*/
00948 
00949   // case 9
00950   if(flag == 9) {
00951 #ifdef DEBUG
00952     std::cout << " -- case 9 --" << std::endl;
00953 #endif
00954 
00955     defcons_size = 12;
00956     prepareVectors (defcons_size);
00957 
00958     // compute the 6 permutations of the 3 variables 
00959     ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3; 
00960     permutation3(ind,ibnd);
00961     int i, flagg=0, idx=0;
00962     i = 0;
00963     while(i < 6 && flagg == 0) {
00964       if(vlb[ind[i][0]] >=0  &&
00965          ((vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00966            <= vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]] &&
00967            vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00968            <= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]) || 
00969           (vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]
00970            <= vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]] &&
00971            vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]
00972            <= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]) ))
00973         {
00974           idx = i;   // store the index of the permutation satisfying the condition
00975           flagg = 1;  // condition is satisfied
00976         }
00977       i++; 
00978     }
00979     if (flagg==0) {
00980       std::cout << "ERROR!!!" << std::endl; exit(0);
00981     }
00982     v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
00983 
00984     double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
00985     double xL2(vlb[v2]); double xU2(vub[v2]);
00986     double xL3(vlb[v3]); double xU3(vub[v3]);
00987 
00988     for (int ii = 0; ii < defcons_size; ii++) {
00989 
00990       cutIndices [ii][0] = v1; 
00991       cutIndices [ii][1] = v2; 
00992       cutIndices [ii][2] = v3; 
00993       cutIndices [ii][3] = v4;     
00994       cutCoeff [ii][3] = 1.;
00995     }
00996     
00997     //if(vlb[v1]>=0) {
00998     if(vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
00999        <= vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3] &&
01000        vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
01001        <= vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3]) {
01002 
01003       double theta1 = xL1*xL2*xU3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3;
01004       double theta2 = xU1*xL2*xU3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xU1*xU2*xL3;
01005 
01006       cutCoeff [0][0] = -xU2*xU3; cutCoeff [0][1] = -xL1*xU3; cutCoeff [0][2] = -xL1*xU2;  bnd[0] = - 2.*xL1*xU2*xU3;
01007       cutCoeff [1][0] = -xL2*xL3; cutCoeff [1][1] = -xU1*xL3; cutCoeff [1][2] = -xU1*xL2;  bnd[1] = - 2.*xU1*xL2*xL3;
01008       cutCoeff [2][0] = -xL2*xU3; cutCoeff [2][1] = -xU1*xU3; cutCoeff [2][2] = -xL1*xL2;  bnd[2] = - xU1*xL2*xU3 - xL1*xL2*xU3;
01009       cutCoeff [3][0] = -xU2*xL3; cutCoeff [3][1] = -xL1*xL3; cutCoeff [3][2] = -xU1*xU2;  bnd[3] = - xU1*xU2*xL3 - xL1*xU2*xL3;
01010       cutCoeff [4][0] = -(theta1/(xL1-xU1)); cutCoeff [4][1] = -xL1*xL3; cutCoeff [4][2] = -xL1*xL2;  
01011       bnd[4] = (-(theta1*xU1)/(xL1-xU1)) - xL1*xL2*xU3 - xL1*xU2*xL3 + xU1*xU2*xU3;
01012       cutCoeff [5][0] = -(theta2/(xU1-xL1)); cutCoeff [5][1] = -xU1*xU3; cutCoeff [5][2] = -xU1*xU2;
01013       bnd[5] = (-(theta2*xL1)/(xU1-xL1)) - xU1*xL2*xU3 - xU1*xU2*xL3 + xL1*xL2*xL3;
01014 
01015     } else if(vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]
01016               <= vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3] &&
01017               vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]
01018               <= vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3]) {
01019 
01020       double theta1 = xU1*xU2*xU3 - xL1*xL2*xU3 - xU1*xU2*xL3 + xL1*xU2*xL3;
01021       double theta2 = xL1*xL2*xL3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xU3;
01022 
01023       cutCoeff [0][0] = -xU2*xU3; cutCoeff [0][1] = -xL1*xU3; cutCoeff [0][2] = -xL1*xU2;  bnd[0] = - 2.*xL1*xU2*xU3;
01024       cutCoeff [1][0] = -xL2*xL3; cutCoeff [1][1] = -xU1*xL3; cutCoeff [1][2] = -xU1*xL2;  bnd[1] = - 2.*xU1*xL2*xL3;
01025       cutCoeff [2][0] = -xL2*xU3; cutCoeff [2][1] = -xU1*xU3; cutCoeff [2][2] = -xU1*xU2;  bnd[2] = - xU1*xL2*xU3 - xU1*xU2*xU3;
01026       cutCoeff [3][0] = -xU2*xL3; cutCoeff [3][1] = -xL1*xL3; cutCoeff [3][2] = -xL1*xL2;  bnd[3] = - xL1*xL2*xL3 - xL1*xU2*xL3;
01027       cutCoeff [4][0] = -xU2*xL3; cutCoeff [4][1] = -(theta1/(xU2-xL2)); cutCoeff [4][2] = -xU1*xU2;  
01028       bnd[4] = (-(theta1*xL2)/(xU2-xL2)) - xU1*xU2*xU3 - xL1*xU2*xL3 + xL1*xL2*xU3;
01029       cutCoeff [5][0] = -xL2*xU3; cutCoeff [5][1] = -(theta2/(xL2-xU2)); cutCoeff [5][2] = -xL1*xL2;
01030       bnd[5] = (-(theta2*xU2)/(xL2-xU2)) - xL1*xL2*xL3 - xU1*xL2*xU3 + xU1*xU2*xL3;
01031     }
01032 
01033     cutCoeff [6][0] = -xL2*xL3; cutCoeff [6][1] = -xL1*xU3; cutCoeff [6][2] = -xL1*xL2;  bnd[6] = - xL1*xL2*xU3 - xL1*xL2*xL3;
01034     cutCoeff [7][0] = -xL2*xL3; cutCoeff [7][1] = -xL1*xL3; cutCoeff [7][2] = -xL1*xU2;  bnd[7] = - xL1*xU2*xL3 - xL1*xL2*xL3;
01035     cutCoeff [8][0] = -xL2*xU3; cutCoeff [8][1] = -xL1*xU3; cutCoeff [8][2] = -xU1*xL2;  bnd[8] = - xL1*xL2*xU3 - xU1*xL2*xU3;
01036     cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xU1*xU3; cutCoeff [9][2] = -xU1*xL2;  bnd[9] = - xU1*xU2*xU3 - xU1*xL2*xU3;
01037     cutCoeff [10][0] = -xU2*xU3; cutCoeff [10][1] = -xU1*xL3; cutCoeff [10][2] = -xU1*xU2;  bnd[10] = - xU1*xU2*xU3 - xU1*xU2*xL3;
01038     cutCoeff [11][0] = -xU2*xL3; cutCoeff [11][1] = -xU1*xL3; cutCoeff [11][2] = -xL1*xU2;  bnd[11] = - xL1*xU2*xL3 - xU1*xU2*xL3;
01039     //}
01040   } // end if case 9
01041 
01042   /*----------------------------------------------------------------------------------------*/
01043 
01044   // case 10
01045   if(flag == 10) {
01046 #ifdef DEBUG
01047     std::cout << " -- case 10 --" << std::endl;
01048 #endif
01049 
01050     defcons_size = 12;
01051     prepareVectors (defcons_size);
01052 
01053     ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3; 
01054     permutation3(ind,ibnd);
01055     int i, flagg=0, idx=0;
01056     i = 0;
01057     while(i < 6 && flagg == 0) {
01058       if(vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
01059          >= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]] &&
01060          vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
01061          >= vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]) 
01062         {
01063           idx = i;   // store the index of the permutation satisfying the condition
01064           flagg = 1;  // condition is satisfied
01065         }
01066       i++; 
01067     }
01068     if (flagg==0) {
01069        std::cout << "ERROR!!!" << std::endl; exit(0);
01070     }
01071     v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
01072 
01073       double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
01074       double xL2(vlb[v2]); double xU2(vub[v2]);
01075       double xL3(vlb[v3]); double xU3(vub[v3]);
01076 
01077       for(int ii = 0; ii < defcons_size; ii++) {
01078  
01079         cutIndices [ii][0] = v1; 
01080         cutIndices [ii][1] = v2; 
01081         cutIndices [ii][2] = v3; 
01082         cutIndices [ii][3] = v4;     
01083         cutCoeff [ii][3] = 1.;
01084       }
01085     
01086     // compute the 6 permutations of the 3 variables 
01087       cutCoeff [0][0] = -xL2*xL3; cutCoeff [0][1] = -xU1*xL3; cutCoeff [0][2] = -xU1*xU2;  bnd[0] = - xU1*xU2*xL3 - xU1*xL2*xL3;
01088       cutCoeff [1][0] = -xU2*xU3; cutCoeff [1][1] = -xL1*xL3; cutCoeff [1][2] = -xL1*xU2;  bnd[1] = - xL1*xU2*xU3 - xL1*xU2*xL3;
01089       cutCoeff [2][0] = -xU2*xL3; cutCoeff [2][1] = -xL1*xL3; cutCoeff [2][2] = -xU1*xU2;  bnd[2] = - xU1*xU2*xL3 - xL1*xU2*xL3;
01090       cutCoeff [3][0] = -xU2*xU3; cutCoeff [3][1] = -xL1*xU3; cutCoeff [3][2] = -xL1*xL2;  bnd[3] = - xL1*xU2*xU3 - xL1*xL2*xU3;
01091       cutCoeff [4][0] = -xL2*xU3; cutCoeff [4][1] = -xU1*xU3; cutCoeff [4][2] = -xL1*xL2;  bnd[4] = - xU1*xL2*xU3 - xL1*xL2*xU3;
01092       cutCoeff [5][0] = -xL2*xL3; cutCoeff [5][1] = -xU1*xU3; cutCoeff [5][2] = -xU1*xL2;  bnd[5] = - xU1*xL2*xU3 - xU1*xL2*xL3;
01093 
01094       //if(vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3]
01095       // >= vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3] &&
01096       // vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3]
01097       // >= vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]) {
01098 
01099       double theta1c = xL1*xU2*xL3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xL2*xU3;
01100       double theta2c = xU1*xU2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xL2*xU3;
01101 
01102       cutCoeff [6][0] = -xL2*xL3; cutCoeff [6][1] = -xL1*xL3; cutCoeff [6][2] = -xL1*xL2;  bnd[6] = - 2.*xL1*xL2*xL3;
01103       cutCoeff [7][0] = -xU2*xU3; cutCoeff [7][1] = -xU1*xU3; cutCoeff [7][2] = -xU1*xU2;  bnd[7] = - 2.*xU1*xU2*xU3;
01104       cutCoeff [8][0] = -xL2*xU3; cutCoeff [8][1] = -xL1*xU3; cutCoeff [8][2] = -xU1*xL2;  bnd[8] = - xL1*xL2*xU3 - xU1*xL2*xU3;
01105       cutCoeff [9][0] = -xU2*xL3; cutCoeff [9][1] = -xU1*xL3; cutCoeff [9][2] = -xL1*xU2;  bnd[9] = - xL1*xU2*xL3 - xU1*xU2*xL3;
01106       cutCoeff [10][0] = -(theta1c/(xL1-xU1)); cutCoeff [10][1] = -xL1*xU3; cutCoeff [10][2] = -xL1*xU2;  
01107       bnd[10] = (-(theta1c*xU1)/(xL1-xU1)) - xL1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xL3;
01108       cutCoeff [11][0] = -(theta2c/(xU1-xL1)); cutCoeff [11][1] = -xU1*xL3; cutCoeff [11][2] = -xU1*xL2;
01109       bnd[11] = (-(theta2c*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xU2*xU3;
01110 
01111   } // end if case 10
01112 
01113   /*----------------------------------------------------------------------------------------*/
01114 
01115   // get lower and upper bound on the constraints to be added to the problem
01116   for (int ii = 0; ii < defcons_size; ii++) {
01117 
01118     //printf ("flag = %d\n", flag);
01119 
01120     if (ii < defcons_size/2) {
01121 
01122       if (cf > 0) {cutLb[ii] =  bnd[ii];           cutUb[ii] = COUENNE_INFINITY;}
01123       if (cf < 0) {cutLb[ii] = -COUENNE_INFINITY;  cutUb[ii] = bnd[ii];}
01124 
01125     } else {
01126 
01127       if (cf > 0) {cutLb[ii] = -COUENNE_INFINITY;  cutUb[ii] = bnd[ii];}
01128       if (cf < 0) {cutLb[ii] =  bnd[ii];           cutUb[ii] = COUENNE_INFINITY;}
01129     }
01130 #ifdef DEBUG
01131 std::cout << ii << ") cutLb =" << cutLb[ii] << " " << "cutUb = " << cutUb[ii] << std::endl;
01132 #endif
01133   }
01134 
01135   for (int i=0; i<6; i++)
01136     delete [] ind [i];
01137 
01138   delete [] ibnd;
01139   delete [] bnd;
01140   delete [] ind;
01141 }
01142 
01143 
01144 // generate cuts for trilinear expressions
01145 void exprTrilinear::generateCuts (expression *w, 
01146                                   OsiCuts &cs, const CouenneCutGenerator *cg,
01147                                   t_chg_bounds *chg, int wind, 
01148                                   CouNumber lbw, CouNumber ubw) {
01149 
01150   expression **args = w -> Image () -> ArgList ();
01151 
01152   int *varInd = new int [4];
01153 
01154   for (int i=0; i<3; i++)
01155     varInd [i] = args [i] -> Index (); 
01156 
01157   varInd [3] = w -> Index ();   
01158 
01159   std::vector <std::vector <int> >    cutIndices;
01160   std::vector <std::vector <double> > cutCoeff;
01161   std::vector <double>                cutLb, cutUb;
01162 
01163   TriLinCuts (cg -> Problem () -> Lb (),
01164               cg -> Problem () -> Ub (),
01165               varInd,
01166               cutIndices, cutCoeff,
01167               cutLb, cutUb);
01168 
01169   // sanity check on returned vectors
01170   assert (cutIndices.size () == cutCoeff.size () && 
01171           cutIndices.size () == cutLb.size    () && 
01172           cutIndices.size () == cutUb.size    ());
01173 
01174   //printf ("trilinear cuts:\n");
01175 
01176   for (int i = (int) cutIndices.size (); i--;) {
01177 
01178     int 
01179        size = (int) cutIndices [i].size (),
01180       *ind  = new int [size];
01181 
01182     double *coe = new double [size];
01183 
01184     std::copy (cutIndices [i].begin (), cutIndices [i].end (), ind);
01185     std::copy (cutCoeff   [i].begin (), cutCoeff   [i].end (), coe);
01186 
01187     OsiRowCut cut (cutLb [i], cutUb [i], 4, 4, ind, coe);
01188     //cut.print ();
01189 
01190     delete [] ind;
01191     delete [] coe;
01192 
01193     if (cg -> Problem () -> bestSol ()) {
01194 
01195       // check validity of cuts by verifying they don't cut the
01196       // optimum in the node
01197 
01198       double *sol = cg -> Problem () -> bestSol ();
01199       const double
01200         *lb = cg -> Problem () -> Lb (), 
01201         *ub = cg -> Problem () -> Ub ();
01202 
01203       int nVars = cg -> Problem () -> nVars ();
01204 
01205       bool optIn = true;
01206 
01207       for (int i=0; i< nVars; i++)
01208         if ((sol [i] < lb [i] - COUENNE_EPS) ||
01209             (sol [i] > ub [i] + COUENNE_EPS)) {
01210           optIn = false;
01211           break;
01212         }
01213 
01214       if (optIn) {
01215 
01216         for (unsigned int i=0; i<cutIndices.size (); i++) {
01217 
01218           double chs = 0.;
01219 
01220           for (unsigned int j=0; j<cutIndices[i].size(); j++)
01221             chs += cutCoeff [i] [j] * sol [cutIndices [i] [j]];
01222 
01223           if ((chs < cutLb [i] - COUENNE_EPS) ||
01224               (chs > cutUb [i] + COUENNE_EPS)) {
01225 
01226             printf ("cut %d violates optimum:\n", i);
01227 
01228             if (cutLb [i] > -COUENNE_INFINITY) printf ("%g <= ", cutLb [i]);
01229             for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g x%d ", cutCoeff [i] [j],       cutIndices [i] [j]);  printf ("\n = ");
01230             for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g *%g ", cutCoeff [i] [j],  sol [cutIndices [i] [j]]); printf ("\n = ");
01231             for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g ",     cutCoeff [i] [j] * sol [cutIndices [i] [j]]); printf (" = %g", chs);
01232             if (cutUb [i] <  COUENNE_INFINITY) printf (" <= %g", cutUb [i]);
01233             printf ("\n");
01234 
01235           }
01236         }
01237       }
01238     }
01239 
01240     cs.insert (cut);
01241   }
01242 
01243   delete [] varInd;
01244 }

Generated on Thu Nov 10 03:05:43 2011 by  doxygen 1.4.7