Changeset 615


Ignore:
Timestamp:
Jun 8, 2011 4:36:24 PM (9 years ago)
Author:
pbelotti
Message:

fixed horrible bug on differentiating products. Added some latex to doxydoc. Completing Hessian and Jacobian

Location:
trunk/Couenne/src
Files:
32 edited

Legend:

Unmodified
Added
Removed
  • trunk/Couenne/src/bound_tightening/obbt.cpp

    r490 r615  
    180180                          t_chg_bounds *chg_bds) {
    181181
     182  // TODO: set up list of hopeless variables and do different OBBT
     183  // saving lps
     184
    182185  // Check if cs contains only one cut and if it is of the form 1 <=
    183186  // x0 <= -1. That means a previous cut generator has determined that
  • trunk/Couenne/src/bound_tightening/obbt_iter.cpp

    r532 r615  
    5757
    5858int CouenneProblem::obbt_iter (OsiSolverInterface *csi,
    59            t_chg_bounds *chg_bds,
    60            const CoinWarmStart *warmstart,
    61            Bonmin::BabInfo *babInfo,
    62            double *objcoe,
    63            int sense,
    64            int index) const {
     59                               t_chg_bounds *chg_bds,
     60                               const CoinWarmStart *warmstart,
     61                               Bonmin::BabInfo *babInfo,
     62                               double *objcoe,
     63                               int sense,
     64                               int index) const {
    6565
    6666  // TODO: do NOT apply OBBT if this is a variable of the form
  • trunk/Couenne/src/branch/CouenneBranchingObject.cpp

    r604 r615  
    159159  }
    160160
    161   if   (l <  -large_bound)
     161  if   (l <  -large_bound) {
    162162    if (u <=  large_bound) // ]-inf,u]
    163163      brpt =  ((brpt < -COUENNE_EPS) ? (AGGR_MUL * (-1. + brpt)) :
    164164               (brpt >  COUENNE_EPS) ? 0.                        : -AGGR_MUL);
    165   else
     165  } else
    166166    if (u >  large_bound) // [l,+inf[
    167167      brpt = ((brpt >  COUENNE_EPS) ? (AGGR_MUL *  ( 1. + brpt)) :
  • trunk/Couenne/src/branch/CouenneComplObject.hpp

    r490 r615  
    1616namespace Couenne {
    1717
    18 /// OsiObject for complementarity constraints $x_1 x_2 <=/>=/= 0$.
     18/// OsiObject for complementarity constraints \f$ x_1 x_2 \ge,\le,= 0 \f$.
    1919///
    20 /// Associated with two variables x_1 and x_2, branches with either x_1=0 or x_2=0
     20/// Associated with two variables \f$ x_1\f$ and \f$x_2\f$, branches with either \f$x_1=0\f$ or \f$x_2=0\f$
    2121
    2222class CouenneComplObject: public CouenneObject {
  • trunk/Couenne/src/expression/CouenneExprAux.hpp

    r607 r615  
    198198
    199199  /// return its sign in the definition constraint
    200   inline enum auxSign sign () const
     200  virtual inline enum auxSign sign () const
    201201  {return sign_;}
    202202};
  • trunk/Couenne/src/expression/operators/CouenneExprAbs.hpp

    r490 r615  
    1919namespace Couenne {
    2020
    21 /// class for \f$w=|f(x)|\f$
     21/// class for \f$ |f(x)| \f$
    2222
    2323class exprAbs: public exprUnary {
  • trunk/Couenne/src/expression/operators/CouenneExprBinProd.hpp

    r490 r615  
    1717namespace Couenne {
    1818
    19   /// class for multiplications
     19  /// class for \f$ \prod_{i=1}^n f_i(x) \f$ with \f$ f_i(x) \f$ all binary
    2020
    2121  class exprBinProd: public exprMul {
  • trunk/Couenne/src/expression/operators/CouenneExprCeil.hpp

    r595 r615  
    1616namespace Couenne {
    1717
    18 /// class ceiling
     18/// class ceiling, \f$ \lceil f(x) \rceil \f$
    1919
    2020class exprCeil: public exprUnary {
  • trunk/Couenne/src/expression/operators/CouenneExprCos.hpp

    r490 r615  
    1616namespace Couenne {
    1717
    18 /// class cosine
     18/// class cosine, \f$ \cos f(x) \f$
    1919
    2020class exprCos: public exprUnary {
  • trunk/Couenne/src/expression/operators/CouenneExprDiv.hpp

    r490 r615  
    2020#define BR_MULT      1e-3
    2121
    22 /// class for divisions
     22/// class for divisions, \f$ \frac{f(x)}{g(x)} \f$
    2323
    2424class exprDiv: public exprOp {
  • trunk/Couenne/src/expression/operators/CouenneExprExp.hpp

    r490 r615  
    1818namespace Couenne {
    1919
    20 /// class for the exponential
     20/// class for the exponential, \f$ e^{f(x)} \f$
    2121
    2222class exprExp: public exprUnary {
  • trunk/Couenne/src/expression/operators/CouenneExprFloor.hpp

    r595 r615  
    1616namespace Couenne {
    1717
    18 /// class floor
     18/// class floor, \f$ \lfloor f(x) \rfloor \f$
    1919
    2020class exprFloor: public exprUnary {
  • trunk/Couenne/src/expression/operators/CouenneExprGroup.hpp

    r490 r615  
    2121class Domain;
    2222
    23 /// class Group, with constant, linear and nonlinear terms
     23/// class Group, with constant, linear and nonlinear terms: \f$ a_0 + \sum_{i=1}^n a_i x_i \f$
    2424
    2525class exprGroup: public exprSum {
  • trunk/Couenne/src/expression/operators/CouenneExprInv.hpp

    r490 r615  
    3131
    3232
    33 /// class inverse (1/f(x))
     33/// class inverse: \f$ 1/f(x) \f$
    3434
    3535class exprInv: public exprUnary {
  • trunk/Couenne/src/expression/operators/CouenneExprLog.hpp

    r490 r615  
    1717namespace Couenne {
    1818
    19 /// class logarithm
     19/// class logarithm, \f$ \log f(x)\f$
    2020
    2121class exprLog: public exprUnary {
  • trunk/Couenne/src/expression/operators/CouenneExprMul.hpp

    r490 r615  
    1818namespace Couenne {
    1919
    20 /// class for multiplications
     20/// class for multiplications, \f$ \prod_{i=1}^n f_i(x) \f$
    2121
    2222class exprMul: public exprOp {
  • trunk/Couenne/src/expression/operators/CouenneExprMultiLin.hpp

    r490 r615  
    1717namespace Couenne {
    1818
    19   /// class for multiplications
     19  /// another class for multiplications, \f$ \prod_{i=1}^n f_i(x) \f$
    2020
    2121  class exprMultiLin: public exprMul {
  • trunk/Couenne/src/expression/operators/CouenneExprNorm.hpp

    r490 r615  
    1616namespace Couenne {
    1717
    18 class exprNorm: public exprOp {
     18  /// Class for \f$ p \f$-norms, \f$ || f(x)||_p = \left(\sum_{i=1}^n f_i(x)^p\right)^{\frac{1}{p}} \f$
    1919
    20 };
     20  class exprNorm: public exprOp {
     21
     22  };
    2123
    2224}
  • trunk/Couenne/src/expression/operators/CouenneExprOpp.hpp

    r490 r615  
    2323
    2424
    25 /// class opposite 
     25/// class opposite, \f$ -f(x) \f$
    2626
    2727class exprOpp: public exprUnary {
  • trunk/Couenne/src/expression/operators/CouenneExprPow.hpp

    r513 r615  
    2424
    2525
    26 /// Power of an expression (binary operator)
     26/// Power of an expression (binary operator), \f$ f(x)^k\f$ with \f$ k\f$ constant
    2727
    2828class exprPow: public exprOp {
  • trunk/Couenne/src/expression/operators/CouenneExprSin.hpp

    r490 r615  
    4444
    4545
    46 /// class for sin f(x)
     46/// class for \f$ \sin f(x)\f$
    4747class exprSin: public exprUnary {
    4848
  • trunk/Couenne/src/expression/operators/CouenneExprSub.hpp

    r490 r615  
    1818namespace Couenne {
    1919
    20 /// class for subtraction
     20/// class for subtraction, \f$ f(x) - g(x) \f$
    2121
    2222class exprSub: public exprOp {
  • trunk/Couenne/src/expression/operators/CouenneExprSum.hpp

    r490 r615  
    1818namespace Couenne {
    1919
    20 /// class sum 
     20/// class sum, \f$ \sum_{i=1}^n f_i(x) \f$
    2121
    2222class exprSum: public exprOp {
  • trunk/Couenne/src/expression/operators/exprMul.cpp

    r560 r615  
    121121      for (int j = 0; j < nargs_; j++)
    122122        if (i!=j)
    123           alm [j] = arglist_ [i] -> clone (); //new exprClone (arglist_ [j]);
     123          alm [j] = arglist_ [j] -> clone (); //new exprClone (arglist_ [j]);
    124124
    125125      als [nonconst++] = new exprMul (alm, nargs_);
  • trunk/Couenne/src/expression/partial/CouenneExprHess.cpp

    r612 r615  
    1717using namespace Couenne;
    1818
    19 //#define DEBUG
     19#define DEBUG
    2020
    2121ExprHess::ExprHess ():
     
    5151}
    5252
     53
    5354/// code for refilling jacobian
    5455void HessElemFill (int i, int level,
     
    5859                   CouenneProblem *);
    5960
     61
    6062/// code for refilling arrays through realloc
    6163static void reAlloc (int nCur, int &nMax, int *&r, int *&c, int *&n, int **&l, expression ***&e);
     
    7476#ifdef DEBUG
    7577  printf ("creating Hessian\n");
     78
     79  p -> print ();
    7680#endif
    7781
     
    8387  int nLevels = 0; // depth of this 3D structure (<= #con + #var + 1)
    8488
    85   // fill in dependence list for the objective...
     89  // fill in dependence list ///////////////////////////////////////////
     90
     91  // 1) for the objective (nLevels = 0). Note: STOP_AT_AUX is
     92  //    sufficient to get to the leaves of this sum of squares
    8693
    8794  p -> Obj (0) -> Body () -> DepList (deplist [nLevels++], STOP_AT_AUX);
    8895
    89   // ... constraints...
     96  // 2) for the constraints
    9097
    9198  for (int i = 0; i < p -> nCons (); i++) {
    9299
    93     CouenneConstraint *c = p -> Con (i);
    94 
    95     if (c -> Body () -> Type () == AUX ||
    96         c -> Body () -> Type () == VAR ||
    97         c -> Body () -> Type () == CONST)
     100    expression *c = p -> Con (i) -> Body ();
     101
     102    enum nodeType ntype = c -> Type ();
     103
     104    if (ntype == AUX ||
     105        ntype == VAR ||
     106        ntype == CONST)
    98107      continue;
    99108
    100     c -> Body () -> DepList (deplist [nLevels++], STOP_AT_AUX);
    101   }
    102 
    103   // ... and auxiliaries
     109    c -> DepList (deplist [nLevels++], STOP_AT_AUX);
     110  }
     111
     112  // 3) and the auxiliaries
    104113
    105114  for (int i = 0; i < p -> nVars (); i++) {
     
    107116    exprVar *e = p -> Var (i);
    108117
    109     if ((e -> Type () != AUX) ||
     118    if ((e -> Type         () != AUX) ||
    110119        (e -> Multiplicity () <= 0))
    111120      continue;
    112121
    113     e -> Image () -> DepList (deplist [nLevels], STOP_AT_AUX);
    114 
    115     // insert auxiliary variable defined here as well
    116     // No, don't. No second order information will be there
    117     //deplist [nLevels].insert (e -> Index ());
    118 
    119     nLevels++;
     122    e -> Image () -> DepList (deplist [nLevels++], STOP_AT_AUX);
    120123  }
    121124
     
    141144    curSize = 0; // used to expand array through realloc
    142145
    143   /// for each variable (filling a row of the hessian)
     146  /// for each variable, fill a row of the hessian
    144147  for (int i=0; i < nVars; i++) {
    145148
    146149    // create dense row. These will become numL later
    147     int          *row = (int          *) malloc (nVars * sizeof (int));
    148     int         **lam = (int         **) malloc (nVars * sizeof (int *));
     150    int          *rnz = (int          *) malloc (nVars * sizeof (int));   // row's number of nonzero
     151    int         **lam = (int         **) malloc (nVars * sizeof (int *)); // row's vectors of indices of nonzeros
    149152    expression ***eee = (expression ***) malloc (nVars * sizeof (expression **));
    150153
    151     CoinFillN (row, nVars, 0);
    152 
    153     for (int j=0; j<nVars; j++) {
    154       lam [j] = NULL;
    155       eee [j] = NULL;
    156     }
     154    CoinFillN (rnz, nVars, 0);
     155
     156    // No CoinFillN for int** and expression***
     157    for (int j=nVars; j--;) {
     158      *lam++ = NULL;
     159      *eee++ = NULL;
     160    }
     161
     162    lam -= nVars;
     163    eee -= nVars;
    157164
    158165    // scan all levels
    159 
    160166    int level = 0;
    161167
    162168    /// fill term for objective
    163     HessElemFill (i, 0, deplist [0], p -> Obj (0) -> Body (), row, lam, eee, p);
     169    HessElemFill (i, 0, deplist [0], p -> Obj (0) -> Body (), rnz, lam, eee, p);
    164170
    165171    level++;
     
    169175      CouenneConstraint *c = p -> Con (j);
    170176      if (c -> Body () -> Type () == AUX ||
    171           c -> Body () -> Type () == VAR) continue;
    172 
    173       HessElemFill (i, level, deplist [level], c -> Body (), row, lam, eee, p);
     177          c -> Body () -> Type () == VAR)
     178        continue;
     179
     180      HessElemFill (i, level, deplist [level], c -> Body (), rnz, lam, eee, p);
    174181
    175182      level++;
     
    183190
    184191      if ((e -> Type         () != AUX) ||
    185           (e -> Multiplicity () <= 0)) continue;
    186 
    187       HessElemFill (i, level, deplist [level], e -> Image (), row, lam, eee, p);
    188       //HessElemFill (i, level, deplist [level], e,             row, lam, eee);
    189       // no need to pass e itself, as its second derivative is zero
     192          (e -> Multiplicity () <= 0))
     193        continue;
     194
     195      HessElemFill (i, level, deplist [level], e -> Image (), rnz, lam, eee, p);
    190196
    191197      level++;
    192198    }
    193199
    194     // sparsify row, eee, lam
     200    // sparsify rnz, eee, lam
    195201
    196202    for (int j=0; j <= i; j++)
    197       if (row [j]) {
     203      if (rnz [j]) {
    198204
    199205        reAlloc (nnz_, curSize, iRow_, jCol_, numL_, lamI_, expr_);
     
    201207        iRow_ [nnz_] = i;
    202208        jCol_ [nnz_] = j;
    203         numL_ [nnz_] = row [j];
    204         lamI_ [nnz_] = (int         *) realloc (lam [j], row [j] * sizeof (int));
    205         expr_ [nnz_] = (expression **) realloc (eee [j], row [j] * sizeof (expression *));
     209        numL_ [nnz_] = rnz [j];
     210        lamI_ [nnz_] = (int         *) realloc (lam [j], rnz [j] * sizeof (int));
     211        expr_ [nnz_] = (expression **) realloc (eee [j], rnz [j] * sizeof (expression *));
    206212
    207213        ++nnz_;
    208214      }
    209215
    210     /// upper triangular matrix has to be discarded
    211     for (int j=i+1; j < nVars; j++)
    212 
    213       if (row [j]) { // lam and eee are non nil too
    214 
    215         free (lam [j]);
    216 
    217         for (int k=0; k<row[j]; k++)
    218           if (eee [j] [k])
    219             delete eee [j] [k];
    220 
    221         free (eee [j]);
    222       }
    223 
    224     free (row);
     216    // no it doesn't -- not filled at all
     217    // /// upper triangular matrix has to be discarded
     218    // for (int j=i+1; j < nVars; j++)
     219    //   if (rnz [j]) { // lam and eee are non NULL too
     220    //  assert (false);
     221    //  free (lam [j]);
     222    //  for (int k=0; k<rnz[j]; k++)
     223    //    if (eee [j] [k])
     224    //      delete eee [j] [k];
     225    //  free (eee [j]);
     226    //   }
     227
     228    free (rnz);
    225229    free (lam);
    226230    free (eee);
     
    252256#define reallocStep 100
    253257
     258// fills a row (?) of the Hessian using expression
     259
    254260void HessElemFill (int i, int level,
    255261                   std::set <int> &list,
    256262                   expression *expr,
    257                    int *row,
     263                   int *nnz,
    258264                   int **lam,
    259265                   expression ***eee,
    260266                   CouenneProblem *p) {
    261267
    262   if (list. find (i) == list.end ())
     268  if (list. find (i) == list.end () ||   // expression does not depend on x_i
     269      (expr -> Linearity () <= LINEAR))  // no second derivative here
    263270    return;
    264271
     272  // only fills lower triangular part (for symmetry)
    265273  for (int k=0; k <= i; k++)
    266274    if ((k==i) || (list. find (k) != list.end ())) {
    267275
    268       /// common case: expr is linear, just skip
    269 
    270       if (expr -> Linearity () <= LINEAR)
    271         continue;
    272 
    273276      // objective depends on k and i. Is its second derivative, w.r.t. k and i, nonzero?
    274277
    275278      expression
    276         *d1  = expr -> differentiate (k),
    277         *sd1 = d1 -> simplify (),
    278         *rd1 = (sd1 ? sd1 : d1),
    279         *d2  = rd1 -> differentiate (i),
    280         *sd2 = d2 -> simplify (),
    281         *rd2 = (sd2 ? sd2 : d2);
    282 
     279        *d1  = expr -> differentiate (k), // derivative w.r.t. k
     280        *sd1 = d1  -> simplify (),        // simplified
     281        *rd1 = (sd1 ? sd1 : d1),          // actually used
     282        *d2  = rd1 -> differentiate (i),  // its derivative w.r.t. i
     283        *sd2 = d2  -> simplify (),        // simplified
     284        *rd2 = (sd2 ? sd2 : d2);          // actually used
    283285
    284286#ifdef DEBUG
    285       printf (" rd2 [x_%d, x_%d]: ", k, i); fflush (stdout);
    286       rd2 -> print (); printf (" -> ");
     287      printf (" rd2 [x_%d, x_%d]: d/d x_%d = ", k, i, k); fflush (stdout);
     288      rd1 -> print (); printf (" -> d/(d x_%d,d x_%d) = ", k, i);
    287289      rd2 -> print (); printf ("\n");
    288290#endif
     291
    289292      delete d1;
    290293      if (sd1) delete sd1;
     
    296299        // we have a nonzero
    297300
    298         if (!(row [k] % reallocStep)) {
    299 
    300           lam [k] = (int         *) realloc (lam[k], (row[k] + reallocStep) * sizeof (int));
    301           eee [k] = (expression **) realloc (eee[k], (row[k] + reallocStep) * sizeof (expression *));
     301        if (!(nnz [k] % reallocStep)) {
     302
     303          lam [k] = (int         *) realloc (lam [k], (nnz [k] + reallocStep) * sizeof (int));
     304          eee [k] = (expression **) realloc (eee [k], (nnz [k] + reallocStep) * sizeof (expression *));
    302305        }
    303306
    304307        rd2 -> realign (p); // fixes variables' domain with the problem.
    305308
    306         lam [k] [row [k]] = level;
    307         eee [k] [row [k]] = rd2;
    308         row [k] ++;
     309        lam [k] [nnz [k]]   = level;
     310        eee [k] [nnz [k]++] = rd2;
    309311
    310312      } else delete rd2;
     
    325327  }
    326328}
    327 
  • trunk/Couenne/src/expression/partial/CouenneExprHess.hpp

    r610 r615  
    2424
    2525    int   nnz_;  ///< number of (symbolic) nonzeroes
    26     int  *iRow_; ///< col indices
    27     int  *jCol_; ///< row starts
     26    int  *iRow_; ///< row indices (read this way by eval_h)
     27    int  *jCol_; ///< col indices
    2828
    29     /// There are m+1 (m constraints + 1 obj).
     29    /// There are m+1 (m constraints + 1 obj) components:
     30    ///
     31    /// \f$ \nabla^2 \mathcal L (x,\lambda) = \nabla^2 f(x) + \lambda^\top \nabla^2 g(x) \f$
    3032    ///
    3133    /// Implementing a FP requires adding one for gg', the gradient
  • trunk/Couenne/src/expression/partial/CouenneExprJac.cpp

    r613 r615  
    1717using namespace Couenne;
    1818
    19 //#define DEBUG
     19#define DEBUG
    2020
    2121ExprJac::ExprJac ():
     
    9595      continue;
    9696
    97     // this is a constraint of the form f(x) <= 0. Find out the
     97    // This is a constraint of the form f(x) <= 0. Find out the
    9898    // variables (original or aux) it depends on, directly
    9999    // (STOP_AT_AUX)
  • trunk/Couenne/src/expression/partial/CouenneExprJac.hpp

    r610 r615  
    2323  private:
    2424
    25     int          nnz_;  ///< number of (symbolic) nonzeroes
    26     int         *iRow_; ///< col indices
    27     int         *jCol_; ///< row starts
    28     expression **expr_; ///< expression
     25    int          nnz_;   ///< number of (symbolic) nonzeroes
     26    int         *iRow_;  ///< row indices (read this way by eval_jac_g)
     27    int         *jCol_;  ///< col indices
     28
     29    expression **expr_;  ///< nonzero expression elements (there are nnz_ of them)
    2930
    3031    int          nRows_; ///< number of actual constraints
  • trunk/Couenne/src/heuristics/CouenneFPSolveMILP.cpp

    r614 r615  
    277277    delete [] nlpSolExp;
    278278
    279     for (int i=0, j=problem_ -> nVars (), k = problem_ -> nVars (); k--; i++)
     279    for (int i = 0, j = problem_ -> nVars (), k = j; k--; i++)
    280280
    281281      if (!compDistInt_ || milp_ -> isInteger (i)) {
    282282
     283        // create vector with single entry of 1 at i-th position
    283284        double val = 1.;
    284285        CoinPackedVector vec (1, &i, val);
     
    287288
    288289          // reserved for non-UnitMatrix hessian (i.e. betaMILP_ > 0)
    289 
    290290        }
    291291
     
    310310
    311311    findSolution ();
    312     iSol = CoinCopyOfArray (milp_ -> getColSolution (),
     312    iSol = CoinCopyOfArray (milp_    -> getColSolution (),
    313313                            problem_ -> nVars ());
    314314
     
    317317    int
    318318      nDeleted = nFinalRows - nInitRows,
    319       *deleted  = new int [nDeleted],
     319     *deleted  = new int [nDeleted],
    320320      nCurRow  = nInitRows;
    321321
  • trunk/Couenne/src/heuristics/CouenneFeasPump.cpp

    r614 r615  
    3131int CouenneFeasPump::solution (double &objVal, double *newSolution) {
    3232
    33   if (problem_ -> nIntVars () <= 0)
     33  if (problem_ -> nIntVars () <= 0 ||
     34      CoinCpuTime () > problem_ -> getMaxCpuTime ())
    3435    return 0;
    3536
    3637  printf ("FP ====================================\n");
    37 
    38   if (CoinCpuTime () > problem_ -> getMaxCpuTime ())
    39     return 0;
    4038
    4139  // This FP works as follows:
     
    360358  milp_ = NULL;
    361359
    362   delete [] iSol;
    363 
    364360  return retval;
    365361}
  • trunk/Couenne/src/heuristics/CouenneFeasPumpConstructors.cpp

    r613 r615  
    3333
    3434  ApplicationReturnStatus status = app_ -> Initialize ();
    35   app_ -> Options () -> SetIntegerValue ("max_iter", 40);
     35  app_ -> Options () -> SetIntegerValue ("max_iter", 50);
    3636  if (status != Solve_Succeeded)
    3737    printf ("FP: Error in initialization\n");
  • trunk/Couenne/src/interfaces/CouenneTNLP.cpp

    r614 r615  
    2222#include "CoinHelperFunctions.hpp"
    2323#include "CoinFinite.hpp"
     24
     25#define DEBUG
    2426
    2527using namespace Couenne;
     
    193195      continue;
    194196
    195     // prevent ipopt from exiting on inconsistent bounds
    196     if (lb <= ub) {*g_l++ = lb; *g_u++ = ub;}
    197     else          {*g_l++ = ub; *g_u++ = lb;}
     197    *g_l = (e -> sign () != expression::AUX_GEQ) ? 0. : -COIN_DBL_MAX;
     198    *g_u = (e -> sign () != expression::AUX_LEQ) ? 0. :  COIN_DBL_MAX;
     199
     200    ++g_l;
     201    ++g_u;
    198202  }
    199203
     
    366370      continue;
    367371
    368     *g++ = (*e) ();
     372    *g++ = (*(e -> Image ())) () - (*e) ();
    369373
    370374    nEntries ++;
     
    404408#endif
    405409
    406   //problem_ -> domain () -> push (n, x, NULL, NULL);
    407 
    408410  if (values == NULL &&
    409411      iRow   != NULL &&
     
    423425    register expression **e = Jac_. expr ();
    424426
    425     for (register int i=0; i<nele_jac; i++)
     427    for (register int i=nele_jac; i--;)
    426428      *values++ = (**(e++)) ();
    427429  }
     
    434436  } else printf ("empty\n");
    435437#endif
    436 
    437   //problem_ -> domain () -> pop ();
    438438
    439439  return true;
Note: See TracChangeset for help on using the changeset viewer.