Skip to content
Snippets Groups Projects
Commit 4f581352 authored by Liam Healy's avatar Liam Healy
Browse files

Redesign polynomial-solve, invert matrix

The function #'polynomial-solve is now directly defined with defmfun
with a single required argument and two optional arguments: the answer
vector-complex-double-float and the workspace, both of which default
to the right thing.  The returned is now the this marray instead of a
list.  Simplify the definition of #'invert-matrix and return the
matrix-double-float instead of the CL array.  All tests pass SBCL and
CCL with trace.
parent 54353bcb
No related branches found
No related tags found
No related merge requests found
;; 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)))))
;; 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))))
......@@ -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)))))))
......@@ -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)))))))
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment