Week 1: phi functions

 

During the first week of GSoC, I wrote four m-files for matrix functions, phi1m,..., phi4m.
I implemented the four functions, based on [BSW 07].
Matrix functions are defined by the recurrence relation

In my files, I use an algorithm based on first a variant of the scaling and squaring approach and after scaling, I use a Padé approximation.

Below the code phi1m.m

function [N, PHI0] = phi1m (A, d)

N is φ₁(A) and PHIO is φ₀(A). d is the order of the Padé approximation.

First, I scale A by power of 2 so that its norm is < ½

l=1;

s= min (ceil (max (0,1+log2 (norm (A, inf)))), 1023);

A = A/2^s;

then I use a (d,d)-Padé approximation

where the polynomials are

I write polynomials in Horner form.

ID = eye (size (A));

i=d;

Ncoeff = sum (cumprod ([1, d-(0:i-1)])./cumprod ([1, 2*d+l-(0:i-1)]).*(-1).^(0:i)./(factorial ((0:i)).*factorial (l+i-(0:i))));

Dcoeff = ((-1)^(i))*prod ([(d-i+1):1:d])/(prod ([(2*d+l-i+1):1:(2*d+l)])*factorial (i));

N= Ncoeff;

D= Dcoeff;

for i = (d-1):-1:0

   Ncoeff = sum (cumprod ([1, d-(0:i-1)])./cumprod ([1, 2*d+l-(0:i-1)]).*(-1).^(0:i)./(factorial ((0:i)).*factorial (l+i-(0:i))));

   N = A*N + Ncoeff * ID ;

   Dcoeff = ((-1)^(i))* prod([(d-i+1):1:d])/(prod ([(2*d+l-i+1):1:(2*d+l)])*factorial (i));

   D = A*D + Dcoeff * ID ;

endfor

N = full (D\N);

and finally I use the scaling relations

PHI0 = A*N+ID;

for i = 1:s

   N = (PHI0+ID)*N/2;

   PHI0 = PHI0*PHI0;

endfor

in the same way I wrote the other three files.  As soon as possible, I will create a repository on Bitbucket and I will put the codes there.