[Gsll-devel] gsl_linalg_exponential_ss

[Liam Healy]

Tamas,

I added this function to the future (ffa) version of GSLL, and
it works.  This is my attempt at a backport to the current version:

(defmfun matrix-exponential (matrix exponential)
  "gsl_linalg_exponential_ss"
  (((mpointer matrix) :pointer) ((mpointer exponential) :pointer)
  :mode)
  :invalidate (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.")

I have not tried it, so try it and see how it works.
If you give it #2A((0.0d0 1.0d0) (-1.0d0 0.0d0)) it
should return
#2A((0.5403023058681384d0 0.841470984807895d0)
    (-0.841470984807895d0 0.5403023058681384d0))
i.e., [[cos(1), sin(1)], [-sin(1), cos(1)]]

Liam



On Wed, Jul 23, 2008 at 11:45 AM, Tamas K Papp <tpapp at princeton.edu> wrote:
> On Wed, Jul 23, 2008 at 11:16:27AM -0400, Liam Healy wrote:
>
> Hi Liam,
>
>> Hi Tamas,
>>
>> No.  I'm trying to figure out the reason for this, and the big one I
>> can see is that I don't see it in the GSL documentation.  It shows up
>> in my /usr/include/gsl/gsl_linalg.h and in /usr/lib/libgsl.so.   It's
>> possible that it was introduced since GSL 1.8, which is where my
>> porting started.  One of my long-term goals is to have an accounting
>> of the port status of every symbol in the library so that I'll catch
>> gaps and new functions like this.  In the meantime, it shouldn't be
>> hard to make a binding, straightforward use of defmfun I think.
>
> Would you please add this binding?  I am reading the source, but still
> haven't figured out defmfun to know how to do this.
>
>> BTW that paper of Moler and Van Loan is a classic and was recently updated
>> Cleve B. Moler and Charles F. Van Loan. Nineteen dubious ways to
>> compute the exponential of a
>>    matrix, twenty-five years later. SIAM Rev., 45(1):3–49, 2003.
>
> Thanks, I just read this, but didn't feel up to programming the whole
> thing :-) There is a Fortran library called Expokit which is said to
> be nice, but I found it impossible to work with CL.
>
> Thanks,
>
> Tamas
>