function M = euclideancomplexfactory(m, n)
% Returns a manifold struct to optimize over complex matrices.
%
% function M = euclideancomplexfactory(m)
% function M = euclideancomplexfactory(m, n)
% function M = euclideancomplexfactory([n1, n2, ...])
%
% Returns M, a structure describing the vector space of complex matrices,
% as a manifold for Manopt.
%
% The complex plane is here viewed as R^2. The inner product between two
% m-by-n matrices A and B is given by: real(trace(A'*B)). This choice
% guides the proper definition of gradient and Hessian for this geometry.
% This is not the classical Euclidean inner product for complex matrices;
% it is a real inner product.
%
% See also: euclideanfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, April 7, 2015.
% Contributors: 
% Change log: 
%
%   Jan. 25, 2017 (NB):
%       Added functionality to handle multidimensional arrays.

    % The size can be defined using both m and n, or simply with m.
    % If m is a scalar, then n is implicitly 1.
    % This mimics the use of built-in Matlab functions such as zeros(...).
    if ~exist('n', 'var') || isempty(n)
        if numel(m) == 1
            n = 1;
        else
            n = [];
        end
    end
    
    dimensions_vec = [m(:)', n(:)']; % We have a row vector.
    
    M.size = @() dimensions_vec;

    M.name = @() sprintf('Euclidean space C^(%s)', num2str(dimensions_vec));
    
    M.dim = @() 2*prod(dimensions_vec);
    
    M.inner = @(x, d1, d2) real(d1(:)'*d2(:));
    
    M.norm = @(x, d) norm(d(:), 'fro');
    
    M.dist = @(x, y) norm(x(:)-y(:), 'fro');
    
    M.typicaldist = @() sqrt(prod(dimensions_vec));
    
    M.proj = @(x, d) d;
    
    M.egrad2rgrad = @(x, g) g;
    
    M.ehess2rhess = @(x, eg, eh, d) eh;
    
    M.tangent = M.proj;
    
    M.exp = @exp;
    function y = exp(x, d, t)
        if nargin == 3
            y = x + t*d;
        else
            y = x + d;
        end
    end
    
    M.retr = M.exp;
    
    M.log = @(x, y) y-x;

    M.hash = @(x) ['z' hashmd5([real(x(:)) ; imag(x(:))])];
    
    M.rand = @() (randn(dimensions_vec) + 1i*randn(dimensions_vec))/sqrt(2);
    
    M.randvec = @randvec;
    function u = randvec(x) %#ok<INUSD>
        u = randn(dimensions_vec) + 1i*randn(dimensions_vec);
        u = u / norm(u(:), 'fro');
    end
    
    M.lincomb = @matrixlincomb;
    
    M.zerovec = @(x) zeros(dimensions_vec);
    
    M.transp = @(x1, x2, d) d;
    
    M.pairmean = @(x1, x2) .5*(x1+x2);
    
    sz = prod(dimensions_vec);
    M.vec = @(x, u_mat) [real(u_mat(:)) ; imag(u_mat(:))];
    M.mat = @(x, u_vec) reshape(u_vec(1:sz), dimensions_vec) ...
                        + 1i*reshape(u_vec((sz+1):end), dimensions_vec);
    M.vecmatareisometries = @() true;

end
