00001
00002
00003
00004
00005
00006
00007
00008
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
00020
00041 bool exprQuad::alphaConvexify (const CouenneProblem *p,
00042 const OsiSolverInterface &si) {
00043
00044 if (matrix_.size () == 0)
00045 return false;
00046
00047
00048
00049
00050 int k=0,
00051 nDiag = bounds_.size (),
00052 *indexmap = new int [si.getNumCols ()],
00053 *indices = new int [nDiag];
00054
00055 CoinFillN (indexmap, si.getNumCols (), -1);
00056
00057
00058 double *diam = new double [nDiag];
00059
00060 bool changed_bounds = false;
00061
00062 for (std::map <exprVar *, std::pair <CouNumber, CouNumber> >::iterator i = bounds_.begin ();
00063 i != bounds_.end (); ++i, k++) {
00064
00065 #ifdef DEBUG
00066 printf ("b%04d. [%20g, %20g]\n", i->first->Index(), i->second.first, i->second.second);
00067 #endif
00068
00069 int index = i -> first -> Index ();
00070 indexmap [index] = k;
00071 indices [k] = index;
00072
00073 CouNumber
00074 lb = i -> first -> lb (),
00075 ub = i -> first -> ub ();
00076
00077
00078 if ((lb < -COUENNE_INFINITY) ||
00079 (ub > COUENNE_INFINITY)) {
00080
00081 delete [] diam;
00082 delete [] indexmap;
00083 delete [] indices;
00084
00085 #ifdef DEBUG
00086 printf ("unbounded, bailing out\n");
00087 #endif
00088
00089 return false;
00090 }
00091
00092
00093 if (fabs (lb - i->second.first) > COUENNE_EPS) {i -> second.first = lb; changed_bounds = true;}
00094 if (fabs (ub - i->second.second) > COUENNE_EPS) {i -> second.second = ub; changed_bounds = true;}
00095
00096 diam [k] = ub - lb;
00097 #ifdef DEBUG
00098 printf ("diam %4d - %4d = %g - %g = %g\n", index, k, ub, lb, diam [k]);
00099 #endif
00100 }
00101
00102 if (!changed_bounds) {
00103
00104 delete [] diam;
00105 delete [] indexmap;
00106 delete [] indices;
00107
00108 return true;
00109 }
00110
00111
00112
00113 double *matrix = new double [nDiag * nDiag];
00114
00115 CoinFillN (matrix, nDiag * nDiag, 0.);
00116
00117 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
00118
00119 int
00120 xind = row -> first -> Index (),
00121 irow = indexmap [xind];
00122
00123 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
00124
00125 int
00126 yind = col -> first -> Index (),
00127 icol = indexmap [yind];
00128
00129 double cell = col -> second * diam [irow] * diam [icol];
00130
00131 matrix [icol * nDiag + irow] = cell;
00132 if (irow != icol)
00133 matrix [irow * nDiag + icol] = cell;
00134 }
00135 }
00136
00137
00138
00139
00140
00141
00142 delete [] indexmap;
00143
00144
00145 double* eigval = new double [nDiag];
00146 int info;
00147
00148 #ifdef DEBUG
00149 printf ("nDiag = %d\n", nDiag);
00150 for (int i=0; i<nDiag; i++) {
00151 for (int j=0; j<nDiag; j++)
00152 printf ("%6.2f ", matrix [i*nDiag + j]);
00153 printf ("\n");
00154 }
00155 #endif
00156
00157 Ipopt::IpLapackDsyev (true,
00158 nDiag,
00159 matrix,
00160 nDiag,
00161 eigval,
00162 info);
00163
00164 if (info != 0) {
00165 printf ("exprQuad::alphaConvexify, warning: problem computing eigenvalue, info=%d\n", info);
00166 return false;
00167
00168 }
00169
00170
00171 eigen_.erase (eigen_.begin (), eigen_.end ());
00172
00173 for (int i=0; i<nDiag; i++) {
00174
00175 std::pair <CouNumber, std::vector <std::pair <exprVar *, CouNumber> > > eigenCoord;
00176
00177 eigenCoord. first = eigval [i];
00178
00179 for (int j=0; j<nDiag; j++) {
00180
00181 CouNumber elem = matrix [i * nDiag + j];
00182
00183 if (fabs (elem) > COUENNE_EPS)
00184 eigenCoord. second. push_back (std::pair <exprVar *, CouNumber>
00185 (p -> Var (indices [j]), elem));
00186 }
00187
00188 eigen_.push_back (eigenCoord);
00189 }
00190
00191 #ifdef DEBUG
00192 for (std::vector <std::pair <CouNumber,
00193 std::vector <std::pair <exprVar *, CouNumber> > > >::iterator i = eigen_.begin ();
00194 i != eigen_.end (); ++i) {
00195 printf (" [%g] -- ", i -> first);
00196 for (std::vector <std::pair <exprVar *, CouNumber> >::iterator j = i -> second. begin();
00197 j != i -> second. end (); ++j)
00198 printf ("(%d,%g) ", j -> first -> Index (), j -> second);
00199 printf ("\n");
00200 }
00201 #endif
00202
00203 delete [] indices;
00204 delete [] matrix;
00205 delete [] diam;
00206 delete [] eigval;
00207
00208 return true;
00209 }