Save previous iteration value for comparison with next iteration

Question

I don't know if this is the right place (at most they will close the post). Anyway, I am trying to implement the Matrix Rotation Count Algorithm explained in this paper (in particular pages 11-12). In a nutshell, this algorithm makes it possible to correctly implement the complex matrix logarithm, but I have a problem with step 8 that requires a comparison between the values of two previous and next iterations.

Now, for all steps (except step 8) I have no problem, but in step 8 I need to compare a value of the previous iteration m_{k-1,i} with the value of the present iteration m_{k,i}. The value I want to save for comparison at each iteration is m_{k-1,i}, however I have no idea how to save the previous iteration to compare with the next one.

In the piece of code (which I attach below) I tried two attempts: in the first I tried working directly with diagonal matrices (since they are consistent with the result I need in the end), in the second attempt I tried working for each value of the diagonal matrix. However, in both cases I could not find a way to compare the two values.

Ideally, for step 8 I should start a comparison from the second argument (since at first I would have nothing to compare it to but itself).

If anyone has any valuable suggestions or a solution sketch I would be very grateful, thanks in advance

clear variables; close all; clc

% -----------------------------
% Parameters
% -----------------------------
S0 = 100; % Initial stock price
r = 0.05; % Risk free rate 
q = 0; % Dividend yield
t = 1.0; % Time to maturity
beta = 3; % Role of Feller condition

M = [-3, 0.0; 0.0, -3];         
R = [-0.7, 0.0; 0.0, -0.7];    
Q = [0.25, 0.0; 0.0, 0.25];   
Sigma0 = [0.01, 0.0; 0.0, 0.01];  % Initial variance 


i = complex(0,1);
gamma = linspace(0,5,10);


% ------------------------------------------------------------------
% Matrix Rotation Count Algorithm ~~~ 1
% ------------------------------------------------------------------
% ---------- step (1) -----------% 
% Initialize the number of rotations
rotIdx = 0; 
% ---------- step (2) -----------% 
for k = 1:length(gamma)
    w = gamma(1,k);
    % ---------- step (3) -----------% 
    % Compute A22 matrix
    MATEXP = expm(t*[M, -2.0*(Q'*Q); ...
             0.5*i*w*(i*w - 1)*eye(2), -1.0*(M' + 2.0*i*w*R*Q)]);
    A22 = [MATEXP(3:4,3:4)];
    % ---------- step (4) -----------% 
    % Compute the PDP decomposition 
    % where D is the diagonal matrix of eigenvalues 
    [Pk,Dk] = eig(A22); 
    % Vector that contain the eigenvalues 
    lambda_eig = A22(sub2ind(size(A22),1:size(A22,1),1:size(A22,2))); 
    % ---------- step (6) -----------%  
    DkLog = logm(diag(lambda_eig)); % diagonal matrix with eigenvalues
    % ---------- step (7) -----------% 
    Mki = mod(imag(DkLog),pi);
    % ---------- step (8) -----------% 
    % TODO
    % if 1 < 0.5*pi
    %     rotIdx = rotIdx + 1; % positive rotation
    % elseif 2 > -0.5*pi
    %     rotIdx = rotIdx - 1; % negative rotation
    % end
    % ---------- step (9) -----------% 
    % compute the correct branch of the imaginary part of complex log
    DkLog = real(DkLog) + i * (Mki + eye(2)*pi*rotIdx);
    % ---------- step (10) -----------% 
    % Compute the correct complex matrix log
    A22log = Pk * DkLog * inv(Pk);
end



%{
% ------------------------------------------------------------------
% Matrix Rotation Count Algorithm ~~~ 2
% ------------------------------------------------------------------
rotIdx = 0;             % step (1)
for k = 1:length(gamma) % step (2)
    w = gamma(1,k);
    MATEXP = expm(t*[M, -2.0*(Q'*Q); ...
             0.5*i*w*(i*w - 1)*eye(2), -1.0*(M' + 2.0*i*w*R*Q)]); % step(3)
    A22 = [MATEXP(3:4,3:4)];
    % ---------- step (4) -----------% 
    % compute the PDP decomposition 
    % where D is the diagonal matrix of eigenvalues 
    [Pk,Dk] = eig(A22); 
    % Vector that contain the eigenvalues 
    lambda_eig = A22(sub2ind(size(A22),1:size(A22,1),1:size(A22,2))); 
    Dklog = zeros(2,2);
    for j = 1:length(lambda_eig) % step (5)
        Dki = diag(lambda_eig); % diagonal matrix with eigenvalues
        % ---------- step (6) -----------% 
        % Evaluate the complex logarithm
        dki = log(val_eig); 
        % ---------- step (7) -----------% 
        % Produce the sawtooth-like function
        mki = mod(imag(dki),pi); 
        % ---------- step (8) -----------% 
        % Check whether rotation has occured
        temp = 1.0; % insert m_k-1_i
        if mki - temp < 0.5*pi
            rotIdx = rotIdx + 1; % positive rotation
        elseif mki - temp < -0.5*pi
            rotIdx = rotIdx - 1; % negative rotation
        end      
        % ---------- step (9) -----------% 
        % compute the correct branch of the imaginary part of complex log
        IMdki = mki + pi * rotIdx;
    end
    % A22log = Pk * ... * inv(Pk);
end
%}


% Work on diagonal
XX = rand(5);
yy = 1:5;
n = size(XX,1);
XX(1:(n+1):end) = yy;

Here I attach a picture of the algorithm