Changes in / [20:30]


Ignore:
Files:
4 added
22 edited

Legend:

Unmodified
Added
Removed
  • /TODO

    r20 r30  
    11TODO
     2o About arrays. Handle them by making wrappers in the matrix
     3  hierarcy. This will be rather slow but saves work.
     4  And in this way it is possible to avoid the non-speialized methods.     
     5o Matrices with general element type.
     6o Some way to just extend the exact input matrix type
     7o Documentation.
    28o Check if there is some non-threadsafe code in the fortran interface,
    3   i.e. array pre-alocation for the workspaces 
     9  i.e. array pre-alocation for the workspaces.
    410o Find out how to dynamically switch between common-lisp specialized
    511  blas-arrays and fortan specialized blas-arrays. Currently this
    612  is a mess.
    7 o Make test code.
     13o Test code.
    814o Error handling.
    9 o Steal special functions from somewhere. What about f2cl on toms?
    10   Or just take it from Maxima? 
    1115o Added spcialized matrix types, in an ordered way.
     16o Package structure.
    1217
    1318Extensions:
    1419o Symbolic maniputions, similar to Ginac in C++.
    15 o Threaded and paralell execution.
     20o Threaded and paralell execution. Use CUDA?
    1621 
  • /src/core/level0-basic.lisp

    r20 r30  
    2222(in-package :lisplab)
    2323
    24 (export '(*lisplab-print-size* in-dir ))
     24(export '(in-dir ))
    2525
    26 (setf *READ-DEFAULT-FLOAT-FORMAT* 'double-float)
     26(setf *READ-DEFAULT-FLOAT-FORMAT* 'double-float) ; TODO make part of pacakge import?
    2727
    2828(defmacro with-gensyms ((&rest names) . body)
     29  ;; TODO remove? Is it used at all?
    2930  `(let ,(loop for n in names collect `(,n (gensym)))
    3031     ,@body))
     
    4445
    4546(defun strcat (&rest args)
     47  ;; TODO move to the part dealing with files
    4648  (apply #'concatenate (append (list 'string) args)))
    4749
    4850(defmacro in-dir (dir &body body)
     51  ;; TODO move to the part dealing with files
    4952  (let ((path (gensym))
    5053        (dir2 (gensym)))
     
    6063         ,@body))))
    6164
     65(defun to-df (x)
     66  "Coerce x to double float."
     67  (coerce x 'double-float))
     68
     69(defun dvec (n)
     70  "Creates a double vector with n elements."
     71  (make-array n :element-type 'double-float :initial-element 0.0))
  • /src/core/level0-functions.lisp

    r20 r30  
    3434
    3535(defmethod ./= ((a number) (b number) &optional (accuracy))
    36   (apply '.= a b accuracy))
     36  (not (.= a b accuracy)))
    3737
    3838(defmethod .< ((a number) (b number))
     
    6060  (- a b))
    6161
     62
     63
    6264(defmethod .expt ((a number) (b number))
    6365  (expt a b))
     66
     67(defmethod .expt ((a real) (b real))
     68  (expt (to-df a) b))
    6469
    6570(defmethod .sin ((x number))
    6671  (sin x))
    6772
     73(defmethod .sin ((x real))
     74  (sin (to-df x)))
     75
    6876(defmethod .cos ((x number))
    6977  (cos x))
    7078
     79(defmethod .cos ((x real))
     80  (cos (to-df x)))
     81
    7182(defmethod .tan ((x number))
    7283  (tan x))
     84
     85(defmethod .tan ((x real))
     86  (tan (to-df x)))
    7387
    7488(defmethod .log ((x number) &optional (base nil))
     
    7791      (log x)))
    7892
     93(defmethod .log ((x real) &optional (base nil))
     94  (if base
     95      (log (to-df x) base)
     96      (log (to-df x))))
     97
    7998(defmethod .exp ((x number))
    8099  (exp x))
     100
     101(defmethod .exp ((x real))
     102  (exp (to-df x)))
    81103
    82104(defmethod .sinh ((x number))
    83105  (sinh x))
    84106
     107(defmethod .sinh ((x real))
     108  (sinh (to-df x)))
     109
    85110(defmethod .cosh ((x number))
    86111  (cosh x))
     112
     113(defmethod .cosh ((x real))
     114  (cosh (to-df x)))
    87115
    88116(defmethod .tanh ((x number))
    89117  (tanh x))
    90118
     119(defmethod .tanh ((x real))
     120  (tanh (to-df x)))
    91121
    92122
    93123
     124
  • /src/core/level0-interface.lisp

    r20 r30  
    11;;; Lisplab, level0-interface.lisp
    2 ;;; Generic function definitions that also contain
    3 ;;; non-matrix methods.
     2;;; Defines a basic algebra.
     3;;;
    44
    55;;; Copyright (C) 2009 Joern Inge Vestgaarden
     
    2323(export '(copy convert
    2424          scalar?
     25          vector? matrix?
    2526          .abs .imagpart .realpart
    2627          .= ./= .< .<= .> .>=
     
    3940          .gamma))       
    4041
     42(defgeneric scalar? (x)
     43  (:documentation "A scalar is a object with ignored internal structure."))
     44
     45(defgeneric vector? (x)
     46  (:documentation "A vector is a object whose elements are accessible with vref."))
     47
     48(defgeneric matrix? (x)
     49  (:documentation "A matrix is a object whose elements are accesible with mref."))
     50
    4151(defgeneric copy (a)
    4252  (:documentation "Copies the elements and structure, but ignore
     
    4555(defgeneric convert (x type)
    4656  (:documentation "Generalized coerce."))
    47 
    48 (defgeneric scalar? (x)
    49   (:documentation "A scalar is a object with ignored internal structure."))
    5057
    5158(defgeneric .abs (a)
     
    5764(defgeneric .imagpart (a)
    5865  (:documentation "Generialized abs"))
     66
     67;;; Binary boolean operators
    5968
    6069(defgeneric .= (a b &optional (accuracy))
     
    7584(defgeneric .>= (a b)
    7685  (:documentation "Generalized >=." ))
     86
     87;;; Binary operators
    7788
    7889(defgeneric .add (a b)
  • /src/fft/level3-fft-blas.lisp

    r20 r30  
    2020;;; TODO should use the normal ref-blas-complex-store
    2121
     22;;; TODO fix the methods so that they use the actual input matrix type, not just
     23;;; the eql spezializer.
     24
    2225(in-package :lisplab)
    2326
    24 (defmethod fft1 ((x blas-real))
    25   (fft1! (convert x 'blas-complex)))
     27;;;; Real matrices
    2628
    27 (defmethod ifft1 ((x blas-real))
    28   (ifft1! (convert x 'blas-complex)))
     29(defmethod fft1 ((x matrix-lisp-dge))
     30  (fft1! (convert x 'matrix-zge)))
    2931
    30 (defmethod ifft2 ((x blas-real))
    31   (ifft2! (convert x 'blas-complex)))
     32(defmethod ifft1 ((x matrix-lisp-dge))
     33  (ifft1! (convert x 'matrix-zge)))
    3234
    33 (defmethod fft2 ((x blas-real))
    34   (fft2! (convert x 'blas-complex)))
     35(defmethod ifft2 ((x matrix-lisp-dge))
     36  (ifft2! (convert x 'matrix-zge)))
    3537
     38(defmethod fft2 ((x matrix-lisp-dge))
     39  (fft2! (convert x 'matrix-zge)))
    3640
    37 (defmethod fft1! ((x blas-complex))
     41;;; Complex matrices
     42
     43(defmethod fft1 ((x matrix-lisp-zge))
     44  (fft1! (copy x)))
     45
     46(defmethod ifft1 ((x matrix-lisp-zge))
     47  (ifft1! (copy x)))
     48
     49(defmethod ifft2 ((x matrix-lisp-zge))
     50  (ifft2! (copy x)))
     51
     52(defmethod fft2 ((x matrix-lisp-zge))
     53  (fft2! (copy x)))
     54
     55(defmethod fft1! ((x matrix-lisp-zge))
    3856  (dotimes (i (cols x))
    39     (fft-radix-2-blas-complex-store! :f (store x) (rows x) (* (rows x) i) 1))
     57    (fft-radix-2-blas-complex-store! :f (matrix-store x) (rows x) (* (rows x) i) 1))
    4058  x)
    4159
    42 (defmethod ifft1! ((x blas-complex))
     60(defmethod ifft1! ((x matrix-lisp-zge))
    4361  (dotimes (i (cols x))
    44     (fft-radix-2-blas-complex-store! :r (store x) (rows x) (* (rows x) i) 1))
     62    (fft-radix-2-blas-complex-store! :r (matrix-store x) (rows x) (* (rows x) i) 1))
    4563  x)
    4664
    47 (defmethod fft2! ((x blas-complex))
     65(defmethod fft2! ((x matrix-lisp-zge))
    4866  (fft1! x)
    4967  (dotimes (i (rows x))
    50     (fft-radix-2-blas-complex-store! :f (store x) (cols x) i (rows x)))
     68    (fft-radix-2-blas-complex-store! :f (matrix-store x) (cols x) i (rows x)))
    5169  x)
    5270
    53 (defmethod ifft2! ((x blas-complex))
     71(defmethod ifft2! ((x matrix-lisp-zge))
    5472  (ifft1! x)
    5573  (dotimes (i (rows x))
    56     (fft-radix-2-blas-complex-store! :r (store x) (cols x) i (rows x)))
     74    (fft-radix-2-blas-complex-store! :r (matrix-store x) (cols x) i (rows x)))
    5775  x)
    5876
  • /src/linalg/level3-linalg-blas-real.lisp

    r20 r30  
    2020(in-package :lisplab)
    2121
    22 ;; TODO move these optimized functions to library
    23 
    24 (defmacro *df (a b) `(truly-the double-float (* ,a ,b)))
    25 (defmacro /df (a b) `(truly-the double-float (/ ,a ,b)))
    26 (defmacro +df (a b) `(truly-the double-float (+ ,a ,b)))
    27 (defmacro -df (a b) `(truly-the double-float (- ,a ,b)))
    28 
    29 
    3022#+todo (defmethod mtr (matrix)
    3123  (let ((ans 0))
     
    3527
    3628#+todo (defmethod mtp (a)
    37   (let ((b (create a 0 (list (cols a) (rows a)))))
     29  (let ((b (mcreate a 0 (list (cols a) (rows a)))))
    3830    (dotimes (i (rows b))
    3931      (dotimes (j (cols b))
     
    4133    b))
    4234
    43 (defmethod mconj ((a blas-real))
     35(defmethod mconj ((a matrix-lisp-dge))
    4436  (copy a))
    4537
    46 (defmethod mct ((a blas-real))
     38(defmethod mct ((a matrix-lisp-dge))
    4739  (mtp a))
    4840
    49 (defmethod m* ((a blas-real) (b blas-real))
     41(defmethod m* ((a matrix-lisp-dge) (b matrix-lisp-dge))
    5042  (let* ((N (rows a))
    5143         (M (cols b))
    5244         (S (rows b))
    53          (c (create a 0 (list N M)))
    54          (a2 (store a))
    55          (b2 (store b))
    56          (c2 (store c)))
     45         (c (mcreate a 0 (list N M)))
     46         (a2 (matrix-store a))
     47         (b2 (matrix-store b))
     48         (c2 (matrix-store c)))
    5749    (declare (type-blas-store a2 b2 c2)
    5850             (type-blas-idx N M S))
     
    6961      c)))
    7062
    71 (defmethod LU-factor! ((A blas-real) p)
     63(defmethod LU-factor! ((A matrix-lisp-dge) p)
    7264  ;; Translation from GSL.
    7365  ;; Destructive LU factorization. The outout is PA=LU,
     
    7769  (let ((N (rows A))
    7870        (sign 1)
    79         (A2 (store A)))
     71        (A2 (matrix-store A)))
    8072    (declare (type-blas-idx N)
    8173             (fixnum sign)
     
    112104(defun L-solve!-blas-real (L x col)
    113105  ;; Solve Lx=b
    114   (declare (blas-real L x))
    115   (let ((L2 (store L))
    116         (x2 (store x))
     106  (declare (matrix-lisp-dge L x))
     107  (let ((L2 (matrix-store L))
     108        (x2 (matrix-store x))
    117109        (N (rows x)))
    118110    (declare (type-blas-store L2 x2)
     
    129121
    130122(defun U-solve!-blas-real (U x col)
    131   (declare (blas-real U x))
    132   (let* ((U2 (store U))
    133          (x2 (store x))
     123  (declare (matrix-lisp-dge U x))
     124  (let* ((U2 (matrix-store U))
     125         (x2 (matrix-store x))
    134126         (N (rows x))
    135127         (N-1 (1- N))
     
    153145  x)
    154146
    155 (defmethod lin-solve ((A blas-real) (b blas-real))
     147(defmethod lin-solve ((A matrix-lisp-dge) (b matrix-lisp-dge))
    156148  (destructuring-bind (LU pvec sign) (LU-factor A)
    157149    (let ((b2 (copy b)))
     
    171163  A)
    172164     
    173 (defmethod minv! ((A blas-real))
     165(defmethod minv! ((A matrix-lisp-dge))
    174166  (minv!-blas-real A))
    175167
    176 (defmethod minv ((A blas-real))
     168(defmethod minv ((A matrix-lisp-dge))
    177169  (minv! (copy A)))
    178170
  • /src/linalg/level3-linalg-generic.lisp

    r20 r30  
    1818;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
    1919
     20
     21;;; TODO clean up. Move read and write out
     22
    2023(in-package :lisplab)
    2124
     
    2932
    3033(defmethod mtp (a)
    31   (let ((b (create a 0 (list (cols a) (rows a)))))
     34  (let ((b (mcreate a 0 (list (cols a) (rows a)))))
    3235    (dotimes (i (rows b))
    3336      (dotimes (j (cols b))
     
    3639
    3740(defmethod mconj (a)
    38   (let ((b (create a #C(0 0) (list (rows a) (cols a)) )))
     41  (let ((b (mcreate a #C(0 0) (list (rows a) (cols a)) )))
    3942    (dotimes (i (size b))
    4043      (setf (vref b i) (conjugate (vref a i))))
     
    4548
    4649(defmethod m* (a b)
    47   (let ((c (create a 0 (list (rows a) (cols b)))))
     50  (let ((c (mcreate a 0 (list (rows a) (cols b)))))
    4851    (dotimes (i (rows c))
    4952      (dotimes (j (cols c))
     
    198201    (let ((L A)
    199202          (U (copy A))
    200           (Pmat (create A 0)))
     203          (Pmat (mcreate A 0)))
    201204      (w/mat L (x i j) (cond ((> i j) x) ((= i j) 1) (t 0)))
    202205      (w/mat U (x i j) (cond ((<= i j) x) (t 0)))
     
    266269
    267270
    268 
    269 ;;; TRASH
    270 
    271 #+nil (defun apply-permutation (a p)
    272         this is inverse
    273   (let ((b (create a 0)))
    274     (dotimes (row (rows a))
    275       (let ((s (vref p row)))
    276         (dotimes (col (cols a))   
    277           (setf (mref b row col) (mref a s col)))))
    278     b))
    279 
    280 #+nil (defun apply-inverse-permutation (a p)
    281         this is forward
    282   (let ((b (create a 0)))
    283     (dotimes (row (rows a))
    284       (let ((s (vref p row)))
    285         (dotimes (col (cols a))   
    286           (setf (mref b s col) (mref a row col)))))
    287     b))
    288 
    289 
    290 #+nil (defmethod LU-factor! (A p)
    291         ;; Versjon fra boka!
    292   (let* ((N (rows A))
    293          (N-1 (1- N))
    294          (det 1))
    295     (dotimes (i N)
    296       (setf (vref p i) i))
    297     (dotimes (i N-1)
    298       (let ((i-pivot i))
    299         (loop for j from (1+ i) below N do
    300              (when (> (abs (mref A j i))
    301                       (abs (mref A i-pivot i)))
    302                (setf i-pivot j)))
    303         (unless (= i-pivot i)
    304           (let ((tmp (vref p i)))
    305             (setf (vref p i) (vref p i-pivot)
    306                   (vref p i-pivot) tmp
    307                   det (- det)))
    308           (dotimes (j N)
    309             (let ((tmp (mref A i j)))
    310               (setf (mref A i j) (mref A i-pivot j)
    311                     (mref A i-pivot j) tmp)))))
    312       ;; Now reduce all elementsbelow the i'th row
    313       (unless (zerop (mref A i i))
    314         (loop for r from (1+ i) below N do
    315              (print (list 'foerr 'i i 'r r A))
    316              (setf (mref A r i) (./ (mref A r i) (mref A i i)))
    317              (loop for c from (1+ i) below N do
    318                   (setf (mref A r c) (.- (mref A r c)
    319                                          (.* (mref A r i) (mref A i c))))
    320                   (print (list 'mellom 'r r 'c c (mref A r c))))
    321              (print (list 'etter 'i i 'r r A)))))
    322     (list A p det)))
    323          
    324 #+nil (defun tmp-LU-mul (A)
    325   ;; Test code. TODO move and make to an automated test
    326   (destructuring-bind (LU p det)
    327       (LU-factor  A)
    328     (let ((L (create LU 0))
    329           (U (create LU 0)))
    330       (dotimes (i (rows A))
    331         (setf (mref L i i) 1)
    332         (loop for j from 0 below i do
    333              (setf (mref L i j)
    334                    (mref LU i j))))
    335       (dotimes (i (rows A))
    336         #+nil (setf (mref U i i) 1)
    337         (loop for j from i below (cols A) do
    338              (setf (mref U i j)
    339                    (mref LU i j))))
    340       (list A
    341             (apply-inverse-permutation (m* L U) p)         
    342             LU
    343             p L U))))
    344            
    345271     
    346272
  • /src/matlisp/geev.lisp

    r20 r30  
    3333(in-package :lisplab)
    3434
    35 (defmethod eigenvectors ((a blas-real))
     35(defmethod eigenvectors ((a matrix-blas-dge))
    3636   (destructuring-bind (evals vl-mat vr-mat)
    3737       (dgeev (copy a) nil (create a 0))
    3838     (list evals vr-mat)))
    3939
    40 (defmethod eigenvalues ((a blas-real))
     40(defmethod eigenvalues ((a matrix-blas-dge))
    4141  (destructuring-bind (evals vl-mat vr-mat)
    4242      (dgeev (copy a) nil nil)
     
    4848  p)
    4949
    50 (defmethod rearrange-eigenvector-matrix ((evals blas-complex) (p blas-real ))
     50(defmethod rearrange-eigenvector-matrix ((evals matrix-blas-zge) (p matrix-blas-dge))
    5151  (let* ((n (size evals))
    5252         (evec (cnew 0 n n)))
     
    6767(defun combine-ri-vectors (wr wi)
    6868  (let* ((n (size wr))
    69          (wr2 (make-instance 'blas-real :rows n :cols 1 :size n :store wr))
    70          (wi2 (make-instance 'blas-real :rows n :cols 1 :size n :store wi)))
     69         (wr2 (make-instance 'matrix-dge :rows n :cols 1 :store wr))
     70         (wi2 (make-instance 'matrix-dge :rows n :cols 1 :store wi)))
    7171    (if (.= wi2 0)
    7272        wr2
    73         (.+ wr2 (.* %i (convert wi2 'blas-complex))))))
     73        (.+ wr2 (.* %i (convert wi2 'matrix-zge))))))
    7474
    7575(defun dgeev-workspace-size (n lv? rv?)
     
    103103         (wr (allocate-real-store n))
    104104         (wi (allocate-real-store n))
    105          (vl (if vl-mat (store vl-mat) xxx))
    106          (vr (if vr-mat (store vr-mat) xxx))
     105         (vl (if vl-mat (matrix-store vl-mat) xxx))
     106         (vr (if vr-mat (matrix-store vr-mat) xxx))
    107107         (lwork (dgeev-workspace-size n (if vl-mat t nil) (if vr-mat t nil) ))
    108108         (work (allocate-real-store lwork)))
     
    113113         (if vr-mat "V" "N")    ; JOBVR
    114114         n                      ; N
    115          (store a)              ; A
     115         (matrix-store a)       ; A
    116116         n                      ; LDA
    117117         wr                     ; WR
     
    129129        (list evec vl-mat2 vr-mat2)))))
    130130
    131 (defmethod eigenvectors ((a blas-complex))
     131(defmethod eigenvectors ((a matrix-zge))
    132132   (destructuring-bind (evals vl-mat vr-mat)
    133133       (zgeev (copy a) nil (create a 0))
    134134     (list evals vr-mat)))
    135135
    136 (defmethod eigenvalues ((a blas-complex))
     136(defmethod eigenvalues ((a matrix-zge))
    137137  (destructuring-bind (evals vl-mat vr-mat)
    138138      (zgeev (copy a) nil nil)
     
    167167         (xxx (allocate-real-store 2))
    168168         (w  (cnew 0 n 1))
    169          (vl (if vl-mat (store vl-mat) xxx))
    170          (vr (if vr-mat (store vr-mat) xxx))
     169         (vl (if vl-mat (matrix-store vl-mat) xxx))
     170         (vr (if vr-mat (matrix-store vr-mat) xxx))
    171171         (lwork (zgeev-workspace-size n (if vl-mat t nil) (if vr-mat t nil)))
    172172         (work (allocate-real-store lwork))
     
    177177     (if vr-mat "V" "N") ;; JOBVR
    178178     n                   ;; N
    179      (store a)           ;; A
     179     (matrix-store a)    ;; A
    180180     n                   ;; LDA
    181      (store w)           ;; W
     181     (matrix-store w)    ;; W
    182182     vl                  ;; VL
    183183     (if vl-mat n 1)     ;; LDVL
  • /src/matlisp/inv.lisp

    r20 r30  
    2020(in-package :lisplab)
    2121
    22 (defmethod minv! ((a blas-real))
     22(defmethod minv! ((a matrix-blas-dge))
    2323  (let* ((N (rows a))
    2424         (ipiv (make-array N :element-type '(unsigned-byte 32)))
    2525         (msg "argument A given to minv is singular to working machine precision"))
    2626    (multiple-value-bind (_ ipiv info)
    27         (f77::dgetrf N N (store a) N ipiv 0)
     27        (f77::dgetrf N N (matrix-store a) N ipiv 0)
    2828      (declare (ignore _))
    2929      (unless (zerop  info)
     
    3131      (let ((work (make-array N :element-type 'double-float)))
    3232      (multiple-value-bind (_ __ info)   
    33           (f77::dgetri N (store a) N ipiv work N 0)
     33          (f77::dgetri N (matrix-store a) N ipiv work N 0)
    3434        (declare (ignore _ __))
    3535        (unless (zerop info)
     
    3737        a)))))
    3838                       
    39 (defmethod minv ((a blas-real))
     39(defmethod minv ((a matrix-blas-dge))
    4040  (minv! (copy a)))
    4141
    42 (defmethod minv! ((a blas-complex))
     42(defmethod minv! ((a matrix-blas-zge))
    4343  (let* ((N (rows a))
    4444         (ipiv (make-array N :element-type '(unsigned-byte 32)))
    4545         (msg "argument A given to mdiv is singular to working machine precision"))
    4646    (multiple-value-bind (_ ipiv info)
    47         (f77::zgetrf N N (store a) N ipiv 0)
     47        (f77::zgetrf N N (matrix-store a) N ipiv 0)
    4848      (declare (ignore _))
    4949      (unless (zerop  info)
     
    5151      (let ((work (make-array (* 2 N) :element-type 'double-float)))
    5252        (multiple-value-bind (_ __ info)
    53             (f77::zgetri N (store a) N ipiv work N 0)
     53            (f77::zgetri N (matrix-store a) N ipiv work N 0)
    5454          (declare (ignore _ __))
    5555          (unless (zerop info)
     
    5757          a)))))
    5858                       
    59 (defmethod minv ((a blas-complex))
     59(defmethod minv ((a matrix-blas-zge))
    6060  (minv! (copy a)))
    6161
  • /src/matlisp/mul.lisp

    r20 r30  
    2020(in-package :lisplab)
    2121
    22 (defmethod m* ((a blas-real) (b blas-real))
     22(defmethod m* ((a matrix-blas-dge) (b matrix-blas-dge))
    2323  (let* ((m (rows a))
    2424         (n (cols b))
    2525         (k (cols a))
    26          (c (new 'blas-real (list m n) t 0.0)))
    27     (f77::dgemm "N" "N" m n k 1.0 (store a) m (store b) k 0.0 (store c) m)
     26         (c (mcreate a 0 (list m n))))
     27    (f77::dgemm "N" "N" m n k 1.0 (matrix-store a) m (matrix-store b) k 0.0 (matrix-store c) m)
    2828    c))
    2929
    30 (defmethod m* ((a blas-complex) (b blas-complex))
     30(defmethod m* ((a matrix-blas-zge) (b matrix-blas-zge))
    3131  (let* ((m (rows a))
    3232         (n (cols b))
    3333         (k (cols a))
    34          (c (new 'blas-complex (list m n) t 0.0)))
    35     (f77::zgemm "N" "N" m n k #C(1.0 0.0) (store a) m (store b) k #C(0.0 0.0) (store c) m)
     34         (c (mcreate a 0 (list m n))))
     35    (f77::zgemm "N" "N" m n k #C(1.0 0.0) (matrix-store a) m (matrix-store b) k #C(0.0 0.0) (matrix-store c) m)
    3636    c))
  • /src/matrix/level1-array.lisp

    r20 r30  
    2525
    2626(defmethod vector? ((a array))
    27   "True for an array of rank 2"
    28   (= (rank a) 1))
     27  "True for any array through row-major-aref"
     28  t)
    2929
    30 (defmethod copy ((a array))
    31   (if (vector? a)
    32       (copy-seq a)
    33       (let ((y (make-array (dim a) :element-type (element-type a))))
    34         (dotimes (i (size a))
    35           (setf (row-major-aref y i)
    36                 (row-major-aref a i)))
    37         y)))
     30(defmethod dim ((a array) &optional axis)
     31  (if axis
     32      (array-dimension a axis)
     33      (array-dimensions a)))
    3834
    39 (defmethod new ((class (eql 'array)) dim &optional (element-type t) (value 0))
    40   (unless (consp dim) (setf dim (list dim 1)))
    41   (make-array dim
    42               :element-type element-type
    43               :initial-element (convert value element-type)))
     35(defmethod size ((a array))
     36  (reduce #'* (dim a)))
    4437
    45 (defmethod new ((class (eql 'simple-array)) dim &optional (element-type t) (value 0))
    46   (unless (consp dim) (setf dim (list dim 1)))
    47   (make-array dim
    48               :element-type element-type
    49               :initial-element (convert value element-type)))
     38(defmethod rank ((a array))
     39  (array-rank a))
    5040
    51 (defmethod new ((class (eql 'vector)) dim &optional (element-type t) (value 0))
    52   (unless (consp dim) (setf dim (list dim 1)))
    53   (make-array dim
    54               :element-type element-type
    55               :initial-element (convert value element-type)))
     41(defmethod rows ((a array))
     42  (array-dimension a 0))
    5643
    57 (defmethod new ((class (eql 'simple-vector)) dim &optional (element-type t) (value 0))
    58   (unless (consp dim) (setf dim (list dim 1)))
    59   (make-array dim
    60               :element-type element-type
    61               :initial-element (convert value element-type)))
     44(defmethod cols ((a array))
     45  (array-dimension a 1))
    6246
    6347(defmethod element-type ((a array))
     
    7963  (setf (row-major-aref a idx) value))
    8064
    81 (defmethod dim ((a array) &optional axis)
    82   (if axis
    83       (array-dimension a axis)
    84       (array-dimensions a)))
     65(defmethod make-matrix-instance ((x (eql 'array)) dim value)
     66  (make-array dim :initial-element value))
    8567
    86 (defmethod vector? ((a array))
    87   "True for an array of rank 1"
    88   (= (rank a) 1))
     68(defmethod copy ((a array))
     69  ;; TODO move to level2
     70  (if (vectorp a)
     71      (copy-seq a)
     72      (let ((y (make-array (dim a) :element-type (element-type a))))
     73        (dotimes (i (size a))
     74          (setf (row-major-aref y i)
     75                (row-major-aref a i)))
     76        y)))
    8977
    90 (defmethod rank ((a array))
    91   (array-rank a))
     78#+nil (defmethod new ((class (eql 'array)) dim &optional (element-type t) (value 0))
     79  (unless (consp dim) (setf dim (list dim 1)))
     80  (make-array dim
     81              :element-type element-type
     82              :initial-element (convert value element-type)))
    9283
    93 (defmethod rows ((a array))
    94   (array-dimension a 0))
     84#+nil (defmethod new ((class (eql 'simple-array)) dim &optional (element-type t) (value 0))
     85  (unless (consp dim) (setf dim (list dim 1)))
     86  (make-array dim
     87              :element-type element-type
     88              :initial-element (convert value element-type)))
    9589
    96 (defmethod cols ((a array))
    97   (array-dimension a 1))
     90#+nil (defmethod new ((class (eql 'vector)) dim &optional (element-type t) (value 0))
     91  (unless (consp dim) (setf dim (list dim 1)))
     92  (make-array dim
     93              :element-type element-type
     94              :initial-element (convert value element-type)))
    9895
     96#+nil (defmethod new ((class (eql 'simple-vector)) dim &optional (element-type t) (value 0))
     97  (unless (consp dim) (setf dim (list dim 1)))
     98  (make-array dim
     99              :element-type element-type
     100              :initial-element (convert value element-type)))
  • /src/matrix/level1-blas-real.lisp

    r20 r30  
    9898  "Creates a new blas-real matrix"
    9999  (new 'blas-real (list rows cols) t value))
    100 
    101 
    102 
  • /src/matrix/level1-classes.lisp

    r20 r30  
    11;;; Lisplab, level1-classes.lisp
    22;;; Level1, matrix classes
    3 ;;;
    43
    54;;; Copyright (C) 2009 Joern Inge Vestgaarden
     
    1817;;; with this program; if not, write to the Free Software Foundation, Inc.,
    1918;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
     19
     20;;; The class structure is inspired by the stream example
     21;;; in Object-Oriented programming in Common Lisp, by Sonja E. Keene.
     22
    2023
    2124(in-package :lisplab)
     
    4346    :reader element-type)))
    4447
    45 ;;; A way to solve conflicts if there is one foreign and one local implementation
     48;;; A way to solve conflicts if there is one foreign and one native implementation
    4649
    4750(defclass matrix-implementation-base () ())
     
    101104(defclass matrix-lisp-dge (matrix-implementation-lisp matrix-base-dge) ())
    102105
    103 (defclass matrix-dge (matrix-implementation-blas matrix-lisp-dge) ())
     106(defclass matrix-blas-dge (matrix-implementation-blas matrix-lisp-dge) ())
    104107
    105108(defclass matrix-dge (matrix-blas-dge) ()
     
    120123    (setf size (* rows cols))
    121124    (unless matrix-store
    122       ;; Todo: fix initialization!
    123       (setf matrix-store (allocate-real-store (* 2 size) value)))))
     125      (setf matrix-store (allocate-complex-store size value)))))
    124126
    125127(defclass matrix-lisp-zge (matrix-implementation-lisp matrix-base-zge) ())
     
    132134;;; Double float diagonal matrices
    133135
     136;;; TODO
     137
    134138(defclass matrix-base-ddi
    135139    (matrix-structure-diagonal matrix-element-double-float matrix-implementation-base)
     
    138142;;; Complex double float diagonal matrices
    139143
     144;;; TODO
     145
    140146(defclass matrix-base-zdi
    141147    (matrix-structure-diagonal matrix-element-complex-double-float matrix-implementation-base)
    142148  ())
     149
     150;;; Function matrices (matrices without a store)
     151
     152(defclass function-matrix
     153    (matrix-structure-general matrix-element-base matrix-implementation-base)
     154  ((mref
     155    :initarg :mref
     156    :initform (constantly 0)
     157    :accessor function-matrix-mref
     158    :type function)
     159   (set-mref
     160    :initarg :set-mref
     161    :initform (constantly nil)
     162    :accessor function-matrix-set-mref
     163    :type function)
     164   (vref
     165    :initarg :vref
     166    :initform (constantly 0)
     167    :accessor function-matrix-vref
     168    :type function)
     169   (set-vref
     170    :initarg :set-vref
     171    :initform (constantly nil)
     172    :accessor function-matrix-set-vref
     173    :type function))
     174  (:documentation "Matrix without a store."))
    143175
    144176
  • /src/matrix/level1-generic.lisp

    r20 r30  
    1818;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
    1919
     20
     21
     22;;; TODO Get rid of this file and have no non-specialized matrix level1 methods
     23
     24
    2025(in-package :lisplab)
    2126
    22 (defmethod scalar? (x)
    23   (numberp x))
    2427
    25 (defmethod vector? (x)
    26   nil)
    27 
    28 (defmethod matrix? (x)
    29   nil)
    30 
    31 (defmethod size (matrix) (reduce '* (dim matrix)))
    32 
    33 (defmethod rank (matrix) (length (dim matrix)))
    34 
    35 (defmethod cols (matrix) (dim matrix 0))
    36 
    37 (defmethod rows (matrix) (dim matrix 1))
    38 
    39 (defmethod convert (obj type)
    40   (if (not (or (vector? obj) (matrix? obj)))
    41       (coerce obj type)
    42       (let ((new (new type (dim obj) (element-type obj))))
    43         (ecase (rank obj)
    44           (1 (dotimes (i (size obj))
    45                (setf (vref new i) (vref obj i))))
    46           (2 (dotimes (i (rows obj))
    47                (dotimes (j (cols obj))
    48                  (setf (mref new i j) (mref obj i j))))))
    49         new)))
    50 
    51 (defmethod copy (a)
    52   (typecase a
    53     (list (copy-list a))
    54     (sequence (copy-seq a))
    55     (t (let ((b (create a)))
    56          (dotimes (i (size a))
    57            (setf (vref b i) (vref a i)))
    58          b))))
    59 
    60 (defmacro w/mat (a args &body body)
    61   (let ((a2 (gensym))
    62         (x (first args))
    63         (i (second args))
    64         (j (third args)))
    65   `(let ((,a2 ,a))
    66      (dotimes (,i (rows ,a2))
    67        (dotimes (,j (cols ,a2))
    68          (let ((,x (mref ,a2 ,i ,j)))
    69            (setf (mref ,a2 ,i ,j)
    70                  ,@body))))
    71      ,a2)))
    72 
    73 (defmacro mat (type &body args)
    74   "Creates a matrics"
    75   `(convert
    76     ,(cons 'list (mapcar (lambda (x)
    77                            (cons 'list x))
    78                          args))
    79     ,type))
    80 
    81 (defun col (type &rest args)
    82   "Creates a column matrix"
    83   (convert (mapcar 'list args) type))
    84 
    85 (defun row (type &rest args)
    86   "Creates a row matrix"
    87   (convert args type))
    88 
    89 (defmethod create (a &optional value dim)
    90   (unless dim (setf dim (dim a)))
    91   (unless (consp dim) (setf dim (list dim 1)))
    92   (if value
    93     (new (class-name (class-of a)) dim (element-type a) value)
    94     (new (class-name (class-of a)) dim)))
    95 
  • /src/matrix/level1-interface.lisp

    r20 r30  
    2121(in-package :lisplab)
    2222
    23 (export '( *lisplab-print-size*
    24           vector? matrix? new ref mref vref
    25           dim element-type create
    26           size rank rows cols ))
     23(export '(*lisplab-print-size*
     24          make-matrix-instance
     25          ref mref vref
     26          dim element-type
     27          size rank rows cols))
    2728
    2829(defvar *lisplab-print-size* 10 "Suggested number of rows and columns printed to standard output. Not all matrices, such as ordinary lisp arrays, will care about the value.")
    2930
    30 (defgeneric vector? (x)
    31   (:documentation "A vector is a object whose elements are accessible with vref."))
    32 
    33 (defgeneric matrix? (x)
    34   (:documentation "A matrix is a object whose elements are accesible with mref."))
    35 
    36 (defgeneric new (class dim &optional element-type value)
    37   (:documentation "Creates a new matrix filled with numeric arguments."))
     31(defgeneric make-matrix-instance (type dim value)
     32  (:documentation "Creates a new matrix instance"))
    3833
    3934(defgeneric ref (matrix &rest subscripts)
     
    6358(defgeneric (setf element-type) (value matrix))
    6459
    65 (defgeneric create (a &optional value dim)
    66   (:documentation "Creates a new matrix of the same type and with the same value as the other,
    67 but with all elements set to value."))
    68 
    6960(defgeneric size (matrix)
    7061  (:documentation "Gives the number of elements in the object."))
  • /src/matrix/level1-matrix.lisp

    r20 r30  
    1919;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
    2020
     21;;; TODO: clean up
     22
    2123(in-package :lisplab)
    2224
    23 ;;; Generic methods
     25(defmethod scalar? ((x matrix-base)) nil)
     26
     27(defmethod vector? ((x matrix-base)) t)
     28
     29(defmethod matrix? ((x matrix-base)) t)
     30
     31(defmethod rank ((matrix matrix-base)) 2)
    2432
    2533(defmethod dim ((matrix matrix-base) &optional direction)
     
    4351      (when (< rows (rows matrix))
    4452        (format stream "...~%")))))
     53
     54;;;; General cration
     55
     56(defmethod make-matrix-instance ((type symbol) dim value)
     57  (make-instance type :rows (car dim) :cols (cadr dim) :value value))
     58
     59(defmethod make-matrix-instance ((type standard-class) dim value)
     60  (make-instance type :rows (car dim) :cols (cadr dim) :value value))
     61
     62(defmethod make-matrix-instance ((description list) dim value)
     63  (make-matrix-instance (find-matrix-class description) dim value))
     64 
    4565
    4666;;; Spcialized for blas-dge
     
    6686        (the double-float (coerce value 'double-float))))
    6787
     88
    6889;;; Spcialized for blas-zge
    6990
     
    81102
    82103(defmethod vref ((matrix  matrix-base-zge) i)
    83   (ref-blas-complex-store (store matrix) i 0 1))
     104  (ref-blas-complex-store (matrix-store matrix) i 0 1))
    84105
    85106(defmethod (setf vref) (value (matrix  matrix-base-zge) i)
     
    88109  value)
    89110
     111;;; Function matrix
    90112
     113(defmethod mref ((f function-matrix) row col)
     114  (funcall (function-matrix-mref f) f row col))
     115
     116(defmethod (setf mref) (value (f function-matrix) row col)
     117  (funcall (function-matrix-set-mref f) value f row col))
     118
     119(defmethod vref ((f function-matrix) idx)
     120  (funcall (function-matrix-vref f) f idx))
     121
     122(defmethod (setf vref) (value (f function-matrix) idx)
     123  (funcall (function-matrix-set-vref f) value f idx))
     124
  • /src/matrix/level1-util.lisp

    r20 r30  
    2020
    2121(in-package :lisplab)
     22
     23(defun fill-matrix-with-list (m x) 
     24  (let* ((rows (rows m))
     25         (cols (cols m)))
     26    (do ((xx x (cdr xx))
     27         (i 0 (1+ i)))
     28        ((= i rows))
     29      (do ((yy (car xx) (cdr yy))
     30           (j 0 (1+ j)))
     31          ((= j cols))
     32        (setf (mref m i j) (car yy))))
     33    m))
    2234
    2335(deftype type-blas-store ()
     
    119131          (aref store (1+ idx)) (imagpart value))
    120132    value))
     133
     134(defun allocate-complex-store (size &optional (value 0.0))
     135  (let* ((2size (* 2 size))
     136         (rv (coerce (realpart value) 'double-float))
     137         (iv (coerce (imagpart value) 'double-float))
     138         (store (allocate-real-store 2size iv)))
     139    (loop for i from 0 below 2size by 2 do
     140         (setf (aref store i) rv))
     141    store))
     142
  • /src/matrix/level2-blas-complex.lisp

    r20 r30  
    2424                 :store (copy (store matrix))
    2525                 :rows (rows matrix)
    26                  :cols (cols matrix)
    27                  :size (size matrix)))
     26                 :cols (cols matrix)))
    2827
    2928(defmethod convert ((a blas-complex) 'blas-real)
  • /src/matrix/level2-generic.lisp

    r20 r30  
    1818;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
    1919
     20;;; TODO clean it up.
     21
     22;;; TOOD introduce an array wrapper matrix type
     23;;;      and spezialize these methods to matrix-base
     24
    2025(in-package :lisplab)
    2126
    22 (defmethod square-matrix? (x)
    23   (and (matrix? x) (= (rows x) (cols x))))
    24 
    25 (defmethod diag (v)
    26   (let* ((n (size v))
    27          (a (create v 0 (list n n))))
    28     (dotimes (i n)
    29       (setf (mref a i i) (vref v i)))
    30     a))
    31 
    32 (defmethod msum (m)
    33   "Sums all elements of m."
    34   (let ((sum 0))
    35     (dotimes (i (size m))
    36       (incf sum (vref m i)))
    37     sum))
    38 
    39 (defmethod mmax (m)
    40   "Retuns the maximum element of m."
    41   (let ((max (vref m 0)))
    42     (dotimes (i (size m))
    43       (when (> (vref m i) max)
    44         (setf max (vref m i))))
    45     max))
    46 
    47 (defmethod mmin (m)
    48   "Retuns the minimum element of m."
    49   (let ((min (vref m 0)))
    50     (dotimes (i (size m))
    51       (when (< (vref m i) min)
    52         (setf min (vref m i))))
    53     min))
    54 
    55 (defmethod mabsmax (m)
    56   "Retuns the element of m with highes absolute value."
    57   (let ((max (vref m 0)))
    58     (dotimes (i (size m))
    59       (when (> (abs (vref m i)) (abs max))
    60         (setf max (vref m i))))
    61     max))
    62 
    63 (defmethod mabsmin (m)
    64   "Retuns the element of m with smallest absolute value."
    65   (let ((min (vref m 0)))
    66     (dotimes (i (size m))
    67       (when (< (abs (vref m i)) (abs min))
    68         (setf min (vref m i))))
    69     min))
    70 
    71 (defmethod fill! (a val)
    72   "Sets all elemnts of a to val."
     27;; Helper function.       
     28#+nil (defun convert-list-to-matrix (list type)
     29  (let* ((rows (length list))
     30         (cols (length (car list)))
     31         (m (make-matrix-instance type (list rows cols) 0)))
     32    (fill-matrix-with-list m list)))
     33
     34;; Helper function.
     35#+nil (defun convert-matrix-to-matrix (m0 type)
     36  (let* ((rows (rows m0))
     37         (cols (cols m0))
     38         (m (make-matrix-instance type (dim m0) 0)))
     39    (dotimes (i rows)
     40      (dotimes (j cols)
     41        (setf (mref m i j) (mref m0 i j))))
     42    m))
     43
     44(defmethod square-matrix? ((x matrix-base))
     45  (= (rows x) (cols x)))
     46
     47
     48;;; This is OK, but could be optimzied!
     49(defmacro w/mat (a args &body body)
     50  (let ((a2 (gensym))
     51        (x (first args))
     52        (i (second args))
     53        (j (third args)))
     54  `(let ((,a2 ,a))
     55     (dotimes (,i (rows ,a2))
     56       (dotimes (,j (cols ,a2))
     57         (let ((,x (mref ,a2 ,i ,j)))
     58           (setf (mref ,a2 ,i ,j)
     59                 ,@body))))
     60     ,a2)))
     61
     62(defmethod copy-contents ((a matrix-base) (b matrix-base) &optional (converter #'identity))
     63  (dotimes (i (rows a))
     64    (dotimes (j (cols a))
     65      (setf (mref b i j) (funcall converter (mref a i j))))
     66    b))
     67
     68(defmethod .some (pred (a matrix-base) &rest args)
    7369  (dotimes (i (size a))
    74     (setf (vref a i) val))
    75   val)
    76 
    77 (defmethod mmap (type f a &rest args) 
    78   (let ((b (new type (dim a) )))
     70    (when (apply pred (mapcar (lambda (x) (vref x i)) (cons a args)))
     71      (return-from .some t))
     72    nil))
     73
     74(defmethod .every (pred (a matrix-base) &rest args)
     75  (dotimes (i (size a))
     76    (unless (apply pred (mapcar (lambda (x) (vref x i)) (cons a args)))
     77      (return-from .every nil))
     78    t))
     79     
     80(defmethod mmap (type f (a matrix-base) &rest args) 
     81  (let ((b (make-matrix-instance type (dim a) 0)))
    7982    (cond ((not args)
    8083           (dotimes (i (size a))
     
    9194    b))
    9295
    93 (defmethod .map (f a &rest args)
     96(defmethod .map (f (a matrix-base) &rest args)
    9497  (apply #'mmap (class-name (class-of a)) f a args))
    9598
    96 (defmethod convert ((x cons) type)
    97   ;; TODO some better way ... some more general guessing routine
    98   ;; like guess-best-element-type
    99   (if (consp (car x))
    100       (let* ((cols (length (car x)))
    101              (rows (length x))
    102              (m (new type (list rows cols))))
    103         (do ((xx x (cdr xx))
    104              (i 0 (1+ i)))
    105             ((= i rows))
    106           (do ((yy (car xx) (cdr yy))
    107                (j 0 (1+ j)))
    108               ((= j cols))
    109             (setf (mref m i j) (car yy))))
    110         m)
    111       ;; else make a row vector
    112       (convert (list x) type)))
    113 
    114 (defmethod circ-shift (A shift)
    115   (let ((B (create A))   
     99#+todo-remove (defmethod diag (v)
     100  (let* ((n (size v))
     101         (a (mcreate v 0 (list n n))))
     102    (dotimes (i n)
     103      (setf (mref a i i) (vref v i)))
     104    a))
     105
     106(defmethod msum ((m matrix-base))
     107  "Sums all elements of m."
     108  (let ((sum 0))
     109    (dotimes (i (size m))
     110      (setf sum (.+ sum (vref m i))))
     111    sum))
     112
     113(defmethod mmax ((m matrix-base))
     114  "Retuns the maximum element of m."
     115  (let ((max (vref m 0)))
     116    (dotimes (i (size m))
     117      (when (.> (vref m i) max)
     118        (setf max (vref m i))))
     119    max))
     120
     121(defmethod mmin ((m matrix-base))
     122  "Retuns the minimum element of m."
     123  (let ((min (vref m 0)))
     124    (dotimes (i (size m))
     125      (when (.< (vref m i) min)
     126        (setf min (vref m i))))
     127    min))
     128
     129(defmethod mabsmax ((m matrix-base))
     130  "Retuns the element of m with highes absolute value."
     131  (let ((max (vref m 0)))
     132    (dotimes (i (size m))
     133      (when (.> (abs (vref m i)) (abs max))
     134        (setf max (vref m i))))
     135    max))
     136
     137(defmethod mabsmin ((m matrix-base))
     138  "Retuns the element of m with smallest absolute value."
     139  (let ((min (vref m 0)))
     140    (dotimes (i (size m))
     141      (when (.< (abs (vref m i)) (abs min))
     142        (setf min (vref m i))))
     143    min))
     144
     145(defmethod fill! ((a matrix-base) val)
     146  "Sets all elemnts of a to val."
     147  (dotimes (i (size a))
     148    (setf (vref a i) val))
     149  val)
     150
     151(defmethod circ-shift ((A matrix-base) shift)
     152  ;; TODO move to level3
     153  (let ((B (mcreate A)) 
    116154        (rows (rows A))
    117155        (cols (cols A))
     
    124162      B))
    125163
    126 (defmethod pad-shift (A shift &optional (value 0))
    127   (let ((B (create A value))
     164(defmethod pad-shift ((A matrix-base) shift &optional (value 0))
     165  ;; TODO move to level3
     166  (let ((B (mcreate A value))
    128167        (rows (rows A))
    129168        (cols (cols A))
     
    136175      B))
    137176
    138 (defmethod reshape (a shape)
    139   (let ((B (create a 0 shape)))
     177(defmethod reshape ((a matrix-base) shape)
     178  (let ((B (mcreate a 0 shape)))
    140179    (dotimes (i (size B))
    141180      (setf (vref B i) (vref A i)))
    142181    B))
    143182
    144 (defmethod to-vector (a)
     183(defmethod to-vector ((a matrix-base))
    145184  (reshape a (list (size a) 1)))
    146185
    147 (defmethod to-matrix (a rows)
     186(defmethod to-matrix ((a matrix-base) rows)
    148187  (reshape a (list rows (/ (size a) rows) 1)))
    149188
    150189
     190;;;; Basic boolean operators
     191
     192
     193;;;; The boolean operators
     194
     195(defmethod .= ((a matrix-base) (b matrix-base) &optional acc)
     196  (if acc
     197      (.every (lambda (a b) (.= a b acc)) a b)
     198      (.every #'.= a b)))
     199
     200(defmethod .= ((a matrix-base) (b number) &optional acc)
     201  (if acc
     202      (.every (lambda (a) (.= a b acc)) a)
     203      (.every (lambda (a) (.= a b)) a)))
     204
     205(defmethod .= ((a number) (b matrix-base) &optional acc)
     206  (if acc
     207      (.every (lambda (b) (.= a b acc)) b)
     208      (.every (lambda (b) (.= a b)) b)))
     209
     210(defmethod ./= ((a matrix-base) (b matrix-base) &optional acc)
     211  (not (.= a b acc)))
     212
     213(defmethod ./= ((a matrix-base) (b number) &optional acc)
     214  (not (.= a b acc)))
     215
     216(defmethod ./= ((a number) (b matrix-base) &optional acc)
     217  (not (.= a b acc)))
     218
     219(defmacro def-matrix-base-boolean-operator (op)
     220  (let ((a (gensym))
     221        (b (gensym)))
     222    `(progn
     223       (defmethod ,op ((,a matrix-base) (,b matrix-base))
     224         (.every #',op ,a ,b))
     225       (defmethod ,op ((,a matrix-base) (,b number))
     226         (.every (lambda (,a) (,op ,a ,b)) ,a))
     227       (defmethod ,op ((,a number) (,b matrix-base))     
     228         (.every (lambda (,b) (,op ,a ,b)) ,b)))))
     229
     230(def-matrix-base-boolean-operator .<)
     231
     232(def-matrix-base-boolean-operator .<=)
     233
     234(def-matrix-base-boolean-operator .>)
     235
     236(def-matrix-base-boolean-operator .>=)
     237
     238
     239 
     240
     241
     242
     243
    151244
    152245;;; TRASH
     246
     247
     248#+todo-remove(defmethod new (class dim &optional (element-type t) (value 0))
     249  ;;; TODO get rid of this default that calls the new constructor
     250  (mnew class value (car dim) (cadr dim)))
     251
     252#+todo-remove(defmethod convert (obj type)
     253  (if (not (or (vector? obj) (matrix? obj)))
     254      (coerce obj type)
     255      (let ((new (new type (dim obj) (element-type obj))))
     256        (ecase (rank obj)
     257          (1 (dotimes (i (size obj))
     258               (setf (vref new i) (vref obj i))))
     259          (2 (dotimes (i (rows obj))
     260               (dotimes (j (cols obj))
     261                 (setf (mref new i j) (mref obj i j))))))
     262        new)))
     263
     264#+todo-remove(defmethod copy (a)
     265  (typecase a
     266    (list (copy-list a))
     267    (sequence (copy-seq a))
     268    (t (let ((b (create a)))
     269         (dotimes (i (size a))
     270           (setf (vref b i) (vref a i)))
     271         b))))
     272
     273#+todo-remove (defmethod create (a &optional value dim)
     274  (mcreate a value dim))
     275
     276;;; TODO move to dge code
     277
     278#+todo-remove(defmethod convert ((x cons) (type (eql 'matrix-dge)))
     279  (convert-list-to-matrix x type))
     280
     281#+todo-remove(defmethod convert ((x matrix-base) (type (eql 'matrix-dge)))
     282  (convert-matrix-to-matrix x type))
     283
     284#+todo-remove(defmethod mnew ((class (eql 'matrix-dge)) value rows &optional (cols 1))
     285  (make-matrix-instance class (list rows cols) value))
     286
     287;;; TODO move to zge code
     288
     289#+todo-remove(defmethod convert ((x cons) (type (eql 'matrix-zge)))
     290  (convert-list-to-matrix x type))
     291
     292#+todo-remove(defmethod convert ((x matrix-base) (type (eql 'matrix-zge)))
     293  (convert-matrix-to-matrix x type))
     294
     295#+todo-remove(defmethod mnew ((class (eql 'matrix-zge)) value rows &optional (cols 1))
     296  (make-matrix-instance class (list rows cols) value))
  • /src/matrix/level2-interface.lisp

    r20 r30  
    2020(in-package :lisplab)
    2121
    22 (export '(square-matrix?
    23           diag
     22;;; TODO sort and possibly move to other levels
     23
     24(export '(
     25          .every .some ; to level0 ?     
     26          square-matrix?
     27          ; new
     28          mnew
     29          ; create
     30          mcreate
     31          copy-contents
     32          ; diag
    2433          .map mmap fill!         
    2534          dlmwrite dlmread
     
    3746          circ-shift
    3847          pad-shift))
     48
     49(defgeneric .some (pred a &rest matrices)
     50  (:documentation "Generalizes some"))
     51
     52(defgeneric .every (pred a &rest matrices)
     53  (:documentation "Generalizes every."))
     54
     55(defgeneric copy-contents (a b &optional converter)
     56  (:documentation "Copies all elements from a to b."))
     57
     58(defgeneric new (class dim &optional element-type value)
     59  (:documentation "Deprecated. Use mnew in stead. Creates a new matrix filled with numeric arguments."))
     60
     61(defgeneric mnew (class value rows &optional cols)
     62  (:documentation "General matrix constructor. Creates a new matrix filled with numeric arguments."))
     63
     64(defgeneric create (a &optional value dim)
     65  (:documentation "Deprecated. Use mcreate in stead. Creates a new matrix of the same type and with the same value as the other,
     66but with all elements set to value."))
     67
     68(defgeneric mcreate (a &optional value dim)
     69  (:documentation "Creates a new matrix of the same type and with the same value as the other,
     70but with all elements set to value."))
     71
     72(defgeneric mmcreate (a b &optional value dim)
     73  (:documentation "Creates a new matrix. The new matrix has a type derived from a and b,
     74and all elements set to value."))
    3975
    4076(defgeneric square-matrix? (x)
  • /src/specfunc/level0-specfunc.lisp

    r20 r30  
    2020
    2121(in-package :lisplab)
    22 
    23 (defun to-df (x)
    24   (coerce x 'double-float))
    25 
    26 (defun dvec (n)
    27   (make-array n :element-type 'double-float))
    2822
    2923(defmethod .besj (n (x number))
  • /system/lisplab.asd

    r20 r30  
    2424
    2525   ;;
     26   ;; Special functions
     27   ;;
     28   (:module :specfunc
     29    :depends-on (:core)
     30    :pathname "../src/specfunc/"
     31    :serial t
     32    :components
     33    (
     34     (:file "level0-specfunc")))
     35
     36   ;;
    2637   ;; All core matrix stuff (level 1 and 2)
    2738   ;;
    2839   (:module :matrix
    29     :depends-on (:core)
     40    :depends-on (:core :specfunc)
    3041    :pathname "../src/matrix/"
    3142    :serial t
    3243    :components
    3344    (
    34      (:file "level1-interface")
    35      (:file "level1-util")
    36      (:file "level1-generic")
    37      (:file "level1-array")
    38      (:file "level1-list")
    39      (:file "level1-blas")
    40      (:file "level1-blas-real")
    41      (:file "level1-blas-complex")
    42      (:file "level1-funmat")
    43 
    44      (:file "level2-interface")
    45      (:file "level2-array-functions")
    46      (:file "level2-generic")
    47      (:file "level2-funmat")
    48      (:file "level2-blas")
    49      (:file "level2-blas-real")
    50      (:file "level2-blas-complex")))
     45      (:file "level1-interface")
     46
     47      (:file "level1-array")
     48
     49      (:file "level1-util")     
     50      (:file "level1-classes")
     51      (:file "level1-constructors")
     52      (:file "level1-matrix")
     53
     54      (:file "level2-interface")
     55      (:file "level2-constructors")
     56      (:file "level2-generic")
     57      (:file "level2-array-functions")
     58      (:file "level2-matrix-dge")
     59      (:file "level2-matrix-zge")
     60
     61
     62
     63;     (:file "level1-interface")
     64;     (:file "level1-util")
     65;     (:file "level1-generic")
     66;     (:file "level1-array")
     67;     (:file "level1-list")
     68
     69
     70;     (:file "level1-blas")
     71;     (:file "level1-blas-real")
     72;     (:file "level1-blas-complex")
     73;     (:file "level1-funmat")
     74
     75;     (:file "level2-interface")
     76;     (:file "level2-array-functions")
     77;     (:file "level2-generic")
     78;     (:file "level2-funmat")
     79;     (:file "level2-blas")
     80;     (:file "level2-blas-real")
     81;     (:file "level2-blas-complex")
     82      ))
    5183
    5284   ;;
     
    83115    (
    84116     (:file "level3-fft-interface")
    85      (:file "level3-fft-generic")
     117     #+nil (:file "level3-fft-generic")
    86118     (:file "level3-fft-blas")))
    87119
Note: See TracChangeset for help on using the changeset viewer.