function [X,iter,SEQ]=bbkarcher(varargin)

% X=BBKARCHER(A1,...,Ap) computes the Karcher mean of positive 
%  definite matrices A1,...,Ap, using a Riemannian Barzilai and
%  Borwein (BB) algorithm where the initial value is the arithmetic mean
% X=BBKARCHER(A{1:p}) for a cell-array input
% Do not use X=BBKARCHER(A) for a cell-array A, use BBKARCHER(A{1:p}) instead
% [X iter SEQ]=BBKARCHER(A1,...,Ap) computes also the sequence of the iterates,
%   in debug mode, with no stopping criterion
%
% varargin: positive definite matrix arguments A1,...,Ap
% X: the Karcher mean of A1,...,Ap
% iter: the number of iterations needed by the iteration
% SEQ: the sequence of the iterates
%
% Code by B. Iannazzo 2016 (Updated on 2019)
% 
% Reference
% B. Iannazzo and M. Porcelli. The Riemannian Barzilai-Borwein method with
%   nonmonotone line-search and the Karcher mean computation. 
%   IMA J. Numer. Anal. 38-1 (2018), pp. 495-517. 

% global settings
maxiter=200; % maximum number of iterations
tol=100*eps; % tolerance for the stopping criterion
ni(1)=Inf; % useful for the stopping criterion

% read the arguments, computes the arithmetic mean as initial value
% and compute the Cholesky factorizations for future usage
p=nargin;
n=size(varargin{1},1);
X=zeros(n);
for h=1:p
  A{h}=varargin{h};
  X=X+A{h}/p;
  R{h}=chol(A{h});
end
if nargout>2 
  SEQ{1}=X; % only for debug
end

% the first step-length is obtained using the Richardson-like strategy
Rx=chol(X);
iRx=inv(Rx);
Sum=zeros(n);
for h=1:p
  Z=R{h}*iRx;   
  [Uz{h} Dz]=schur(Z'*Z);
  V{h}=diag(Dz);
  Sum=Sum+Uz{h}*diag(log(diag(Dz)))*Uz{h}';
end

% the next decomposition is for future use
[Q Ds]=schur(Sum);

% computes the gradient explicitly
Gold=-2*Rx'*Sum*Rx;

% computes the first steplength
beta=0;gamma=0;
for h=1:p
  ch=max(V{h})/min(V{h});
  if (abs(ch-1)<0.5)
    dh=log1p(ch-1)/(ch-1);
  else
    dh=log(ch)/(ch-1);
  end
  beta=beta+dh;gamma=gamma+ch*dh;
end
alfa=1/(gamma+beta);

% start the main iteration
for k=1:maxiter
  
  % compute the new value of X using the old gradient
  Z=Q'*Rx;
  X=Z'*diag(exp(diag(2*alfa*Ds)))*Z;
  Rx=chol(X);
  if nargout>2
    SEQ{k+1}=X; % just for debug
  end

  % computes S and trasport the old gradient
  S=2*alfa*Z'*Ds*diag(exp(diag(2*alfa*Ds)))*Z;  
  Gold_trasp=-1/alfa*S;
  
  % computes the gradient of f at the new value of X
  iRx=inv(Rx);
  Sum=zeros(n);
  for h=1:p
	Z=R{h}*iRx;
	[U D]=schur(Z'*Z);
	Sum=Sum+U*diag(log(diag(D)))*U';
  end
  Gnew=-2*Rx'*Sum*Rx;

  [Q Ds]=schur(Sum);

  % compute Y
  Y=Gnew-Gold_trasp;
  
  % compute the tbb and the new alfa
  TMP=iRx'*S*iRx;
  tbb=-2*trace(TMP*Sum)/trace(TMP^2)+1/alfa;
  alfa=1/tbb;
    
  % the stopping criterion is based on the norm of Gnew
  % stopping criterion is not applied if in debug mode
  if nargout<3
    nGnew=norm(Gnew);
    ni=sort([ni nGnew]);
    if ( (nGnew<norm(X)*tol) || (nGnew>ni(2)) )
      iter=k;break
    end
  end

  if (k==maxiter)
    disp('Max number of iterations reached');
    iter=k; break;
  end
  
end