[lisplab-cvs] r2 - src src/matlisp
Jørn Inge Vestgården
jivestgarden at common-lisp.net
Sun Mar 1 10:44:25 UTC 2009
Author: jivestgarden
Date: Sun Mar 1 10:44:24 2009
New Revision: 2
Log:
bugix in geev. Cleaning.
Modified:
src/level3-linalg-interface.lisp
src/matlisp/geev.lisp
start.lisp
Modified: src/level3-linalg-interface.lisp
==============================================================================
--- src/level3-linalg-interface.lisp (original)
+++ src/level3-linalg-interface.lisp Sun Mar 1 10:44:24 2009
@@ -21,67 +21,66 @@
(export '(mtp mtp! mconj mconj! mct mct!
mtr mdet minv! minv
- m* m*! m/ m/!))
+ m* m*! m/ m/!
+ LU-factor LU-factor!
+ lin-solve
+ eigenvalues
+ eigenvectors))
(defgeneric mtp (matrix)
- (:documentation "Matrix transpose"))
+ (:documentation "Matrix transpose."))
(defgeneric mtp! (matrix)
- (:documentation "Matrix transpose. Destructive"))
+ (:documentation "Matrix transpose. Destructive."))
(defgeneric mconj (matrix)
- (:documentation "Matrix conjugate"))
+ (:documentation "Matrix conjugate."))
(defgeneric mconj! (matrix)
- (:documentation "Matrix conjugate. Destructive"))
+ (:documentation "Matrix conjugate. Destructive."))
(defgeneric mct (matrix)
- (:documentation "Matrix conjugate transpose"))
+ (:documentation "Matrix conjugate transpose."))
(defgeneric mct! (matrix)
- (:documentation "Matrix conjugate transpose. Destructive"))
+ (:documentation "Matrix conjugate transpose. Destructive."))
(defgeneric mtr (matrix)
- (:documentation "Matrix trace (sum of diagonal elements)"))
+ (:documentation "Matrix trace (sum of diagonal elements)."))
(defgeneric mdet (matrix)
- (:documentation "Matrix determinant"))
+ (:documentation "Matrix determinant."))
(defgeneric minv! (a)
- (:documentation "Matrix inverse. Destructive"))
+ (:documentation "Matrix inverse. Destructive."))
(defgeneric minv (a)
- (:documentation "Matrix inverse"))
+ (:documentation "Matrix inverse."))
(defgeneric m* (a b)
- (:documentation "Matrix product"))
+ (:documentation "Matrix product."))
(defgeneric m*! (a b)
- (:documentation "Matrix product. Destructive"))
+ (:documentation "Matrix product. Destructive."))
(defgeneric m/ (a b)
- (:documentation "Short for (m* a (minv b))"))
+ (:documentation "Short for (m* a (minv b))."))
(defgeneric m/! (a b)
- (:documentation "Short for (m*! a (minv b)). Destructive"))
+ (:documentation "Short for (m*! a (minv b)). Destructive."))
(defgeneric LU-factor! (matrix pivotes)
(:documentation "LU-factorization with pivoting. Destructive."))
(defgeneric LU-factor (matrix)
(:documentation "LU-factorization with pivoting. Outputs (m p) where
-m is the LU-matrix and p is the pivot permutations"))
-
-#+nil (defgeneric L-solve (L b w/diag)
- (:documentation "Solves the system Lx=b, when L is an upper triangular matrix"))
-
-#+nil (defgeneric U-solve (U b w/diag)
- (:documentation "Solves the system Ux=b, when U is an upper triangular matrix"))
-
-#+nil (defgeneric LU-solve (LU b)
- (:documentation "Solves the system LUx=Pb, when L,U,and P are the outputs from
-the LU-factorization"))
+m is the LU-matrix and p is the pivot permutations."))
(defgeneric lin-solve (A b)
- (:documentation "Solves the linear system of equations Ax=b"))
+ (:documentation "Solves the linear system of equations Ax=b."))
+
+(defgeneric eigenvectors (a)
+ (:documentation "Returns (P d) where P is matrix of right eigenvector and d is a vector of eigenvalues."))
+(defgeneric eigenvalues (a)
+ (:documentation "Returns the vector of eigenvalues."))
Modified: src/matlisp/geev.lisp
==============================================================================
--- src/matlisp/geev.lisp (original)
+++ src/matlisp/geev.lisp Sun Mar 1 10:44:24 2009
@@ -32,20 +32,6 @@
(in-package :lisplab)
-(defgeneric eigenvectors (a)
- (:documentation "Returns (P d) where P is matrix of right eigenvector and d is a vector of eigenvalues"))
-
-(defgeneric eigenvalues (a)
- (:documentation "Returns the vector of eigenvalues"))
-
-(defgeneric diagonalize (a)
- (:documentation "Return (U D V) where A=UDV, and D is diagnoal matrix."))
-
-(defmethod diagonalize ((a blas-real))
- (destructuring-bind (evals vl-mat vr-mat)
- (dgeev (copy a) (create a 0) (create a 0))
- (list vl-mat (diag evals) vr-mat)))
-
(defmethod eigenvectors ((a blas-real))
(destructuring-bind (evals vl-mat vr-mat)
(dgeev (copy a) nil (create a 0))
@@ -57,7 +43,8 @@
evals))
(defgeneric rearrange-eigenvector-matrix (v p))
-(defmethod rearrange-eigenvector-matrix ((v blas-real) (p blas-real))
+
+(defmethod rearrange-eigenvector-matrix (v p)
p)
(defmethod rearrange-eigenvector-matrix ((evals blas-complex) (p blas-real ))
@@ -110,6 +97,7 @@
(ceiling (realpart (aref work 0))))))
(defun dgeev (a vl-mat vr-mat)
+ "Wrapper for f77:dgeev. Potentially destructive."
(let* ((n (rows a))
(xxx (allocate-real-store 1))
(wr (allocate-real-store n))
@@ -140,502 +128,65 @@
(vr-mat2 (rearrange-eigenvector-matrix evec vr-mat)))
(list evec vl-mat2 vr-mat2)))))
-(defmethod eigenvectors ((a blas-complex))
- (destructuring-bind (w p)
- (zgeev-right! (copy a) (create a 0))
- (list p w)))
-
+(defmethod eigenvectors ((a blas-complex))
+ (destructuring-bind (evals vl-mat vr-mat)
+ (zgeev (copy a) nil (create a 0))
+ (list evals vr-mat)))
+
(defmethod eigenvalues ((a blas-complex))
- (car (zgeev-right! (copy a) nil)))
+ (destructuring-bind (evals vl-mat vr-mat)
+ (zgeev (copy a) nil nil)
+ evals))
-(defun zgeev-workspace-size (n vectors?)
+(defun zgeev-workspace-size (n lv? rv?)
;; Ask geev how much space it wants for the work array
(let* ((work (allocate-real-store 1)))
- (multiple-value-bind (store-a store-w store-vl store-vr work info)
- (f77::zgeev
- "N"
- (if vectors? "V" "N")
- n ; N
- work ; A
- n ; LDA
- work ; W
- work ; VL
- 1 ; LDVL
- work ; VR
- (if vectors? n 1) ; LDVR
- work ; WORK
- -1 ; LWORK
- work ; RWORK
- 0 ) ; INFO
- (declare (ignore store-a store-w store-vl store-vr info))
- ;; The desired size in in work[0], which we convert to an
- ;; integer.
- (ceiling (aref work 0)))))
+ (f77::zgeev
+ (if lv? "V" "N")
+ (if rv? "V" "N")
+ n ; N
+ work ; A
+ n ; LDA
+ work ; W
+ work ; VL
+ (if lv? n 1) ; LDVL
+ work ; VR
+ (if rv? n 1) ; LDVR
+ work ; WORK
+ -1 ; LWORK
+ work ; RWORK
+ 0 ) ; INFO
+ ;; The desired size in in work[0], which we convert to an
+ ;; integer.
+ (ceiling (aref work 0))))
-(defun zgeev-right! (a vr-mat)
+(defun zgeev (a vl-mat vr-mat)
+ "Wrapper for f77:zgeev. Potentially destructive."
(let* ((n (rows a))
(2n (* 2 n))
(xxx (allocate-real-store 2))
(w (cnew 0 n 1))
+ (vl (if vl-mat (store vl-mat) xxx))
(vr (if vr-mat (store vr-mat) xxx))
- (lwork (zgeev-workspace-size n (if vr-mat t nil)))
+ (lwork (zgeev-workspace-size n (if vl-mat t nil) (if vr-mat t nil)))
(work (allocate-real-store lwork))
(rwork (allocate-real-store 2n)))
- (multiple-value-bind (store-a store-w store-vl store-vr work info)
- (f77::zgeev
- "N" ;; JOBVL
- (if vr-mat "V" "N") ;; JOBVR
- n ;; N
- (store a) ;; A
- n ;; LDA
- (store w) ;; W
- xxx ;; VL
- 1 ;; LDVL
- vr ;; VR
- (if vr-mat n 1) ;; LDVR
- work ;; WORK
- lwork ;; LWORK
- rwork ;; RWORK
- 0 ) ;; INFO
- (declare (ignore store-a store-w store-vl store-vr work))
- (list w vr-mat)
- )))
-
-#+nil (defun geev-values (a)
- (destructuring-bind (wr wi vr)
- (dgeev-right! (copy a) nil)
- (let* ((n (rows a))
- (wr2 (convert (make-instance 'blas-real :rows n :cols 1 :size n :store wr)
- 'blas-complex))
- (wi2 (convert (make-instance 'blas-real :rows n :cols 1 :size n :store wi)
- 'blas-complex)))
- (.+ wr2 (.* %i wi2)))))
-
-#+nil (defun geev-vectors (a)
- (destructuring-bind (wr wi vr)
- (dgeev-right! (copy a) (create a 0))
- (let* ((n (rows a))
- (wr2 (convert (make-instance 'blas-real :rows n :cols 1 :size n :store wr)
- 'blas-complex))
- (wi2 (convert (make-instance 'blas-real :rows n :cols 1 :size n :store wi)
- 'blas-complex))
- (evals (.+ wr2 (.* %i wi2)))
- (vr2 (make-instance 'blas-real :rows n :cols n :size (* n n) :store vr))
- (evec (cnew 0 n n)))
- (do ((col 0 (incf col)))
- ((>= col n))
- (if (zerop (imagpart (vref evals col )))
- (dotimes (row n)
- (setf (mref evec row col)
- (mref vr2 row col)))
- (progn
- (dotimes (row n)
- (let ((c (complex (mref vr2 row col) (mref vr2 row (1+ col)))))
- (setf (mref evec row col) c
- (mref evec row (1+ col)) (conjugate c))))
- (incf col))))
- (list evals evec))))
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-#+nil (defgeneric geev (a)
- (:documentation ""))
-
-#+nil (defun dgeev-workspace-inquiry (n job)
- ;; Ask geev how much space it wants for the work array
- (multiple-value-bind (jobvl jobvr)
- (case job
- (:nn (values "N" "N"))
- ((:vn t) (values "N" "V"))
- (:nv (values "V" "N"))
- (:vv (values "V" "V")))
- (let* ((work (allocate-real-store 1))
- (ldvr (if (equal jobvr "V") n 1))
- (ldvl (if (equal jobvl "V") n 1)))
- (multiple-value-bind (store-a store-wr store-wi store-vl store-vr
- work info)
- (f77::dgeev jobvl
- jobvr
- n ; N
- work ; A
- n ; LDA
- work ; WR
- work ; WI
- work ; VL
- ldvl ; LDVL
- work ; VR
- ldvr ; LDVR
- work ; WORK
- -1 ; LWORK
- 0 ) ; INFO
- (declare (ignore store-a store-wr store-wi store-vl store-vr))
- (assert (zerop info))
- (ceiling (realpart (aref work 0)))))))
-
-#+nil (defmethod geev ((a blas-real))
- (let* ((n (rows a))
- (a (copy a))
- (xxx (allocate-real-store 1))
- (wr (allocate-real-store n))
- (wi (allocate-real-store n))
- (lwork (dgeev-workspace-inquiry n :nn))
- (work (allocate-real-store lwork)))
- (format t "Now call fortran. Lwork = ~A " lwork)
- (multiple-value-bind (a wr wi vl vr work info)
- (f77::dgeev "N" ; JOBVL
- "N" ; JOBVR
- n ; N
- (store a) ; A
- n ; LDA
- wr ; WR
- wi ; WI
- xxx ; VL
- 1 ; LDVL
- xxx ; VR
- 1 ; LDVR
- work ; WORK
- lwork ; LWORK
- 0 ) ; INFO
- (declare (ignore a work))
- (if (zerop info)
- (values n wr wi vr vl nil nil)
- (values nil nil)))))
-
-
-
-#|
-
-
-(defun geev-fix-up-eigen (n wr wi vr vl left-eig right-eig)
- (let ((res nil)
- ;; Eigenvalues are real unless the max of wi is not zero.
- (real-eig (zerop (aref wi (1- (blas::idamax n wi 1))))))
-
- (when right-eig
- (push (geev-fix-up-eigvec n real-eig wi vr) res))
-
- (if real-eig
- (if (or right-eig left-eig)
- (let ((eigval (make-real-matrix n n)))
- (dotimes (k n)
- (setf (matrix-ref eigval k k) (aref wr k)))
- (push eigval res))
- (let ((eigval (make-real-matrix n 1)))
- (dotimes (k n)
- (setf (matrix-ref eigval k) (aref wr k)))
- (push eigval res)))
- (if (or right-eig left-eig)
- (let ((eigval (make-complex-matrix n n)))
- (dotimes (k n)
- (setf (matrix-ref eigval k k) (complex (aref wr k) (aref wi k))))
- (push eigval res))
- (let ((eigval (make-complex-matrix n 1)))
- (dotimes (k n)
- (setf (matrix-ref eigval k) (complex (aref wr k) (aref wi k))))
- (push eigval res))))
-
- (when left-eig
- (push (geev-fix-up-eigvec n real-eig wi vl) res))
-
- (push t res)
- (values-list (nreverse res))))
-
-
-(let ((work (allocate-real-store 1)))
- (defun dgeev-workspace-inquiry (n job)
- ;; Ask geev how much space it wants for the work array
- (multiple-value-bind (jobvl jobvr)
- (case job
- (:nn (values "N" "N"))
- ((:vn t) (values "N" "V"))
- (:nv (values "V" "N"))
- (:vv (values "V" "V")))
-
- (let* ((ldvr (if (equal jobvr "V") n 1))
- (ldvl (if (equal jobvl "V") n 1)))
-
- (multiple-value-bind (store-a store-wr store-wi store-vl store-vr
- work info)
- (dgeev jobvl
- jobvr
- n ; N
- work ; A
- n ; LDA
- work ; WR
- work ; WI
- work ; VL
- ldvl ; LDVL
- work ; VR
- ldvr ; LDVR
- work ; WORK
- -1 ; LWORK
- 0 ) ; INFO
- (declare (ignore store-a store-wr store-wi store-vl store-vr))
- (assert (zerop info))
- (ceiling (realpart (aref work 0))))))))
-
-
-(defmethod geev ((a real-matrix) &optional (job :NN))
- (let* ((n (nrows a))
- (a (copy a))
- (xxx (allocate-real-store 1))
- (wr (allocate-real-store n))
- (wi (allocate-real-store n))
- (lwork (dgeev-workspace-inquiry n job))
- (work (allocate-real-store lwork)))
-
- (declare (type fixnum n)
- (type (simple-array real-matrix-element-type (*)) xxx wr wi))
-
- (case job
- (:nn
- (multiple-value-bind (a wr wi vl vr work info)
- (dgeev "N" ; JOBVL
- "N" ; JOBVR
- n ; N
- (store a) ; A
- n ; LDA
- wr ; WR
- wi ; WI
- xxx ; VL
- 1 ; LDVL
- xxx ; VR
- 1 ; LDVR
- work ; WORK
- lwork ; LWORK
- 0 ) ; INFO
- (declare (ignore a work))
- (if (zerop info)
- (geev-fix-up-eigen n wr wi vr vl nil nil)
- (values nil nil))))
-
- ((:vn t)
- (let* ((vr (allocate-real-store (* n n))))
- (multiple-value-bind (a wr wi vl vr work info)
- (dgeev "N" ;; JOBVL
- "V" ;; JOBVR
- n ;; N
- (store a) ;; A
- n ;; LDA
- wr ;; WR
- wi ;; WI
- xxx ;; VL
- 1 ;; LDVL
- vr ;; VR
- n ;; LDVR
- work ;; WORK
- lwork ;; LWORK
- 0 ) ;; INFO
- (declare (ignore a work))
- (if (zerop info)
- (geev-fix-up-eigen n wr wi vr vl nil t)
- (values nil nil)))))
-
- (:nv
- (let* ((vl (allocate-real-store (* n n))))
-
- (multiple-value-bind (a wr wi vl vr work info)
- (dgeev "V" ;; JOBVL
- "N" ;; JOBVR
- n ;; N
- (store a) ;; A
- n ;; LDA
- wr ;; WR
- wi ;; WI
- vl ;; VL
- n ;; LDVL
- xxx ;; VR
- 1 ;; LDVR
- work ;; WORK
- lwork ;; LWORK
- 0 ) ;; INFO
- (declare (ignore a work))
- (if (zerop info)
- (geev-fix-up-eigen n wr wi vr vl t nil)
- (values nil nil)))))
-
- (:vv
- (let* ((vl (allocate-real-store (* n n)))
- (vr (allocate-real-store (* n n))))
-
- (multiple-value-bind (a wr wi vl vr work info)
- (dgeev "V" ;; JOBVL
- "V" ;; JOBVR
- n ;; N
- (store a) ;; A
- n ;; LDA
- wr ;; WR
- wi ;; WI
- vl ;; VL
- n ;; LDVL
- vr ;; VR
- n ;; LDVR
- work ;; WORK
- lwork ;; LWORK
- 0 ) ;; INFO
- (declare (ignore a work))
- (if (zerop info)
- (geev-fix-up-eigen n wr wi vr vl t t)
- (values nil nil)))))
-
-
- )))
-
-
-(let ((work (allocate-complex-store 1)))
- (defun zgeev-workspace-inquiry (n job)
- ;; Ask geev how much space it wants for the work array
- (multiple-value-bind (jobvl jobvr)
- (case job
- (:nn (values "N" "N"))
- ((:vn t) (values "N" "V"))
- (:nv (values "V" "N"))
- (:vv (values "V" "V")))
-
- (let* ((ldvr (if (equal jobvr "V") n 1))
- (ldvl (if (equal jobvl "V") n 1)))
-
- (multiple-value-bind (store-a store-w store-vl store-vr work info)
- (zgeev jobvl
- jobvr
- n ; N
- work ; A
- n ; LDA
- work ; W
- work ; VL
- 1 ; LDVL
- work ; VR
- ldvr ; LDVR
- work ; WORK
- -1 ; LWORK
- work ; RWORK
- 0 ) ; INFO
- (declare (ignore store-a store-w store-vl store-vr info))
- ;; The desired size in in work[0], which we convert to an
- ;; integer.
- (ceiling (aref work 0)))))))
-
-;; Hmm, should this really be 4 (5) different methods, one for each
-;; possible value of job?
-
-(defmethod geev ((a complex-matrix) &optional (job :NN))
- (let* ((n (nrows a))
- (a (copy a))
- (w (make-complex-matrix-dim n 1))
- (xxx (allocate-complex-store 1))
- (lwork (zgeev-workspace-inquiry n job))
- (work (allocate-complex-store lwork))
- (rwork (allocate-complex-store n)))
-
- (declare (type fixnum lwork n)
- (type (simple-array complex-matrix-element-type (*)) xxx work)
- (type (simple-array real-matrix-element-type (*)) rwork))
-
- (case job
- (:nn
- (multiple-value-bind (store-a store-w store-vl store-vr work info)
- (zgeev "N" ; JOBVL
- "N" ; JOBVR
- n ; N
- (store a) ; A
- n ; LDA
- (store w) ; W
- xxx ; VL
- 1 ; LDVL
- xxx ; VR
- 1 ; LDVR
- work ; WORK
- lwork ; LWORK
- rwork ; RWORK
- 0 ) ; INFO
- (declare (ignore store-a store-w store-vl store-vr work))
- (if (zerop info)
- (values w t)
- (values nil nil))))
-
- ((:vn t)
- (let* ((vr (make-complex-matrix-dim n n)))
-
- (multiple-value-bind (store-a store-w store-vl store-vr work info)
- (zgeev "N" ;; JOBVL
- "V" ;; JOBVR
- n ;; N
- (store a) ;; A
- n ;; LDA
- (store w) ;; W
- xxx ;; VL
- 1 ;; LDVL
- (store vr) ;; VR
- n ;; LDVR
- work ;; WORK
- lwork ;; LWORK
- rwork ;; RWORK
- 0 ) ;; INFO
- (declare (ignore store-a store-w store-vl store-vr work))
- (if (zerop info)
- (values vr (diag w) t)
- (values nil nil)))))
-
- (:nv
- (let* ((vl (make-complex-matrix-dim n n)))
-
- (multiple-value-bind (store-a store-w store-vl store-vr work info)
- (zgeev "V" ;; JOBVL
- "N" ;; JOBVR
- n ;; N
- (store a) ;; A
- n ;; LDA
- (store w) ;; W
- (store vl) ;; VL
- n ;; LDVL
- xxx ;; VR
- 1 ;; LDVR
- work ;; WORK
- lwork ;; LWORK
- rwork ;; RWORK
- 0 ) ;; INFO
- (declare (ignore store-a store-w store-vl store-vr work))
- (if (zerop info)
- (values (diag w) vl t)
- (values nil nil)))))
-
- (:vv
- (let* ((vr (make-complex-matrix-dim n n))
- (vl (make-complex-matrix-dim n n)))
-
-
- (multiple-value-bind (store-a store-w store-vl store-vr work info)
- (zgeev "V" ;; JOBVL
- "V" ;; JOBVR
- n ;; N
- (store a) ;; A
- n ;; LDA
- (store w) ;; W
- (store vl) ;; VL
- n ;; LDVL
- (store vr) ;; VR
- n ;; LDVR
- work ;; WORK
- lwork ;; LWORK
- rwork ;; RWORK
- 0 ) ;; INFO
- (declare (ignore store-a store-w store-vl store-vr work))
- (if (zerop info)
- (values vr (diag w) vl t)
- (values nil nil)))))
-
-
- )))
-
-|#
\ No newline at end of file
+ ;; TODO get info for error
+ (f77::zgeev
+ (if vl-mat "V" "N") ;; JOBVL
+ (if vr-mat "V" "N") ;; JOBVR
+ n ;; N
+ (store a) ;; A
+ n ;; LDA
+ (store w) ;; W
+ vl ;; VL
+ (if vl-mat n 1) ;; LDVL
+ vr ;; VR
+ (if vr-mat n 1) ;; LDVR
+ work ;; WORK
+ lwork ;; LWORK
+ rwork ;; RWORK
+ 0 ) ;; INFO
+ (let* ((vl-mat2 (rearrange-eigenvector-matrix w vl-mat))
+ (vr-mat2 (rearrange-eigenvector-matrix w vr-mat)))
+ (list w vl-mat2 vr-mat2))))
Modified: start.lisp
==============================================================================
--- start.lisp (original)
+++ start.lisp Sun Mar 1 10:44:24 2009
@@ -5,25 +5,8 @@
*default-pathname-defaults*)
*central-registry*)
-(defparameter *lisplab-external-libraries* '("/usr/lib/atlas/libblas.so.3.0"
- "/usr/lib/atlas/liblapack.so.3.0"))
-
-#+nil (defmethod asdf:perform :after ((op asdf:load-op) (c (eql :lisplab-matlisp)))
- (print "hei")
- (sb-alien:load-shared-object "/usr/lib/atlas/libblas.so.3.0")
- (sb-alien:load-shared-object "/usr/lib/atlas/liblapack.so.3.0"))
+(defparameter *lisplab-external-libraries*
+ '("/usr/lib/atlas/libblas.so.3.0"
+ "/usr/lib/atlas/liblapack.so.3.0"))
(asdf:operate 'asdf:load-op 'lisplab)
-
-#+nil (in-package :ll)
-
-#+nil (run-program "/usr/bin/make" '("-C" "lib-src"))
-
-#+nil (load-libs "lib/")
-
-#+nil (asdf:operate 'asdf:load-op 'brandt)
-
-#+nil (in-dir "fft-sim/"
- (load "region.lisp")
- (load "fft-sim.lisp")
- (load "script.lisp"))
More information about the lisplab-cvs
mailing list