From e48ba0cbb3b7febc2cc1b7c2becb6f06b49cc5a2 Mon Sep 17 00:00:00 2001
From: Liam Healy <liam@thinkpad.local>
Date: Tue, 22 Dec 2009 23:24:49 -0500
Subject: [PATCH] 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.
---
 data/marray.lisp                      |  6 +-
 linear-algebra/matrix-generation.lisp | 85 +++++----------------------
 linear-algebra/qr.lisp                |  8 +--
 3 files changed, 21 insertions(+), 78 deletions(-)

diff --git a/data/marray.lisp b/data/marray.lisp
index 316ea621..21553c26 100644
--- a/data/marray.lisp
+++ b/data/marray.lisp
@@ -1,6 +1,6 @@
 ;; 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))
 
diff --git a/linear-algebra/matrix-generation.lisp b/linear-algebra/matrix-generation.lisp
index a87459de..1f2d3152 100644
--- a/linear-algebra/matrix-generation.lisp
+++ b/linear-algebra/matrix-generation.lisp
@@ -1,6 +1,6 @@
 ;; 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))
-	    row-index-max)
-     do
-     (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)))
-	do
-	(setf (maref matrix i j)
-	      (coerce (if pass-element
-			  (funcall function i j (maref matrix i j))
-			  (funcall function i j))
-		      (element-type matrix)))))
-  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."
-  (set-matrix
-   (make-marray (c-array:cl-single element-type) :dimensions (list dim0 dim1))
-   function
-   pass-element))
-
-(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."
-  (create-matrix
-   (lambda (i j) (if (= i j) (maref vector i) 0))
-   (dim0 vector)))
+  (grid:map-grid
+   :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)))
+  (grid:make-grid
+   `((marray ,dim0 ,dim1) ,element-type)
+   :initial-element constant))
 
 ;;;;****************************************************************************
 ;;;; Specific arrays used in linear algebra tests
@@ -127,13 +70,13 @@
   (create-matrix
    (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)
-      (create-vector
-       (lambda (i) (complex (1+ (* 2 i)) (+ 2 (* 2 i)))) 7 element-type)
-      (create-vector '1+ dim element-type)))
+      (create-matrix
+       (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))
diff --git a/linear-algebra/qr.lisp b/linear-algebra/qr.lisp
index bf463a1f..e2aad0ca 100644
--- a/linear-algebra/qr.lisp
+++ b/linear-algebra/qr.lisp
@@ -1,6 +1,6 @@
 ;; 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))
+	     dim1 nil))
 	 (qr1
 	  (create-matrix
 	   (lambda (i j) (+ (maref matrix i j) (* (maref u i) (maref v j))))
-- 
GitLab