Commit 2a6a743f authored by Liam Healy's avatar Liam Healy
Browse files

Be strict on initial marray values; default marrays

Removed some marrays marked as :input when their values were not used.
This failed in CCL because the arrays in those cases were created with
no initial values, and no values were set, so that they were marked
cl-invalid.  Created #'make-marray-or-default to facilitate defaulting
of marrays in defmfun, and marking validity appropriately.
Test results:
SBCL 64:
TOTAL: 1627 assertions passed, 6 failed, 0 execution errors.
CCL 64:
TOTAL: 1620 assertions passed, 13 failed, 0 execution errors.
parent a13502c2
;; Functions for both vectors and matrices.
;; Liam Healy 2008-04-26 20:48:44EDT both.lisp
;; Time-stamp: <2009-12-21 10:47:55EST both.lisp>
;; Time-stamp: <2009-12-26 10:18:59EST both.lisp>
(in-package :gsl)
......@@ -39,7 +39,6 @@
(((mpointer object) :pointer) (value :element-c-type))
:definition :generic
:element-types #+fsbv t #-fsbv :no-complex
:inputs (object)
:outputs (object)
:c-return :void
:documentation "Set all elements to the value.")
......
;; Combinations
;; Liam Healy, Sun Mar 26 2006 - 11:51
;; Time-stamp: <2009-12-21 10:21:30EST combination.lisp>
;; Time-stamp: <2009-12-26 10:38:37EST combination.lisp>
;; $Id$
(in-package :gsl)
......@@ -60,7 +60,6 @@
"gsl_combination_init_first"
(((mpointer combination) :pointer))
:c-return :void
:inputs (combination)
:outputs (combination)
:documentation ; FDL
"Initialize the combination c to the lexicographically
......@@ -70,7 +69,6 @@
"gsl_combination_init_last"
(((mpointer combination) :pointer))
:c-return :void
:inputs (combination)
:outputs (combination)
:documentation ; FDL
"Initialize the combination c to the lexicographically
......@@ -106,7 +104,6 @@
(((mpointer c) :pointer))
:definition :method
:c-return sizet
:inputs (c)
:documentation ; FDL
"The number of elements (k) in the combination c.")
......
;; A "marray" is an array in both GSL and CL
;; Liam Healy 2008-04-06 21:23:41EDT
;; Time-stamp: <2009-12-25 23:04:11EST marray.lisp>
;; Time-stamp: <2009-12-26 12:11:01EST marray.lisp>
(in-package :gsl)
......@@ -234,3 +234,24 @@
;; foreign-friendly array. There is no choice but to copy over the
;; data even on native implementations; because GSL is doing the
;; mallocing, the data are not CL-accessible.
;;;;****************************************************************************
;;;; Convenient defaulting of defmfun arguments
;;;;****************************************************************************
(defun make-marray-or-default
(default dimensions
&optional dont-make-if-nil (element-type 'double-float) initial-element)
"If default is T or dont-make-if-nil, make a marray and
returned with the dimensions and element-type. If default is a
marray, it is returned after being synced for non-native."
(if (member default (list t dont-make-if-nil))
(if initial-element
(make-marray
element-type :dimensions dimensions :initial-element initial-element)
(make-marray element-type :dimensions dimensions))
(progn #-native
(when (typep default 'marray)
(c-array:copy-cl-to-c default))
default)))
;; BLAS level 2, Matrix-vector operations
;; Liam Healy, Wed Apr 26 2006 - 21:08
;; Time-stamp: <2009-09-26 12:52:44EDT blas2.lisp>
;; Time-stamp: <2009-12-26 11:54:21EST blas2.lisp>
;; $Id$
(in-package :gsl)
......@@ -36,16 +36,15 @@
y
(alpha 1) (beta 1) (TransA :notrans) TransB
&aux
(yarr (or y
(make-marray
element-type :dimensions (matrix-product-dimensions A x)
:initial-element 0))))
(yarr
(make-marray-or-default
y (matrix-product-dimensions A x) nil element-type 0)))
("gsl_blas_" :type "gemv")
((transa cblas-transpose) (alpha :element-c-type) ((mpointer A) :pointer)
((mpointer x) :pointer) (beta :element-c-type) ((mpointer yarr) :pointer))
:definition :generic
:element-types #+fsbv :float-complex #-fsbv :float
:inputs (A x yarr)
:inputs (A x)
:outputs (yarr)
:documentation ; FDL
"If the second and third arguments are vectors, compute
......
;; Cholesky Decomposition
;; Liam Healy, Wed May 3 2006 - 16:38
;; Time-stamp: <2009-09-26 12:36:23EDT cholesky.lisp>
;; Time-stamp: <2009-12-26 12:38:34EST cholesky.lisp>
;; $Id$
(in-package :gsl)
......@@ -38,9 +38,7 @@
(defmfun cholesky-solve
((A matrix) (b vector) &optional x-spec
&aux
(x (if (eq x-spec t)
(make-marray element-type :dimensions (dimensions b))
x-spec)))
(x (make-marray-or-default x-spec (dimensions b) t)))
(("gsl_linalg" :complex "_cholesky_svx")
("gsl_linalg" :complex "_cholesky_solve"))
((((mpointer A) :pointer) ((mpointer b) :pointer))
......
;; Householder Transformations
;; Liam Healy, Wed May 10 2006 - 10:03
;; Time-stamp: <2009-09-24 22:18:00EDT householder.lisp>
;; Time-stamp: <2009-12-26 12:38:35EST householder.lisp>
;; $Id$
;;; For householder-transform, it would be nice to be able to pick the
......@@ -72,9 +72,7 @@
(defmfun householder-solve
(A b &optional x-spec
&aux
(x (if (eq x-spec t)
(make-marray 'double-float :dimensions (dimensions b))
x-spec)))
(x (make-marray-or-default x-spec (dimensions b) t)))
("gsl_linalg_HH_svx" "gsl_linalg_HH_solve")
((((mpointer A) :pointer) ((mpointer b) :pointer))
(((mpointer A) :pointer) ((mpointer b) :pointer)
......
;; QR decomposition
;; Liam Healy 2008-02-17 11:05:20EST qr.lisp
;; Time-stamp: <2009-12-22 22:34:28EST qr.lisp>
;; Time-stamp: <2009-12-26 12:15:35EST qr.lisp>
;; $Id$
(in-package :gsl)
......@@ -46,16 +46,14 @@
(defmfun QR-solve
(QR tau b &optional x-spec
&aux
(x (if (eq x-spec t)
(make-marray 'double-float :dimensions (dimensions b))
x-spec)))
&aux
(x (make-marray-or-default x-spec (dimensions b) t)))
("gsl_linalg_QR_svx" "gsl_linalg_QR_solve")
((((mpointer QR) :pointer) ((mpointer tau) :pointer)
((mpointer b) :pointer))
(((mpointer QR) :pointer) ((mpointer tau) :pointer)
((mpointer b) :pointer) ((mpointer x) :pointer)))
:inputs (QR tau b x)
:inputs (QR tau b)
:outputs (x)
:return ((or x b))
:documentation ; FDL
......@@ -116,9 +114,7 @@
(defmfun QR-Rsolve
(QR b &optional x-spec
&aux
(x (if (eq x-spec t)
(make-marray 'double-float :dimensions (dimensions b))
x-spec)))
(x (make-marray-or-default x-spec (dimensions b) t)))
("gsl_linalg_QR_Rsvx" "gsl_linalg_QR_Rsolve")
((((mpointer QR) :pointer) ((mpointer b) :pointer))
(((mpointer QR) :pointer) ((mpointer b) :pointer) ((mpointer x) :pointer)))
......@@ -178,9 +174,7 @@
(defmfun R-solve
(R b &optional x-spec
&aux
(x (if (eq x-spec t)
(make-marray 'double-float :dimensions (dimensions b))
x-spec)))
(x (make-marray-or-default x-spec (dimensions b) t)))
("gsl_linalg_R_svx" "gsl_linalg_R_solve")
((((mpointer R) :pointer) ((mpointer b) :pointer))
(((mpointer R) :pointer) ((mpointer b) :pointer) ((mpointer x) :pointer)))
......@@ -252,7 +246,7 @@
(lambda (i j) (+ (maref matrix i j) (* (maref u i) (maref v j))))
dim0 dim1))
(qr2 (copy matrix))
(w (make-marray 'double-float :dimensions dim0)))
(w (make-marray 'double-float :dimensions dim0 :initial-element 0)))
(multiple-value-bind (QR2 tau)
(QR-decomposition qr2)
(multiple-value-bind (Q2 R2)
......
;; QR with column pivoting
;; Liam Healy, Fri Apr 28 2006 - 16:53
;; Time-stamp: <2009-09-26 16:46:01EDT qrpt.lisp>
;; Time-stamp: <2009-12-26 12:33:04EST qrpt.lisp>
;; $Id$
(in-package :gsl)
......@@ -65,16 +65,14 @@
(defmfun QRPT-solve
(QR tau permutation b &optional x-spec
&aux
(x (if (eq x-spec t)
(make-marray 'double-float :dimensions (dimensions b))
x-spec)))
(x (make-marray-or-default x-spec (dimensions b) t)))
("gsl_linalg_QRPT_svx" "gsl_linalg_QRPT_solve")
((((mpointer QR) :pointer) ((mpointer tau) :pointer)
((mpointer permutation) :pointer) ((mpointer b) :pointer))
(((mpointer QR) :pointer) ((mpointer tau) :pointer)
((mpointer permutation) :pointer)
((mpointer b) :pointer) ((mpointer x) :pointer)))
:inputs (QR tau permutation b x)
:inputs (QR tau permutation b)
:outputs (x)
:return ((or x b))
:documentation ; FDL
......@@ -118,15 +116,13 @@
(defmfun QRPT-Rsolve
(QR permutation b &optional x-spec
&aux
(x (if (eq x-spec t)
(make-marray 'double-float :dimensions (dimensions b))
x-spec)))
(x (make-marray-or-default x-spec (dimensions b) t)))
("gsl_linalg_QRPT_Rsvx" "gsl_linalg_QRPT_Rsolve")
((((mpointer QR) :pointer) ((mpointer permutation) :pointer)
((mpointer b) :pointer))
(((mpointer QR) :pointer) ((mpointer permutation) :pointer)
((mpointer b) :pointer) ((mpointer x) :pointer)))
:inputs (QR permutation b x)
:inputs (QR permutation b)
:outputs (x b)
:return ((or x b))
:documentation ; FDL
......@@ -168,7 +164,12 @@
(multiple-value-bind (Q R)
(QR-unpack QRPT tau)
(let* ((qr (matrix-product Q R))
(ans (make-marray 'double-float :dimensions (dimensions qr))))
(ans (make-marray
'double-float :dimensions (dimensions qr)
;; It shouldn't be necessary to initialize values
;; because we are going to setf every row,
;; but it is for #-native. See Sat Dec 26 2009.
:initial-element 0)))
(dotimes (i (dim0 qr) ans)
(setf (row ans i) (permute-inverse permutation (row qr i))))))))
......
;; Shuffling and sampling
;; Liam Healy, Sat Dec 2 2006 - 18:40
;; Time-stamp: <2009-12-21 09:04:36EST shuffling-sampling.lisp>
;; Time-stamp: <2009-12-26 10:53:46EST shuffling-sampling.lisp>
;; $Id$
(in-package :gsl)
......@@ -38,18 +38,18 @@
((c-pointer destarr) :pointer) ((dim0 destarr) sizet)
((c-pointer src) :pointer) ((dim0 src) sizet) ((c-array:element-size src) sizet))
:definition :method
:inputs (destarr src)
:inputs (src)
:outputs (destarr)
:documentation ; FDL
"Fill the array dest[k] with k objects taken randomly from the n
"Fill the array destarr[k] with k objects taken randomly from the n
elements of the array src[0...n-1]. The output of the random
number generator r is used to make the selection. The algorithm
ensures all possible samples are equally likely, assuming a perfect
source of randomness.
The objects are sampled without replacement, thus each object can
only appear once in dest[k]. It is required that k be less
than or equal to n. The objects in dest will be in the
only appear once in destarr[k]. It is required that k be less
than or equal to n. The objects in destarr will be in the
same relative order as those in src. You will need to call
with 'shuffle if you want to randomize the order.")
......@@ -66,7 +66,7 @@
((c-pointer destarr) :pointer) ((dim0 destarr) sizet)
((c-pointer src) :pointer) ((dim0 src) sizet) ((c-array:element-size src) sizet))
:definition :method
:inputs (destarr src)
:inputs (src)
:outputs (destarr)
:c-return :void
:documentation
......
......@@ -7,17 +7,17 @@
(LIST #(4 3 6 1 5 7 2 8))
(MULTIPLE-VALUE-LIST
(let ((rng (make-random-number-generator +mt19937+ 0))
(v1 #31m(1 2 3 4 5 6 7 8)))
(cl-array (sample rng 'shuffle :base v1)))))
(v1 #31m(1 2 3 4 5 6 7 8)))
(cl-array (sample rng 'shuffle :base v1)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(2 3 5 8))
(MULTIPLE-VALUE-LIST
(let ((rng (make-random-number-generator +mt19937+ 0))
(v1 #31m(1 2 3 4 5 6 7 8)))
(cl-array (sample rng 'choose-random :src v1 :dest 4)))))
(v1 #31m(1 2 3 4 5 6 7 8)))
(cl-array (sample rng 'choose-random :src v1 :dest 4)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(8 2 3 8 2 4 8 6 5 6))
(MULTIPLE-VALUE-LIST
(let ((rng (make-random-number-generator +mt19937+ 0))
(v1 #31m(1 2 3 4 5 6 7 8)))
(cl-array (sample rng 'random-sample :src v1 :dest 10))))))
(v1 #31m(1 2 3 4 5 6 7 8)))
(cl-array (sample rng 'random-sample :src v1 :dest 10))))))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment