function F = dgeexfhrp(A,k)
%   DGEEXFHRP  Matrix exponential function based on Hermites series.
%
%   dgeexfhrp(A, k) computes the matrix exponential of the 
%   matrix A, where A is a double precision class matrix, using
%   a scaling and squaring algorithm based on a mixed Taylor and 
%   Hermite series approximation, with matrix product reduction 
%   based on relative accuracy bounds on floating point arithmetic
%   
%   Input parameters:
%       A matrix
%       k is the maximum order approximation 'm' to be used:    
%                            k: 1  2  3  4  5  6  7  8
%      Approximation order (m): 4  6  9 12 16 20 25 30
%
%   Output parameter:
%   	F=exponential of matrix A
%
%   Reference: J. Sastre, J. Ibañez, E. Defez and P. Ruiz, 
%   "Eficient orthogonal matrix polynomial based method for computing
%   matrix exponential", Applied Mathematics and Computation, 217 (2011),
%   pp. 6451–6463.

%   Authors: Jorge Sastre, Javier Ibáñez, Pedro A. Ruiz, Emilio Defez.
%   Revised version: 2015/11/25.
%
%   Group of High Performance Scientific Computation (HiPerSC)
%   Universitat Politecnica de Valencia (Spain)
%   http://hipersc.blogs.upv.es

% Check arguments
 if nargin<1 || nargin>2
    error('Incorrect number of input arguments.');
 end
 if nargin==1
    k = -1;
 end
 if k<1||k>8
     k=8;
 end
theta=[3.397168839976962e-4   %m=4  Taylor approximation
       9.065656407595101e-3   %m=6  Taylor approximation
       8.957760203223343e-2   %m=9  Taylor approximation
       2.996158913811581e-1   %m=12 Taylor approximation
       7.802874256626574e-1   %m=16 Taylor approximation
       1.438252596804337      %m=20 Taylor approximation 
       2.441356829252848      %N=25 Hermite aproximation, optimal parameter lambda = 16.66121324200387
       3.578700513755017];    %N=30 Hermite aproximation, optimal parameter lambda = 7.596210771817034
   
u=eps/2;%unit roundoff   
% Compute the Hermite approximation order (m) 
% and the scaling factor (s) that will be used
normA = norm(A,1);
for j=1:k
	if normA<=theta(j)
        break
    end
end

if (j==k)&&(normA>theta(k))
    [t0 s] = log2(normA/theta(j));
    s = s - (t0 == 0.5); % Adjust s if normA/theta(end) is a power of 2.
    if k>1
        if normA/2^s<=theta(k-1) % This condition might occur because theta{k}/2<theta{k-1}
            j=k-1;
        end
    end
    A = A/2^s;    % Scaling the matrix A
 else
    s=0;
 end

%Compute the low bound used in the product reduction;
Low_bound_expAs_u=exp(-(normA/2^s))*2^(-53);

n = size(A,1); 
A2 = A*A;
normA2s = norm(A2,1);
Cp=1;  %Number of matrix products (used only for performance evaluation)
I = eye(n);

% Compute F=Hm(A/2^s,lambda) using the Paterson-Stockmeyer method
% If it is possible, product reductions are applied
switch j 
	case 1 %m=4    %Taylor method for order m<=20
        p=[1.0, 0.5, 0.166666666666666666666666666667, 0.0416666666666666666666666666667];
        F= p(4)*A2+p(3)*A;
        auxnormF=norm(F,1);
        F=F+p(2)*I;
        if (auxnormF<=p(2)*u)||(min(auxnormF,norm(F,1))*normA2s <= Low_bound_expAs_u)
            F= p(2)*A2+p(1)*A+I;   %One matrix product is saved
        else
            F=    F*A2+p(1)*A+I;
            Cp=Cp+1;
        end
	case 2 %m=6 
        p=[1.0, 0.5, 0.166666666666666666666666666667, 0.0416666666666666666666666666667, 0.00833333333333333333333333333333, 0.00138888888888888888888888888889];        
        F= p(6)*A2+p(5)*A;
        auxnormF=norm(F,1);
        F=F+p(4)*I;
        if (auxnormF<=p(4)*u)||(min(auxnormF,norm(F,1))*normA2s^2 <= Low_bound_expAs_u)
            F= p(4)*A2+p(3)*A;   %One matrix product is saved
        else
            F=    F*A2+p(3)*A;
           Cp=Cp+1;
        end
        auxnormF=norm(F,1);
        F=F+p(2)*I;        
        if (auxnormF<=p(2)*u)||(min(auxnormF,norm(F,1))*normA2s <= Low_bound_expAs_u)
            F= p(2)*A2+ p(1)*A+I;   %One matrix product is saved
        else
            
            F=    F*A2+ p(1)*A+I;
           Cp=Cp+1;
        end
	case 3 %m=9
        p=[1.0, 0.5, 0.166666666666666666666666666667, 0.0416666666666666666666666666667, 0.00833333333333333333333333333333, 0.00138888888888888888888888888889, 0.000198412698412698412698412698413, 0.0000248015873015873015873015873016, 0.00000275573192239858925109505932705];
        A3=A*A2;
        Cp=Cp+1;
        normA3s = norm(A3,1);
        F= p(9)*A3+p(8)*A2+p(7)*A;
        auxnormF=norm(F,1);
        F=F+p(6)*I;
        if (auxnormF<=p(6)*u)||(min(auxnormF,norm(F,1))*normA3s^2<= Low_bound_expAs_u)
            F= p(6)*A3+p(5)*A2+p(4)*A;   %One matrix product is saved
        else
            F=    F*A3+p(5)*A2+p(4)*A;
            Cp=Cp+1;
        end
        auxnormF=norm(F,1);
        F=F+p(3)*I;
        if (auxnormF<=p(3)*u)||(min(auxnormF,norm(F,1))*normA3s <= Low_bound_expAs_u)
            F=     p(3)*A3+p(2)*A2+p(1)*A+I;   %One matrix product is saved
        else
            F=        F*A3+p(2)*A2+p(1)*A+I;
            Cp=Cp+1;
        end

	case 4 %m=12
        p=[1.0, 0.5, 0.166666666666666666666666666667, 0.0416666666666666666666666666667, 0.00833333333333333333333333333333, 0.00138888888888888888888888888889, 0.000198412698412698412698412698413, 0.0000248015873015873015873015873016, 0.00000275573192239858925109505932705, 0.00000027557319223985888275785856999, 0.0000000250521083854417202238661793214, 0.00000000208767569878681001865551494345];
        A3=A*A2;
        Cp=Cp+1;
        normA3s = norm(A3,1);
        F= p(12)*A3+p(11)*A2+p(10)*A;
        auxnormF=norm(F,1);
        F=F+p(9)*I;
        if (auxnormF<=p(9)*u)||(min(auxnormF,norm(F,1))*normA3s^3 <= Low_bound_expAs_u)  
            F= p(9)*A3+ p(8)*A2+ p(7)*A;   %One matrix product is saved
        else
            F=    F*A3+ p(8)*A2+ p(7)*A;
            Cp=Cp+1;
        end
        auxnormF=norm(F,1);
        F=F+p(6)*I;       
        if (auxnormF<=p(6)*u)||(min(auxnormF,norm(F,1))*normA3s^2 <= Low_bound_expAs_u) 
            F=  p(6)*A3+ p(5)*A2+ p(4)*A;   %One matrix product is saved
        else
            F=     F*A3+ p(5)*A2+ p(4)*A;
            Cp=Cp+1;
        end  
        auxnormF=norm(F,1);
        F=F+p(3)*I;       
        if (auxnormF<=p(3)*u)||(min(auxnormF,norm(F,1))*normA3s <= Low_bound_expAs_u)
            F=      p(3)*A3+ p(2)*A2+ p(1)*A+I;   %One matrix product is saved
        else
            F=         F*A3+ p(2)*A2+ p(1)*A+I;
            Cp=Cp+1;
        end

	case 5 %m=16
        p=[1.0, 0.5, 0.166666666666666666666666666667, 0.0416666666666666666666666666667, 0.00833333333333333333333333333333, 0.00138888888888888888888888888889, 0.000198412698412698412698412698413, 0.0000248015873015873015873015873016, 0.00000275573192239858925109505932705, 0.00000027557319223985888275785856999, 0.0000000250521083854417202238661793214, 0.00000000208767569878681001865551494345, 0.000000000160590438368216133408629182949, 0.0000000000114707455977297245072965704504, 0.000000000000764716373181981640551384421244, 0.0000000000000477947733238738525344615263278];
        A3=A*A2;A4=A2*A2;
        Cp=Cp+2;
        normA4s = norm(A4,1);
        F=p(16)*A4+p(15)*A3+p(14)*A2+p(13)*A;
        auxnormF=norm(F,1);
        F=F+p(12)*I;
        if (auxnormF<=p(12)*u)||(min(auxnormF,norm(F,1))*normA4s^3 <= Low_bound_expAs_u)
            F= p(12)*A4+p(11)*A3+p(10)*A2 +p(9)*A;   %One matrix product is saved
        else
            F=     F*A4+p(11)*A3+p(10)*A2 +p(9)*A;
            Cp=Cp+1;
        end
        auxnormF=norm(F,1);
        F=F+p(8)*I;        
        if (auxnormF<=p(8)*u)||(min(auxnormF,norm(F,1))*normA4s^2 <= Low_bound_expAs_u)
            F= p(8)* A4+ p(7)*A3+ p(6)*A2+ p(5)*A;   %One matrix product is saved
        else
            F=     F*A4+ p(7)*A3+ p(6)*A2+ p(5)*A;
            Cp=Cp+1;
        end
        auxnormF=norm(F,1);
        F=F+p(4)*I;       
        if (auxnormF<=p(4)*u)||(min(auxnormF,norm(F,1))*normA4s <= Low_bound_expAs_u)
            F=   p(4)*A4+ p(3)*A3+ p(2)*A2+ p(1)*A+ I;   %One matrix product is saved
        else
            F=      F*A4+ p(3)*A3+ p(2)*A2+ p(1)*A+ I;    
            Cp=Cp+1;
        end

	case 6 %m=20
        p=[1.0, 0.5, 0.166666666666666666666666666667, 0.0416666666666666666666666666667, 0.00833333333333333333333333333333, 0.00138888888888888888888888888889, 0.000198412698412698412698412698413, 0.0000248015873015873015873015873016, 0.00000275573192239858925109505932705, 0.00000027557319223985888275785856999, 0.0000000250521083854417202238661793214, 0.00000000208767569878681001865551494345, 0.000000000160590438368216133408629182949, 0.0000000000114707455977297245072965704504, 0.000000000000764716373181981640551384421244, 0.0000000000000477947733238738525344615263278, 0.0000000000000028114572543455205981105182744, 0.000000000000000156192069685862252711483983227, 0.00000000000000000822063524662401518377011604394, 0.000000000000000000411031762331200788077454968005];
        A3=A*A2;A4=A2*A2;
        Cp=Cp+2;
        normA4s = norm(A4,1);
        F= p(20)*A4+p(19)*A3+p(18)*A2+p(17)*A;
        auxnormF=norm(F,1);
        F=F+p(16)*I;        
        if (auxnormF<=p(16)*u)||(min(auxnormF,norm(F,1))*normA4s^4 <= Low_bound_expAs_u)
            F= p(16)*A4+p(15)*A3+p(14)*A2+p(13)*A;   %One matrix product is saved
        else
            F=     F*A4+p(15)*A3+p(14)*A2+p(13)*A;
            Cp=Cp+1;
        end
        auxnormF=norm(F,1);
        F=F+p(12)*I;        
        if (auxnormF<=p(12)*u)||(min(auxnormF,norm(F,1))*normA4s^3<= Low_bound_expAs_u)
            F= p(12)*A4+p(11)*A3+p(10)*A2+ p(9)*A;   %One matrix product is saved
        else
            F=     F*A4+p(11)*A3+p(10)*A2+ p(9)*A;
            Cp=Cp+1;
        end
        auxnormF=norm(F,1);
        F=F+p(8)*I;
        if (auxnormF<=p(8)*u)||(min(auxnormF,norm(F,1))*normA4s^2 <= Low_bound_expAs_u)
            F=  p(8)*A4+ p(7)*A3+ p(6)*A2+ p(5)*A;   %One matrix product is saved
        else
            F=     F*A4+ p(7)*A3+ p(6)*A2+ p(5)*A;
            Cp=Cp+1;
        end
        auxnormF=norm(F,1);
        F=F+p(4)*I;        
        if (auxnormF<=p(4)*u)||(min(auxnormF,norm(F,1))*normA4s <= Low_bound_expAs_u)
            F=  p(4)*A4+ p(3)*A3+ p(2)*A2+ p(1)*A+ I;   %One matrix product is saved
        else
            F=     F*A4+ p(3)*A3+ p(2)*A2+ p(1)*A+ I;
            Cp=Cp+1;
        end

            
	case 7 %m=25   %Hermite method for order m > 20
        p=[1.0, 0.5, 0.166666666666666666666666666667, 0.0416666666666666666666666666667, 0.00833333333333333333333333333333, 0.00138888888888888888888888888889, 0.000198412698412698412698412698413, 0.0000248015873015873015873015873016, 0.00000275573192239858925109505932705, 0.00000027557319223985888275785856999, 0.0000000250521083854417202238661793214, 0.00000000208767569878681001865551494345, 0.000000000160590438368216133408629182949, 0.0000000000114707455977297245072965704504, 0.000000000000764716373181981640551384421244, 0.0000000000000477947733238740986590639552835, 0.00000000000000281145725434553519203726486311, 0.000000000000000156192069684763147603381520372, 0.00000000000000000822063524656648226372742502475, 0.000000000000000000411031765542346347406064263718, 0.000000000000000000019572941216302206876038034553, 0.000000000000000000000889673352690311238723066279236, 0.0000000000000000000000386814501169700482355951169037, 0.00000000000000000000000161173757109611838589528686931, 0.0000000000000000000000000644695028438447358949889555624];
        A3=A*A2;A4=A*A3;A5=A*A4;
        Cp=Cp+3;
        normA5s = norm(A5,1);
        F= p(25)*A5+p(24)*A4+p(23)*A3+p(22)*A2+p(21)*A;
        auxnormF=norm(F,1);
        F=F+p(20)*I;      
        if (auxnormF<=p(20)*u)||(min(auxnormF,norm(F,1))*normA5s^4<= Low_bound_expAs_u)
            F= p(20)*A5+p(19)*A4+p(18)*A3+p(17)*A2+p(16)*A;   %One matrix product is saved
        else
            F=     F*A5+p(19)*A4+p(18)*A3+p(17)*A2+p(16)*A;
            Cp=Cp+1;
        end
        auxnormF=norm(F,1);
        F=F+p(15)*I;
        if (auxnormF<=p(15)*u)||(min(auxnormF,norm(F,1))*normA5s^3 <= Low_bound_expAs_u)  
            F= p(15)*A5+p(14)*A4+p(13)*A3+p(12)*A2+p(11)*A;   %One matrix product is saved
        else
            F=     F*A5+p(14)*A4+p(13)*A3+p(12)*A2+p(11)*A;
            Cp=Cp+1;
        end
        auxnormF=norm(F,1);
        F=F+p(10)*I;        
        if (auxnormF<=p(10)*u)||(min(auxnormF,norm(F,1))*normA5s^2 <= Low_bound_expAs_u)
            F=  p(10)*A5+ p(9)*A4+ p(8)*A3+ p(7)*A2+ p(6)*A;   %One matrix product is saved
        else          
            F=      F*A5+ p(9)*A4+ p(8)*A3+ p(7)*A2+ p(6)*A;
            Cp=Cp+1;
        end
        auxnormF=norm(F,1);
        F=F+p(5)*I;        
        if (auxnormF<=p(5)*u)||(min(auxnormF,norm(F,1))*normA5s <= Low_bound_expAs_u)
            F=   p(5)*A5+ p(4)*A4+ p(3)*A3+ p(2)*A2+ p(1)*A+  I;   %One matrix product is saved
        else      
            F=      F*A5+ p(4)*A4+ p(3)*A3+ p(2)*A2+ p(1)*A+  I;
            Cp=Cp+1;
        end

         
	case 8 %m=30
        p=[1.0, 0.5, 0.166666666666666666666666666667, 0.0416666666666666666666666666667, 0.00833333333333333333333333333333, 0.00138888888888888888888888888889, 0.000198412698412698412698412698413, 0.0000248015873015873015873015873016, 0.00000275573192239858925109505932705, 0.00000027557319223985888275785856999, 0.0000000250521083854417202238661793214, 0.00000000208767569878681001865551494345, 0.000000000160590438368216133408629182949, 0.0000000000114707455977297245072965704504, 0.000000000000764716373181981640551384421244, 0.0000000000000477947733238738525344615263278, 0.0000000000000028114572543455205981105182744, 0.000000000000000156192069685862252711483983227, 0.00000000000000000822063524662401518377011604394, 0.000000000000000000411031762331200788077454968005, 0.000000000000000000019572941063649948088650082313, 0.000000000000000000000889679139256815736720869136242, 0.0000000000000000000000386817015588920091210520467065, 0.00000000000000000000000161173756495383383582449682343, 0.000000000000000000000000064469559502767552664955459304, 0.00000000000000000000000000247959844241413675133902747082, 0.0000000000000000000000000000918229471527745410777697068915, 0.00000000000000000000000000000327939096974194829600562219607, 1.13099628864477180712860737885e-31, 3.76998762881590606937720898174e-33];
        A3=A*A2;A4=A*A3;A5=A*A4;
        Cp=Cp+3;
        normA5s = norm(A5,1);
        F= p(30)*A5+p(29)*A4+p(28)*A3+p(27)*A2+p(26)*A;
        auxnormF=norm(F,1);
        F=F+p(25)*I;
        if (auxnormF<=p(25)*u)||(min(auxnormF,norm(F,1))*normA5s^5 <= Low_bound_expAs_u)
            F= p(25)*A5+p(24)*A4+p(23)*A3+p(22)*A2+p(21)*A;   %One matrix product is saved
        else
            F=     F*A5+p(24)*A4+p(23)*A3+p(22)*A2+p(21)*A;
            Cp=Cp+1;
        end
        auxnormF=norm(F,1);
        F=F+p(20)*I;      
        if (auxnormF<=p(20)*u)||(min(auxnormF,norm(F,1))*normA5s^4 <= Low_bound_expAs_u)
            F= p(20)*A5+p(19)*A4+p(18)*A3+p(17)*A2+p(16)*A;   %One matrix product is saved
        else
            F=     F*A5+p(19)*A4+p(18)*A3+p(17)*A2+p(16)*A;
            Cp=Cp+1;
        end
        auxnormF=norm(F,1);
        F=F+p(15)*I;
        if (auxnormF<=p(15)*u)||(min(auxnormF,norm(F,1))*normA5s^3 <= Low_bound_expAs_u)
            F= p(15)*A5+p(14)*A4+p(13)*A3+p(12)*A2+p(11)*A;   %One matrix product is saved
        else
            F=     F*A5+p(14)*A4+p(13)*A3+p(12)*A2+p(11)*A;
            Cp=Cp+1;
        end
        auxnormF=norm(F,1);
        F=F+p(10)*I;        
        if (auxnormF<=p(10)*u)||(min(auxnormF,norm(F,1))*normA5s^2 <= Low_bound_expAs_u)
            F= p(10)*A5+ p(9)*A4+ p(8)*A3+ p(7)*A2+ p(6)*A;   %One matrix product is saved
        else
            F=     F*A5+ p(9)*A4+ p(8)*A3+ p(7)*A2+ p(6)*A;
            Cp=Cp+1;
        end
        auxnormF=norm(F,1);
        F=F+p(5)*I;       
        if (auxnormF<=p(5)*u)||(min(auxnormF,norm(F,1))*normA5s <= Low_bound_expAs_u)
            F=      p(5)*A5+ p(4)*A4+ p(3)*A3+ p(2)*A2+ p(1)*A+ I;   %One matrix product is saved
        else
            F=         F*A5+ p(4)*A4+ p(3)*A3+ p(2)*A2+ p(1)*A+ I;
            Cp=Cp+1;
        end
        
        
end
% Compute F=(Hm(A/2^s))^(2^s) (Squaring phase)
for i = 1:s 
    F = F*F;  
    Cp=Cp+1;
end
end



