Re: [Gsl-shell-info] Is there a function to compute matrix powers or the matrix exponential?

[Francesco Abbate]

2013/8/15 Omar Antolín Camarena <address@hidden>:
> Thanks for the reply, Francesco. I was just surprised that GSL Shell
> did have A^n to compute integral powers of a matrix and started
> browsing the documentation to see if powers weren't available under a
> differentsyntax. That's when I realized that it didn't have the
> exponential function either and I wanted to confirm I hadn't somehow
> missed those things. (So I never actually needed to compute exp(M),
> just A^n.)

Hmmm, it seems that I actually misunderstood your question, sorry :-)

Of course you are right, the matrix powers, at least for integer
exponents, are basic functions that I ... cough cough... "forgot" to
implement a long time ago :-)

Now it is done. In the github master branch I've implemented the
matrix powers for integer exponents for both real or complex matrices.
The nice things is that you can also give negative integers at the
exponent to obtain the inverse matrix and its powers.

I've used the "exponentiation by squaring" algorithm as this seems to
be the "way to go" when the exponent is an integer. The algorithm is
actually very simple, it is described for example here:

http://stackoverflow.com/questions/101439/the-most-efficient-way-to-implement-an-integer-based-power-function-powint-int

the real challenge was to implement it in an optimal way using a
limited amount of memory. I believe I've been quite successful since
only 2*n^2 addition space is required where n is the size of the
matrix. Otherwise the algorithm is executed in O(log(P)) where P is
the exponent. Of course I've used the BLAS dgemm/zgemm functions under
the hood to perform efficiently the multiplications.

Thank you for pointing out this missing feature in GSL Shell.

Best regards
Francesco