function d = drift(A)
% On input: the 3-dimensional array A(:,:,:) associated with an M/G/1 type
% Markov chain; the third index is the block index.
% On output: the drift d ( d<0 --> positive recurrent, d>0 --> transient) 
m = size(A,1);
k = size(A,3);
B= sum(A, 3);
x = pvgth(B);
v = zeros(m,1);
for j=1:k
    v=v+(j-2)*sum(A(:,:,j),2);
end
d = v'*x;
