matrix.lisp 23.9 KB
Newer Older
1 2
;; Matrices
;; Liam Healy, Sun Mar 26 2006 - 11:51
3
;; Time-stamp: <2008-01-21 11:47:04EST matrix.lisp>
4
;; $Id: $
liam's avatar
liam committed
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

(in-package :gsl)

;;;;****************************************************************************
;;;; Matrix object definition, allocation, reading & writing
;;;;****************************************************************************

;;; GSL-matrix definition
(cffi:defcstruct gsl-matrix-c
  (size1 :size)
  (size2 :size)
  (tda :size)
  (data :pointer)
  (block :pointer)
  (owner :int))

21 22
(defclass gsl-matrix (gsl-data) ())

liam's avatar
liam committed
23
;;; Allocation, freeing, reading and writing
24 25 26 27 28 29 30 31
(defdata "matrix" matrix-double double-float gsl-matrix 2)
(defdata "matrix_float" matrix-single single-float gsl-matrix 2)
(defdata "matrix_int" matrix-fixnum fixnum gsl-matrix 2)
(defdata "matrix_complex" matrix-complex complex gsl-matrix 2)

(defmacro defun-gsl-mdsfc (&rest args)
  "A defun-gsl for matrices of double, single, fixnum, and complex."
  (defun-gsl-all
32
      '(double single fixnum complex)
33 34 35 36 37 38 39 40
      '(:double :float :int gsl-complex)
    "matrix"
    'gsl-matrix
    args))

(defmacro defun-gsl-mdsf (&rest args)
  "A defun-gsl for matrices of double, single, and fixnum."
  (defun-gsl-all
41
      '(double single fixnum)
42 43 44 45
      '(:double :float :int)
    "matrix"
    'gsl-matrix
    args))
liam's avatar
liam committed
46

47 48 49
(defmethod gsl-array ((object gsl-matrix))
  (foreign-slot-value (pointer object) 'gsl-matrix-c 'data))

