function [y,mvec_count,H] = phiRT(A,v,g,t,toler,m)
%
% [y,mvec,H] = phiRT(A,v,g,t,toler,m)
%
% solves approximately to the given tolerance an IVP
% y'(t) = -A*y(t) + g,  y(0) = v
% using Arnoldi method with the RT (residual time) restarting
% as proposed in
%
% M. A. Botchev, L. Knizhnerman, and E. E. Tyrtyshnikov. Residual
% and restarting in Krylov subspace evaluation of the phi function. 
% SIAM J. Sci. Comput., 43(6):A3733–A3759, 2021.
% https://doi.org/10.1137/20M1375383 (see also
% https://arxiv.org/abs/2010.08494)
% 
% input:  A       (n x n)-matrix
%         v       n-vector
%         g       n-vector
%         t>0     length of the time interval
%         toler>0 tolerance
%         m       maximal Krylov dimension
% output: y       the approximate solution
%         mvec    number of matrix-vector multiplications done to
%                 compute y
%         H       last Krylov projection matrix
%   
% Copyright (c) 2020 by M.A. Botchev
% Permission to copy all or part of this work is granted,
% provided that the copies are not made or distributed
% for resale, and that the copyright notice and this
% notice are retained.
%
% THIS WORK IS PROVIDED ON AN "AS IS" BASIS.  THE AUTHOR
% PROVIDES NO WARRANTY WHATSOEVER, EITHER EXPRESSED OR IMPLIED,
% REGARDING THE WORK, INCLUDING WARRANTIES WITH RESPECT TO ITS
% MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE.
%
n = size (g,1);
V = zeros(n  ,m+1);
H = zeros(m+1,m);

symm_A = ( norm(A-A',1)/norm(A,1)<1e-12 );
if symm_A 
   disp('A is symmetric => Lanczos will be used')
end

convergence = false;
mvec_count  = 0;
while ~convergence
   g_rst = g;     % source term for this restart
   if norm(v)
       % change of variables y(t) -> y(t)-v :
      g_rst = g_rst - A*v;
   end
   % beta = norm(g);  % wrong
   beta = norm(g_rst);  % correct
   V(:,1) = g_rst/beta;

   for j=1:m
       w = A*V(:,j);
       mvec_count = mvec_count + 1;
       j_start = 1;
       if symm_A
          j_start = max([j-1,1]);
       end
       for i=j_start:j
           H(i,j) = w'*V(:,i);
           w      = w - H(i,j)*V(:,i);
       end
       H(j+1,j) = norm(w);
       e1 = zeros(j,1); e1(1) = 1;
       ej = zeros(j,1); ej(j) = 1;
       s  = [1 : 6]/6*t;
       % mb:
       src_term   = s(1)*phim(-s(1)*H(1:j,1:j),e1);
       u_true      = src_term;
       beta_j1(1) = -H(j+1,j)* (ej'*u_true);
       expmH = expm( -s(1)*H(1:j,1:j) );
       for q=2:length(s)
           u_true      = expmH*u_true + src_term;
           beta_j1(q) = -H(j+1,j)* (ej'*u_true);
       end
       resnorm = norm(beta_j1,'inf');
       if resnorm<=toler
           fprintf('\nj = %d, resnorm = %.2e - convergence!\n',j,resnorm);
           convergence = true;
           break
       elseif j==m
           n_tsteps    = 100;   
           resid_value = 2*toler;
           while resid_value>toler
               
               src_term = (t/n_tsteps)*phim(-(t/n_tsteps)*H(1:j,1:j),e1);
               u_true   = src_term;
               resid_value = abs(-H(j+1,j)* (ej'*u_true));

               if resid_value<=toler
                  break
               end
               n_tsteps = 2*n_tsteps;
           end
           % keyboard     % to plot residual vs time - MB
           % Compute residual for intermediate time points until its
           % value exceeds tolerance
           u_old = zeros(j,1);
           expmH = expm( -(t/n_tsteps)*H(1:j,1:j) );
           for k=1:n_tsteps
               u_true = expmH*u_old + src_term;
               
               resid_value = abs(-H(j+1,j)* (ej'*u_true));
               if resid_value>toler
                   u_true = u_old;
                   t_ok   = (k-1)/n_tsteps*t;        
                   y_ok   = V(:,1:j)*(beta*u_true);
                   % change of variables y(t) -> y(t) + v :
                   y_ok = y_ok + v;
                   fprintf('\rtime interval reduced by %6.2f%%',t_ok/t*100);
                   t = t - t_ok;
                   
                   v = y_ok;
             
                   break          % restart
               end 
               u_old  = u_true;
               
           end
           break                  % restart
       end
       V(:,j+1) = w/H(j+1,j);
   end
end
   
H = H(1:j,1:j);
y = V(:,1:j)*(beta*u_true);
y = y + v;

% ------------------------------------------------------------------------------------
function y = phim(A,x)
%
% y = phim(A,x)
% 
% computes y = phi(A)*x
%
use_schur    = true;

n            = size(A,1);
Aug          = zeros(n+1,n+1);
Aug(1:n,1:n) = A;
Aug(1:n,n+1) = x;

if use_schur
   [Ua,Sa] = schur(Aug);
   expmAug = Ua*expm(Sa)*Ua';
else
   expmAug = expm(Aug);
end

y = expmAug(1:n,n+1);
