Commit e48ba0cb authored by Liam Healy's avatar Liam Healy
Browse files

Replace create-matrix and constant-matrix using grid functions

Replace create-matrix and constant-matrix by using the grid functions
grid:map-grid and grid:make-grid; eliminate create-vector because
vectors can be generated with create-matrix and using nil for the
second dimension argument.  In the grid:make-grid-data method for
make-marray, the first argument, which is the element type, needs to
be extracted from the rest spec with a newly-changed
grid:spec-scalar-p because compound types like (complex double-float)
were not being recognized as element types.
parent 1a4cb94f
;; A "marray" is an array in both GSL and CL
;; Liam Healy 2008-04-06 21:23:41EDT
;; Time-stamp: <2009-12-21 22:26:21EST marray.lisp>
;; Time-stamp: <2009-12-22 23:17:22EST marray.lisp>
(in-package :gsl)
......@@ -138,8 +138,8 @@
(defmethod grid:make-grid-data
((type (eql 'marray)) dimensions rest-spec
&key (initial-element nil initial-element-p))
(make-marray rest-spec
&key initial-element)
(make-marray (grid:spec-scalar-p rest-spec)
:dimensions dimensions
:initial-element initial-element))
;; Generate matrices used in tests of linear algebra functions
;; Liam Healy 2009-09-19 18:28:31EDT matrix-generation.lisp
;; Time-stamp: <2009-12-06 19:50:10EST matrix-generation.lisp>
;; Time-stamp: <2009-12-22 22:37:22EST matrix-generation.lisp>
(in-package :gsl)
......@@ -16,77 +16,20 @@
;;;; General array creation from indices
(defun set-matrix
(matrix function
&optional pass-element
(row-index-min 0) (row-index-max t) (col-index-min 0) (col-index-max t))
"Set the matrix elements to the values defined by the function.
The function will be called with two arguments if pass-element is nil:
the row index, the column index,
or three arguments if it is non-nil:
the row index, the column index, and the current value of the element.
The range of the indices for rows and columns are defined by the
optional arguments
row-index-min, row-index-max, col-index-min, col-index-max;
if the -max values are T, then the maximum row or column index possible
will be substituted. Either col-index-min or col-index-max or both
can be functions, they will be called with the current row index."
(loop for i
from row-index-min
to (if (eq row-index-max t)
(1- (dim0 matrix))
(loop for j
from (max 0
(typecase col-index-min
(function (funcall col-index-min i))
(integer col-index-min)
(t 0)))
to (min (1- (dim1 matrix))
(typecase col-index-max
(function (funcall col-index-max i))
(integer col-index-max)
(t array-dimension-limit)))
(setf (maref matrix i j)
(coerce (if pass-element
(funcall function i j (maref matrix i j))
(funcall function i j))
(element-type matrix)))))
(defun create-matrix
(function dim0
&optional (dim1 dim0) pass-element (element-type 'double-float))
"Make a matrix of the specified dimensions, with contents
(function dim0 &optional (dim1 dim0) (element-type 'double-float))
"Make a matrix or vector of the specified dimensions, with contents
based on a function of the element indices i, j."
(make-marray (c-array:cl-single element-type) :dimensions (list dim0 dim1))
(defun create-vector
(function dim &optional (element-type 'double-float))
"Make a vector of the specified dimension, with contents
based on a function of the element index."
(let ((vector
(make-marray (c-array:cl-single element-type) :dimensions dim)))
(dotimes (i dim vector)
(setf (maref vector i)
(coerce (funcall function i) element-type)))))
(defun create-diagonal-matrix (vector)
"Place the vector along the diagonal of square matrix."
(lambda (i j) (if (= i j) (maref vector i) 0))
(dim0 vector)))
:source function
:destination-gtype `((marray ,(if dim1 2 1)) ,element-type)
:source-dims (if dim1 (list dim0 dim1) (list dim0))))
(defun constant-matrix
(constant dim0 &optional (dim1 dim0) (element-type 'double-float))
(let ((cst (coerce constant element-type)))
(create-matrix (lambda (i j) (declare (ignore i j)) cst)
dim0 dim1 nil element-type)))
`((marray ,dim0 ,dim1) ,element-type)
:initial-element constant))
;;;; Specific arrays used in linear algebra tests
......@@ -127,13 +70,13 @@
(lambda (i j)
(complex (/ (+ 1 i j)) (+ 1/2 (expt i 2) (expt j 2))))
dim dim nil '(complex double-float)))
dim dim '(complex double-float)))
(defun create-rhs-vector (dim &optional (element-type 'double-float))
(if (subtypep element-type 'complex)
(lambda (i) (complex (1+ (* 2 i)) (+ 2 (* 2 i)))) 7 element-type)
(create-vector '1+ dim element-type)))
(lambda (i) (complex (1+ (* 2 i)) (+ 2 (* 2 i)))) 7 nil element-type)
(create-matrix '1+ dim nil element-type)))
(defparameter *hilb2* (create-hilbert-matrix 2))
(defparameter *hilb3* (create-hilbert-matrix 3))
;; QR decomposition
;; Liam Healy 2008-02-17 11:05:20EST qr.lisp
;; Time-stamp: <2009-09-20 22:26:14EDT qr.lisp>
;; Time-stamp: <2009-12-22 22:34:28EST qr.lisp>
;; $Id$
(in-package :gsl)
......@@ -243,10 +243,10 @@
"Test QR rank-1 update; this should return a matrix with all
elements near zero."
(let* ((dim0 (dim0 matrix)) (dim1 (dim1 matrix))
(u (create-vector (lambda (i) (sin (1+ i))) dim0))
(v (create-vector
(u (create-matrix (lambda (i) (sin (1+ i))) dim0 nil))
(v (create-matrix
(lambda (i) (+ (cos (+ 2 i)) (sin (+ 3 (expt i 2)))))
dim1 nil))
(lambda (i j) (+ (maref matrix i j) (* (maref u i) (maref v j))))
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