From 4a3fb0c5b9ea32ffdfdba252593c165c03986399 Mon Sep 17 00:00:00 2001 From: Liam Healy <liam@thinkpad.local> Date: Sun, 10 Aug 2008 18:44:53 -0400 Subject: [PATCH] Exponential of a matrix Function #'matrix-exponential to compute the exponential of a matrix using the methods of Moler & Van Loan. This function is undocumented in GSL but seems to work. --- gsll.asd | 3 ++- linear-algebra/exponential.lisp | 33 +++++++++++++++++++++++++++++++++ 2 files changed, 35 insertions(+), 1 deletion(-) create mode 100644 linear-algebra/exponential.lisp diff --git a/gsll.asd b/gsll.asd index 1fc4d032..6fe0ffac 100644 --- a/gsll.asd +++ b/gsll.asd @@ -1,6 +1,6 @@ ;; Definition of GSLL system ;; Liam Healy -;; Time-stamp: <2008-08-10 17:01:08EDT gsll.asd> +;; Time-stamp: <2008-08-10 18:27:45EDT gsll.asd> ;; $Id$ (asdf:defsystem "gsll" @@ -82,6 +82,7 @@ ((:file "blas1") (:file "blas2") (:file "blas3" :depends-on (blas2)) + (:file "exponential") ;(:file "lu") ;(:file "qr") ;(:file "qrpt") diff --git a/linear-algebra/exponential.lisp b/linear-algebra/exponential.lisp new file mode 100644 index 00000000..ecd54695 --- /dev/null +++ b/linear-algebra/exponential.lisp @@ -0,0 +1,33 @@ +;; Exponential of a matrix +;; Liam Healy 2008-08-10 17:25:35EDT exponential.lisp +;; Time-stamp: <2008-08-10 18:27:04EDT exponential.lisp> +;; $Id: $ + +(in-package :gsl) + +;;; No documentation +(defmfun matrix-exponential (matrix exponential &optional (mode :double)) + "gsl_linalg_exponential_ss" + (((mpointer matrix) :pointer) ((mpointer exponential) :pointer) + (mode sf-mode)) + :inputs (matrix) + :outputs (exponential) + :return (exponential) + :documentation + "Calculate the matrix exponential by the scaling and + squaring method described in Moler + Van Loan, + SIAM Rev 20, 801 (1978). The matrix to be exponentiated + is matrix, the returned exponential is exponential. + The mode argument allows + choosing an optimal strategy, from the table + given in the paper, for a given precision.") + +#| +;;; A matrix of the form [[0, 1], [-1, 0]] +;;; when exponentiated gives [[cos x, sin x], [-sin x, cos x]] +(letm ((mat (matrix-double-float (a (0.0d0 1.0d0) (-1.0d0 0.0d0)))) + (exp (matrix-double-float '(2 2)))) + (cl-array (matrix-exponential mat exp))) +#2A((0.5403023058681384d0 0.841470984807895d0) + (-0.841470984807895d0 0.5403023058681384d0)) +|# -- GitLab