/home/coin/SVN-release/OS-2.2.0/Couenne/src/convex/operators/alphaConvexify.cpp

Go to the documentation of this file.
00001 /* $Id: alphaConvexify.cpp 217 2009-07-08 17:02:07Z pbelotti $
00002  *
00003  * Name:    alphaConvexify.cpp
00004  * Author:  Stefan Vigerske
00005  * Purpose: create alpha-convexification of a quadratic expression
00006  *
00007  * (C) Carnegie-Mellon University, 2007. 
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 
00012 #include "CoinHelperFunctions.hpp"
00013 #include "OsiSolverInterface.hpp"
00014 #include "IpLapack.hpp"
00015 
00016 #include "exprQuad.hpp"
00017 #include "CouenneProblem.hpp"
00018 
00019 //#define DEBUG
00020 
00041 bool exprQuad::alphaConvexify (const CouenneProblem *p) {
00042   //const OsiSolverInterface &si) {
00043 
00044   if (matrix_.size () == 0)
00045     return false;
00046 
00047   // inverse of dIndex_ mapping: for each variable tell me the index
00048   // that it will have in dIndex_, or -1 if not there
00049 
00050   int k=0,
00051      nDiag    = bounds_.size (),
00052     *indexmap = new int [p -> nVars ()],
00053     *indices  = new int [nDiag];
00054 
00055   //CoinFillN (indexmap, si.getNumCols (), -1);
00056   CoinFillN (indexmap, p -> nVars (), -1);
00057 
00058   // box diameter
00059   double *diam = new double [nDiag];
00060 
00061   bool changed_bounds = false;
00062 
00063   for (std::map <exprVar *, std::pair <CouNumber, CouNumber> >::iterator i = bounds_.begin ();
00064        i != bounds_.end (); ++i, k++) {
00065 
00066 #ifdef DEBUG
00067   printf ("b%04d. [%20g, %20g]\n", i->first->Index(), i->second.first, i->second.second);
00068 #endif
00069 
00070     int index = i -> first -> Index ();
00071     indexmap [index] = k;
00072     indices [k] = index;
00073 
00074     CouNumber
00075       lb = i -> first -> lb (),
00076       ub = i -> first -> ub ();
00077 
00078     // if one variable unbounded, bail out
00079     if ((lb < -COUENNE_INFINITY) ||
00080         (ub >  COUENNE_INFINITY)) {
00081 
00082       delete [] diam;
00083       delete [] indexmap;
00084       delete [] indices;
00085 
00086 #ifdef DEBUG
00087       printf ("unbounded, bailing out\n");
00088 #endif
00089 
00090       return false;
00091     }
00092 
00093     // if no variable has changed bounds, no need to convexify
00094     if (fabs (lb - i->second.first)  > COUENNE_EPS) {i -> second.first  = lb; changed_bounds = true;}
00095     if (fabs (ub - i->second.second) > COUENNE_EPS) {i -> second.second = ub; changed_bounds = true;}
00096 
00097     diam [k] = ub - lb;
00098 #ifdef DEBUG
00099     printf ("diam %4d - %4d = %g - %g = %g\n", index, k, ub, lb, diam [k]);
00100 #endif
00101   }
00102 
00103   if (!changed_bounds) {
00104 
00105     delete [] diam;
00106     delete [] indexmap;
00107     delete [] indices;
00108 
00109     return true;
00110   }
00111 
00112   // lower triangular of quadratic term matrix, scaled by box diameter
00113 
00114   double *matrix = new double [nDiag * nDiag];
00115 
00116   CoinFillN (matrix, nDiag * nDiag, 0.);
00117 
00118   for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
00119 
00120     int 
00121       xind = row -> first -> Index (),
00122       irow = indexmap [xind];
00123 
00124     for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
00125 
00126       int 
00127         yind = col -> first -> Index (),
00128         icol = indexmap [yind];
00129 
00130       double cell = col -> second * diam [irow] * diam [icol];
00131 
00132       matrix          [icol * nDiag + irow] = cell;
00133       if (irow != icol) 
00134         matrix        [irow * nDiag + icol] = cell;
00135     }
00136   }
00137 
00138   // compute value of matrix entry = q_ij * (u_i-l_i) * (u_j-l_j)
00139   // I (Stefan) do not understand the Lapack docu; it says it needs
00140   // only the lower triangular but it seem to need both parts to
00141   // work correct
00142 
00143   delete [] indexmap;
00144 
00145   // compute minimum and maximum eigenvalue of matrix
00146   double* eigval = new double [nDiag];
00147   int info;
00148 
00149 #ifdef DEBUG
00150   printf ("nDiag = %d\n", nDiag);
00151   for (int i=0; i<nDiag; i++) {
00152     for (int j=0; j<nDiag; j++)
00153       printf ("%6.2f ", matrix [i*nDiag + j]);
00154     printf ("\n");
00155   }
00156 #endif 
00157 
00158   Ipopt::IpLapackDsyev (true,   // compute eigenvector
00159                         nDiag,  // dimension
00160                         matrix, // matrix
00161                         nDiag,  // "leading dimension" (number of columns, I think)
00162                         eigval, // output vector to store eigenvalues
00163                         info);  // output status variable
00164 
00165   if (info != 0) {
00166     printf ("exprQuad::alphaConvexify, warning: problem computing eigenvalue, info=%d\n", info);
00167     return false;
00168     //TODO error handling
00169   }
00170 
00171   // clean eigenvector structure
00172   eigen_.erase (eigen_.begin (), eigen_.end ());
00173 
00174   for (int i=0; i<nDiag; i++) {
00175 
00176     std::pair <CouNumber, std::vector <std::pair <exprVar *, CouNumber> > > eigenCoord;
00177 
00178     eigenCoord. first = eigval [i];
00179 
00180     for (int j=0; j<nDiag; j++) {
00181 
00182       CouNumber elem = matrix [i * nDiag + j];
00183 
00184       if (fabs (elem) > COUENNE_EPS) 
00185         eigenCoord. second. push_back (std::pair <exprVar *, CouNumber> 
00186                                        (p -> Var (indices [j]), elem));
00187     }
00188 
00189     eigen_.push_back (eigenCoord);
00190   }
00191 
00192 #ifdef DEBUG
00193   for (std::vector <std::pair <CouNumber, 
00194          std::vector <std::pair <exprVar *, CouNumber> > > >::iterator i = eigen_.begin ();
00195        i != eigen_.end (); ++i) {
00196     printf (" [%g] -- ", i -> first);
00197     for (std::vector <std::pair <exprVar *, CouNumber> >::iterator j = i -> second. begin();
00198          j != i -> second. end (); ++j)
00199       printf ("(%d,%g) ", j -> first -> Index (), j -> second);
00200     printf ("\n");
00201   }
00202 #endif
00203 
00204   delete [] indices;
00205   delete [] matrix;
00206   delete [] diam;
00207   delete [] eigval;
00208 
00209   return true;
00210 }

Generated on Thu Aug 5 03:02:56 2010 by  doxygen 1.4.7