50 51 52 53 54 55
(export 'matrix-data)
(defun matrix-data (pointer)
  "A pointer to the GSL array with the data contents, from the
   sruct pointer."
  (cffi:foreign-slot-value pointer 'gsl-matrix-c 'data))

liam's avatar
liam committed
56
;;;;****************************************************************************
57
;;;; Getting values
liam's avatar
liam committed
58 59
;;;;****************************************************************************

60
(defun-gsl-mdsfc gsl-aref ((matrix gsl-matrix) &rest indices)
liam's avatar
liam committed
61
  "gsl_matrix_get"
62 63 64 65
  (((pointer matrix) :pointer)
   ((first indices) :size)
   ((second indices) :size))
  :c-return :c-base-type
66
  :documentation "The (i,j)-th element of the matrix.")	;FDL
liam's avatar
liam committed
67

68 69 70 71 72
(defun-gsl mref (pointer index0 index1)
  "gsl_matrix_get"
  ((pointer :pointer) (index0 :size) (index1 :size))
  :c-return :double
  :index nil
73 74 75 76 77 78 79
  :documentation			; FDL
  "An element of the matrix of doubles, computed from the pointer.")

(export 'gsl-matrix-ptr)
(defgeneric gsl-matrix-ptr (matrix i j)
  (:documentation
   "A pointer to the i,j-th element of a matrix."))
80

81 82
(defun-gsl-mdsfc gsl-matrix-ptr ((matrix gsl-matrix) i j)
  "gsl_matrix_ptr" (((pointer matrix) :pointer) (i :size) (j :size))
83
  :c-return :pointer)
84 85 86 87

(defmethod data ((object gsl-matrix) &optional array)
  (let ((arr (or array
		 (make-array (storage-size object)
88
			     :element-type (cl-base-type object)))))
89
    (loop for i from 0
90 91 92 93 94
       below (min (array-dimension arr 0) (first (storage-size object)))
       do
       (loop for j from 0
	  below
	  (min (array-dimension arr 1) (second (storage-size object)))
95
	  do
96
	  (setf (aref arr i j) (gsl-aref object i j))))
97 98 99 100 101 102
    arr))

;;;;****************************************************************************
;;;; Setting values
;;;;****************************************************************************

103
(defun-gsl-mdsfc (setf gsl-aref) (value (matrix gsl-matrix) &rest indices)
liam's avatar
liam committed
104
  "gsl_matrix_set"
105 106 107 108 109
  (((pointer matrix) :pointer)
   ((first indices) :size)
   ((second indices) :size)
   (value :c-base-type))
  :c-return :void
110 111
  :documentation			; FDL
  "Set the (i,j)-th element of the matrix.")
liam's avatar
liam committed
112

113 114 115 116 117
(defun-gsl (setf mref) (value pointer index0 index1)
  "gsl_matrix_set"
  ((pointer :pointer) (index0 :size) (index1 :size) (value :double))
  :c-return :void
  :index nil
118 119
  :documentation			; FDL
  "Set an element of the matrix of doubles, using its pointer.")
120

121 122 123 124 125 126 127 128 129 130
(defmethod (setf data) (array (object gsl-matrix))
  (loop for i from 0
	below
	(min (array-dimension array 0) (first (storage-size object)))
	do
	(loop for j from 0
	      below
	      (min (array-dimension array 1) (second (storage-size object)))
	      do
	      (setf (gsl-aref object i j) (aref array i j)))))
liam's avatar
liam committed
131

132
(defun-gsl-mdsfc set-all ((object gsl-matrix) value)
133
  "gsl_matrix_set_all"
134 135
  (((pointer object) :pointer) (value :c-base-type))
  :c-return :void)
136

137
(defun-gsl-mdsfc set-zero ((object gsl-matrix))
138
  "gsl_matrix_set_zero"
139 140
  (((pointer object) :pointer))
  :c-return :void)
141

142 143 144
(defun-gsl-mdsfc set-identity ((matrix gsl-matrix))
  "gsl_matrix_set_identity" (((pointer matrix) gsl-matrix-c))
  :c-return :void
145
  :documentation			; FDL
146
  "Set the elements of the matrix @var{m} to the
147 148 149
  corresponding elements of the identity matrix, @math{m(i,j) =
  \delta(i,j)}, i.e. a unit diagonal with all off-diagonal elements zero.
  This applies to both square and rectangular matrices.")
liam's avatar
liam committed
150 151

;;;;****************************************************************************
152
;;;; Matrix Views
liam's avatar
liam committed
153 154 155 156 157
;;;;****************************************************************************

(cffi:defcstruct gsl-matrix-view
  (matrix gsl-matrix-c))

158 159 160 161 162 163 164 165 166
(export 'submatrix)
(defgeneric submatrix (matrix k1 k2 n1 n2)
  (:documentation			; FDL
   "A matrix view of a submatrix of the matrix.
   The upper-left element of the submatrix is the element
   (k1, k2) of the original matrix.  The submatrix has n1
   rows and n2 columns.  The physical number of columns in memory
   is unchanged."))

167
(defun-gsl-mdsfc submatrix ((matrix gsl-matrix) k1 k2 n1 n2)
liam's avatar
liam committed
168
  "gsl_matrix_submatrix"
169
  (((pointer matrix) gsl-matrix-c) (k1 :size) (k2 :size) (n1 :size) (n2 :size))
170 171 172 173 174 175 176 177
  :c-return gsl-matrix-view)

(export 'matrix-array)
(defgeneric matrix-array (matrix n1 n2)
  (:documentation			; FDL
   "A matrix view of the array.  The
  matrix has n1 rows and n2 columns.  The physical number of
  columns in memory is also given by n2."))
liam's avatar
liam committed
178

179
(defun-gsl-mdsfc matrix-array ((matrix gsl-matrix) n1 n2)
liam's avatar
liam committed
180
  "gsl_matrix_view_array"
181
  (((pointer matrix) gsl-matrix-c) (n1 :size) (n2 :size))
182 183 184 185 186 187 188 189 190 191
  :c-return gsl-matrix-view)

(export 'matrix-array-tda)
(defgeneric matrix-array-tda (matrix i j tda)
  (:documentation			; FDL
   "A matrix view of the array with a
  physical number of columns tda which may differ from the corresponding
  dimension of the matrix.  The matrix has n1 rows and n2
  columns, and the physical number of columns in memory is given by
  tda."))
liam's avatar
liam committed
192

193
(defun-gsl-mdsfc matrix-array-tda ((matrix gsl-matrix) n1 n2 tda)
liam's avatar
liam committed
194
  "gsl_matrix_view_array_with_tda"
195
  (((pointer matrix) gsl-matrix-c) (n1 :size) (n2 :size) (tda :size))
196
  :c-return gsl-matrix-view)
197 198 199 200 201

(defmacro defun-gsl-mvdsfc (&rest args)
  "A defun-gsl for vectors of double, single, fixnum, and complex,
  translating to a GSL function named matrix_*."
  (defun-gsl-all
202
      '(double single fixnum complex)
203 204 205 206 207
      '(:double :float :int gsl-complex)
    "matrix"
    'gsl-vector
    args))

208 209 210 211 212
(export 'matrix-vector)
(defgeneric matrix-vector (vector n1 n2)
  (:documentation			; FDL
   "A matrix view of the vector.  The matrix
  has n1 rows and n2 columns. The vector must have unit
213
  stride. The physical number of columns in memory is also given by
214
  n2.  Mathematically, the (i,j)-th element of the new
215
  matrix is given by m'(i,j) = v->data[i*n2 + j]
216 217 218 219
  where the index i runs from 0 to n1-1 and the index j
  runs from 0 to n2-1.
  The new matrix is only a view of the vector.  When the view
  goes out of scope the original vector will continue to exist.
220 221
  The original memory can only be deallocated by freeing the original
  vector.  Of course, the original vector should not be deallocated while
222
  the view is still in use."))
223

224 225 226 227 228 229 230 231 232 233
(defun-gsl-mvdsfc matrix-vector ((v gsl-vector) n1 n2)
  "gsl_matrix_view_vector"
  (((pointer v) gsl-vector-c) (n1 :size) (n2 :size))
  :c-return gsl-matrix-view)

(export 'matrix-vector-tda)
(defgeneric matrix-vector-tda (vector n1 n2 tda)
  (:documentation			; FDL
  "A matrix view of the vector with a
  physical number of columns tda which may differ from the
234
  corresponding matrix dimension.  The vector must have unit stride. The
235 236 237
  matrix has n1 rows and n2 columns, and the physical number
  of columns in memory is given by tda.  Mathematically, the
  (i,j)-th element of the new matrix is given by
238
  m'(i,j) = v->data[i*tda + j]
239 240 241 242
  where the index i runs from 0 to n1-1 and the index j
  runs from 0 to n2-1.
  The new matrix is only a view of the vector.  When the view
  goes out of scope the original vector will continue to exist.
243 244
  The original memory can only be deallocated by freeing the original
  vector.  Of course, the original vector should not be deallocated while
245 246 247 248 249 250
  the view is still in use."))

(defun-gsl-mvdsfc matrix-vector-tda ((v gsl-vector) n1 n2 tda)
    "gsl_matrix_view_vector_with_tda"
  (((pointer v) gsl-vector-c) (n1 :size) (n2 :size) (tda :size))
  :c-return gsl-matrix-view)
liam's avatar
liam committed
251

252 253 254
;;;;****************************************************************************
;;;; Row and Column Views
;;;;****************************************************************************
liam's avatar
liam committed
255

256 257 258 259 260
(export 'row-view)
(defgeneric row-view (matrix i)
  (:documentation			; FDL
   "A vector view of the ith row of the matrix."))

261 262 263
(defun-gsl-mdsfc row-view ((matrix gsl-matrix) i) 
  "gsl_matrix_row" (((pointer matrix) gsl-matrix-c) (i :size))
  :c-return gsl-vector-view
264 265 266 267 268 269
  :null-pointer-info (:EFAULT (format nil "index ~d out of range" i)))

(export 'column-view)
(defgeneric column-view (matrix j)
  (:documentation			; FDL
   "A vector view of the jth column of the matrix."))
liam's avatar
liam committed
270

271 272 273
(defun-gsl-mdsfc column-view ((matrix gsl-matrix) j)
  "gsl_matrix_column" ((matrix gsl-matrix-c) (j :size))
  :c-return gsl-matrix-view
274
  :null-pointer-info (:EFAULT (format nil "index ~d out of range" j)))
liam's avatar
liam committed
275

276 277 278 279
(export 'diagonal-view)
(defgeneric diagonal-view (matrix)
  (:documentation			; FDL
   "A vector view of the diagonal of the matrix.
280
   The matrix is not required to be square.
281
   For a rectangular matrix the length of the diagonal is the same as the smaller
282 283 284 285 286 287 288 289 290 291 292 293
   dimension of the matrix."))

(defun-gsl-mdsfc diagonal-view ((matrix gsl-matrix))
  "gsl_matrix_diagonal" ((matrix gsl-matrix-c))
  :c-return gsl-matrix-view)

(export 'subdiagonal-view)
(defgeneric subdiagonal-view (matrix k)
  (:documentation			; FDL
   "A vector view of the kth subdiagonal of the matrix; it is not
   required to be square.  The diagonal of the matrix corresponds to
   k = 0."))
294

295 296
(defun-gsl-mdsfc subdiagonal-view ((matrix gsl-matrix) k)
  "gsl_matrix_subdiagonal" ((matrix gsl-matrix-c) (k :size))
297 298 299 300 301 302 303
  :c-return gsl-matrix-view)

(export 'superdiagonal-view)
(defgeneric superdiagonal-view (matrix k)
  (:documentation			; FDL
  "A vector view of the kth superdiagonal of the matrix; it is not
   required to be square. The diagonal of the matrix corresponds to k = 0."))
304

305
(defun-gsl-mdsfc superdiagonal-view ((matrix gsl-matrix) k)
306
  "gsl_matrix_superdiagonal"
307
  ((matrix gsl-matrix-c) (k :size))
308
  :c-return gsl-matrix-view)
liam's avatar
liam committed
309 310 311 312 313

;;;;****************************************************************************
;;;; Copying
;;;;****************************************************************************

314
(defun-gsl-mdsfc copy ((destination gsl-matrix) (source gsl-matrix))
liam's avatar
liam committed
315
  "gsl_matrix_memcpy"
316
  (((pointer destination) gsl-matrix-c) ((pointer source) gsl-matrix-c))
317
  :invalidate (destination)
318 319 320
  :documentation			; FDL
  "Copy the elements of the matrix source into the
   matrix destination.  The two matrices must have the same size.")
liam's avatar
liam committed
321

322 323 324
(defun-gsl-mdsfc swap ((m1 gsl-matrix) (m2 gsl-matrix))
  "gsl_matrix_swap"
  (((pointer m1) gsl-matrix-c) ((pointer m2) gsl-matrix-c))
325
  :invalidate (m1 m2)
326 327 328
  :documentation			; FDL
  "Exchange the elements of the matrices m1 and
   m2 by copying.  The two matrices must have the same size.")
liam's avatar
liam committed
329 330

;;;;****************************************************************************
331
;;;; Copying rows and columns
liam's avatar
liam committed
332 333
;;;;****************************************************************************

334 335 336 337 338 339 340
(export 'row)
(defgeneric row (vector matrix i)
  (:documentation			; FDL
   "Copy the elements of the ith row of the matrix
   into the vector.  The length of the vector must be the
   same as the length of the row."))

341
(defun-gsl-mdsfc row ((vector gsl-vector) (matrix gsl-matrix) i)
342
  "gsl_matrix_get_row"
343
  (((pointer vector) gsl-vector-c) ((pointer matrix) gsl-matrix-c) (i :size))
344 345 346 347 348 349
  :invalidate (vector))

(export 'column)
(defgeneric column (vector matrix j)
  (:documentation			; FDL
  "Copy the elements of the jth column of the matrix
350
   into the vector.  The length of the vector must be the
351
   same as the length of the column."))
352

353
(defun-gsl-mdsfc column ((vector gsl-vector) (matrix gsl-matrix) j)
354
  "gsl_matrix_get_col"
355
  (((pointer vector) gsl-vector-c) ((pointer matrix) gsl-matrix-c) (j :size))
356 357 358 359 360 361 362 363
  :invalidate (vector))

(export 'set-row)
(defgeneric set-row (vector matrix i)
  (:documentation			; FDL
  "Copy the elements of the vector into the
   ith row of the matrix.  The length of the vector must be
   the same as the length of the row."))
liam's avatar
liam committed
364

365
(defun-gsl-mdsfc set-row ((matrix gsl-matrix) i (vector gsl-vector))
366
  "gsl_matrix_set_row"
367
  (((pointer matrix) gsl-matrix-c) (i :size) ((pointer vector) gsl-vector-c))
368 369 370 371 372 373 374
  :invalidate (matrix))

(export 'set-column)
(defgeneric set-column (vector matrix j)
  (:documentation 			; FDL
   "Copy the elements of the vector into the jth column of the matrix.
  The length of the vector must be the same as the length of the column."))
375

376
(defun-gsl-mdsfc set-column ((matrix gsl-matrix) j (vector gsl-vector))
377
  "gsl_matrix_set_col"
378
  (((pointer matrix) gsl-matrix-c) (j :size) ((pointer vector) gsl-vector-c))
379
  :invalidate (matrix))
380 381 382 383 384 385 386 387

(defun (setf row) (vector matrix i) (set-row matrix i vector))
(defun (setf column) (vector matrix j) (set-column matrix j vector))

;;;;****************************************************************************
;;;; Exchanging rows and columns
;;;;****************************************************************************

388 389 390 391 392
(export 'swap-rows)
(defgeneric swap-rows (matrix i j)
  (:documentation 			; FDL
  "Exchange the ith and jth rows of the matrix in-place."))

393
(defun-gsl-mdsfc swap-rows ((matrix gsl-matrix) i j)
394
  "gsl_matrix_swap_rows"
395
  (((pointer matrix) gsl-matrix-c) (i :size) (j :size))
396 397 398 399 400 401
  :invalidate (matrix))

(export 'swap-columns)
(defgeneric swap-columns (matrix i j)
  (:documentation 			; FDL
  "Exchange the ith and jth columns of the matrix in-place."))
402

403
(defun-gsl-mdsfc swap-columns ((matrix gsl-matrix) i j)
404 405
  "gsl_matrix_swap_columns"
  (((pointer matrix) gsl-matrix-c) (i :size) (j :size))
406 407 408 409 410 411 412 413
  :invalidate (matrix))

(export 'swap-rowcol)
(defgeneric swap-rowcol (matrix i j)
  (:documentation			; FDL
  "Exchange the ith row and jth column of the
   matrix in-place.  The matrix must be square for this operation to
   be possible."))
414

415
(defun-gsl-mdsfc swap-rowcol ((matrix gsl-matrix) i j)
416
  "gsl_matrix_swap_rowcol"
417
  (((pointer matrix) gsl-matrix-c) (i :size) (j :size))
418 419 420 421 422 423 424 425
  :invalidate (matrix))

(export 'matrix-transpose-copy)
(defgeneric matrix-transpose-copy (destination source)
  (:documentation			; FDL
   "Make the destination matrix the transpose of the source matrix
   by copying the elements.  The dimensions of the destination
   matrix must match the transposed dimensions of the source."))
426

427 428
(defun-gsl-mdsfc matrix-transpose-copy
    ((destination gsl-matrix) (source gsl-matrix))
429
  "gsl_matrix_transpose_memcpy"
430
  (((pointer destination) gsl-matrix-c) ((pointer source) gsl-matrix-c))
431 432 433 434 435 436 437 438
  :invalidate (destination))

(export 'matrix-transpose)
(defgeneric matrix-transpose (matrix)
  (:documentation			; FDL
   "Replace the matrix by its transpose by copying the elements
   of the matrix in-place.  The matrix must be square for this
   operation to be possible."))
439

440
(defun-gsl-mdsfc matrix-transpose ((matrix gsl-matrix))
441
  "gsl_matrix_transpose"
442
  (((pointer matrix) gsl-matrix-c))
443
  :invalidate (matrix))
liam's avatar
liam committed
444 445 446 447 448

;;;;****************************************************************************
;;;; Arithmetic operations
;;;;****************************************************************************

449 450
(defun-gsl-mdsfc gsl+ ((a gsl-matrix) (b gsl-matrix))
    "gsl_matrix_add"
451
  (((pointer a) gsl-matrix-c) ((pointer b) gsl-matrix-c))
452
  :invalidate (a)
453 454 455
  :documentation			; FDL
  "Add the elements of b to the elements of a,
   a'_i = a_i + b_i. The two matrices must have the
456
   same dimensions.")
liam's avatar
liam committed
457

458
(defun-gsl-mdsfc gsl- ((a gsl-matrix) (b gsl-matrix))
459
  "gsl_matrix_sub" (((pointer a) gsl-matrix-c) ((pointer b) gsl-matrix-c))
460
  :invalidate (a)
461 462 463
  :documentation			; FDL
  "Subtract the elements of matrix b from the elements of matrix
   a, a'_i = a_i - b_i. The two matrices must have the
464
   same dimensions.")
liam's avatar
liam committed
465

466
(defun-gsl-mdsfc gsl* ((a gsl-matrix) (b gsl-matrix))
467
  "gsl_matrix_mul_elements"
468
  (((pointer a) gsl-matrix-c) ((pointer b) gsl-matrix-c))
469
  :invalidate (a)
470 471 472
  :documentation			; FDL
  "Multiply the elements of matrix a by the elements of
  matrix b, a'(i,j) = a(i,j) * b(i,j). The two matrices must have the
473
  same dimensions.")
liam's avatar
liam committed
474

475
(defun-gsl-mdsfc gsl/ ((a gsl-matrix) (b gsl-matrix))
476 477
  "gsl_matrix_div_elements"
  (((pointer a) gsl-matrix-c) ((pointer b) gsl-matrix-c))
478
  :invalidate (a)
479 480 481
  :documentation			; FDL
  "Divide the elements of matrix a by the elements of
   matrix b, a'(i,j) = a(i,j) / b(i,j). The two matrices must have the
482
   same dimensions.")
liam's avatar
liam committed
483

484
(defun-gsl-mdsfc gsl*c ((a gsl-matrix) x)
485
  "gsl_matrix_scale" (((pointer a) gsl-matrix-c) (x :c-base-type))
486
  :invalidate (a)
487 488 489
  :documentation			; FDL
  "Multiply the elements of matrix a by the constant
  factor x, a'(i,j) = x a(i,j).")
liam's avatar
liam committed
490

491
(defun-gsl-mdsfc gsl+c ((a gsl-matrix) x)
492
  "gsl_matrix_add_constant" (((pointer a) gsl-matrix-c) (x :c-base-type))
493
  :invalidate (a)
494 495 496
  :documentation			; FDL
  "Add the constant value x to the elements of the
  matrix a, a'(i,j) = a(i,j) + x.")
liam's avatar
liam committed
497 498 499 500 501

;;;;****************************************************************************
;;;; Maximum and minimum elements
;;;;****************************************************************************

502
(defun-gsl-mdsf gsl-max ((m gsl-matrix))
503
  "gsl_matrix_max" (((pointer m) gsl-matrix-c))
504 505
  :documentation			; FDL
  "The maximum value in the matrix m."
506
  :c-return :c-base-type)
liam's avatar
liam committed
507

508
(defun-gsl-mdsf gsl-min ((m gsl-matrix))
509
  "gsl_matrix_min" (((pointer m) gsl-matrix-c))
510 511
  :documentation			; FDL
  "The minimum value in the matrix m."
512
  :c-return :c-base-type)
liam's avatar
liam committed
513

514
(defun-gsl-mdsf gsl-minmax ((m gsl-matrix))
liam's avatar
liam committed
515
  "gsl_matrix_minmax"
516
  (((pointer m) gsl-matrix-c) (min :c-base-type) (max :c-base-type))
517 518
  :documentation			; FDL
  "The minimum and maximum values in the matrix m."
519
  :c-return :void)
liam's avatar
liam committed
520

521
(defun-gsl-mdsf gsl-max-index ((m gsl-matrix))
522 523
  "gsl_matrix_max_index"
  (((pointer m) gsl-matrix-c) (imax :size) (jmax :size))
524 525
  :documentation			; FDL
  "The index of the maximum value in the matrix m
526
   When there are several equal maximum elements then the lowest index is
liam's avatar
liam committed
527
   returned."
528 529
  :c-return :void
  :return ((list (size-to-cl imax) (size-to-cl jmax))))
liam's avatar
liam committed
530

531
(defun-gsl-mdsf gsl-min-index ((m gsl-matrix))
532 533
  "gsl_matrix_min_index"
  (((pointer m) gsl-matrix-c) (imin :size) (jmin :size))
534 535
  :documentation			; FDL
  "The index of the minimum value in the matrix m
536 537
  When there are several equal minimum elements then the
  lowest index is returned."
538 539
  :c-return :void
  :return ((list (size-to-cl imin) (size-to-cl jmin))))
liam's avatar
liam committed
540

541
(defun-gsl-mdsf gsl-minmax-index ((m gsl-matrix))
542 543 544
  "gsl_matrix_minmax_index"
  (((pointer m) gsl-matrix-c)
   (imin :size) (jmin :size) (imax :size) (jmax :size))
545 546
  :documentation			; FDL
  "The indices of the minimum and maximum values in the matrix m.
liam's avatar
liam committed
547 548
  When there are several equal minimum elements then the lowest index is
  returned."
549
  :c-return :void
550 551
  :return ((list (size-to-cl imin) (size-to-cl jmin))
	   (list (size-to-cl imax) (size-to-cl jmax))))
liam's avatar
liam committed
552 553 554 555 556

;;;;****************************************************************************
;;;; Properties
;;;;****************************************************************************

557 558
(defun-gsl-mdsfc gsl-zerop ((m gsl-matrix))
  "gsl_matrix_isnull" ((m gsl-matrix-c))
559 560
  :documentation			; FDL
  "All elements of matrix m are zero."
561
  :c-return :boolean)
liam's avatar
liam committed
562

563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************

(defparameter *intmat-1* (make-data 'matrix-fixnum nil 2 2))
(defparameter *intmat-2* (make-data 'matrix-fixnum nil 2 2))
(defparameter *intmatvec* (make-data 'vector-fixnum nil 2))

(lisp-unit:define-test matrix-fixnum
  (lisp-unit:assert-eql			;(setf gsl-aref), gsl-aref
   77
   (progn
     (setf (gsl-aref *intmat-1* 0 1) 77)
     (gsl-aref *intmat-1* 0 1)))
  (lisp-unit:assert-equalp		;(setf data)
   #2A((4 6) (8 2))
   (progn (setf (data *intmat-1*) #2A((4 6) (8 2))) (data *intmat-1*)))
  (lisp-unit:assert-equalp		;set-zero
   #2A((0 0) (0 0))
   (progn (set-zero *intmat-1*) (data *intmat-1*)))
  (lisp-unit:assert-equalp		;set-all
   #2A((44 44)(44 44))
   (progn (set-all *intmat-1* 44) (data *intmat-1*)))
  (lisp-unit:assert-equalp		;set-basis
   #2A((1 0)(0 1))
   (progn (set-identity *intmat-1*) (data *intmat-1*)))
  (lisp-unit:assert-equalp		;row
   #(4 6)
   (progn
     (setf (data *intmat-1*) #2A((4 6) (8 2)))
     (row *intmatvec* *intmat-1* 0)
     (data *intmatvec*)))
  (lisp-unit:assert-equalp		;column
   #(6 2)
   (progn
     (setf (data *intmat-1*) #2A((4 6) (8 2)))
     (column *intmatvec* *intmat-1* 1)
     (data *intmatvec*)))
  (lisp-unit:assert-eql			;gsl-min
   -12
   (progn
     (setf (data *intmat-1*) #2A((-1 -12) (8 3)))
     (gsl-min *intmat-1*)))
  (lisp-unit:assert-eql			;gsl-max
   8
   (progn
     (setf (data *intmat-1*) #2A((-1 -12) (8 3)))
     (gsl-max *intmat-1*)))
  (lisp-unit:assert-equal		;gsl-minmax
   '(-12 8)
   (progn
     (setf (data *intmat-1*) #2A((-1 -12) (8 3)))
     (multiple-value-list (gsl-minmax *intmat-1*))))
  (lisp-unit:assert-equal		;gsl-min-index
   '(0 1)
   (progn
     (setf (data *intmat-1*) #2A((-1 -12) (8 3)))
     (gsl-min-index *intmat-1*)))
  (lisp-unit:assert-equal		;gsl-max-index
   '(1 0)
   (progn
     (setf (data *intmat-1*) #2A((-1 -12) (8 3)))
     (gsl-max-index *intmat-1*)))
  (lisp-unit:assert-equal		;gsl-minmax-index
   '((0 1) (1 0))
   (progn
     (setf (data *intmat-1*) #2A((-1 -12) (8 3)))
     (multiple-value-list (gsl-minmax-index *intmat-1*))))
  (lisp-unit:assert-equalp		;copy
   #2A((1 2)(3 4))
   (progn
     (setf (data *intmat-1*) #2A((1 2)(3 4)))
     (copy *intmat-2* *intmat-1*) (data *intmat-2*)))
  (lisp-unit:assert-equalp		;swap
   #2A((5 6) (7 8))
   (progn
     (setf (data *intmat-1*) #2A((1 2) (3 4))
	   (data *intmat-2*) #2A((5 6) (7 8)))
     (swap *intmat-1* *intmat-2*)
     (data *intmat-1*)))
  (lisp-unit:assert-equalp		;swap-rows
   #2A((3 4) (1 2))
   (progn
     (setf (data *intmat-1*) #2A((1 2) (3 4)))
     (swap-rows *intmat-1* 0 1)
     (data *intmat-1*)))
  (lisp-unit:assert-equalp		;swap-columns
   #2A((2 1) (4 3))
   (progn
     (setf (data *intmat-1*) #2A((1 2) (3 4)))
     (swap-columns *intmat-1* 0 1)
     (data *intmat-1*)))
  (lisp-unit:assert-equalp		;swap-rowcol
   #2A((2 4) (3 1))
   (progn
     (setf (data *intmat-1*) #2A((1 2) (3 4)))
     (swap-rowcol *intmat-1* 0 1)
     (data *intmat-1*))))

liam's avatar
liam committed
662 663 664 665 666 667
;;;;****************************************************************************
;;;; Examples
;;;;****************************************************************************

#|

668
(with-data (mat matrix (10 3))
liam's avatar
liam committed
669 670 671 672 673 674 675
  (loop for i from 0 below 10
	do
	(loop for j from 0 below 3
	      do (setf (gsl-aref mat i j) (+ 0.23d0 j (* 100 i)))))
  (loop for i from 0 below 10
	do
	(loop for j from 0 below 3 do (print (gsl-aref mat i j)))))
liam's avatar
liam committed
676

677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707
(with-data (mat matrix (10 3))
  (loop for i from 0 below 10
	do
	(loop for j from 0 below 3
	      do (setf (gsl-aref mat i j) (+ 0.23d0 j (* 100 i)))))
  (data mat))

#2A((0.23d0 1.23d0 2.23d0)
    (100.23d0 101.23d0 102.23d0)
    (200.23d0 201.23d0 202.23d0)
    (300.23d0 301.23d0 302.23d0)
    (400.23d0 401.23d0 402.23d0)
    (500.23d0 501.23d0 502.23d0)
    (600.23d0 601.23d0 602.23d0)
    (700.23d0 701.23d0 702.23d0)
    (800.23d0 801.23d0 802.23d0)
    (900.23d0 901.23d0 902.23d0))

(with-data (mat matrix (2 2))
  (setf (data mat) #2A((1.0d0 2.0d0) (3.0d0 4.0d0)))
  (with-data (ans matrix (2 2))
    (matrix-copy ans mat)
    (data ans)))
;;; #2A((1.0d0 2.0d0) (3.0d0 4.0d0))

(with-data (mat matrix (2 2))
  (setf (data mat) #2A((1.0d0 2.0d0) (3.0d0 4.0d0)))
  (matrix* mat mat)
  (data mat)))
;;; #2A((1.0d0 4.0d0) (9.0d0 16.0d0))

liam's avatar
liam committed
708
|#