Commit 2597eaec authored by Liam M. Healy's avatar Liam M. Healy

Use grid:copy-to instead of cl-array

Do not use cl-array in tests because it isn't always available.  The
safer way is by doing a copy out to a new CL array.
parent 99257936
;; Numerical integration
;; Liam Healy, Wed Jul 5 2006 - 23:14
;; Time-stamp: <2010-07-06 23:26:39EDT numerical-integration.lisp>
;; Time-stamp: <2010-07-07 14:21:55EDT numerical-integration.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -367,6 +367,6 @@
(integration-QAGil (integration-test-myfn2 1.0d0) 1.0d0 1.0d-7 0.0d0 1000)
(integration-QAGp
'integration-test-f454
(copy (vector 0.0d0 1.0d0 (sqrt 2.0d0) 3.0d0) :grid-type 'foreign-array)
(grid:copy-to (vector 0.0d0 1.0d0 (sqrt 2.0d0) 3.0d0))
0.0d0 1.0d-3 1000)
(integration-QAWc 'integration-test-f459 -1.0d0 5.0d0 0.0d0 0.0d0 1.0d-3 1000))
;; Tests of array functions
;; Liam Healy 2008-10-20 22:41:48EDT array-tests.lisp
;; Time-stamp: <2010-07-06 23:51:52EDT array-tests.lisp>
;; Time-stamp: <2010-07-07 14:25:02EDT array-tests.lisp>
;;
;; Copyright 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -32,33 +32,33 @@
(generate-all-array-tests vector-set-all #+fsbv t #-fsbv :no-complex
(let ((v1 (array-default 3 t)))
(cl-array (set-all v1 (scalar-default)))))
(grid:copy-to (set-all v1 (scalar-default)))))
(generate-all-array-tests matrix-set-all #+fsbv t #-fsbv :no-complex
(let ((m1 (array-default '(3 3) t)))
(cl-array (set-all m1 (scalar-default)))))
(grid:copy-to (set-all m1 (scalar-default)))))
(generate-all-array-tests vector-set-zero t
(let ((v1 (array-default 3)))
(set-zero v1)
(cl-array v1)))
(grid:copy-to v1)))
(generate-all-array-tests matrix-set-zero t
(let ((m1 (array-default '(3 3))))
(set-zero m1)
(cl-array m1)))
(grid:copy-to m1)))
(generate-all-array-tests vector-swap t
(let ((v1 (array-default 3))
(v2 (array-default 3)))
(swap v2 v1)
(list (cl-array v1) (cl-array v2))))
(list (grid:copy-to v1) (grid:copy-to v2))))
(generate-all-array-tests matrix-swap t
(let ((m1 (array-default '(3 3)))
(m2 (array-default '(3 3))))
(swap m2 m1)
(list (cl-array m1) (cl-array m2))))
(list (grid:copy-to m1) (grid:copy-to m2))))
;;;;****************************************************************************
;;;; Arithmetic operations
......@@ -74,53 +74,53 @@
(generate-all-array-tests matrix-add :no-complex
(let ((m1 (array-default '(3 3)))
(m2 (array-default '(3 3))))
(cl-array (elt+ m1 m2))))
(grid:copy-to (elt+ m1 m2))))
(generate-all-array-tests vector-sub :no-complex
(let ((v1 (array-default 3))
(v2 (array-default 3)))
(cl-array (elt- v1 v2))))
(grid:copy-to (elt- v1 v2))))
(generate-all-array-tests matrix-sub :no-complex
(let ((m1 (array-default '(3 3)))
(m2 (array-default '(3 3))))
(cl-array (elt- m1 m2))))
(grid:copy-to (elt- m1 m2))))
(generate-all-array-tests vector-mult :no-complex
(let ((v1 (array-default 3))
(v2 (array-default 3)))
(cl-array (elt* v1 v2))))
(grid:copy-to (elt* v1 v2))))
(generate-all-array-tests matrix-mult :no-complex
(let ((m1 (array-default '(3 3)))
(m2 (array-default '(3 3))))
(cl-array (elt* m1 m2))))
(grid:copy-to (elt* m1 m2))))
(generate-all-array-tests vector-div :no-complex
(let ((v1 (array-default 3))
(v2 (array-default 3)))
(cl-array (elt/ v1 v2))))
(grid:copy-to (elt/ v1 v2))))
(generate-all-array-tests matrix-div :no-complex
(let ((m1 (array-default '(3 3)))
(m2 (array-default '(3 3))))
(cl-array (elt/ m1 m2))))
(grid:copy-to (elt/ m1 m2))))
(generate-all-array-tests vector-mult-scalar :no-complex
(let ((v1 (array-default 3)))
(cl-array (elt* v1 1.39d0))))
(grid:copy-to (elt* v1 1.39d0))))
(generate-all-array-tests matrix-mult-scalar :no-complex
(let ((m1 (array-default '(3 3))))
(cl-array (elt* m1 1.39d0))))
(grid:copy-to (elt* m1 1.39d0))))
(generate-all-array-tests vector-add-scalar :no-complex
(let ((v1 (array-default 3)))
(cl-array (elt+ v1 18.19d0))))
(grid:copy-to (elt+ v1 18.19d0))))
(generate-all-array-tests matrix-add-scalar :no-complex
(let ((m1 (array-default '(3 3))))
(cl-array (elt+ m1 18.19d0))))
(grid:copy-to (elt+ m1 18.19d0))))
;;;;****************************************************************************
;;;; Maximum and minimum elements
......@@ -183,16 +183,16 @@
(generate-all-array-tests set-basis t
(let ((v1 (array-default 8)))
(set-basis v1 2)
(cl-array v1)))
(grid:copy-to v1)))
(generate-all-array-tests swap-elements t
(let ((v1 (array-default 8)))
(swap-elements v1 2 5)
(cl-array v1)))
(grid:copy-to v1)))
(generate-all-array-tests vector-reverse t
(let ((v1 (array-default 8)))
(cl-array (vector-reverse v1))))
(grid:copy-to (vector-reverse v1))))
;;;;****************************************************************************
;;;; Matrix-only function definitions
......@@ -200,47 +200,47 @@
(generate-all-array-tests set-identity t
(let ((m1 (array-default '(3 3))))
(cl-array (set-identity m1))))
(grid:copy-to (set-identity m1))))
(generate-all-array-tests row t
(let ((m1 (array-default '(3 3)))
(row (array-default 3 t)))
(cl-array (row m1 1 row))))
(grid:copy-to (row m1 1 row))))
(generate-all-array-tests setf-row t
(let ((m1 (array-default '(3 3)))
(row (array-default 3)))
(setf (row m1 2) row)
(cl-array m1)))
(grid:copy-to m1)))
(generate-all-array-tests column t
(let ((m1 (array-default '(3 3)))
(col (array-default 3 t)))
(cl-array (column m1 1 col))))
(grid:copy-to (column m1 1 col))))
(generate-all-array-tests setf-column t
(let ((m1 (array-default '(3 3)))
(col (array-default 3)))
(setf (column m1 2) col)
(cl-array m1)))
(grid:copy-to m1)))
(generate-all-array-tests swap-rows t
(let ((m1 (array-default '(3 3))))
(cl-array (swap-rows m1 0 1))))
(grid:copy-to (swap-rows m1 0 1))))
(generate-all-array-tests swap-columns t
(let ((m1 (array-default '(3 3))))
(cl-array (swap-columns m1 1 2))))
(grid:copy-to (swap-columns m1 1 2))))
(generate-all-array-tests swap-row-column t
(let ((m1 (array-default '(3 3))))
(cl-array (swap-row-column m1 0 2))))
(grid:copy-to (swap-row-column m1 0 2))))
(generate-all-array-tests matrix-transpose* t
(let ((m1 (array-default '(3 3))))
(cl-array (matrix-transpose* m1))))
(grid:copy-to (matrix-transpose* m1))))
(generate-all-array-tests matrix-transpose t
(let ((m1 (array-default '(3 3)))
(m2 (array-default '(3 3) t)))
(cl-array (matrix-transpose m1 m2))))
(grid:copy-to (matrix-transpose m1 m2))))
;; Combinations
;; Liam Healy, Sun Mar 26 2006 - 11:51
;; Time-stamp: <2010-07-01 21:13:06EDT combination.lisp>
;; Time-stamp: <2010-07-07 14:30:07EDT combination.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -94,7 +94,7 @@
"Initialize the combination c to the lexicographically
last combination, i.e. (n-k,n-k+1,...,n-1).")
(defmfun grid:copy-to-destination ((source combination) (destination combination))
(defmfun grid::copy-to-destination ((source combination) (destination combination))
"gsl_combination_memcpy"
(((mpointer destination) :pointer)
((mpointer source) :pointer))
......
;; Permutations
;; Liam Healy, Sun Mar 26 2006 - 11:51
;; Time-stamp: <2010-07-01 21:12:59EDT permutation.lisp>
;; Time-stamp: <2010-07-07 14:30:06EDT permutation.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -79,7 +79,7 @@
"Initialize the permutation p to the identity, i.e.
(0,1,2,...,n-1).")
(defmfun grid:copy-to-destination ((source permutation) (destination permutation))
(defmfun grid::copy-to-destination ((source permutation) (destination permutation))
"gsl_permutation_memcpy"
(((mpointer destination) :pointer)
((mpointer source) :pointer))
......
;; Eigenvectors and eigenvalues
;; Liam Healy, Sun May 21 2006 - 19:52
;; Time-stamp: <2010-06-29 22:15:24EDT symmetric-hermitian.lisp>
;; Time-stamp: <2010-07-07 14:25:01EDT symmetric-hermitian.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -155,7 +155,7 @@
(let ((m1 (array-default '(5 5))))
(multiple-value-bind (eval evec)
(eigenvalues-eigenvectors (elt+ m1 (matrix-transpose-copy m1)))
(list (cl-array eval) (cl-array evec)))))
(list (grid:copy-to eval) (grid:copy-to evec)))))
|#
(defun eigenvalue-eigenvectors-example ()
......@@ -165,7 +165,7 @@
-10.0d0 30.0d0 0.0d0 ^
0.0d0 0.0d0 40.0d0)))
(eigenvalues-eigenvectors mat evals evecs)
(values (cl-array evals) (cl-array evecs))))
(values (grid:copy-to evals) (grid:copy-to evecs))))
(save-test eigensystems
(eigenvalue-eigenvectors-example))
......
;; Example FFT: transform a pulse (using the "clean" fft interface)
;; Sumant Oemrawsingh, Sat Oct 31 2009 - 00:24
;; Time-stamp: <2010-07-06 23:16:03EDT example.lisp>
;; Time-stamp: <2010-07-07 14:21:55EDT example.lisp>
;;
;; Copyright 2009 Sumant Oemrawsingh, Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -147,7 +147,7 @@
:initial-contents
(map 'list
(lambda (x) (/ x (floor (size backward) stride)))
(copy backward :grid-type 'array))))))
(grid:copy-to backward))))))
(multiple-value-bind (forward inverse)
(test-real-fft-noise random-vector :stride stride)
(values dft-random-vector forward random-vector inverse)))))
......
;; Discrete Hankel Transforms.
;; Liam Healy, Sat Dec 8 2007 - 16:50
;; Time-stamp: <2010-07-06 23:13:16EDT hankel.lisp>
;; Time-stamp: <2010-07-07 14:19:28EDT hankel.lisp>
;;
;; Copyright 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -69,8 +69,7 @@
(save-test hankel ; tests from dht/test.c
(copy (apply-hankel (make-hankel 3 1.0d0 1.0d0) #m(1.0d0 2.0d0 3.0d0))
:grid-type 'array)
(grid:copy-to (apply-hankel (make-hankel 3 1.0d0 1.0d0) #m(1.0d0 2.0d0 3.0d0)))
;; Exact, forward-inverse transform should be accurate to 2e-5
(let ((hank (make-hankel 3 1.0d0 1.0d0)))
(copy
......@@ -84,16 +83,14 @@
(loop for n from 0 below 128
do (setf (grid:gref in n)
(/ (1+ (expt (sample-x-hankel hank n) 2)))))
(copy (apply-hankel hank in)
:grid-type 'array))
(grid:copy-to (apply-hankel hank in)))
;; Integrate[ x exp(-x) J_1(a x), {x,0,Inf}] = a F(3/2, 2; 2; -a^2)
;; expected accuracy only about 2%
(let ((hank (make-hankel 128 1.0d0 20.0d0))
(in (grid:make-foreign-array 'double-float :dimensions 128)))
(loop for n from 0 below 128
do (setf (grid:gref in n) (exp (- (sample-x-hankel hank n)))))
(copy (apply-hankel hank in)
:grid-type 'array))
(grid:copy-to (apply-hankel hank in)))
;; Integrate[ x^2 (1-x^2) J_1(a x), {x,0,1}] = 2/a^2 J_3(a)
(let ((hank (make-hankel 128 1.0d0 1.0d0))
(in (grid:make-foreign-array 'double-float :dimensions 128)))
......@@ -101,4 +98,4 @@
do (setf (grid:gref in n)
(let ((x (sample-x-hankel hank n)))
(* x (- 1 (expt x 2))))))
(copy (apply-hankel hank in) :grid-type 'array)))
(grid:copy-to (apply-hankel hank in))))
;; The histogram structure
;; Liam Healy, Mon Jan 1 2007 - 11:32
;; Time-stamp: <2010-07-01 21:14:48EDT histogram.lisp>
;; Time-stamp: <2010-07-07 14:30:30EDT histogram.lisp>
;;
;; Copyright 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -79,7 +79,7 @@
:definition :method
:return (histogram))
(defmfun grid:copy-to-destination ((source histogram) (destination histogram))
(defmfun grid::copy-to-destination ((source histogram) (destination histogram))
"gsl_histogram_memcpy"
(((mpointer destination) :pointer) ((mpointer source) :pointer))
:definition :method
......@@ -91,7 +91,7 @@
an exact copy of the former.
The two histograms must be of the same size.")
(defmfun grid:copy-to-destination ((source histogram2d) (destination histogram2d))
(defmfun grid::copy-to-destination ((source histogram2d) (destination histogram2d))
"gsl_histogram2d_memcpy"
(((mpointer destination) :pointer) ((mpointer source) :pointer))
:definition :method
......@@ -103,14 +103,14 @@
an exact copy of the former.
The two histograms must be of the same size.")
(defmfun grid:copy-making-destination ((source histogram))
(defmfun grid::copy-making-destination ((source histogram))
"gsl_histogram_clone"
(((mpointer source) :pointer))
:definition :method
:c-return :pointer
:index copy)
(defmfun grid:copy-making-destination ((source histogram2d))
(defmfun grid::copy-making-destination ((source histogram2d))
"gsl_histogram2d_clone"
(((mpointer source) :pointer))
:definition :method
......
;; Definition of GSL objects and ways to use them.
;; Liam Healy, Sun Dec 3 2006 - 10:21
;; Time-stamp: <2010-07-01 19:37:25EDT mobject.lisp>
;; Time-stamp: <2010-07-07 14:30:30EDT mobject.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -351,7 +351,7 @@
;;;;****************************************************************************
#|
(defmethod grid:copy-making-destination :around ((object mobject))
(defmethod grid::copy-making-destination :around ((object mobject))
(if (next-method-p)
;; The subclass method should only return the malloced
;; mpointer (as from a "_clone" function); it will be put into
......
;; BLAS level 1, Vector operations
;; Liam Healy, Wed Apr 26 2006 - 15:23
;; Time-stamp: <2010-06-30 19:57:28EDT blas1.lisp>
;; Time-stamp: <2010-07-07 14:25:01EDT blas1.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -233,29 +233,29 @@
(let ((v1 (array-default 3))
(v2 (array-default 3)))
(blas-swap v2 v1)
(list (cl-array v1) (cl-array v2))))
(list (grid:copy-to v1) (grid:copy-to v2))))
(generate-all-array-tests blas-copy :float-complex
(let ((v1 (array-default 3))
(v2 (array-default 3 t)))
(blas-copy v1 v2)
(cl-array v2)))
(grid:copy-to v2)))
(generate-all-array-tests axpy :float-complex
(let ((v1 (array-default 8))
(v2 (array-default 8))
(scalar (scalar-default)))
(cl-array (axpy scalar v1 v2))))
(grid:copy-to (axpy scalar v1 v2))))
(generate-all-array-tests scale #+fsbv :float-complex #-fsbv :float
(let ((v1 (array-default 8))
(scalar (scalar-default)))
(cl-array (scale scalar v1))))
(grid:copy-to (scale scalar v1))))
(generate-all-array-tests scale :complex
(let ((v1 (array-default 8))
(scalar (scalar-default t)))
(cl-array (scale scalar v1))))
(grid:copy-to (scale scalar v1))))
(generate-all-array-tests givens :float
(let ((v1 (array-default 8))
......@@ -267,4 +267,4 @@
(setf (grid:gref sines i) (sin (grid:gref angles i)))
(setf (grid:gref cosines i) (cos (grid:gref angles i))))
(givens-rotation v1 v2 cosines sines)
(list (cl-array v1) (cl-array v2))))
(list (grid:copy-to v1) (grid:copy-to v2))))
;; BLAS level 2, Matrix-vector operations
;; Liam Healy, Wed Apr 26 2006 - 21:08
;; Time-stamp: <2010-06-29 22:15:24EDT blas2.lisp>
;; Time-stamp: <2010-07-07 14:25:00EDT blas2.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -340,21 +340,21 @@
(v2 (array-default 3))
(s1 (scalar-default))
(s2 (scalar-default)))
(cl-array (matrix-product m1 v1 v2 s1 s2))))
(grid:copy-to (matrix-product m1 v1 v2 s1 s2))))
(generate-all-array-tests matrix-product-triangular
#+fsbv :float-complex #-fsbv :float
(let ((m1 (array-default '(3 3)))
(v1 (array-default 3))
(s1 (scalar-default)))
(cl-array (matrix-product-triangular m1 v1 s1))))
(grid:copy-to (matrix-product-triangular m1 v1 s1))))
(generate-all-array-tests inverse-matrix-product
#+fsbv :float-complex #-fsbv :float
(let ((m1 (array-default '(3 3)))
(v1 (array-default 3))
(s1 (scalar-default)))
(cl-array (inverse-matrix-product m1 v1 s1))))
(grid:copy-to (inverse-matrix-product m1 v1 s1))))
(generate-all-array-tests matrix-product-symmetric :float
(let ((m1 (array-default '(3 3)))
......@@ -362,7 +362,7 @@
(v3 (array-default 3))
(s1 (scalar-default))
(s2 (scalar-default)))
(cl-array (matrix-product-symmetric m1 v1 v3 s1 s2))))
(grid:copy-to (matrix-product-symmetric m1 v1 v3 s1 s2))))
#+fsbv
(generate-all-array-tests matrix-product-hermitian :complex
......@@ -371,7 +371,7 @@
(v2 (array-default 3))
(s1 (scalar-default))
(s2 (scalar-default)))
(cl-array (matrix-product-hermitian m1 v1 v2 s1 s2))))
(grid:copy-to (matrix-product-hermitian m1 v1 v2 s1 s2))))
#|
;;; Error, needs to be tracked down
......@@ -383,5 +383,5 @@
(v1 (array-default 3))
(v2 (array-default 3))
(s1 (scalar-default)))
(cl-array (rank-1-update s1 v1 v2 m1))))
(grid:copy-to (rank-1-update s1 v1 v2 m1))))
|#
;; BLAS level 3, Matrix-matrix operations
;; Liam Healy, Wed Apr 26 2006 - 21:08
;; Time-stamp: <2010-06-29 22:15:23EDT blas3.lisp>
;; Time-stamp: <2010-07-07 14:25:00EDT blas3.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -177,26 +177,26 @@
(m3 (array-default '(3 3)))
(s1 (scalar-default))
(s2 (scalar-default)))
(cl-array (matrix-product m1 m2 m3 s1 s2))))
(grid:copy-to (matrix-product m1 m2 m3 s1 s2))))
(generate-all-array-tests matrix-product-nonsquare :float-complex
(let ((m1 (array-default '(2 3)))
(m2 (array-default '(3 2))))
(cl-array (matrix-product m1 m2))))
(grid:copy-to (matrix-product m1 m2))))
(generate-all-array-tests matrix-product-triangular
#+fsbv :float-complex #-fsbv :float
(let ((m1 (array-default '(3 3)))
(m2 (array-default '(3 3)))
(s1 (scalar-default)))
(cl-array (matrix-product-triangular m1 m2 s1))))
(grid:copy-to (matrix-product-triangular m1 m2 s1))))
(generate-all-array-tests inverse-matrix-product
#+fsbv :float-complex #-fsbv :float
(let ((m1 (array-default '(3 3)))
(m2 (array-default '(3 3)))
(s1 (scalar-default)))
(cl-array (inverse-matrix-product m1 m2 s1))))
(grid:copy-to (inverse-matrix-product m1 m2 s1))))
(generate-all-array-tests matrix-product-symmetric :float
(let ((m1 (array-default '(3 3)))
......@@ -204,7 +204,7 @@
(m3 (array-default '(3 3)))
(s1 (scalar-default))
(s2 (scalar-default)))
(cl-array (matrix-product-symmetric m1 m2 m3 s1 s2))))
(grid:copy-to (matrix-product-symmetric m1 m2 m3 s1 s2))))
#+fsbv
(generate-all-array-tests matrix-product-hermitian :complex
......@@ -213,4 +213,4 @@
(m3 (array-default '(3 3)))
(s1 (scalar-default))
(s2 (scalar-default)))
(cl-array (matrix-product-hermitian m1 m2 m3 s1 s2))))
(grid:copy-to (matrix-product-hermitian m1 m2 m3 s1 s2))))
;; Exponential of a matrix
;; Liam Healy 2008-08-10 17:25:35EDT exponential.lisp
;; Time-stamp: <2010-06-29 22:15:23EDT exponential.lisp>
;; Time-stamp: <2010-07-07 14:25:00EDT exponential.lisp>
;;
;; Copyright 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -41,7 +41,7 @@
;;; when exponentiated gives [[cos x, sin x], [-sin x, cos x]]
(let ((mat #m(0.0d0 1.0d0 ^ -1.0d0 0.0d0))
(exp (grid:make-foreign-array 'double-float :dimensions '(2 2))))
(cl-array (matrix-exponential mat exp)))
(grid:copy-to (matrix-exponential mat exp)))
#2A((0.5403023058681384d0 0.841470984807895d0)
(-0.841470984807895d0 0.5403023058681384d0))
|#
;; LU decomposition
;; Liam Healy, Thu Apr 27 2006 - 12:42
;; Time-stamp: <2010-06-29 22:15:23EDT lu.lisp>
;; Time-stamp: <2010-07-07 14:24:59EDT lu.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -164,7 +164,7 @@
(save-test
lu
(cl-array
(grid:copy-to
(invert-matrix
(grid:make-foreign-array 'double-float
:dimensions '(2 2)
......@@ -176,7 +176,7 @@
(multiple-value-bind (matrix perm)
(lu-decomposition matrix)
(let ((x (lu-solve matrix vec perm)))
(cl-array
(grid:copy-to
(permute-inverse
perm
(matrix-product-triangular
......
;; Polynomials
;; Liam Healy, Tue Mar 21 2006 - 18:33
;; Time-stamp: <2010-07-06 23:17:06EDT polynomial.lisp>
;; Time-stamp: <2010-07-07 14:18:55EDT polynomial.lisp>
;;
;; Copyright 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -211,8 +211,7 @@
(solve-cubic -6.0d0 -13.0d0 42.0d0)
(solve-cubic-complex -1.0d0 1.0d0 -1.0d0)
;; Example from GSL manual
(copy (polynomial-solve #m(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0))
:grid-type 'array)
(grid:copy-to (polynomial-solve #m(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0)))
;; tests from gsl-1.11/poly/test.c
(evaluate #m(1 0.5 0.3) 0.5d0)
(evaluate #m(1 -1 1 -1 1 -1 1 -1 1 -1 1) 1.0d0)
......@@ -243,18 +242,15 @@
(solve-cubic-complex -57.0d0 1071.0d0 -6647.0d0)
(solve-cubic-complex -11.0d0 -493.0d0 +6647.0d0)
(solve-cubic-complex -143.0d0 5087.0d0 -50065.0d0)
(copy (polynomial-solve #m(-120 274 -225 85 -15 1.0))
:grid-type 'array)
(copy (polynomial-solve #m(1 0 0 0 1 0 0 0 1))
:grid-type 'array)
(grid:copy-to (polynomial-solve #m(-120 274 -225 85 -15 1.0)))
(grid:copy-to (polynomial-solve #m(1 0 0 0 1 0 0 0 1)))
(let* ((xa #m(0.16 0.97 1.94 2.74 3.58 3.73 4.70))
(ya #m(0.73 1.11 1.49 1.84 2.30 2.41 3.07))
(dd (divided-difference xa ya)))
(list
(copy dd :grid-type 'array)
(grid:copy-to dd)
(map 'vector (lambda (x) (evaluate xa x :divided-difference dd))
(copy xa :grid-type 'array))
(grid:copy-to xa))
(map 'vector
(lambda (x) (evaluate (taylor-divided-difference 1.5d0 dd xa) (- x 1.5d0)))
(copy xa :grid-type 'array)))))
(grid:copy-to xa)))))
;; Dirichlet distribution
;; Liam Healy, Sun Oct 29 2006
;; Time-stamp: <2010-06-27 18:13:44EDT dirichlet.lisp>
;; Time-stamp: <2010-07-07 14:24:59EDT dirichlet.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -83,7 +83,7 @@
(save-test dirichlet
(let ((rng (make-random-number-generator +mt19937+ 0))
(alpha #m(1.0d0 2.0d0 3.0d0 4.0d0)))
(cl-array (sample rng :dirichlet :alpha alpha)))
(grid:copy-to (sample rng :dirichlet :alpha alpha)))
(let ((alpha #m(1.0d0 2.0d0 3.0d0 4.0d0))
(theta #m(1.0d0 2.0d0 3.0d0 4.0d0)))
(dirichlet-pdf alpha theta))
......
;; Generators of random numbers.
;; Liam Healy, Sat Jul 15 2006 - 14:43
;; Time-stamp: <2010-07-01 21:15:11EDT generators.lisp>
;; Time-stamp: <2010-07-07 14:30:30EDT generators.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -176,7 +176,7 @@
;;;; Copying state
;;;;****************************************************************************
(defmfun grid:copy-to-destination
(defmfun grid::copy-to-destination
((source random-number-generator) (destination random-number-generator))
"gsl_rng_memcpy"
(((mpointer destination) :pointer) ((mpointer source) :pointer))
......@@ -188,7 +188,7 @@
making destination into an exact copy
of source. The two generators must be of the same type.")
(defmfun grid:copy-making-destination ((instance random-number-generator))
(defmfun grid::copy-making-destination ((instance random-number-generator))
"gsl_rng_clone" (((mpointer instance) :pointer))
:definition :method
:c-return :pointer
......
;; Multinomial distribution
;; Liam Healy, Sat Nov 25 2006 - 16:00
;; Time-stamp: <2010-06-29 22:15:17EDT multinomial.lisp>
;; Time-stamp: <2010-07-07 14:24:58EDT multinomial.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -78,7 +78,7 @@
(save-test multinomial
(let ((rng (make-random-number-generator +mt19937+ 0))
(p #m(0.1d0 0.2d0 0.3d0 0.4d0)))
(cl-array (sample rng :multinomial :sum 8 :probabilities p)))
(grid:copy-to (sample rng :multinomial :sum 8 :probabilities p)))
(let ((p #m(0.1d0 0.2d0 0.3d0 0.4d0))
(n #31m(5 0 1 2)))
(multinomial-pdf p N))
......
;; Quasi-random sequences in arbitrary dimensions.
;; Liam Healy, Sun Jul 16 2006 - 15:54
;; Time-stamp: <2010-07-01 21:15:03EDT quasi.lisp>
;; Time-stamp: <2010-07-07 14:30:06EDT quasi.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -64,7 +64,7 @@
:export nil
:index gsl-random-state)
(defmfun grid:copy-to-destination
(defmfun grid::copy-to-destination
((source quasi-random-number-generator)
(destination quasi-random-number-generator))
"gsl_qrng_memcpy"
......@@ -76,7 +76,7 @@