added (untested) interface for specifying preconditioners

Ignore-this: 246b0790dff127b54fbf72d34d2da402

darcs-hash:20111115221929-c8de0-a27654751310b4251bdf6961ca7f1f3fc13164f4.gz

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) */

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,

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;
      }
 

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)						\

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;
 }

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 */

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;