From 5aafab6619c03a20e654753c817bda7b0706bcdf Mon Sep 17 00:00:00 2001
From: liam <liam@a3d8a0fb-c1db-0310-ace7-a616afeb9e30>
Date: Wed, 16 Jan 2008 03:57:28 +0000
Subject: [PATCH] Multivariate root solving with derivative now functions. 
 Symbol macro 'success returns the GSL_SUCCESS code, required for many
 callback definitions.  Introduced macros vref and mref (setfable) to get
 components of GSL vector doubles and matrix doubles directly from the
 pointers.

git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3265 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
---
 data/data.lisp                                |  20 +--
 data/matrix.lisp                              |  32 +++-
 data/vector.lisp                              |  39 +++--
 general/conditions.lisp                       |  19 ++-
 general/functions.lisp                        |  14 +-
 index.html                                    |   4 +-
 init/interface.lisp                           | 144 +++++++++---------
 .../ode-system.lisp                           |  16 +-
 roots-multi.lisp                              |  81 +++++++++-
 roots-one.lisp                                |  67 ++++----
 10 files changed, 267 insertions(+), 169 deletions(-)

diff --git a/data/data.lisp b/data/data.lisp
index 8db3fc8d..5a7be85e 100644
--- a/data/data.lisp
+++ b/data/data.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        data.lisp                                 
-; description: Using GSL storage.                        
-; date:        Sun Mar 26 2006 - 16:32                   
-; author:      Liam M. Healy                             
-; modified:    Sun Jan 13 2008 - 22:49
-;********************************************************
-;;; $Id: $
+;; Using GSL storage.
+;; Liam Healy, Sun Mar 26 2006 - 16:32
+;; Time-stamp: <2008-01-14 21:09:15 liam data.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -234,14 +230,6 @@
 	 (progn ,@body)
       (free ,symbol))))
 
-(defun make-data-from-pointer (pointer &optional (class 'gsl-vector-double) size)
-  "Given a C pointer to a GSL data type, make the CL object."
-  (make-instance
-   class
-   :pointer pointer
-   :storage-size
-   (or size (cffi:foreign-slot-value pointer 'gsl-vector-c 'size))))
-
 ;;;;****************************************************************************
 ;;;; Getting values into CL
 ;;;;****************************************************************************
diff --git a/data/matrix.lisp b/data/matrix.lisp
index e9c516b0..c886e790 100644
--- a/data/matrix.lisp
+++ b/data/matrix.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        matrix.lisp                        
-; description: Matrices
-; date:        Sun Mar 26 2006 - 11:51                   
-; author:      Liam M. Healy                             
-; modified:    Sun Dec  3 2006 - 18:53
-;********************************************************
-;;; $Id: $
+;; Matrices
+;; Liam Healy, Sun Mar 26 2006 - 11:51
+;; Time-stamp: <2008-01-15 19:34:56 liam matrix.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -51,6 +47,12 @@
 (defmethod gsl-array ((object gsl-matrix))
   (foreign-slot-value (pointer object) 'gsl-matrix-c 'data))
 
+(export 'matrix-data)
+(defun matrix-data (pointer)
+  "A pointer to the GSL array with the data contents, from the
+   sruct pointer."
+  (cffi:foreign-slot-value pointer 'gsl-matrix-c 'data))
+
 ;;;;****************************************************************************
 ;;;; Getting values
 ;;;;****************************************************************************
@@ -63,6 +65,13 @@
   :c-return :c-base-type
   :documentation "The (i,j)-th element of the matrix.")
 
+(defun-gsl mref (pointer index0 index1)
+  "gsl_matrix_get"
+  ((pointer :pointer) (index0 :size) (index1 :size))
+  :c-return :double
+  :index nil
+  :documentation "An element of the matrix of doubles, computed from the pointer.")
+
 (defun-gsl-mdsfc gsl-matrix-ptr ((matrix gsl-matrix) i j)
   "gsl_matrix_ptr" (((pointer matrix) :pointer) (i :size) (j :size))
   :c-return :pointer
@@ -96,6 +105,13 @@
   :c-return :void
   :documentation "Set the (i,j)-th element of the matrix.")
 
+(defun-gsl (setf mref) (value pointer index0 index1)
+  "gsl_matrix_set"
+  ((pointer :pointer) (index0 :size) (index1 :size) (value :double))
+  :c-return :void
+  :index nil
+  :documentation "Set an element of the matrix of doubles, using its pointer.")
+
 (defmethod (setf data) (array (object gsl-matrix))
   (loop for i from 0
 	below
diff --git a/data/vector.lisp b/data/vector.lisp
index c35867cf..c5159706 100644
--- a/data/vector.lisp
+++ b/data/vector.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        vector.lisp                        
-; description: Vectors
-; date:        Sun Mar 26 2006 - 11:51                   
-; author:      Liam M. Healy                             
-; modified:    Sun Jan 13 2008 - 22:40
-;********************************************************
-;;; $Id: $
+;; Vectors
+;; Liam Healy, Sun Mar 26 2006 - 11:51
+;; Time-stamp: <2008-01-15 19:01:37 liam vector.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -79,7 +75,16 @@ deallocated with the vector.
 (defmethod gsl-array ((object gsl-vector))
   (cffi:foreign-slot-value (pointer object) 'gsl-vector-c 'data))
 
-(defun gsl-array-p (pointer)
+(defun make-data-from-pointer (pointer &optional (class 'gsl-vector-double) size)
+  "Given a C pointer to a GSL data type, make the CL object."
+  (make-instance
+   class
+   :pointer pointer
+   :storage-size
+   (or size (cffi:foreign-slot-value pointer 'gsl-vector-c 'size))))
+
+(export 'vector-data)
+(defun vector-data (pointer)
   "A pointer to the GSL array with the data contents, from the
    sruct pointer."
   (cffi:foreign-slot-value pointer 'gsl-vector-c 'data))
@@ -94,6 +99,13 @@ deallocated with the vector.
   :c-return :c-base-type
   :documentation "The ith element of the vector.")
 
+(defun-gsl vref (pointer index)
+  "gsl_vector_get"
+  ((pointer :pointer) (index :size))
+  :c-return :double
+  :index nil
+  :documentation "An element of the vector of doubles, computed from the pointer.")
+
 (defun-gsl gsl-vector-ptr (vector i)
   "gsl_vector_ptr" (((pointer vector) :pointer) (i :size))
   :c-return :pointer
@@ -107,7 +119,14 @@ deallocated with the vector.
   "gsl_vector_set"
   (((pointer vector) :pointer) ((first indices) :size) (value :c-base-type))
   :c-return :void
-  :documentation "Set the ith element of the vector.")
+  :documentation "Set an element of the vector.")
+
+(defun-gsl (setf vref) (value pointer index)
+  "gsl_vector_set"
+  ((pointer :pointer) (index :size) (value :double))
+  :c-return :void
+  :index nil
+  :documentation "Set an element of the vector of doubles, using its pointer.")
 
 (defun-gsl-vdsfc set-all ((object gsl-vector) value)
   "gsl_vector_set_all"
diff --git a/general/conditions.lisp b/general/conditions.lisp
index 8373a1d0..3d1e0bcb 100644
--- a/general/conditions.lisp
+++ b/general/conditions.lisp
@@ -1,10 +1,7 @@
-;********************************************************
-; file:        conditions.lisp
-; description: GSL errors                                
-; date:        Sat Mar  4 2006 - 18:33                   
-; author:      Liam M. Healy
-; modified:    Sun Apr  2 2006 - 15:56
-;********************************************************
+;; GSL errors                                
+;; Liam Healy Sat Mar  4 2006 - 18:33
+;; Time-stamp: <2008-01-14 22:22:55 liam conditions.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -98,3 +95,11 @@
 (cffi:foreign-funcall
  "gsl_set_error_handler"
  :pointer (cffi:callback gsl-error))
+
+
+(defmacro gsl-errorno-sm (keyword)
+  `(define-symbol-macro
+    ,(intern (symbol-name keyword) (find-package :gsl))
+    (cffi:foreign-enum-value 'gsl-errorno ,keyword)))
+
+(gsl-errorno-sm :success)
diff --git a/general/functions.lisp b/general/functions.lisp
index 991000d4..8c0c259c 100644
--- a/general/functions.lisp
+++ b/general/functions.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        functions.lisp                            
-; description: Foreign callback functions.               
-; date:        Sun Dec  9 2007 - 22:08                   
-; author:      Liam Healy                                
-; modified:    Sun Jan 13 2008 - 13:26
-;********************************************************
-;;; $Id: $
+;; Foreign callback functions.               
+;; Liam Healy 
+;; Time-stamp: <2008-01-14 22:56:45 liam functions.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -90,7 +86,7 @@
 	,@(when (eq return-type :success-failure)
 		;; We always return success, because if there was a
 		;; problem, a CL error would be signalled.
-		'((cffi:foreign-enum-value 'gsl-errorno :SUCCESS))))
+		'(success)))
       ,@(when
 	 structure
 	 ;; Assume that defcallback does not bind the variable 'name.
diff --git a/index.html b/index.html
index c4b886f9..4f927e3c 100644
--- a/index.html
+++ b/index.html
@@ -313,7 +313,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>
 	<tr>
 	<td><a href="Multidimensional-Root_Finding.html">
 	Multidimensional Root Finding</a></td>
-	<td>Started</td>
+	<td>Done</td>
 	</tr>
 	<tr>
 	<td><a href="Multidimensional-Minimization.html">
@@ -357,7 +357,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>
 <!-- Created: Feb 25 2005 -->
 <!-- hhmts start -->
     <small>
-       Time-stamp: <2008-01-13 16:51:23 liam index.html>
+       Time-stamp: <2008-01-15 22:37:31 liam index.html>
        </small>
 <!-- hhmts end -->
  </div>
diff --git a/init/interface.lisp b/init/interface.lisp
index f0bf70e8..b2011ecf 100644
--- a/init/interface.lisp
+++ b/init/interface.lisp
@@ -1,10 +1,7 @@
-;********************************************************
-; file:        interface.lisp                            
-; description: Macros to interface GSL functions.
-; date:        Mon Mar  6 2006 - 22:35                   
-; author:      Liam M. Healy
-; modified:    Sun Jan  6 2008 - 13:48
-;********************************************************
+;; Macros to interface GSL functions.
+;; Liam Healy 
+;; Time-stamp: <2008-01-15 18:55:45 liam interface.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -101,7 +98,7 @@
 
 (defun success-continue (value)
   "If status indicates success, return T, othewise return NIL."
-  (eql value (cffi:foreign-enum-value 'gsl-errorno :SUCCESS)))
+  (eql value success))
 
 ;;;;****************************************************************************
 ;;;; Macro defun-gsl 
@@ -125,7 +122,8 @@
 ;;;            or its derivatives.
 ;;;            Default are allocated quantities in c-arguments, or :c-return if none.
 ;;; type       :function or :method
-;;; index      Name under which this function should be cross-referenced
+;;; index      Name under which this function should be cross-referenced; t
+;;;            to use name, nil to not index.
 ;;; export     Whether to export the symbol.
 ;;; null-pointer-info Return value if C function returns a null pointer.
 ;;; documentation
@@ -136,7 +134,7 @@
     (name arglist gsl-name c-arguments
      &key (c-return :error-code)
      (return nil return-supplied-p)
-     (type :function) index (export (not (eq type :method)))
+     (type :function) (index t) (export (not (eq type :method)))
      null-pointer-info documentation invalidate after enumeration)
   (let* ((cargs (substitute '(mode sf-mode) :mode c-arguments))
 	 (carg-symbs
@@ -164,69 +162,69 @@
 		    (unless (eq c-return :void)
 		      (list cret-name)))))
     `(progn
-       (,(if (eq type :function) 'defun 'defmethod)
-	 ,name
-	 ,(if (member :mode c-arguments)
-	      `(,@clargs &optional (mode :double-prec))
-	      `(,@clargs))
-	 ,@(when documentation (list documentation))
-	 (,@(if allocated
-		`(cffi:with-foreign-objects
-		     ,(mapcar #'wfo-declare allocated-decl))
-		'(let ()))
-	    (let ((,cret-name
-		   (cffi:foreign-funcall
-		    ,gsl-name
-		    ,@(mapcan
-		       (lambda (arg)
-			 (list (if (member (st-symbol arg) allocated)
-				   :pointer
-				   (st-type arg))
-			       (st-symbol arg)))
-		       cargs)
-		    ,cret-type)))
-	      ,@(case c-return
-		      (:void `((declare (ignore ,cret-name))))
-		      (:error-code	; fill in arguments
-		       `((check-gsl-status ,cret-name ',name))))
-	      ,@(when invalidate `((cl-invalidate ,@invalidate)))
-	      ,@(when (or null-pointer-info (eq c-return :pointer))
-		      `((check-null-pointer ,cret-name
-					    ,@(or null-pointer-info
-						  '(:ENOMEM "No memory allocated")))))
-	      ,@after
-	      (values
-	       ,@(case c-return
-		       (:number-of-answers
-			(mapcan
-			 (lambda (vbl seq)
-			   `((when (> ,cret-name ,seq) ,vbl)))
-			 clret
-			 (loop for i below (length clret) collect i)))
-		       (:success-failure
-			(if (equal clret invalidate)
-			    ;; success-failure more important than passed-in
-			    `(success-failure ,@clret)
-			    (remove cret-name ; don't return c-return itself
-				    `(,@clret (success-failure ,cret-name)))))
-		       (:success-continue
-			(if (equal clret invalidate)
-			    ;; success-failure more important than passed-in
-			    `(success-continue ,@clret)
-			    (remove cret-name ; don't return c-return itself
-				    `(,@clret (success-continue ,cret-name)))))
-		       (:true-false
-			`((not (zerop ,cret-name))))
-		       (:enumerate
-			`((cffi:foreign-enum-keyword ',enumeration ,cret-name)))
-		       (t (unless
-			      (or
-			       (and (eq c-return :error-code)
-				    (not allocated))
-			       (and (null return) return-supplied-p))
-			    clret)))))))
-       (map-name ',(or index name) ,gsl-name)
-       ,@(when export `((export ',name))))))
+      (,(if (eq type :function) 'defun 'defmethod)
+       ,name
+       ,(if (member :mode c-arguments)
+	    `(,@clargs &optional (mode :double-prec))
+	    `(,@clargs))
+       ,@(when documentation (list documentation))
+       (,@(if allocated
+	      `(cffi:with-foreign-objects
+		,(mapcar #'wfo-declare allocated-decl))
+	      '(let ()))
+	(let ((,cret-name
+	       (cffi:foreign-funcall
+		,gsl-name
+		,@(mapcan
+		   (lambda (arg)
+		     (list (if (member (st-symbol arg) allocated)
+			       :pointer
+			       (st-type arg))
+			   (st-symbol arg)))
+		   cargs)
+		,cret-type)))
+	  ,@(case c-return
+		  (:void `((declare (ignore ,cret-name))))
+		  (:error-code		; fill in arguments
+		   `((check-gsl-status ,cret-name ',name))))
+	  ,@(when invalidate `((cl-invalidate ,@invalidate)))
+	  ,@(when (or null-pointer-info (eq c-return :pointer))
+		  `((check-null-pointer ,cret-name
+		     ,@(or null-pointer-info
+			   '(:ENOMEM "No memory allocated")))))
+	  ,@after
+	  (values
+	   ,@(case c-return
+		   (:number-of-answers
+		    (mapcan
+		     (lambda (vbl seq)
+		       `((when (> ,cret-name ,seq) ,vbl)))
+		     clret
+		     (loop for i below (length clret) collect i)))
+		   (:success-failure
+		    (if (equal clret invalidate)
+			;; success-failure more important than passed-in
+			`(success-failure ,@clret)
+			(remove cret-name ; don't return c-return itself
+				`(,@clret (success-failure ,cret-name)))))
+		   (:success-continue
+		    (if (equal clret invalidate)
+			;; success-failure more important than passed-in
+			`(success-continue ,@clret)
+			(remove cret-name ; don't return c-return itself
+				`(,@clret (success-continue ,cret-name)))))
+		   (:true-false
+		    `((not (zerop ,cret-name))))
+		   (:enumerate
+		    `((cffi:foreign-enum-keyword ',enumeration ,cret-name)))
+		   (t (unless
+			  (or
+			   (and (eq c-return :error-code)
+				(not allocated))
+			   (and (null return) return-supplied-p))
+			clret)))))))
+      ,@(when index `((map-name ',(if (eql index t) name index) ,gsl-name)))
+      ,@(when export `((export ',name))))))
 
 (defmacro defun-optionals
     (name arglist no-optional optionals &optional documentation)
diff --git a/ordinary-differential-equations/ode-system.lisp b/ordinary-differential-equations/ode-system.lisp
index bba290ad..5d11509d 100644
--- a/ordinary-differential-equations/ode-system.lisp
+++ b/ordinary-differential-equations/ode-system.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        ode-system.lisp      
-; description: ODE system setup
-; date:        Sun Apr 15 2007 - 14:19                   
-; author:      Liam Healy                                
-; modified:    Sat Jan  5 2008 - 21:39
-;********************************************************
-;;; $Id: $
+;; ODE system setup
+;; Liam Healy, Sun Apr 15 2007 - 14:19
+;; Time-stamp: <2008-01-14 22:58:19 liam ode-system.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -42,7 +38,7 @@
 	   (,params :pointer))
 	(declare (ignore ,params))
 	(,name ,time ,dependent ,derivatives)
-	(cffi:foreign-enum-value 'gsl-errorno :SUCCESS))
+	success)
       ;; The function should take four arguments: time, dependent, dfdy, dfdt
       ;; The last three will be arrays.
       (cffi:defcallback ,jacobian :int
@@ -53,7 +49,7 @@
 	   (,params :pointer))	
 	(declare (ignore ,params))
 	(,name ,time ,dependent ,dfdy ,dfdt)
-	(cffi:foreign-enum-value 'gsl-errorno :SUCCESS))
+	success)
       ;; Assume that defcallback does not bind the variable 'name.
       (defparameter ,name (cffi:foreign-alloc 'ode-system))
       (set-slot-function ,name 'ode-system 'function ',name)
diff --git a/roots-multi.lisp b/roots-multi.lisp
index 68c34f63..f5c4cf1a 100644
--- a/roots-multi.lisp
+++ b/roots-multi.lisp
@@ -1,11 +1,15 @@
 ;;; Multivariate roots.                
 ;;; Liam Healy 2008-01-12 12:49:08
-;;; Time-stamp: <2008-01-13 22:51:02 liam roots-multi.lisp>
+;;; Time-stamp: <2008-01-15 22:44:27 liam roots-multi.lisp>
 ;;; $Id: $
 
 (in-package :gsl)
 
-;;; I don't like using make-data-from-pointer.
+;;; I don't like using make-data-from-pointer, but it's the only way
+;;; to have access to the GSL functions when given a pointer.
+;;; Alternatively, I could provide the GSL pointer and then the only
+;;; thing the user could use is vref, or, of course
+;;; make-data-from-pointer.
 
 ;;;;****************************************************************************
 ;;;; Function definition
@@ -25,6 +29,16 @@
     ((dimensions ,dimensions))
     ((returned-value gsl-vector-c))))
 
+(cffi:defcstruct gsl-mfunction-fdf
+  ;; See /usr/include/gsl/gsl_multiroots.h
+  "The definition of a function and its derivatives for multiroot
+   finding in GSL."
+  (function :pointer)
+  (df :pointer)
+  (fdf :pointer)
+  (dimensions :size)
+  (parameters :pointer))
+
 ;;;;****************************************************************************
 ;;;; Initialization
 ;;;;****************************************************************************
@@ -55,7 +69,7 @@
 (defun-gsl set-mfdfsolver (solver function-derivative initial)
   "gsl_multiroot_fdfsolver_set"
   ((solver :pointer) (function-derivative :pointer)
-   ((gsl-array initial) :pointer))
+   ((pointer initial) :pointer))
   :documentation
   "Set or reset an existing solver to use the function and derivative
    (fdf) and the initial guess.")
@@ -383,13 +397,22 @@
 (defparameter *gsl-vector*
   (make-instance 'gsl-vector-double :pointer nil :storage-size nil))
 
+;;; One alternative way of writing the function, not recommended.
 (defun rosenbrock (argument return)
   "Rosenbrock test function."
-  (with-c-doubles (((gsl-array-p argument) x0 x1)
-		   ((gsl-array-p return) f0 f1))
+  (with-c-doubles (((vector-data argument) x0 x1)
+		   ((vector-data return) f0 f1))
     (setf f0 (* *rosenbrock-a* (- 1 x0))
 	  f1 (* *rosenbrock-b* (- x1 (expt x0 2))))))
 
+;;; The recommended alternative
+(defun rosenbrock (argument return)
+  "Rosenbrock test function."
+  (setf (vref return 0)
+	(* *rosenbrock-a* (- 1 (vref argument 0)))
+	(vref return 1)
+	(* *rosenbrock-b* (- (vref argument 1) (expt (vref argument 0) 2)))))
+
 (def-mfunction rosenbrock 2)
 
 (defun roots-multi-example ()
@@ -416,3 +439,51 @@
 				  (gsl-aref argval 1)
 				  (gsl-aref fnval 0)
 				  (gsl-aref fnval 1)))))))))
+
+(defun rosenbrock-df (argument jacobian)
+  "The partial derivatives of the Rosenbrock functions."
+  (setf (mref jacobian 0 0) (- *rosenbrock-a*)
+	(mref jacobian 0 1) 0.0d0
+	(mref jacobian 1 0) (* -2 *rosenbrock-b* (vref argument 0))
+	(mref jacobian 1 1) *rosenbrock-b*))
+
+(defun rosenbrock-fdf (argument value jacobian)
+  (rosenbrock argument value)
+  (rosenbrock-df argument jacobian))
+
+;;; Because def-solver-functions and def-scalar-function bind a symbol
+;;; of the same name as the first function, and we want both to run,
+;;; we'll make an alias function so we can use both.  
+(eval-when (:load-toplevel :execute)
+  (setf (fdefinition 'rosenbrock-f) #'rosenbrock))
+
+(def-solver-functions rosenbrock-f rosenbrock-df rosenbrock-fdf 2)
+
+(defun roots-multi-example-df ()
+  "Solving Rosenbrock with derivatives, the example given in Sec. 34.8
+   of the GSL manual."
+  (flet ((print-state (iter argval fnval)
+	   (format t "~&iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g"
+		   iter
+		   (gsl-aref argval 0)
+		   (gsl-aref argval 1)
+		   (gsl-aref fnval 0)
+		   (gsl-aref fnval 1))))
+    (let ((max-iter 1000))
+      (with-data (vect vector-double 2)
+	(setf (data vect) #(-10.0d0 -5.0d0))
+	(with-mfdfsolver (solver *gnewton-mfdfsolver* rosenbrock-f vect)
+	  (let ((fnval (mfdfsolver-f solver))
+		(argval (mfdfsolver-root solver)))
+	    (loop for iter from 0
+		  while (and (< iter max-iter)
+			     (not (multiroot-test-residual solver 1.0d-7)))
+		  initially (print-state iter argval fnval)
+		  do
+		  (iterate-mfdfsolver solver)
+		  (print-state iter argval fnval)
+		  finally (return
+			    (values (gsl-aref argval 0)
+				    (gsl-aref argval 1)
+				    (gsl-aref fnval 0)
+				    (gsl-aref fnval 1))))))))))
diff --git a/roots-one.lisp b/roots-one.lisp
index 0b0f833a..1958edbd 100644
--- a/roots-one.lisp
+++ b/roots-one.lisp
@@ -1,12 +1,7 @@
-;********************************************************
-; file:        roots-one.lisp                           
-; description: One-dimensional root solver.              
-; date:        Sun Dec  9 2007 - 15:47                   
-; author:      Liam Healy                                
-; modified:    Sun Jan 13 2008 - 17:00
-;********************************************************
-;;; Time-stamp: <2008-01-13 17:00:49 liam roots-one.lisp>
-;;; $Id: $
+;; One-dimensional root solver.
+;; Liam Healy 
+;; Time-stamp: <2008-01-15 22:57:11 liam roots-one.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -25,38 +20,51 @@
   (parameters :pointer))
 
 (export 'def-solver-functions)
-(defmacro def-solver-functions (function df fdf)
+(defmacro def-solver-functions (function df fdf &optional dimensions)
   "Setup functions for solvers.
    The CL functions name and derivative should be defined previously
-   with defuns."
+   with defuns.  If dimensions is non-nil, set multiroot solver
+   functions."
   (let ((arg (make-symbol "ARG"))
 	(params (make-symbol "PARAMS"))
 	(func (make-symbol "F"))
-	(deriv (make-symbol "DF")))
+	(deriv (make-symbol "DF"))
+	(struct (if dimensions 'gsl-mfunction-fdf 'gsl-function-fdf))
+	(argtype (if dimensions :pointer :double))
+	(rettype (if dimensions :int :double)))
     `(progn
-      (cffi:defcallback ,function :double
-	  ((,arg :double)
-	   (,params :pointer))
+      (cffi:defcallback ,function ,rettype
+	  ((,arg ,argtype)
+	   (,params :pointer)
+	   ,@(when dimensions `((,func :pointer))))
 	(declare (ignore ,params))
-	(,function ,arg))
-      (cffi:defcallback ,df :double
-	  ((,arg :double)
-	   (,params :pointer))
+	(,function ,arg ,@(when dimensions `(,func)))
+	,@(when dimensions '(success)))
+      (cffi:defcallback ,df ,rettype
+	  ((,arg ,argtype)
+	   (,params :pointer)
+	   ,@(when dimensions `((,deriv :pointer))))
 	(declare (ignore ,params))
-	(,df ,arg))
-      (cffi:defcallback ,fdf :pointer
-	  ((,arg :double)
+	(,df ,arg ,@(when dimensions `(,deriv)))
+	,@(when dimensions '(success)))
+      (cffi:defcallback ,fdf ,(if dimensions :int :pointer)
+	  ((,arg ,argtype)
 	   (,params :pointer)
 	   (,func :pointer)
 	   (,deriv :pointer))
 	(declare (ignore ,params))
 	(,fdf ,arg ,func ,deriv)
-	(cffi:null-pointer))
-      (defparameter ,function (cffi:foreign-alloc 'gsl-function-fdf))
-      (set-slot-function ,function 'gsl-function-fdf 'function ',function)
-      (set-slot-function ,function 'gsl-function-fdf 'df ',df)
-      (set-slot-function ,function 'gsl-function-fdf 'fdf ',fdf)
-      (set-parameters ,function 'gsl-function-fdf))))
+	;; fdf function returns sucess code in multivariate case, void
+	;; pointer in univariate
+	,(if dimensions 'success
+	     '(cffi:null-pointer)))
+      (defparameter ,function (cffi:foreign-alloc ',struct))
+      (set-slot-function ,function ',struct 'function ',function)
+      (set-slot-function ,function ',struct 'df ',df)
+      (set-slot-function ,function ',struct 'fdf ',fdf)
+      ,@(when dimensions
+	      `((set-structure-slot ,function ',struct 'dimensions ,dimensions)))
+      (set-parameters ,function ',struct))))
 
 ;;;;****************************************************************************
 ;;;; Initialization
@@ -395,7 +403,7 @@
 		    root (- root (sqrt 5.0d0))
 		    (- upper lower))))))
 
-;;; Because def-solver-functions and def-scalar-fucntion bind a symbol
+;;; Because def-solver-functions and def-scalar-function bind a symbol
 ;;; of the same name as the first function, and we want both to run,
 ;;; we'll make an alias function so we can use both.  
 (eval-when (:load-toplevel :execute)
@@ -420,3 +428,4 @@
 	    do
 	    (format t "~&~d~6t~10,8g ~18t~10,6g~34t~10,6g"
 		    iter root (- root (sqrt 5.0d0)) (- root oldroot))))))
+
-- 
GitLab