matrix-generation.lisp 4.33 KB
Newer Older
1
2
;; Generate matrices used in tests of linear algebra functions
;; Liam Healy 2009-09-19 18:28:31EDT matrix-generation.lisp
3
;; Time-stamp: <2009-12-22 22:37:22EST matrix-generation.lisp>
4
5
6
7
8
9
10

(in-package :gsl)

;;; These are general to all the linear solver techniques, so more
;;; tests need to be made.  The symbols are not exported because it
;;; assumed they will only be used internally by the test functions.

11
12
;;; When sufficiently stable and featureful, these will be exported.

13
14
;;; See linalg/test.c.

Liam Healy's avatar
Liam Healy committed
15
16
17
18
;;;;****************************************************************************
;;;; General array creation from indices
;;;;****************************************************************************

19
(defun create-matrix
20
21
    (function dim0 &optional (dim1 dim0) (element-type 'double-float))
  "Make a matrix or vector of the specified dimensions, with contents
22
   based on a function of the element indices i, j."
23
24
25
26
  (grid:map-grid
   :source function
   :destination-gtype `((marray ,(if dim1 2 1)) ,element-type)
   :source-dims (if dim1 (list dim0 dim1) (list dim0))))
Liam Healy's avatar
Liam Healy committed
27

28
29
(defun constant-matrix
    (constant dim0 &optional (dim1 dim0) (element-type 'double-float))
30
31
32
  (grid:make-grid
   `((marray ,dim0 ,dim1) ,element-type)
   :initial-element constant))
33

Liam Healy's avatar
Liam Healy committed
34
35
36
;;;;****************************************************************************
;;;; Specific arrays used in linear algebra tests
;;;;****************************************************************************
37
38

(defun create-general-matrix (dim0 dim1)
39
40
41
  (create-matrix
   (lambda (i j) (/ (+ 1 i j)))
   dim0 dim1))
42

43
(defun create-singular-matrix (dim0 dim1)
44
45
46
  (create-matrix
   (lambda (i j) (if (zerop i) 0 (/ (+ 1 i j))))
   dim0 dim1))
47

48
49
50
51
52
53
(defun create-hilbert-matrix (dim)
  "Make Hilbert matrix used to test linear algebra functions."
  (create-general-matrix dim dim))

(defun create-vandermonde-matrix (dim)
  "Make Van der Monde matrix used to test linear algebra functions."
54
55
56
  (create-matrix
   (lambda (i j) (expt (1+ i) (- dim j 1)))
   dim))
57
58

(defun create-moler-matrix (dim)
59
60
61
  (create-matrix
   (lambda (i j) (- (min (1+ i) (1+ j)) 2))
   dim))
62
63
64

(defun create-row-matrix (dim0 dim1)
  ;; This would be better named a column matrix, but they call it a row.
65
66
67
  (create-matrix
   (lambda (i j) (if (zerop j) (/ (1+ i)) 0))
   dim0 dim1))
68
69
70
71
72

(defun create-complex-matrix (dim)
  (create-matrix
   (lambda (i j)
     (complex (/ (+ 1 i j)) (+ 1/2 (expt i 2) (expt j 2))))
73
   dim dim '(complex double-float)))
74

Liam Healy's avatar
Liam Healy committed
75
76
(defun create-rhs-vector (dim &optional (element-type 'double-float))
  (if (subtypep element-type 'complex)
77
78
79
      (create-matrix
       (lambda (i) (complex (1+ (* 2 i)) (+ 2 (* 2 i)))) 7 nil element-type)
      (create-matrix '1+ dim nil element-type)))
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119

(defparameter *hilb2* (create-hilbert-matrix 2))
(defparameter *hilb3* (create-hilbert-matrix 3))
(defparameter *hilb4* (create-hilbert-matrix 4))
(defparameter *hilb12* (create-hilbert-matrix 12))

(defparameter *hilb2-soln*
  (make-marray 'double-float :initial-contents '(-8.0d0 18.0d0)))
(defparameter *hilb3-soln*
  (make-marray 'double-float :initial-contents '(27.0d0 -192.0d0 210.0d0)))
(defparameter *hilb4-soln*
  (make-marray 'double-float
	       :initial-contents '(-64.0d0 900.0d0 -2520.0d0 1820.0d0)))
(defparameter *hilb12-soln*
  (make-marray 'double-float :initial-contents
	       '(-1728.0d0 245388.0d0 -8528520.0d0
		 127026900.0d0 -1009008000.0d0 4768571808.0d0
		 -14202796608.0d0 27336497760.0d0 -33921201600.0d0
		 26189163000.0d0 -11437874448.0d0 2157916488.0d0)))

(defparameter *vander2* (create-vandermonde-matrix 2))
(defparameter *vander3* (create-vandermonde-matrix 3))
(defparameter *vander4* (create-vandermonde-matrix 4))
(defparameter *vander12* (create-vandermonde-matrix 12))

(defparameter *vander2-soln*
  (make-marray 'double-float :initial-contents '(1.0d0 0.0d0)))
(defparameter *vander3-soln*
  (make-marray 'double-float :initial-contents
	       '(0.0d0 1.0d0 0.0d0)))
(defparameter *vander4-soln*
  (make-marray 'double-float :initial-contents
	       '(0.0d0 0.0d0 1.0d0 0.0d0)))
(defparameter *vander12-soln*
  (make-marray 'double-float :initial-contents
	       '(0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0
		 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0)))

(defparameter *m35* (create-general-matrix 3 5))
(defparameter *m53* (create-general-matrix 5 3))
120
121
(defparameter *s35* (create-singular-matrix 3 5))
(defparameter *s53* (create-singular-matrix 5 3))