diff --git a/linear-algebra/lu.lisp b/linear-algebra/lu.lisp index 808cf66927d67ae02f788f59d3569c8b19dc3c63..b2fa5a49317b0a524a3825eb20a5f83eee0bbcea 100644 --- a/linear-algebra/lu.lisp +++ b/linear-algebra/lu.lisp @@ -1,6 +1,6 @@ ;; LU decomposition ;; Liam Healy, Thu Apr 27 2006 - 12:42 -;; Time-stamp: <2008-12-27 10:57:53EST lu.lisp> +;; Time-stamp: <2008-12-30 22:32:38EST lu.lisp> ;; $Id$ (in-package :gsl) @@ -133,16 +133,16 @@ (export 'invert-matrix) (defun invert-matrix (mat) "Invert the matrix." - (let* ((mmat mat) ; BROKEN, requires copy of mat - (dim (first (dimensions mat))) + (let* ((dim (first (dimensions mat))) (per (make-permutation dim)) (inv (make-marray 'double-float :dimensions (list dim dim)))) - (LU-decomposition mmat per) - (lu-invert mmat per inv) - (cl-array inv))) + (LU-decomposition mat per) + (lu-invert mat per inv))) -(save-test lu - (invert-matrix - (make-marray 'double-float - :dimensions '(2 2) - :initial-contents '(1.0d0 2.0d0 3.0d0 4.0d0)))) +(save-test + lu + (cl-array + (invert-matrix + (make-marray 'double-float + :dimensions '(2 2) + :initial-contents '(1.0d0 2.0d0 3.0d0 4.0d0))))) diff --git a/polynomial.lisp b/polynomial.lisp index 83c8e4386fe16fee290266feb5661bf31cb1cd58..1063aecce13c7f0e515e990232469c67a30445f1 100644 --- a/polynomial.lisp +++ b/polynomial.lisp @@ -1,6 +1,6 @@ ;; Polynomials ;; Liam Healy, Tue Mar 21 2006 - 18:33 -;; Time-stamp: <2008-12-30 18:38:17EST polynomial.lisp> +;; Time-stamp: <2008-12-30 21:33:37EST polynomial.lisp> ;; $Id$ (in-package :gsl) @@ -127,39 +127,26 @@ ((n sizet)) "complex workspace for polynomials") -(export 'polynomial-solve) -(defun polynomial-solve (coefficients) - ;; FDL - "The roots of the general polynomial - P(x) = a_0 + a_1 x + a_2 x^2 + ... + a_{n-1} x^{n-1} using - balanced-QR reduction of the companion matrix. The parameter n - specifies the length of the coefficient array. The coefficient of the - highest order term must be non-zero. The function requires a workspace - w of the appropriate size. The n-1 roots are returned in - the packed complex array z of length 2(n-1), alternating - real and imaginary parts." - (let ((len (total-size coefficients))) - ;; Should this be making a complex array? - (let ((answer (make-marray 'double-float :dimensions (* 2 (1- len)))) - (ws (make-polynomial-complex-workspace len))) - (values-list (polynomial-solve-ws coefficients ws answer))))) - -(defmfun polynomial-solve-ws (coefficients workspace answer-pd) +(defmfun polynomial-solve + (coefficients + &optional + (answer (make-marray '(complex double-float) + :dimensions (1- (total-size coefficients)))) + (workspace (make-polynomial-complex-workspace (total-size coefficients)))) "gsl_poly_complex_solve" (((c-pointer coefficients) :pointer) ((dim0 coefficients) sizet) - ((mpointer workspace) :pointer) ((c-pointer answer-pd) :pointer)) + ((mpointer workspace) :pointer) ((c-pointer answer) :pointer)) :inputs (coefficients) - :outputs (answer-pd) - :return - ((loop for i from 0 below (dim0 answer-pd) by 2 - collect (complex (maref answer-pd i) - (maref answer-pd (1+ i))))) - :export nil - :index polynomial-solve + :outputs (answer) + :return (answer) :documentation ; FDL - "Arguments are: - a GSL array of coefficients, a workspace, a gsl-array of doubles.") - + "Arguments are: a vector-double-float of coefficients, a complex + vector of length one less than coefficients that will hold the + answer, and a workspace made by make-polynomial-complex-workspace. + The roots of the general polynomial + P(x) = a_0 + a_1 x + a_2 x^2 + ... + a_{n-1} x^{n-1} using + balanced-QR reduction of the companion matrix. The coefficient of the + highest order term must be non-zero.") ;;;;**************************************************************************** ;;;; Examples and unit test @@ -183,5 +170,5 @@ (solve-cubic -6.0d0 -13.0d0 42.0d0) (solve-cubic-complex -1.0d0 1.0d0 -1.0d0) ;; Example from GSL manual - (polynomial-solve #m(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0))) + (cl-array (polynomial-solve #m(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0)))) diff --git a/tests/lu.lisp b/tests/lu.lisp index 78116612fa1a57cebe0c43d787f408b6d2874b81..1556be227fb829d9d48537823fe89cf6d28338c5 100644 --- a/tests/lu.lisp +++ b/tests/lu.lisp @@ -8,8 +8,9 @@ #2A((-1.9999999999999998d0 1.0d0) (1.4999999999999998d0 -0.49999999999999994d0))) (MULTIPLE-VALUE-LIST - (INVERT-MATRIX - (MAKE-MARRAY 'DOUBLE-FLOAT :DIMENSIONS '(2 2) - :INITIAL-CONTENTS - '(1.0d0 2.0d0 3.0d0 4.0d0)))))) + (CL-ARRAY + (INVERT-MATRIX + (MAKE-MARRAY 'DOUBLE-FLOAT :DIMENSIONS '(2 2) + :INITIAL-CONTENTS + '(1.0d0 2.0d0 3.0d0 4.0d0))))))) diff --git a/tests/polynomial.lisp b/tests/polynomial.lisp index 577fa5d5540a34e36509dc262b023e8fbc5034fe..d71a547a6e0b3b252975eb510d51a3e80ee23be8 100644 --- a/tests/polynomial.lisp +++ b/tests/polynomial.lisp @@ -57,14 +57,16 @@ (MULTIPLE-VALUE-LIST (SOLVE-CUBIC-COMPLEX -1.0d0 1.0d0 -1.0d0))) (LISP-UNIT::ASSERT-NUMERICAL-EQUAL - (LIST #C(-0.8090169943749477d0 0.5877852522924734d0) - #C(-0.8090169943749477d0 -0.5877852522924734d0) - #C(0.3090169943749475d0 0.951056516295153d0) - #C(0.3090169943749475d0 -0.951056516295153d0) - #C(0.9999999999999999d0 0.0d0)) + (LIST + #(#C(-0.8090169943749477d0 0.5877852522924734d0) + #C(-0.8090169943749477d0 -0.5877852522924734d0) + #C(0.3090169943749475d0 0.951056516295153d0) + #C(0.3090169943749475d0 -0.951056516295153d0) + #C(0.9999999999999999d0 0.0d0))) (MULTIPLE-VALUE-LIST - (POLYNOMIAL-SOLVE - (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS - '(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 - 1.0d0)))))) + (CL-ARRAY + (POLYNOMIAL-SOLVE + (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS + '(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 + 1.0d0)))))))