Commit 74aed234 authored by Liam M. Healy's avatar Liam M. Healy

Use grid:aref instead of grid:gref

parent fcc1723d
;; Basis splines. ;; Basis splines.
;; Liam Healy 2008-02-18 14:43:20EST basis-splines.lisp ;; Liam Healy 2008-02-18 14:43:20EST basis-splines.lisp
;; Time-stamp: <2011-01-13 10:24:59EST basis-splines.lisp> ;; Time-stamp: <2011-05-26 12:37:29EDT basis-splines.lisp>
;; ;;
;; Copyright 2008, 2009, 2011 Liam M. Healy ;; Copyright 2008, 2009, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License ;; Distributed under the terms of the GNU General Public License
...@@ -125,18 +125,18 @@ ...@@ -125,18 +125,18 @@
(let* ((xi (coerce (* i (/ 15 (1- ndata))) 'double-float)) (let* ((xi (coerce (* i (/ 15 (1- ndata))) 'double-float))
(yi (+ (* (cos xi) (exp (* -0.1d0 xi))) (yi (+ (* (cos xi) (exp (* -0.1d0 xi)))
(sample rng :gaussian :sigma sigma)))) (sample rng :gaussian :sigma sigma))))
(setf (grid:gref x i) xi (setf (grid:aref x i) xi
(grid:gref y i) yi (grid:aref y i) yi
(grid:gref w i) (/ (expt sigma 2))))) (grid:aref w i) (/ (expt sigma 2)))))
;; Uniform breakpoints [0, 15] ;; Uniform breakpoints [0, 15]
(uniform-knots 0.0d0 15.0d0 bw) (uniform-knots 0.0d0 15.0d0 bw)
;; Fit matrix ;; Fit matrix
(dotimes (i ndata) (dotimes (i ndata)
;; Compute B_j ;; Compute B_j
(evaluate bw (grid:gref x i) :b B) (evaluate bw (grid:aref x i) :b B)
;; Fill in row i of X ;; Fill in row i of X
(dotimes (j ncoeffs) (dotimes (j ncoeffs)
(setf (grid:gref Xmatrix i j) (grid:gref B j)))) (setf (grid:aref Xmatrix i j) (grid:aref B j))))
;; Do the fit ;; Do the fit
(linear-mfit Xmatrix y c w nil cov mw) (linear-mfit Xmatrix y c w nil cov mw)
;; Return the smoothed curve ;; Return the smoothed curve
......
;; Functions for both vectors and matrices. ;; Functions for both vectors and matrices.
;; Liam Healy 2008-04-26 20:48:44EDT both.lisp ;; Liam Healy 2008-04-26 20:48:44EDT both.lisp
;; Time-stamp: <2011-01-12 00:38:30EST both.lisp> ;; Time-stamp: <2011-05-26 12:37:36EDT both.lisp>
;; ;;
;; Copyright 2008, 2009, 2010, 2011 Liam M. Healy ;; Copyright 2008, 2009, 2010, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License ;; Distributed under the terms of the GNU General Public License
...@@ -81,7 +81,7 @@ ...@@ -81,7 +81,7 @@
;;;;**************************************************************************** ;;;;****************************************************************************
;;;; Array elements; used in callbacks scalarsp=T only ;;;; Array elements; used in callbacks scalarsp=T only
;;;;**************************************************************************** ;;;;****************************************************************************
;;; Normal foreign array access is with grid:gref, but in order to ;;; Normal foreign array access is with grid:aref, but in order to
;;; avoid the overhead of instantiating a foreign-array object to ;;; avoid the overhead of instantiating a foreign-array object to
;;; access components, we use these macros which expand to gsl_*_get ;;; access components, we use these macros which expand to gsl_*_get
;;; and gsl_*_set. ;;; and gsl_*_set.
......
;; Permutations ;; Permutations
;; Liam Healy, Sun Mar 26 2006 - 11:51 ;; Liam Healy, Sun Mar 26 2006 - 11:51
;; Time-stamp: <2011-01-10 18:16:26EST permutation.lisp> ;; Time-stamp: <2011-05-26 12:37:36EDT permutation.lisp>
;; ;;
;; Copyright 2006, 2007, 2008, 2009, 2010, 2011 Liam M. Healy ;; Copyright 2006, 2007, 2008, 2009, 2010, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License ;; Distributed under the terms of the GNU General Public License
...@@ -339,7 +339,7 @@ ...@@ -339,7 +339,7 @@
(save-test permutation (save-test permutation
(let ((perm-1 (make-permutation 4 t))) (let ((perm-1 (make-permutation 4 t)))
(grid:gref perm-1 2)) (grid:aref perm-1 2))
(let ((perm-1 (make-permutation 4 t))) ;grid:contents (let ((perm-1 (make-permutation 4 t))) ;grid:contents
(grid:contents perm-1)) (grid:contents perm-1))
(let ((perm-1 (make-permutation 4 t))) ;permutation-reverse (let ((perm-1 (make-permutation 4 t))) ;permutation-reverse
......
;; Example FFT: transform a pulse (using the "clean" fft interface) ;; Example FFT: transform a pulse (using the "clean" fft interface)
;; Sumant Oemrawsingh, Sat Oct 31 2009 - 00:24 ;; Sumant Oemrawsingh, Sat Oct 31 2009 - 00:24
;; Time-stamp: <2011-01-11 23:32:17EST example.lisp> ;; Time-stamp: <2011-05-26 12:37:35EDT example.lisp>
;; ;;
;; Copyright 2009, 2010, 2011 Sumant Oemrawsingh, Liam M. Healy ;; Copyright 2009, 2010, 2011 Sumant Oemrawsingh, Liam M. Healy
;; Distributed under the terms of the GNU General Public License ;; Distributed under the terms of the GNU General Public License
...@@ -42,10 +42,10 @@ ...@@ -42,10 +42,10 @@
(assert (and (integerp dimension) (> dimension 20))) (assert (and (integerp dimension) (> dimension 20)))
(let ((pulse (grid:make-foreign-array element-type :dimensions dimension)) (let ((pulse (grid:make-foreign-array element-type :dimensions dimension))
(init-value (coerce 1 element-type))) (init-value (coerce 1 element-type)))
(setf (grid:gref pulse 0) init-value) (setf (grid:aref pulse 0) init-value)
(loop for i from 1 to 10 (loop for i from 1 to 10
do (setf (grid:gref pulse i) init-value do (setf (grid:aref pulse i) init-value
(grid:gref pulse (- dimension i)) init-value)) (grid:aref pulse (- dimension i)) init-value))
(forward-fourier-transform pulse))) (forward-fourier-transform pulse)))
(save-test (save-test
...@@ -82,7 +82,7 @@ ...@@ -82,7 +82,7 @@
(when init-offset (when init-offset
(loop for i from 0 below size (loop for i from 0 below size
do do
(setf (grid:gref vec i) (setf (grid:aref vec i)
(if (subtypep element-type 'complex) (if (subtypep element-type 'complex)
(coerce (complex (+ (* 2 i) init-offset) (coerce (complex (+ (* 2 i) init-offset)
(+ (* 2 i) init-offset 1)) (+ (* 2 i) init-offset 1))
...@@ -97,7 +97,7 @@ ...@@ -97,7 +97,7 @@
(* stride dimension) :init-offset init-offset))) (* stride dimension) :init-offset init-offset)))
(loop for i from 0 below (* stride dimension) by stride (loop for i from 0 below (* stride dimension) by stride
do do
(setf (grid:gref vec i) (setf (grid:aref vec i)
(if (subtypep element-type 'complex) (if (subtypep element-type 'complex)
(coerce (complex (urand) (urand)) element-type) (coerce (complex (urand) (urand)) element-type)
(complex (coerce (urand) element-type))))) (complex (coerce (urand) element-type)))))
...@@ -112,8 +112,8 @@ ...@@ -112,8 +112,8 @@
:init-offset init-offset))) :init-offset init-offset)))
(loop for i below (size complex-vector) by stride (loop for i below (size complex-vector) by stride
do do
(setf (grid:gref real-vector i) (setf (grid:aref real-vector i)
(realpart (grid:gref complex-vector i)))) (realpart (grid:aref complex-vector i))))
real-vector)) real-vector))
(defun copy-with-stride (vector &key (stride 1) init-offset) (defun copy-with-stride (vector &key (stride 1) init-offset)
...@@ -125,7 +125,7 @@ ...@@ -125,7 +125,7 @@
:init-offset init-offset))) :init-offset init-offset)))
(loop for i below (size vector) by stride (loop for i below (size vector) by stride
do do
(setf (grid:gref vec i) (grid:gref vector i))) (setf (grid:aref vec i) (grid:aref vector i)))
vec)) vec))
(defun size-vector-scalar (vector &key (stride 1)) (defun size-vector-scalar (vector &key (stride 1))
...@@ -144,7 +144,7 @@ ...@@ -144,7 +144,7 @@
(loop with length = (size-vector-scalar vector :stride stride) (loop with length = (size-vector-scalar vector :stride stride)
for i from 0 below (size vector) by stride for i from 0 below (size vector) by stride
do do
(setf (grid:gref vector i) (coerce (/ (grid:gref vector i) length) element-type))) (setf (grid:aref vector i) (coerce (/ (grid:aref vector i) length) element-type)))
vector)) vector))
(defun test-real-fft-noise (vector &key (stride 1) non-radix-2) (defun test-real-fft-noise (vector &key (stride 1) non-radix-2)
......
...@@ -76,12 +76,12 @@ ...@@ -76,12 +76,12 @@
(zero-pos (if (evenp n) (- split 1) split))) (zero-pos (if (evenp n) (- split 1) split)))
;; Positive frequencies ;; Positive frequencies
(loop for i from zero-pos to (- n 1) (loop for i from zero-pos to (- n 1)
do (setf (grid:gref shifted (* i stride)) do (setf (grid:aref shifted (* i stride))
(grid:gref vector (* (- i zero-pos) stride)))) (grid:aref vector (* (- i zero-pos) stride))))
;; Negative frequencies ;; Negative frequencies
(loop for i from (+ split 1) to (- n 1) (loop for i from (+ split 1) to (- n 1)
do (setf (grid:gref shifted (* (- i split 1) stride)) do (setf (grid:aref shifted (* (- i split 1) stride))
(grid:gref vector (* i stride)))) (grid:aref vector (* i stride))))
shifted)) shifted))
(export 'fft-inverse-shift) (export 'fft-inverse-shift)
...@@ -96,10 +96,10 @@ ...@@ -96,10 +96,10 @@
(zero-pos (if (evenp n) (- split 1) split))) (zero-pos (if (evenp n) (- split 1) split)))
;; Positive frequencies ;; Positive frequencies
(loop for i from zero-pos to (- n 1) (loop for i from zero-pos to (- n 1)
do (setf (grid:gref inv-shifted (* (- i zero-pos) stride)) do (setf (grid:aref inv-shifted (* (- i zero-pos) stride))
(grid:gref vector (* i stride)))) (grid:aref vector (* i stride))))
;; Negative frequencies ;; Negative frequencies
(loop for i from (+ split 1) to (- n 1) (loop for i from (+ split 1) to (- n 1)
do (setf (grid:gref inv-shifted (* i stride)) do (setf (grid:aref inv-shifted (* i stride))
(grid:gref vector (* (- i split 1) stride)))) (grid:aref vector (* (- i split 1) stride))))
inv-shifted)) inv-shifted))
;; Discrete Hankel Transforms. ;; Discrete Hankel Transforms.
;; Liam Healy, Sat Dec 8 2007 - 16:50 ;; Liam Healy, Sat Dec 8 2007 - 16:50
;; Time-stamp: <2011-01-10 18:16:23EST hankel.lisp> ;; Time-stamp: <2011-05-26 12:37:29EDT hankel.lisp>
;; ;;
;; Copyright 2007, 2008, 2009, 2011 Liam M. Healy ;; Copyright 2007, 2008, 2009, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License ;; Distributed under the terms of the GNU General Public License
...@@ -81,7 +81,7 @@ ...@@ -81,7 +81,7 @@
(let ((hank (make-hankel 128 0.0d0 100.0d0)) (let ((hank (make-hankel 128 0.0d0 100.0d0))
(in (grid:make-foreign-array 'double-float :dimensions 128))) (in (grid:make-foreign-array 'double-float :dimensions 128)))
(loop for n from 0 below 128 (loop for n from 0 below 128
do (setf (grid:gref in n) do (setf (grid:aref in n)
(/ (1+ (expt (sample-x-hankel hank n) 2))))) (/ (1+ (expt (sample-x-hankel hank n) 2)))))
(grid:copy-to (apply-hankel hank in))) (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) ;; Integrate[ x exp(-x) J_1(a x), {x,0,Inf}] = a F(3/2, 2; 2; -a^2)
...@@ -89,13 +89,13 @@ ...@@ -89,13 +89,13 @@
(let ((hank (make-hankel 128 1.0d0 20.0d0)) (let ((hank (make-hankel 128 1.0d0 20.0d0))
(in (grid:make-foreign-array 'double-float :dimensions 128))) (in (grid:make-foreign-array 'double-float :dimensions 128)))
(loop for n from 0 below 128 (loop for n from 0 below 128
do (setf (grid:gref in n) (exp (- (sample-x-hankel hank n))))) do (setf (grid:aref in n) (exp (- (sample-x-hankel hank n)))))
(grid:copy-to (apply-hankel hank in))) (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) ;; 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)) (let ((hank (make-hankel 128 1.0d0 1.0d0))
(in (grid:make-foreign-array 'double-float :dimensions 128))) (in (grid:make-foreign-array 'double-float :dimensions 128)))
(loop for n from 0 below 128 (loop for n from 0 below 128
do (setf (grid:gref in n) do (setf (grid:aref in n)
(let ((x (sample-x-hankel hank n))) (let ((x (sample-x-hankel hank n)))
(* x (- 1 (expt x 2)))))) (* x (- 1 (expt x 2))))))
(grid:copy-to (apply-hankel hank in)))) (grid:copy-to (apply-hankel hank in))))
;; N-tuples ;; N-tuples
;; Liam Healy Sat Feb 3 2007 - 12:53 ;; Liam Healy Sat Feb 3 2007 - 12:53
;; Time-stamp: <2010-06-30 19:57:28EDT ntuple.lisp> ;; Time-stamp: <2011-05-26 12:37:35EDT ntuple.lisp>
;; ;;
;; Copyright 2007, 2008, 2009 Liam M. Healy ;; Copyright 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License ;; Distributed under the terms of the GNU General Public License
...@@ -235,7 +235,7 @@ ...@@ -235,7 +235,7 @@
(close-ntuple ntuple) (close-ntuple ntuple)
;; Check the histogram to see if it matches ntuple-example-values ;; Check the histogram to see if it matches ntuple-example-values
(dotimes (row 100 t) (dotimes (row 100 t)
(unless (= (grid:gref histo row) (aref answer row)) (unless (= (grid:aref histo row) (aref answer row))
(return nil)))))) (return nil))))))
(save-test ntuple (save-test ntuple
......
;; Updating and accessing histogram elements. ;; Updating and accessing histogram elements.
;; Liam Healy, Mon Jan 1 2007 - 14:43 ;; Liam Healy, Mon Jan 1 2007 - 14:43
;; Time-stamp: <2010-11-14 17:46:45EST updating-accessing.lisp> ;; Time-stamp: <2011-05-26 12:37:35EDT updating-accessing.lisp>
;; ;;
;; Copyright 2007, 2008, 2009 Liam M. Healy ;; Copyright 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License ;; Distributed under the terms of the GNU General Public License
...@@ -42,12 +42,12 @@ ...@@ -42,12 +42,12 @@
histograms for a small range of a larger dataset, ignoring the values histograms for a small range of a larger dataset, ignoring the values
outside the range of interest.") outside the range of interest.")
(defmfun grid:gref ((histogram histogram) &rest indices) (defmfun grid:aref ((histogram histogram) &rest indices)
"gsl_histogram_get" "gsl_histogram_get"
(((mpointer histogram) :pointer) ((first indices) sizet)) (((mpointer histogram) :pointer) ((first indices) sizet))
:definition :method :definition :method
:c-return :double :c-return :double
:index grid:gref :index grid:aref
:documentation ; FDL :documentation ; FDL
"Return the contents of the i-th bin of the histogram. "Return the contents of the i-th bin of the histogram.
If i lies outside the valid range of index for the If i lies outside the valid range of index for the
...@@ -128,22 +128,22 @@ ...@@ -128,22 +128,22 @@
(set-ranges-uniform histo 0.0d0 10.0d0) (set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo 2.7d0) (increment histo 2.7d0)
(increment histo 6.9d0 2.0d0) (increment histo 6.9d0 2.0d0)
(grid:gref histo 1)) (grid:aref histo 1))
(let ((histo (make-histogram 10))) (let ((histo (make-histogram 10)))
(set-ranges-uniform histo 0.0d0 10.0d0) (set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo 2.7d0) (increment histo 2.7d0)
(increment histo 6.9d0 2.0d0) (increment histo 6.9d0 2.0d0)
(grid:gref histo 2)) (grid:aref histo 2))
(let ((histo (make-histogram 10))) (let ((histo (make-histogram 10)))
(set-ranges-uniform histo 0.0d0 10.0d0) (set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo 2.7d0) (increment histo 2.7d0)
(increment histo 6.9d0 2.0d0) (increment histo 6.9d0 2.0d0)
(grid:gref histo 6)) (grid:aref histo 6))
(let ((histo (make-histogram 10))) (let ((histo (make-histogram 10)))
(set-ranges-uniform histo 0.0d0 10.0d0) (set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo 2.7d0) (increment histo 2.7d0)
(increment histo 6.9d0 2.0d0) (increment histo 6.9d0 2.0d0)
(grid:gref histo 16)) (grid:aref histo 16))
(let ((histo (make-histogram 10))) (let ((histo (make-histogram 10)))
(set-ranges-uniform histo 0.0d0 10.0d0) (set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo 2.7d0) (increment histo 2.7d0)
......
;; Example spline ;; Example spline
;; Liam Healy, Sat Nov 10 2007 - 21:18 ;; Liam Healy, Sat Nov 10 2007 - 21:18
;; Time-stamp: <2010-08-07 21:37:07EDT spline-example.lisp> ;; Time-stamp: <2011-05-26 12:37:34EDT spline-example.lisp>
;; ;;
;; Copyright 2007, 2008, 2009 Liam M. Healy ;; Copyright 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License ;; Distributed under the terms of the GNU General Public License
...@@ -35,7 +35,7 @@ ...@@ -35,7 +35,7 @@
(loop for i from 0.0d0 below 10.0d0 (loop for i from 0.0d0 below 10.0d0
collect (+ i (cos (expt i 2)))))) collect (+ i (cos (expt i 2))))))
(spline (make-spline +cubic-spline-interpolation+ xarr yarr))) (spline (make-spline +cubic-spline-interpolation+ xarr yarr)))
(loop for xi from (grid:gref xarr 0) below (grid:gref xarr 9) by step (loop for xi from (grid:aref xarr 0) below (grid:aref xarr 9) by step
collect (list xi (evaluate spline xi))))) collect (list xi (evaluate spline xi)))))
(defun evaluate-integral-example (&optional (intervals 4)) (defun evaluate-integral-example (&optional (intervals 4))
......
;; BLAS level 1, Vector operations ;; BLAS level 1, Vector operations
;; Liam Healy, Wed Apr 26 2006 - 15:23 ;; Liam Healy, Wed Apr 26 2006 - 15:23
;; Time-stamp: <2011-04-23 17:04:06EDT blas1.lisp> ;; Time-stamp: <2011-05-26 12:37:34EDT blas1.lisp>
;; ;;
;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy ;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License ;; Distributed under the terms of the GNU General Public License
...@@ -264,7 +264,7 @@ ...@@ -264,7 +264,7 @@
(sines (array-default 8 t)) (sines (array-default 8 t))
(cosines (array-default 8 t))) (cosines (array-default 8 t)))
(loop for i below 8 do (loop for i below 8 do
(setf (grid:gref sines i) (sin (grid:gref angles i))) (setf (grid:aref sines i) (sin (grid:aref angles i)))
(setf (grid:gref cosines i) (cos (grid:gref angles i)))) (setf (grid:aref cosines i) (cos (grid:aref angles i))))
(givens-rotation v1 v2 cosines sines) (givens-rotation v1 v2 cosines sines)
(list (grid:copy-to v1) (grid:copy-to v2)))) (list (grid:copy-to v1) (grid:copy-to v2))))
;; Cholesky Decomposition ;; Cholesky Decomposition
;; Liam Healy, Wed May 3 2006 - 16:38 ;; Liam Healy, Wed May 3 2006 - 16:38
;; Time-stamp: <2011-04-23 17:13:59EDT cholesky.lisp> ;; Time-stamp: <2011-05-26 12:37:33EDT cholesky.lisp>
;; ;;
;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy ;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License ;; Distributed under the terms of the GNU General Public License
...@@ -98,7 +98,7 @@ ...@@ -98,7 +98,7 @@
(let ((decomp (cholesky-decomposition (copy matrix)))) (let ((decomp (cholesky-decomposition (copy matrix))))
(dotimes (row (dim0 matrix) decomp) (dotimes (row (dim0 matrix) decomp)
(loop for col from (1+ row) below (dim1 matrix) do (loop for col from (1+ row) below (dim1 matrix) do
(setf (grid:gref decomp row col) 0.0d0))) (setf (grid:aref decomp row col) 0.0d0)))
(matrix-product decomp decomp nil 1.0d0 0.0d0 :notrans :trans))) (matrix-product decomp decomp nil 1.0d0 0.0d0 :notrans :trans)))
(defun test-cholesky-invert-dim (matrix) (defun test-cholesky-invert-dim (matrix)
......
;; Tridiagonal and Bidiagonal matrices ;; Tridiagonal and Bidiagonal matrices
;; Liam Healy, Thu May 4 2006 - 15:43 ;; Liam Healy, Thu May 4 2006 - 15:43
;; Time-stamp: <2011-01-12 00:49:45EST diagonal.lisp> ;; Time-stamp: <2011-05-26 12:37:34EDT diagonal.lisp>
;; ;;
;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy ;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License ;; Distributed under the terms of the GNU General Public License
...@@ -255,11 +255,11 @@ ...@@ -255,11 +255,11 @@
(subdiag (mvec (1- n) :initial-element 0.5d0))) (subdiag (mvec (1- n) :initial-element 0.5d0)))
(setf (setf
;; i=0 matrix elements ;; i=0 matrix elements
(grid:gref diag 0) 1d0 (grid:aref diag 0) 1d0
(grid:gref superdiag 0) 0d0 (grid:aref superdiag 0) 0d0
(grid:gref b 0) (coerce (expt (1- n) 2) 'double-float) (grid:aref b 0) (coerce (expt (1- n) 2) 'double-float)
;; i=n-1 matrix elements ;; i=n-1 matrix elements
(grid:gref diag (1- n)) 1d0 (grid:aref diag (1- n)) 1d0
(grid:gref subdiag (- n 2)) 0d0 (grid:aref subdiag (- n 2)) 0d0
(grid:gref b (1- n)) 0d0) (grid:aref b (1- n)) 0d0)
(solve-tridiagonal diag superdiag subdiag b x)))) (solve-tridiagonal diag superdiag subdiag b x))))
;; QR decomposition ;; QR decomposition
;; Liam Healy 2008-02-17 11:05:20EST qr.lisp ;; Liam Healy 2008-02-17 11:05:20EST qr.lisp
;; Time-stamp: <2011-04-23 17:14:30EDT qr.lisp> ;; Time-stamp: <2011-05-26 12:37:33EDT qr.lisp>
;; ;;
;; Copyright 2008, 2009, 2011 Liam M. Healy ;; Copyright 2008, 2009, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License ;; Distributed under the terms of the GNU General Public License
...@@ -258,7 +258,7 @@ ...@@ -258,7 +258,7 @@
dim1 nil)) dim1 nil))
(qr1 (qr1
(create-matrix (create-matrix
(lambda (i j) (+ (grid:gref matrix i j) (* (grid:gref u i) (grid:gref v j)))) (lambda (i j) (+ (grid:aref matrix i j) (* (grid:aref u i) (grid:aref v j))))
dim0 dim1)) dim0 dim1))
(qr2 (copy matrix)) (qr2 (copy matrix))
(w (grid:make-foreign-array 'double-float :dimensions dim0 :initial-element 0))) (w (grid:make-foreign-array 'double-float :dimensions dim0 :initial-element 0)))
......
;; ODE system setup ;; ODE system setup
;; Liam Healy, Sun Apr 15 2007 - 14:19 ;; Liam Healy, Sun Apr 15 2007 - 14:19
;; Time-stamp: <2011-05-04 10:25:19EDT ode-system.lisp> ;; Time-stamp: <2011-05-26 12:37:33EDT ode-system.lisp>
;; ;;
;; Copyright 2007, 2008, 2009, 2010, 2011 Liam M. Healy ;; Copyright 2007, 2008, 2009, 2010, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License ;; Distributed under the terms of the GNU General Public License
...@@ -47,11 +47,11 @@ ...@@ -47,11 +47,11 @@
(,ctime (grid:make-foreign-array 'double-float :dimensions 1)) (,ctime (grid:make-foreign-array 'double-float :dimensions 1))
(,cstep (grid:make-foreign-array 'double-float :dimensions 1))) (,cstep (grid:make-foreign-array 'double-float :dimensions 1)))
(symbol-macrolet (symbol-macrolet
((,time (grid:gref ,ctime 0)) ((,time (grid:aref ,ctime 0))
(,step-size (grid:gref ,cstep 0)) (,step-size (grid:aref ,cstep 0))
,@(loop for symb in dependent ,@(loop for symb in dependent
for i from 0 for i from 0
collect `(,symb (grid:gref ,dep ,i)))) collect `(,symb (grid:aref ,dep ,i))))
(flet ((next-step () (flet ((next-step ()
(apply-evolution (apply-evolution
evolve ,ctime ,dep ,cstep control stepperobj ,max-time))) evolve ,ctime ,dep ,cstep control stepperobj ,max-time)))
......
;; Series acceleration. ;; Series acceleration.
;; Liam Healy, Wed Nov 21 2007 - 18:41 ;; Liam Healy, Wed Nov 21 2007 - 18:41
;; Time-stamp: <2011-01-10 17:59:30EST series-acceleration.lisp> ;; Time-stamp: <2011-05-26 12:37:29EDT series-acceleration.lisp>
;; ;;
;; Copyright 2007, 2008, 2009, 2011 Liam M. Healy ;; Copyright 2007, 2008, 2009, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License ;; Distributed under the terms of the GNU General Public License
...@@ -95,8 +95,8 @@ ...@@ -95,8 +95,8 @@
(let ((levin (make-levin maxterms)) (let ((levin (make-levin maxterms))
(array (grid:make-foreign-array 'double-float :dimensions maxterms))) (array (grid:make-foreign-array 'double-float :dimensions maxterms)))
(dotimes (n maxterms) (dotimes (n maxterms)
(setf (grid:gref array n) (coerce (/ (expt (1+ n) 2)) 'double-float)) (setf (grid:aref array n) (coerce (/ (expt (1+ n) 2)) 'double-float))
(incf sum (grid:gref array n))) (incf sum (grid:aref array n)))
(multiple-value-bind (accelerated-sum error) (multiple-value-bind (accelerated-sum error)
(accelerate array levin) (accelerate array levin)
(when print-explanation (when print-explanation
......
;; Linear least squares, or linear regression ;; Linear least squares, or linear regression
;; Liam Healy <2008-01-21 12:41:46EST linear-least-squares.lisp> ;; Liam Healy <2008-01-21 12:41:46EST linear-least-squares.lisp>
;; Time-stamp: <2011-01-10 18:16:18EST linear-least-squares.lisp> ;; Time-stamp: <2011-05-26 12:37:31EDT linear-least-squares.lisp>
;; ;;
;; Copyright 2008, 2009, 2011 Liam M. Healy ;; Copyright 2008, 2009, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License ;; Distributed under the terms of the GNU General Public License
...@@ -315,15 +315,15 @@ ...@@ -315,15 +315,15 @@
(loop for i from 0 below (dim0 x) (loop for i from 0 below (dim0 x)
do do
(format t "data: ~12,5f ~12,5f ~12,5f~&" (format t "data: ~12,5f ~12,5f ~12,5f~&"
(grid:gref x i) (grid:aref x i)
(grid:gref y i) (grid:aref y i)
(/ (grid:gref w i)))) (/ (grid:aref w i))))
(loop for i from -30 below 130 by 10 ; don't print everything (loop for i from -30 below 130 by 10 ; don't print everything
for for
xf = (+ (grid:gref x 0) xf = (+ (grid:aref x 0)
(* (/ i 100) (* (/ i 100)
(- (grid:gref x (1- (dim0 x))) (- (grid:aref x (1- (dim0 x)))