00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <stdio.h>
00011
00012 #include "CoinHelperFunctions.hpp"
00013
00014 #include "CouenneExprJac.hpp"
00015 #include "CouenneProblem.hpp"
00016 #include "CouenneProblemElem.hpp"
00017 #include "CouenneExprAux.hpp"
00018
00019 using namespace Couenne;
00020
00021
00022
00023
00024 ExprJac::ExprJac ():
00025 nnz_ (0),
00026 iRow_ (NULL),
00027 jCol_ (NULL),
00028 expr_ (NULL),
00029 nRows_ (0) {}
00030
00031
00032
00033 ExprJac::~ExprJac () {
00034
00035 if (nnz_) {
00036
00037 free (iRow_);
00038 free (jCol_);
00039
00040 for (int i=0; i<nnz_; i++)
00041 delete expr_ [i];
00042
00043 free (expr_);
00044 }
00045 }
00046
00047
00048 ExprJac::ExprJac (const ExprJac &rhs)
00049 {operator= (rhs);}
00050
00051
00052
00053 ExprJac *ExprJac::clone ()
00054 {return new ExprJac (*this);}
00055
00056
00057
00058 ExprJac &ExprJac::operator= (const ExprJac &rhs) {
00059
00060 nnz_ = rhs. nnz_;
00061 nRows_ = rhs. nRows_;
00062
00063 iRow_ = (nnz_ && rhs.iRow_) ? (int *) malloc (nnz_ * sizeof (int)) : NULL;
00064 jCol_ = (nnz_ && rhs.jCol_) ? (int *) malloc (nnz_ * sizeof (int)) : NULL;
00065
00066 CoinCopyN (rhs.iRow_, nnz_, iRow_);
00067 CoinCopyN (rhs.jCol_, nnz_, jCol_);
00068
00069 if (nnz_) {
00070
00071 expr_ = (expression **) malloc (nnz_ * sizeof (expression *));
00072
00073 for (int i=0; i<nnz_; i++)
00074 expr_ [i] = expr_ [i] -> clone ();
00075
00076 } else expr_ = NULL;
00077
00078 return *this;
00079 }
00080
00081
00083
00084 #define reallocStep 100
00085 static void reAlloc (int nCur, int &nMax, int *&r, int *&c, expression **&e) {
00086
00087 if (nCur >= nMax) {
00088
00089 nMax += reallocStep;
00090
00091 r = (int *) realloc (r, nMax * sizeof (int));
00092 c = (int *) realloc (c, nMax * sizeof (int));
00093 e = (expression **) realloc (e, nMax * sizeof (expression *));
00094 }
00095 }
00096
00097
00098 ExprJac::ExprJac (CouenneProblem *p):
00099
00100 nnz_ (0),
00101 iRow_ (NULL),
00102 jCol_ (NULL),
00103 expr_ (NULL),
00104 nRows_ (0) {
00105
00113
00115 int
00116 cursize = 0,
00117 nRealCons = 0;
00118
00119 reAlloc (nnz_, cursize, iRow_, jCol_, expr_);
00120
00121
00122
00123 for (int i = 0; i < p -> nCons (); i++) {
00124
00125 CouenneConstraint *c = p -> Con (i);
00126
00127 if (c -> Body () -> Type () == AUX ||
00128 c -> Body () -> Type () == VAR)
00129 continue;
00130
00131
00132
00133
00134
00135 std::set <int> deplist;
00136
00137 c -> Body () -> DepList (deplist, STOP_AT_AUX);
00138
00139 int nTerms = 0;
00140
00141 for (std::set <int>::iterator k = deplist.begin (); k != deplist.end (); ++k) {
00142
00143 expression
00144 *J = c -> Body () -> differentiate (*k),
00145
00146
00147
00148 *sJ = J -> simplify (),
00149 *rJ = sJ ? sJ : J;
00150
00151 if (sJ)
00152 delete J;
00153
00154 if ((rJ -> Type () == CONST) &&
00155 (rJ -> Value () == 0.))
00156 continue;
00157
00158
00159
00160 reAlloc (nnz_ + 1, cursize, iRow_, jCol_, expr_);
00161
00162 rJ -> realign (p);
00163
00164 iRow_ [nnz_] = nRealCons;
00165 jCol_ [nnz_] = *k;
00166 expr_ [nnz_] = rJ;
00167
00168 nnz_++;
00169 nTerms++;
00170 }
00171
00172 if (nTerms) {
00173 ++nRealCons;
00174 ++nRows_;
00175 }
00176 }
00177
00178
00179
00181
00182 for (int i = 0; i < p -> nVars (); i++) {
00183
00184 exprVar *e = p -> Var (i);
00185
00186 if ((e -> Type () != AUX) ||
00187 (e -> Multiplicity () <= 0))
00188 continue;
00189
00190
00191
00192
00193
00194 std::set <int> deplist;
00195
00196 e -> Image () -> DepList (deplist, STOP_AT_AUX);
00197
00198 deplist.insert (e -> Index ());
00199
00200 int nTerms = 0;
00201
00202 for (std::set <int>::iterator k = deplist.begin (); k != deplist.end (); ++k) {
00203
00204 expression
00205 *J = (*k == e -> Index ()) ?
00206 new exprConst (-1.) :
00207 e -> Image () -> differentiate (*k),
00208
00209
00210
00211 *sJ = J -> simplify (),
00212 *rJ = sJ ? sJ : J;
00213
00214 if (sJ)
00215 delete J;
00216
00217 if ((rJ -> Type () == CONST) &&
00218 (rJ -> Value () == 0.))
00219 continue;
00220
00221 rJ -> realign (p);
00222
00223
00224
00225 reAlloc (nnz_ + 1, cursize, iRow_, jCol_, expr_);
00226
00227 iRow_ [nnz_] = nRealCons;
00228 jCol_ [nnz_] = *k;
00229 expr_ [nnz_] = rJ;
00230
00231 ++nnz_;
00232 ++nTerms;
00233 }
00234
00235 if (nTerms) {
00236 ++nRealCons;
00237 ++nRows_;
00238 }
00239 }
00240
00241 #ifdef DEBUG
00242 printf ("jacobian: %d nonzeros, %d rows\n", nnz_, nRows_);
00243
00244 for (int i=0; i<nnz_; i++) {
00245
00246 printf ("[%d,%d]: ", iRow_ [i], jCol_ [i]);
00247
00248 fflush (stdout);
00249 expr_ [i] -> print ();
00250 printf ("\n");
00251 }
00252 #endif
00253 }