Commit 4a3fb0c5 authored by Liam Healy's avatar Liam Healy

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.
parent 6f90daf7
;; 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")
;; 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))
(((mpointer matrix) :pointer) ((mpointer exponential) :pointer)
(mode sf-mode))
:inputs (matrix)
:outputs (exponential)
:return (exponential)
"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))
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment