Menu
▾
▴
matlisp-commit
|
From: Akshay S. <ak...@us...> - 2013-09-10 20:52:35
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, classy has been updated
via f8b87a620796e228cadb86996b85f4298409ed75 (commit)
via d5f56f654435a06be255e30e2f660360e6920ced (commit)
via 739477d2ea4ae8e582b2355220d502443cd722a3 (commit)
from c6c440e0043ee6633cb729a0bc590e9ca97d5eff (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit f8b87a620796e228cadb86996b85f4298409ed75
Author: Akshay Srinivasan <aks...@gm...>
Date: Tue Sep 10 13:44:24 2013 -0700
Made changes to build on a Mac.
diff --git a/configure b/configure
index 097bf32..5a562fa 100755
--- a/configure
+++ b/configure
@@ -15464,6 +15464,14 @@ else
fi
+#Create libdir if missing (seems to be an issue on the Mac)
+if test -d $libdir; then
+ echo "libdir: Directory present"
+else
+ echo "libdir: Directory not present, creating."
+ mkdir $libdir
+fi;
+
# Check to see if the BLAS library is compatible with matlisp's
# ffi. Basically the same test as above that checks to see if -ff2c
# is needed. We call zdotu which is a Fortran function returning a
diff --git a/configure.ac b/configure.ac
index 6b0aeba..bce467b 100644
--- a/configure.ac
+++ b/configure.ac
@@ -313,6 +313,14 @@ AC_HELP_STRING([--with-external-blas-lapack=libpath], [Location of the BLAS/LAPA
])
AM_CONDITIONAL([EXT_BLAPACK], [test x$ext_blapack = xtrue])
+#Create libdir if missing (seems to be an issue on the Mac)
+if test -d $libdir; then
+ echo "libdir: Directory present"
+else
+ echo "libdir: Directory not present, creating."
+ mkdir $libdir
+fi;
+
# Check to see if the BLAS library is compatible with matlisp's
# ffi. Basically the same test as above that checks to see if -ff2c
# is needed. We call zdotu which is a Fortran function returning a
diff --git a/lib-src/matlisp/Makefile.am b/lib-src/matlisp/Makefile.am
index 63da032..37703b4 100644
--- a/lib-src/matlisp/Makefile.am
+++ b/lib-src/matlisp/Makefile.am
@@ -13,5 +13,4 @@ compat.f\
descal.f\
zescal.f\
dediv.f\
-zediv.f\
-xerbla.f
+zediv.f
diff --git a/lib-src/matlisp/xerbla.f b/lib-src/matlisp/xerbla.f
deleted file mode 100644
index 5de4b64..0000000
--- a/lib-src/matlisp/xerbla.f
+++ /dev/null
@@ -1,41 +0,0 @@
- SUBROUTINE XERBLA(SRNAME,INFO)
-*
-* -- LAPACK auxiliary routine (preliminary version) --
-* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-* November 2006
-*
-* .. Scalar Arguments ..
- INTEGER INFO
- CHARACTER*6 SRNAME
-* ..
-*
-* Purpose
-* =======
-*
-* XERBLA is an error handler for the LAPACK routines.
-* It is called by an LAPACK routine if an input parameter has an
-* invalid value. A message is printed and execution stops.
-*
-* Installers may consider modifying the STOP statement in order to
-* call system-specific exception-handling facilities.
-*
-* Arguments
-* =========
-*
-* SRNAME (input) CHARACTER*6
-* The name of the routine which called XERBLA.
-*
-* INFO (input) INTEGER
-* The position of the invalid parameter in the parameter list
-* of the calling routine.
-*
-*
- WRITE (*,FMT=9999) SRNAME,INFO
-*
-*
- 9999 FORMAT (' ** On entry to ',A6,' parameter number ',I2,' had ',
- + 'an illegal value')
-*
-* End of XERBLA
-*
- END
diff --git a/src/foreign-core/lazy-loader.lisp.in b/src/foreign-core/lazy-loader.lisp.in
index f7e5b4e..025c65e 100644
--- a/src/foreign-core/lazy-loader.lisp.in
+++ b/src/foreign-core/lazy-loader.lisp.in
@@ -72,10 +72,10 @@
(progn
(push "@BLAS_LAPACK_DIR@" cffi:*foreign-library-directories*)
(cffi:define-foreign-library blas
- (:darwin "libblas.dylib")
+ (:darwin (:or "libBLAS.dylib" "libblas.dylib"))
(t (:default "@BLAS_LAPACK_DIR@/libblas")))
(cffi:define-foreign-library lapack
- (:darwin "liblapack.dylib")
+ (:darwin (:or "libLAPACK.dylib" "liblapack.dylib"))
(t (:default "@BLAS_LAPACK_DIR@/liblapack"))))
(progn
;; Use our blas and lapack libraries
@@ -88,8 +88,8 @@
(defun load-blas-&-lapack-libraries ()
;; Load the additional matlisp libraries
- (cffi:use-foreign-library matlisp)
(cffi:use-foreign-library blas)
+ (cffi:use-foreign-library matlisp)
(cffi:use-foreign-library lapack)
(cffi:use-foreign-library dfftpack)
(cffi:use-foreign-library toms715)
diff --git a/src/lapack/eig.lisp b/src/lapack/eig.lisp
index 94bee6a..c603342 100644
--- a/src/lapack/eig.lisp
+++ b/src/lapack/eig.lisp
@@ -191,7 +191,7 @@
(values-list (remove-if #'null
(list
(let ((*check-after-initializing?* nil))
- (make-instance 'complex-tensor ;',(complexified-type cla)
+ (make-instance ',(complexified-type cla)
:dimensions (make-index-store (list (nrows A)))
:strides (make-index-store (list 1))
:head 0
@@ -213,10 +213,10 @@
(vr (when revec? (zeros (list n n) (class-of matrix)))))
(geev! (copy matrix) vl vr)))
+
(defun geev-fix-up-eigvec (n eigval eigvec)
- (declare (type complex-tensor eigval)
- (type real-tensor eigvec))
- (let* ((evec (copy! eigvec (zeros (list n n) 'complex-tensor)))
+ (let* ((evec (copy! eigvec (zeros (list n n) (complexified-type (class-of eigvec)))))
+ (tmp (zeros n (complexified-type (class-of eigvec))))
(cviewa (col-slice~ evec 0))
(cviewb (col-slice~ evec 0))
(cst (aref (strides evec) 1)))
@@ -228,9 +228,10 @@
(progn
(setf (slot-value cviewa 'head) (* i cst)
(slot-value cviewb 'head) (* (1+ i) cst))
- (axpy! #c(0d0 1d0) cviewb cviewa)
- (scal! #c(0d0 -2d0) cviewb)
- (axpy! #c(1d0 0d0) cviewa cviewb)
+ (copy! cviewb tmp)
+ (copy! cviewa cviewb)
+ (axpy! #c(0 1) tmp cviewa)
+ (axpy! #c(0 -1) tmp cviewb)
(incf i 2)))
(return nil)))
evec))
diff --git a/src/lapack/geqr.lisp b/src/lapack/geqr.lisp
index 3bdf501..3d129a1 100644
--- a/src/lapack/geqr.lisp
+++ b/src/lapack/geqr.lisp
@@ -32,7 +32,7 @@
;;; Initial revision for QR routines.
;;;
-(in-package "MATLISP")
+(in-package #:matlisp)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Set up the methods required to handle general matricies of Real
commit d5f56f654435a06be255e30e2f660360e6920ced
Author: Akshay Srinivasan <aks...@gm...>
Date: Sun Sep 8 01:02:54 2013 -0700
Added eig.
diff --git a/lib-src/matlisp/xerbla.f b/lib-src/matlisp/xerbla.f
new file mode 100644
index 0000000..5de4b64
--- /dev/null
+++ b/lib-src/matlisp/xerbla.f
@@ -0,0 +1,41 @@
+ SUBROUTINE XERBLA(SRNAME,INFO)
+*
+* -- LAPACK auxiliary routine (preliminary version) --
+* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+* November 2006
+*
+* .. Scalar Arguments ..
+ INTEGER INFO
+ CHARACTER*6 SRNAME
+* ..
+*
+* Purpose
+* =======
+*
+* XERBLA is an error handler for the LAPACK routines.
+* It is called by an LAPACK routine if an input parameter has an
+* invalid value. A message is printed and execution stops.
+*
+* Installers may consider modifying the STOP statement in order to
+* call system-specific exception-handling facilities.
+*
+* Arguments
+* =========
+*
+* SRNAME (input) CHARACTER*6
+* The name of the routine which called XERBLA.
+*
+* INFO (input) INTEGER
+* The position of the invalid parameter in the parameter list
+* of the calling routine.
+*
+*
+ WRITE (*,FMT=9999) SRNAME,INFO
+*
+*
+ 9999 FORMAT (' ** On entry to ',A6,' parameter number ',I2,' had ',
+ + 'an illegal value')
+*
+* End of XERBLA
+*
+ END
diff --git a/src/classes/numeric.lisp b/src/classes/numeric.lisp
index e266cfc..f360bf9 100644
--- a/src/classes/numeric.lisp
+++ b/src/classes/numeric.lisp
@@ -97,7 +97,7 @@
(imagpart (imagpart element)))
(format stream (if (zerop imagpart)
"~11,5,,,,,'Eg"
- "#C(~11,4,,,,,'Ee ~11,4,,,,,'Ee)")
+ "#C(~0,4,,,,,'Ee, ~0,4,,,,,'Ee)")
realpart imagpart)))
;;
(defleaf complex-tensor (complex-numeric-tensor) ())
diff --git a/src/lapack/eig.lisp b/src/lapack/eig.lisp
index 923aee7..94bee6a 100644
--- a/src/lapack/eig.lisp
+++ b/src/lapack/eig.lisp
@@ -198,3 +198,56 @@
:store (t/geev-output-fix ,cla wr wi)))
vl vr)))))))
(geev! A vl vr))
+;;
+
+
+(defgeneric eig (matrix &optional job)
+ (:method :before ((matrix standard-tensor) &optional (job :nn))
+ (assert (tensor-matrixp matrix) nil 'tensor-dimension-mismatch)
+ (assert (member job '(:nn :nv :vn :vv)) nil 'invalid-arguments)))
+
+(defmethod eig ((matrix complex-numeric-tensor) &optional (job :nn))
+ (mlet* ((n (nrows matrix))
+ ((levec? revec?) (values-list (mapcar #'(lambda (x) (char= x #\V)) (split-job job))))
+ (vl (when levec? (zeros (list n n) (class-of matrix))))
+ (vr (when revec? (zeros (list n n) (class-of matrix)))))
+ (geev! (copy matrix) vl vr)))
+
+(defun geev-fix-up-eigvec (n eigval eigvec)
+ (declare (type complex-tensor eigval)
+ (type real-tensor eigvec))
+ (let* ((evec (copy! eigvec (zeros (list n n) 'complex-tensor)))
+ (cviewa (col-slice~ evec 0))
+ (cviewb (col-slice~ evec 0))
+ (cst (aref (strides evec) 1)))
+ (loop
+ :with i := 0
+ :do (if (< i n)
+ (if (zerop (imagpart (ref eigval i)))
+ (incf i)
+ (progn
+ (setf (slot-value cviewa 'head) (* i cst)
+ (slot-value cviewb 'head) (* (1+ i) cst))
+ (axpy! #c(0d0 1d0) cviewb cviewa)
+ (scal! #c(0d0 -2d0) cviewb)
+ (axpy! #c(1d0 0d0) cviewa cviewb)
+ (incf i 2)))
+ (return nil)))
+ evec))
+
+(defmethod eig ((matrix real-tensor) &optional (job :nn))
+ (mlet* ((n (nrows matrix))
+ ((levec? revec?) (values-list (mapcar #'(lambda (x) (char= x #\V)) (split-job job))))
+ (ret (multiple-value-list
+ (geev! (copy matrix)
+ (when levec? (zeros (list n n) 'real-tensor))
+ (when revec? (zeros (list n n) 'real-tensor)))))
+ (eig (car ret)))
+ (if (let ((stoe (store eig)))
+ (loop :for i :from 0 :below n
+ :do (unless (zerop (t/fc (t/field-type complex-tensor) (t/store-ref complex-tensor stoe i)))
+ (return nil))
+ :finally (return t)))
+ (values-list ret)
+ (values-list (cons eig (mapcar #'(lambda (mat) (geev-fix-up-eigvec n eig mat)) (cdr ret)))))))
+
diff --git a/src/level-2/gemv.lisp b/src/level-2/gemv.lisp
index f793ccd..3dfb285 100644
--- a/src/level-2/gemv.lisp
+++ b/src/level-2/gemv.lisp
@@ -157,3 +157,8 @@
(defmethod gemv (alpha (A standard-tensor) (x standard-tensor)
beta (y standard-tensor) &optional (job :n))
(gemv! alpha A x beta (copy y) job))
+
+(defmethod gemv (alpha (A standard-tensor) (x standard-tensor)
+ (beta (eql nil)) (y (eql nil)) &optional (job :n))
+ (let ((ret (zeros (nrows A) (class-of A))))
+ (gemv! alpha A x 1 ret job)))
diff --git a/src/old/foreign-real-matrix.lisp b/src/old/foreign-real-matrix.lisp
index a0c0248..6d933b1 100644
--- a/src/old/foreign-real-matrix.lisp
+++ b/src/old/foreign-real-matrix.lisp
@@ -34,4 +34,4 @@
(declare (type fixnum n m))
(make-instance 'foreign-real-matrix :nrows n :ncols m
- :store store))
\ No newline at end of file
+ :store store))
diff --git a/src/old/help.lisp b/src/old/help.lisp
index 38f0e68..dbebe14 100644
--- a/src/old/help.lisp
+++ b/src/old/help.lisp
@@ -236,4 +236,3 @@ For example, (HELP matlisp) (HELP mapcar)")
(format stream "~&~%No documentation available for symbol ~a"
item)
(values)))))
-
commit 739477d2ea4ae8e582b2355220d502443cd722a3
Author: Akshay Srinivasan <aks...@gm...>
Date: Tue Aug 27 15:41:13 2013 -0700
Added a nicer dlsode interface. Added stuff to the gnuplot interface.
diff --git a/lib-src/gnuplot/gnuplot.lisp b/lib-src/gnuplot/gnuplot.lisp
index d34a360..f9b48dc 100644
--- a/lib-src/gnuplot/gnuplot.lisp
+++ b/lib-src/gnuplot/gnuplot.lisp
@@ -2,11 +2,20 @@
(defvar *current-gnuplot-process* nil)
(defun open-gnuplot-stream (&optional (gnuplot-binary (pathname "/usr/bin/gnuplot")))
- (#+:sbcl
- sb-ext:run-program
- #+:ccl
- ccl:run-program
- gnuplot-binary nil :input :stream :wait nil :output t))
+ (setf *current-gnuplot-process* (#+:sbcl
+ sb-ext:run-program
+ #+:ccl
+ ccl:run-program
+ gnuplot-binary nil :input :stream :wait nil :output t))
+ (gnuplot-send "
+set datafile fortran
+")
+ *current-gnuplot-process*)
+
+(defun close-gnuplot-stream ()
+ (when *current-gnuplot-process*
+ (gnuplot-send "quit~%")
+ (setf *current-gnuplot-process* nil)))
(defun gnuplot-send (str &rest args)
(unless *current-gnuplot-process*
@@ -19,14 +28,35 @@
(apply #'format (append (list stream str) args))
(finish-output stream)))
-(defun plot2d (data &key (lines t) (color (list "#FF0000")))
- (with-open-file (s "/tmp/matlisp-gnuplot.out" :direction :output :if-exists :supersede :if-does-not-exist :create)
- (loop :for i :from 0 :below (loop :for x :in data :minimizing (size x))
- :do (loop :for x :in data :do (format s "~a " (coerce (ref x i) 'single-float)) :finally (format s "~%"))))
- (if lines
- (gnuplot-send "plot '/tmp/matlisp-gnuplot.out' with lines linecolor rgb ~s~%" color)
- (gnuplot-send "plot '/tmp/matlisp-gnuplot.out'~%")))
+(defun splitcol (num)
+ (multiple-value-bind (a b0) (floor num 256)
+ (multiple-value-bind (b2 b1) (floor a 256)
+ (list b2 b1 b0))))
+(defun plot2d (data &key (lines t) (color nil))
+ (let ((fname "/tmp/matlisp-gnuplot.out"))
+ (with-open-file (s fname :direction :output :if-exists :supersede :if-does-not-exist :create)
+ (loop :for i :from 0 :below (loop :for x :in data :minimizing (size x))
+ :do (loop :for x :in data :do (format s "~a " (coerce (ref x i) 'single-float)) :finally (format s "~%"))))
+ (let ((col (if (listp color) color
+ (let ((lst (list color)))
+ (setf (cdr lst) lst)
+ lst))))
+ (let ((cmd (apply #'string+ (cons "plot " (loop :for x :in (cdr data)
+ :for i := 2 :then (1+ i)
+ :for clist := col :then (cdr clist)
+ :collect (string+ "'" fname "' using 1:" (format nil "~a " i)
+ "with " (if lines "lines" "points") " "
+ (if (car clist)
+ (apply #'(lambda (r g b) (format nil "linecolor rgb(~a, ~a, ~a)" r g b))
+ (splitcol (car clist)))
+ "")
+ (format nil "title \"~a\"" (1- i))
+ ", "))))))
+ (setf (aref cmd (- (length cmd) 2)) #\;
+ (aref cmd (- (length cmd) 1)) #\Newline)
+ (gnuplot-send cmd)))))
+
;; (defclass gnuplot-plot-info ()
;; ((title
;; :initform "GNU PLOT"
diff --git a/packages.lisp b/packages.lisp
index 04096db..f3d665d 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -61,6 +61,7 @@
#:tensor-cannot-find-counter-class
#:tensor-cannot-find-optimization
#:tensor-dimension-mismatch
+ #:tensor-type-mismatch
#:tensor-store-not-consecutive
#:tensor-method-does-not-exist
#:tensor-abstract-class
diff --git a/src/base/print.lisp b/src/base/print.lisp
index d058298..3ad6c7f 100644
--- a/src/base/print.lisp
+++ b/src/base/print.lisp
@@ -29,13 +29,13 @@
(in-package #:matlisp)
;; Routines for printing a tensors/matrices nicely.
-(defparameter *print-max-len* 5
+(defparameter *print-max-len* 10
"
Maximum number of elements in any particular argument to print.
Set this to T to print all the elements.
")
-(defparameter *print-max-args* 2
+(defparameter *print-max-args* 5
"
Maximum number of arguments of the tensor to print.
Set this to T to print all the arguments.
diff --git a/src/base/standard-tensor.lisp b/src/base/standard-tensor.lisp
index 426536f..cccd99b 100644
--- a/src/base/standard-tensor.lisp
+++ b/src/base/standard-tensor.lisp
@@ -397,46 +397,54 @@
(declare (type standard-tensor tensor)
(type list subscripts)
(type boolean preserve-rank))
- (let-typed ((dims (dimensions tensor) :type index-store-vector)
- (stds (strides tensor) :type index-store-vector)
- (rank (rank tensor) :type index-type))
- (loop :for (start step end) :in subscripts
- :for i :of-type index-type := 0 :then (1+ i)
- :with ndims :of-type index-store-vector := (allocate-index-store rank)
- :with nstds :of-type index-store-vector := (allocate-index-store rank)
- :with nrank :of-type index-type := 0
- :with nhd :of-type index-type := (head tensor)
- :do (assert (< i rank) nil 'tensor-index-rank-mismatch :index-rank (1+ i) :rank rank)
- :do (let* ((start (if (eq start '*) 0
- (progn
- (assert (and (typep start 'index-type) (< -1 start (aref dims i))) nil 'tensor-index-out-of-bounds :argument i :index start :dimension (aref dims i))
- start)))
- (step (if (eq step '*) 1
- (progn
- (assert (and (typep step 'index-type) (< 0 step)) nil 'invalid-value :given step :expected '(< 0 step) :message "STEP cannot be <= 0.")
- step)))
- (end (if (eq end '*) (aref dims i)
- (progn
- (assert (and (typep end 'index-type) (<= 0 end (aref dims i))) nil 'tensor-index-out-of-bounds :argument i :index start :dimension (aref dims i))
- end))))
- (declare (type index-type start step end))
- ;;
- (let-typed ((dim (ceiling (the index-type (- end start)) step) :type index-type))
- (unless (and (= dim 1) (not preserve-rank))
- (setf (aref ndims nrank) dim
- (aref nstds nrank) (* step (aref stds i)))
- (incf nrank))
- (when (/= start 0)
- (incf nhd (the index-type (* start (aref stds i)))))))
- :finally (return
- (if (= nrank 0) (store-ref tensor nhd)
- (let ((*check-after-initializing?* nil))
- (make-instance (class-of tensor)
- :head nhd
- :dimensions (very-quickly (vectorify (the index-store-vector ndims) nrank 'index-type))
- :strides (very-quickly (vectorify (the index-store-vector nstds) nrank 'index-type))
- :store (store tensor)
- :parent-tensor tensor)))))))
+ (if (null subscripts)
+ (let ((*check-after-initializing?* nil))
+ (make-instance (class-of tensor)
+ :head (head tensor)
+ :dimensions (copy-seq (dimensions tensor))
+ :strides (copy-seq (strides tensor))
+ :store (store tensor)
+ :parent-tensor tensor))
+ (let-typed ((dims (dimensions tensor) :type index-store-vector)
+ (stds (strides tensor) :type index-store-vector)
+ (rank (rank tensor) :type index-type))
+ (loop :for (start step end) :in subscripts
+ :for i :of-type index-type := 0 :then (1+ i)
+ :with ndims :of-type index-store-vector := (allocate-index-store rank)
+ :with nstds :of-type index-store-vector := (allocate-index-store rank)
+ :with nrank :of-type index-type := 0
+ :with nhd :of-type index-type := (head tensor)
+ :do (assert (< i rank) nil 'tensor-index-rank-mismatch :index-rank (1+ i) :rank rank)
+ :do (let* ((start (if (eq start '*) 0
+ (progn
+ (assert (and (typep start 'index-type) (< -1 start (aref dims i))) nil 'tensor-index-out-of-bounds :argument i :index start :dimension (aref dims i))
+ start)))
+ (step (if (eq step '*) 1
+ (progn
+ (assert (and (typep step 'index-type) (< 0 step)) nil 'invalid-value :given step :expected '(< 0 step) :message "STEP cannot be <= 0.")
+ step)))
+ (end (if (eq end '*) (aref dims i)
+ (progn
+ (assert (and (typep end 'index-type) (<= 0 end (aref dims i))) nil 'tensor-index-out-of-bounds :argument i :index start :dimension (aref dims i))
+ end))))
+ (declare (type index-type start step end))
+ ;;
+ (let-typed ((dim (ceiling (the index-type (- end start)) step) :type index-type))
+ (unless (and (= dim 1) (not preserve-rank))
+ (setf (aref ndims nrank) dim
+ (aref nstds nrank) (* step (aref stds i)))
+ (incf nrank))
+ (when (/= start 0)
+ (incf nhd (the index-type (* start (aref stds i)))))))
+ :finally (return
+ (if (= nrank 0) (store-ref tensor nhd)
+ (let ((*check-after-initializing?* nil))
+ (make-instance (class-of tensor)
+ :head nhd
+ :dimensions (very-quickly (vectorify (the index-store-vector ndims) nrank 'index-type))
+ :strides (very-quickly (vectorify (the index-store-vector nstds) nrank 'index-type))
+ :store (store tensor)
+ :parent-tensor tensor))))))))
(definline slice~ (x axis &optional (idx 0))
(let ((slst (make-list (rank x) :initial-element '(* * *))))
diff --git a/src/conditions.lisp b/src/conditions.lisp
index d396359..dc7efc0 100644
--- a/src/conditions.lisp
+++ b/src/conditions.lisp
@@ -218,6 +218,13 @@
(declare (ignore c))
(format stream "The dimensions of the given tensors are not suitable for continuing with the operation."))))
+(define-condition tensor-type-mismatch (tensor-error)
+ ()
+ (:documentation "The types of the given tensors are not suitable for continuing with the operation.")
+ (:report (lambda (c stream)
+ (declare (ignore c))
+ (format stream "The types of the given tensors are not suitable for continuing with the operation."))))
+
(define-condition tensor-store-not-consecutive (tensor-error)
()
(:documentation "The strides of the store, of the given tensor are not conscutive.")
diff --git a/src/packages/odepack/dlsode.lisp b/src/packages/odepack/dlsode.lisp
index 466603c..15206e6 100644
--- a/src/packages/odepack/dlsode.lisp
+++ b/src/packages/odepack/dlsode.lisp
@@ -15,7 +15,7 @@
(c-y (* :double-float :size c-neq) :input)
(c-ydot (* :double-float :size c-neq) :output)))
(neq :integer :input)
- (y (* :double-float) :input-output)
+ (y (* :double-float :inc head-y) :input-output)
(ts :double-float :input-output)
(tout :double-float :input)
(itol :integer :input)
@@ -40,6 +40,55 @@
(mf :integer :input))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+(defun ode-evolve (field y0 t-array)
+ (declare (type real-tensor y0 t-array))
+ (assert (and (tensor-vectorp t-array) (tensor-vectorp y0)) nil 'tensor-dimension-mismatch)
+ (let* ((neq (size y0))
+ (nt (size t-array))
+ (t0 (ref t-array 0))
+ (ts t0)
+ (ret (zeros (list neq nt) 'real-tensor))
+ ;;
+ (lrw (+ 22 (* 9 neq) (* neq neq) 5))
+ (liw (+ 20 neq 5))
+ (itol 1)
+ (atol (make-array 1 :element-type 'double-float :initial-element 1d-12))
+ (rtol (make-array 1 :element-type 'double-float :initial-element 1d-12))
+ (itask 1)
+ (istate 1)
+ (iopt 0)
+ (mf 22)
+ (rwork (make-array lrw :element-type 'double-float :initial-element 0d0))
+ (iwork (make-array liw :element-type '(signed-byte 32) :initial-element 0))
+ ;;
+ (view (slice~ ret 1))
+ (stv (aref (strides ret) 1))
+ ;;
+ (y-tmp (zeros neq 'real-tensor))
+ (stoy (store y-tmp)))
+ (copy! y0 view)
+ (incf (slot-value view 'head) stv)
+ (labels ((field-sugar (neq time yf ydotf)
+ (loop :for i :from 0 :below neq
+ :do (t/store-set real-tensor (fv-ref yf i) stoy i))
+ ;;Because of some black magic, this does not seem to
+ ;;affect the amount of memory allocated!
+ (let ((ydot (funcall field time y-tmp)))
+ (loop :for i :from 0 :below neq
+ :do (setf (fv-ref ydotf i) (ref ydot i))))
+ nil))
+ (loop :for i :from 1 :below nt
+ :do (let ((tout (ref t-array i)))
+ (multiple-value-bind (y-out ts-out istate-out rwork-out iwork-out)
+ (dlsode #'field-sugar neq (store y0) ts tout itol rtol atol itask istate iopt rwork lrw iwork liw #'(lambda (&rest th) (declare (ignore th))) mf (head y0))
+ (declare (ignore y-out rwork-out iwork-out))
+ (setq ts ts-out)
+ (setq istate istate-out))
+ (copy! y0 view)
+ (incf (slot-value view 'head) stv))))
+ ret))
+;;
+
(defun lsode-evolve (field y t-array report)
;;
(let* ((neq (length y))
@@ -61,7 +110,7 @@
do (progn
(setq tout (aref t-array i))
(multiple-value-bind (y-out ts-out istate-out rwork-out iwork-out)
- (dlsode field neq y ts tout itol rtol atol itask istate iopt rwork lrw iwork liw #'(lambda (&rest th) (declare (ignore th))) mf)
+ (dlsode field neq (store y) ts tout itol rtol atol itask istate iopt rwork lrw iwork liw #'(lambda (&rest th) (declare (ignore th))) mf)
(setq ts ts-out)
(setq istate istate-out))
(funcall report ts y)))))
@@ -74,10 +123,32 @@
(defun pend-report (ts y)
(format t "~A ~A ~A ~%" ts (aref y 0) (aref y 1)))
-(defun pcart-field (neq time y ydot)
- (declare (ignore neq time))
- (very-quickly
- (destructuring-bind (x theta xdot thetadot) (mapcar #'(lambda (n) (fv-ref y n)) '(0 1 2 3))
+(defun pcart-field (time y ydot)
+ (declare (ignore time))
+ (destructuring-bind (x theta xdot thetadot) (mapslice #'id y)
+ (setf (ref ydot 0) xdot
+ (ref ydot 1) thetadot
+ (ref ydot 2) (/ (+ (* (cos theta) (sin theta)) (* (sin theta) (expt thetadot 2))) (- 2 (expt (cos theta) 2)))
+ (ref ydot 3) (/ (+ (* 2 (sin theta)) (* (cos theta) (sin theta) (expt thetadot 2))) (- (expt (cos theta) 2) 2)))))
+
+(defun pcart-field (time y)
+ (declare (ignore time))
+ (let ((ydot (zeros 4 'real-tensor)))
+ (destructuring-bind (x theta xdot thetadot) (mapslice #'id y)
+ (setf (ref ydot 0) xdot
+ (ref ydot 1) thetadot
+ (ref ydot 2) (/ (+ (* (cos theta) (sin theta)) (* (sin theta) (expt thetadot 2))) (- 2 (expt (cos theta) 2)))
+ (ref ydot 3) (/ (+ (* 2 (sin theta)) (* (cos theta) (sin theta) (expt thetadot 2))) (- (expt (cos theta) 2) 2))))
+ ydot))
+
+
+(defun pcart-report (ts y)
+ (format t "~A ~A ~A ~A ~A ~%" ts (aref y 0) (aref y 1) (aref y 2) (aref y 3)))
+
+
+(defun pcart-field (time y ydot)
+ (declare (ignore time))
+ (destructuring-bind (x theta xdot thetadot) (mapcar #'(lambda (n) (fv-ref y n)) '(0 1 2 3))
(declare (type double-float x theta xdot thetadot))
(setf (fv-ref ydot 0) xdot
(fv-ref ydot 1) thetadot
@@ -87,6 +158,7 @@
(defun pcart-report (ts y)
(format t "~A ~A ~A ~A ~A ~%" ts (aref y 0) (aref y 1) (aref y 2) (aref y 3)))
+
#+nil
(let ((y (make-array 2 :element-type 'double-float :initial-contents `(,(/ pi 2) 0d0))))
(lsode-evolve #'pend-field y #(0d0 1d0 2d0) #'pend-report))
@@ -94,6 +166,8 @@
;; 1.0d0 1.074911802207049d0 -0.975509986605856d0
;; 2.0d0 -0.20563950412081608d0 -1.3992359518735706d0
+(ode-evolve #'pcart-field (copy! (list 0 (/ pi 3) 0 0) (zeros 4)) (range 0 10))
+
#+nil
(let ((y (make-array 4 :element-type 'double-float :initial-contents `(0d0 ,(/ pi 3) 0d0 0d0)))
diff --git a/src/special/map.lisp b/src/special/map.lisp
index 9cb083a..98b8443 100644
--- a/src/special/map.lisp
+++ b/src/special/map.lisp
@@ -40,13 +40,28 @@
y)))
(mapsor! func x y))
;;
-(defun mapslice (func x &optional (axis 1))
+
+(defun mapslice (func x &optional (axis 0))
+ (declare (type standard-tensor x))
+ (if (tensor-vectorp x)
+ (loop :for i :from 0 :below (aref (dimensions x) axis)
+ :collect (funcall func (ref x i)))
+ (let* ((v-x (slice~ x axis))
+ (st-x (aref (strides x) axis)))
+ (loop :for i :from 0 :below (aref (the index-store-vector (dimensions x)) axis)
+ :collect (prog1 (funcall func (copy v-x))
+ (incf (slot-value v-x 'head) st-x))))))
+
+(defun mapslice~ (func x &optional (axis 0))
(declare (type standard-tensor x))
- (let* ((v-x (slice~ x axis))
- (st-x (aref (strides x) axis)))
- (loop :for i :from 0 :below (aref (the index-store-vector (dimensions x)) axis)
- :collect (prog1 (funcall func (copy v-x))
- (incf (slot-value v-x 'head) st-x)))))
+ (if (tensor-vectorp x)
+ (loop :for i :from 0 :below (aref (dimensions x) axis)
+ :collect (funcall func (ref x i)))
+ (let* ((v-x (slice~ x axis))
+ (st-x (aref (strides x) axis)))
+ (loop :for i :from 0 :below (aref (the index-store-vector (dimensions x)) axis)
+ :collect (prog1 (funcall func (sub-tensor~ v-x nil))
+ (incf (slot-value v-x 'head) st-x))))))
(defmacro tensor-foldl (type func ten init &optional (init-type (field-type type)))
(using-gensyms (decl (ten init))
-----------------------------------------------------------------------
Summary of changes:
configure | 8 +++
configure.ac | 8 +++
lib-src/gnuplot/gnuplot.lisp | 54 ++++++++++++++++-----
lib-src/matlisp/Makefile.am | 3 +-
packages.lisp | 1 +
src/base/print.lisp | 4 +-
src/base/standard-tensor.lisp | 88 ++++++++++++++++++---------------
src/classes/numeric.lisp | 2 +-
src/conditions.lisp | 7 +++
src/foreign-core/lazy-loader.lisp.in | 6 +-
src/lapack/eig.lisp | 56 +++++++++++++++++++++-
src/lapack/geqr.lisp | 2 +-
src/level-2/gemv.lisp | 5 ++
src/old/foreign-real-matrix.lisp | 2 +-
src/old/help.lisp | 1 -
src/packages/odepack/dlsode.lisp | 86 ++++++++++++++++++++++++++++++--
src/special/map.lisp | 27 ++++++++--
17 files changed, 284 insertions(+), 76 deletions(-)
hooks/post-receive
--
matlisp
|
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X