added (untested) interface for specifying preconditioners

Ignore-this: 246b0790dff127b54fbf72d34d2da402

darcs-hash:20111115221929-c8de0-a27654751310b4251bdf6961ca7f1f3fc13164f4.gz

diff --git a/api/nlopt-internal.h b/api/nlopt-internal.h
index 3b621d1393b242d160c1b0d3aaad5a722a9fbdd4..d41099fa01d134d2728a3951087e0a3dcedc2f6a 100644 (file)
--- a/api/nlopt-internal.h
+++ b/api/nlopt-internal.h
@@ -38,6 +38,7 @@ struct nlopt_opt_s {
      unsigned n; /* the dimension of the problem (immutable) */
 
      nlopt_func f; void *f_data; /* objective function to minimize */
+     nlopt_precond pre; /* optional preconditioner for f (NULL if none) */
      int maximize; /* nonzero if we are maximizing, not minimizing */
 
      double *lb, *ub; /* lower and upper bounds (length n) */

diff --git a/api/nlopt.h b/api/nlopt.h
index 2ab6ab6534335e487325cc4bc9709be96d3eee90..46584e9a745aaf3873bd1d84fc6cf7698498d9fa 100644 (file)
--- a/api/nlopt.h
+++ b/api/nlopt.h
@@ -68,10 +68,10 @@ typedef void (*nlopt_mfunc)(unsigned m, double *result,
                             double *gradient, /* NULL if not needed */
                             void *func_data);
 
-/* a preconditioner, which computes the approximate Hessian H(x) at x.
-   In particular, it returns Hdx = H(x) * dx.    [Array lengths = n.] */
-typedef void (*nlopt_precond)(unsigned n, const double *x, const double *dx,
-                             double *Hdx, void *data);
+/* A preconditioner, which preconditions v at x to return vpre. 
+   (The meaning of "preconditioning" is algorithm-dependent.) */
+typedef void (*nlopt_precond)(unsigned n, const double *x, const double *v,
+                             double *vpre, void *data);
 
 typedef enum {
      /* Naming conventions:
@@ -202,6 +202,9 @@ NLOPT_EXTERN(nlopt_result) nlopt_set_min_objective(nlopt_opt opt, nlopt_func f,
 NLOPT_EXTERN(nlopt_result) nlopt_set_max_objective(nlopt_opt opt, nlopt_func f, 
                                                  void *f_data);
 
+NLOPT_EXTERN(nlopt_result) nlopt_set_precond_min_objective(nlopt_opt opt, nlopt_func f, nlopt_precond pre, void *f_data);
+NLOPT_EXTERN(nlopt_result) nlopt_set_precond_max_objective(nlopt_opt opt, nlopt_func f, nlopt_precond pre, void *f_data);
+
 NLOPT_EXTERN(nlopt_algorithm) nlopt_get_algorithm(const nlopt_opt opt);
 NLOPT_EXTERN(unsigned) nlopt_get_dimension(const nlopt_opt opt);
 
@@ -223,6 +226,9 @@ NLOPT_EXTERN(nlopt_result) nlopt_add_inequality_constraint(nlopt_opt opt,
                                                          nlopt_func fc,
                                                          void *fc_data,
                                                          double tol);
+NLOPT_EXTERN(nlopt_result) nlopt_add_precond_inequality_constraint(
+     nlopt_opt opt, nlopt_func fc, nlopt_precond pre, void *fc_data,
+     double tol);
 NLOPT_EXTERN(nlopt_result) nlopt_add_inequality_mconstraint(nlopt_opt opt,
                                                            unsigned m,
                                                            nlopt_mfunc fc,
@@ -234,6 +240,9 @@ NLOPT_EXTERN(nlopt_result) nlopt_add_equality_constraint(nlopt_opt opt,
                                                        nlopt_func h,
                                                        void *h_data,
                                                        double tol);
+NLOPT_EXTERN(nlopt_result) nlopt_add_precond_equality_constraint(
+     nlopt_opt opt, nlopt_func h, nlopt_precond pre, void *h_data,
+     double tol);
 NLOPT_EXTERN(nlopt_result) nlopt_add_equality_mconstraint(nlopt_opt opt,
                                                          unsigned m,
                                                          nlopt_mfunc h,

diff --git a/api/optimize.c b/api/optimize.c
index 236519fb8ea8a7ed225a23a8279f747e6c55adf4..75aa48750eee32acd5c31ceca58bffc85513797c 100644 (file)
--- a/api/optimize.c
+++ b/api/optimize.c
@@ -661,8 +661,7 @@ static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf)
                                      lb, ub, x, minf, &stop, dual_opt);
              else
                   ret = ccsa_quadratic_minimize(
-                       n, f, f_data, opt->m, opt->fc,
-                       NULL, NULL, NULL, NULL,
+                       n, f, f_data, opt->m, opt->fc, opt->pre,
                        lb, ub, x, minf, &stop, dual_opt);
              nlopt_destroy(dual_opt);
              return ret;
@@ -797,6 +796,7 @@ static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf)
 
 typedef struct {
      nlopt_func f;
+     nlopt_precond pre;
      void *f_data;
 } f_max_data;
 
@@ -813,22 +813,31 @@ static double f_max(unsigned n, const double *x, double *grad, void *data)
      return -val;
 }
 
+static void pre_max(unsigned n, const double *x, const double *v,
+                   double *vpre, void *data)
+{
+     f_max_data *d = (f_max_data *) data;
+     unsigned i;
+     d->pre(n, x, v, vpre, data);
+     for (i = 0; i < n; ++i) vpre[i] = -vpre[i];
+}
+
 nlopt_result 
 NLOPT_STDCALL nlopt_optimize(nlopt_opt opt, double *x, double *opt_f)
 {
-     nlopt_func f; void *f_data;
+     nlopt_func f; void *f_data; nlopt_precond pre;
      f_max_data fmd;
      int maximize;
      nlopt_result ret;
 
      if (!opt || !opt_f || !opt->f) return NLOPT_INVALID_ARGS;
-     f = opt->f; f_data = opt->f_data;
+     f = opt->f; f_data = opt->f_data; pre = opt->pre;
 
      /* for maximizing, just minimize the f_max wrapper, which 
        flips the sign of everything */
      if ((maximize = opt->maximize)) {
-         fmd.f = f; fmd.f_data = f_data;
-         opt->f = f_max; opt->f_data = &fmd;
+         fmd.f = f; fmd.f_data = f_data; fmd.pre = pre;
+         opt->f = f_max; opt->f_data = &fmd; opt->pre = pre_max;
          opt->stopval = -opt->stopval;
          opt->maximize = 0;
      }
@@ -853,7 +862,7 @@ done:
      if (maximize) { /* restore original signs */
          opt->maximize = maximize;
          opt->stopval = -opt->stopval;
-         opt->f = f; opt->f_data = f_data;
+         opt->f = f; opt->f_data = f_data; opt->pre = pre;
          *opt_f = -*opt_f;
      }
 

diff --git a/api/options.c b/api/options.c
index 2e7ae0367a32317397e80a49e6e9f50ce407b9fe..f3a432030709a13a8bda174a6691c03a8e1ddf47 100644 (file)
--- a/api/options.c
+++ b/api/options.c
@@ -68,7 +68,7 @@ nlopt_opt NLOPT_STDCALL nlopt_create(nlopt_algorithm algorithm, unsigned n)
      if (opt) {
          opt->algorithm = algorithm;
          opt->n = n;
-         opt->f = NULL; opt->f_data = NULL;
+         opt->f = NULL; opt->f_data = NULL; opt->pre = NULL;
          opt->maximize = 0;
          opt->munge_on_destroy = opt->munge_on_copy = NULL;
 
@@ -212,12 +212,14 @@ oom:
 
 /*************************************************************************/
 
-nlopt_result NLOPT_STDCALL nlopt_set_min_objective(nlopt_opt opt,
-                                                  nlopt_func f, void *f_data)
+nlopt_result NLOPT_STDCALL nlopt_set_precond_min_objective(nlopt_opt opt,
+                                                          nlopt_func f, 
+                                                          nlopt_precond pre,
+                                                          void *f_data)
 {
      if (opt) {
          if (opt->munge_on_destroy) opt->munge_on_destroy(opt->f_data);
-         opt->f = f; opt->f_data = f_data;
+         opt->f = f; opt->f_data = f_data; opt->pre = pre;
          opt->maximize = 0;
          if (nlopt_isinf(opt->stopval) && opt->stopval > 0)
               opt->stopval = -HUGE_VAL; /* switch default from max to min */
@@ -226,12 +228,20 @@ nlopt_result NLOPT_STDCALL nlopt_set_min_objective(nlopt_opt opt,
      return NLOPT_INVALID_ARGS;
 }
 
-nlopt_result NLOPT_STDCALL nlopt_set_max_objective(nlopt_opt opt, 
+nlopt_result NLOPT_STDCALL nlopt_set_min_objective(nlopt_opt opt,
                                                   nlopt_func f, void *f_data)
+{
+     nlopt_set_precond_min_objective(opt, f, NULL, f_data);
+}
+
+nlopt_result NLOPT_STDCALL nlopt_set_precond_max_objective(nlopt_opt opt, 
+                                                          nlopt_func f, 
+                                                          nlopt_precond pre,
+                                                          void *f_data)
 {
      if (opt) {
          if (opt->munge_on_destroy) opt->munge_on_destroy(opt->f_data);
-         opt->f = f; opt->f_data = f_data;
+         opt->f = f; opt->f_data = f_data; opt->pre = pre;
          opt->maximize = 1;
          if (nlopt_isinf(opt->stopval) && opt->stopval < 0)
               opt->stopval = +HUGE_VAL; /* switch default from min to max */
@@ -240,6 +250,12 @@ nlopt_result NLOPT_STDCALL nlopt_set_max_objective(nlopt_opt opt,
      return NLOPT_INVALID_ARGS;
 }
 
+nlopt_result NLOPT_STDCALL nlopt_set_max_objective(nlopt_opt opt,
+                                                  nlopt_func f, void *f_data)
+{
+     nlopt_set_precond_max_objective(opt, f, NULL, f_data);
+}
+
 /*************************************************************************/
 
 nlopt_result
@@ -336,6 +352,7 @@ NLOPT_STDCALL nlopt_remove_inequality_constraints(nlopt_opt opt)
 static nlopt_result add_constraint(unsigned *m, unsigned *m_alloc,
                                   nlopt_constraint **c,
                                   unsigned fm, nlopt_func fc, nlopt_mfunc mfc,
+                                  nlopt_precond pre,
                                   void *fc_data,
                                   const double *tol)
 {
@@ -371,6 +388,7 @@ static nlopt_result add_constraint(unsigned *m, unsigned *m_alloc,
      
      (*c)[*m - 1].m = fm;
      (*c)[*m - 1].f = fc;
+     (*c)[*m - 1].pre = pre;
      (*c)[*m - 1].mf = mfc;
      (*c)[*m - 1].f_data = fc_data;
      (*c)[*m - 1].tol = tolcopy;
@@ -400,26 +418,37 @@ NLOPT_STDCALL nlopt_add_inequality_mconstraint(nlopt_opt opt, unsigned m,
      }
      if (!opt || !inequality_ok(opt->algorithm)) ret = NLOPT_INVALID_ARGS;
      else ret = add_constraint(&opt->m, &opt->m_alloc, &opt->fc,
-                              m, NULL, fc, fc_data, tol);
+                              m, NULL, fc, NULL, fc_data, tol);
      if (ret < 0 && opt && opt->munge_on_destroy)
          opt->munge_on_destroy(fc_data);
      return ret;
 }
 
 nlopt_result
-NLOPT_STDCALL nlopt_add_inequality_constraint(nlopt_opt opt,
-                                             nlopt_func fc, void *fc_data,
-                                             double tol)
+NLOPT_STDCALL nlopt_add_precond_inequality_constraint(nlopt_opt opt,
+                                                     nlopt_func fc, 
+                                                     nlopt_precond pre,
+                                                     void *fc_data,
+                                                     double tol)
 {
      nlopt_result ret;
      if (!opt || !inequality_ok(opt->algorithm)) ret = NLOPT_INVALID_ARGS;
      else ret = add_constraint(&opt->m, &opt->m_alloc, &opt->fc,
-                              1, fc, NULL, fc_data, &tol);
+                              1, fc, NULL, pre, fc_data, &tol);
      if (ret < 0 && opt && opt->munge_on_destroy)
          opt->munge_on_destroy(fc_data);
      return ret;
 }
 
+nlopt_result
+NLOPT_STDCALL nlopt_add_inequality_constraint(nlopt_opt opt,
+                                             nlopt_func fc, void *fc_data,
+                                             double tol)
+{
+     return nlopt_add_precond_inequality_constraint(opt, fc, NULL, fc_data,
+                                                   tol);
+}
+
 nlopt_result
 NLOPT_STDCALL nlopt_remove_equality_constraints(nlopt_opt opt)
 {
@@ -460,28 +489,38 @@ NLOPT_STDCALL nlopt_add_equality_mconstraint(nlopt_opt opt, unsigned m,
         || nlopt_count_constraints(opt->p, opt->h) + m > opt->n) 
          ret = NLOPT_INVALID_ARGS;
      else ret = add_constraint(&opt->p, &opt->p_alloc, &opt->h,
-                              m, NULL, fc, fc_data, tol);
+                              m, NULL, fc, NULL, fc_data, tol);
      if (ret < 0 && opt && opt->munge_on_destroy)
          opt->munge_on_destroy(fc_data);
      return ret;
 }
 
 nlopt_result
-NLOPT_STDCALL nlopt_add_equality_constraint(nlopt_opt opt,
-                                           nlopt_func fc, void *fc_data,
-                                           double tol)
+NLOPT_STDCALL nlopt_add_precond_equality_constraint(nlopt_opt opt,
+                                                   nlopt_func fc,
+                                                   nlopt_precond pre,
+                                                   void *fc_data,
+                                                   double tol)
 {
      nlopt_result ret;
      if (!opt || !equality_ok(opt->algorithm)
         || nlopt_count_constraints(opt->p, opt->h) + 1 > opt->n)
          ret = NLOPT_INVALID_ARGS;
      else ret = add_constraint(&opt->p, &opt->p_alloc, &opt->h,
-                              1, fc, NULL, fc_data, &tol);
+                              1, fc, NULL, pre, fc_data, &tol);
      if (ret < 0 && opt && opt->munge_on_destroy)
          opt->munge_on_destroy(fc_data);
      return ret;
 }
 
+nlopt_result
+NLOPT_STDCALL nlopt_add_equality_constraint(nlopt_opt opt,
+                                           nlopt_func fc, void *fc_data,
+                                           double tol)
+{
+     return nlopt_add_precond_equality_constraint(opt, fc, NULL, fc_data, tol);
+}
+
 /*************************************************************************/
 
 #define SET(param, T, arg)                                             \

diff --git a/mma/ccsa_quadratic.c b/mma/ccsa_quadratic.c
index dbca8fad8005e3381bd5464915d756f771c62a42..eb3715b2248bcacb902423b519748872cb9268a2 100644 (file)
--- a/mma/ccsa_quadratic.c
+++ b/mma/ccsa_quadratic.c
@@ -212,8 +212,7 @@ nlopt_result ccsa_quadratic_minimize(
      unsigned n, nlopt_func f, void *f_data,
      unsigned m, nlopt_constraint *fc,
 
-     nlopt_precond pre, void *pre_data, /* for f */
-     nlopt_precond *prec, void **prec_data, /* length = # constraints */
+     nlopt_precond pre, 
 
      const double *lb, const double *ub, /* bounds */
      double *x, /* in: initial guess, out: minimizer */
@@ -264,14 +263,30 @@ nlopt_result ccsa_quadratic_minimize(
      dd.rhoc = rhoc;
      dd.xcur = xcur;
      dd.gcval = gcval;
-     dd.pre = pre; dd.pre_data = pre_data;
-     dd.prec = prec; dd.prec_data = prec_data;
+     dd.pre = pre; dd.pre_data = f_data;
+     dd.prec = NULL; dd.prec_data = NULL;
      dd.scratch = NULL;
 
+     if (m) {
+         dd.prec = (nlopt_precond *) malloc(sizeof(nlopt_precond) * m);
+         dd.prec_data = (void **) malloc(sizeof(void *) * m);
+         if (!dd.prec || !dd.prec_data) {
+              ret = NLOPT_OUT_OF_MEMORY;
+              goto done;
+         }
+         for (i = ifc = 0; ifc < mfc; ++ifc) {
+              unsigned inext = i + fc[ifc].m;
+              for (; i < inext; ++i) {
+                   dd.prec[i] = fc[ifc].pre;
+                   dd.prec_data[i] = fc[ifc].f_data;
+              }
+         }
+     }
+
      no_precond = pre == NULL;
-     if (prec)
+     if (dd.prec)
          for (i = 0; i < m; ++i)
-              no_precond = no_precond && prec[i] == NULL;
+              no_precond = no_precond && dd.prec[i] == NULL;
 
      if (!no_precond) {
          dd.scratch = (double*) malloc(sizeof(double) * (2*n));
@@ -430,15 +445,14 @@ nlopt_result ccsa_quadratic_minimize(
               }
               for (i = ifc = 0; ifc < mfc; ++ifc) {
                    unsigned i0 = i, inext = i + fc[ifc].m;
-                   for (; i < inext; ++i)
-                        if (!isnan(fcval_cur[i])) {
-                             feasible_cur = feasible_cur 
-                                  && fcval_cur[i] <= fc[ifc].tol[i-i0];
-                             inner_done = inner_done && 
-                                  (dd.gcval[i] >= fcval_cur[i]);
-                             if (fcval_cur[i] > infeasibility_cur)
-                                  infeasibility_cur = fcval_cur[i];
-                        }
+                   for (; i < inext; ++i) {
+                        feasible_cur = feasible_cur 
+                             && fcval_cur[i] <= fc[ifc].tol[i-i0];
+                        inner_done = inner_done && 
+                             (dd.gcval[i] >= fcval_cur[i]);
+                        if (fcval_cur[i] > infeasibility_cur)
+                             infeasibility_cur = fcval_cur[i];
+                   }
               }
 
               if ((fcur < *minf && (inner_done || feasible_cur || !feasible))
@@ -525,6 +539,10 @@ nlopt_result ccsa_quadratic_minimize(
  done:
      nlopt_destroy(pre_opt);
      if (dd.scratch) free(dd.scratch);
+     if (m) {
+         free(dd.prec_data);
+         free(dd.prec);
+     }
      free(sigma);
      return ret;
 }

diff --git a/mma/mma.h b/mma/mma.h
index a08b043bd3954a4906c3934c3cc002f288091340..7968dc86bc99a41930dad7194986bba2db434e74 100644 (file)
--- a/mma/mma.h
+++ b/mma/mma.h
@@ -45,8 +45,7 @@ nlopt_result ccsa_quadratic_minimize(
      unsigned n, nlopt_func f, void *f_data,
      unsigned m, nlopt_constraint *fc,
 
-     nlopt_precond pre, void *pre_data, /* for f */
-     nlopt_precond *prec, void **prec_data, /* length = # constraints */
+     nlopt_precond pre,
 
      const double *lb, const double *ub, /* bounds */
      double *x, /* in: initial guess, out: minimizer */

diff --git a/util/nlopt-util.h b/util/nlopt-util.h
index ddbb801ff1ee44dce0baead559c894c020d40e30..ca33df4b9fb3831f9aa17f52936e86c772422728 100644 (file)
--- a/util/nlopt-util.h
+++ b/util/nlopt-util.h
@@ -111,6 +111,7 @@ typedef struct {
      unsigned m; /* dimensional of constraint: mf maps R^n -> R^m */
      nlopt_func f; /* one-dimensional constraint, requires m == 1 */
      nlopt_mfunc mf;
+     nlopt_precond pre; /* preconditioner for f (NULL if none or if mf) */
      void *f_data;
      double *tol;
 } nlopt_constraint;