function [X, obj, y, info] = runSDP(C, A, b, varargin)

% RUNSDP  Solve a semi-definite program given in primal form
% requires: SeDuMi
% author: Toby Cubitt
% license: GPL2
%
%     [X, OBJ, y] = runSDP(C, A, b) solves the following SDP
%     (where >= denotes matrix inequality):
%
%       minimize     trace(C*X)
%       subject to   X >= 0
%                    trace(A{i}*X) == b(i)
%
%     The matrix C defines the objective function, A is a 1-dimensional cell
%     array containing the matrices appearing in the equality constraints, b
%     is a vector containing the corresponding constants.
%
%     X is the optimal solution matrix, OBJ is the value of the objective
%     function at that point, y is the dual solution vector. INFO is a
%     SeDuMi structure containing information about the outcome of the
%     numerical optimisation.
%
%     [X, OBJ, y] = runSDP(C, A, b, 'real') restricts X to be a real
%     symmetric matrix. Otherwise, X can be Hermitian.
%
%     If your problem only requires a single constraint, A is allowed to be
%     a single matrix rather than a cell array (b will be a scalar in this
%     case).
%
%
%     See also RUNLMI, SEDUMI.
%
%     Note that the optimisation problems solved by RUNSDP and RUNLMI are
%     not strictly dual to each other: if [X,f,y] = runSDP(C,A,b), then
%     [-y,-f,X] = runLMI(b,A,-C).


%% Copyright (C) 2011 Toby Cubitt
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License
%% as published by the Free Software Foundation; either version 2
%% of the License, or (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; if not, write to the Free Software
%% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
%% MA 02110-1301, USA.


% sort out input arguments
if isnumeric(A)
  A = {A};
end
if ~ length(b) == length(A);
  error('dimensions of A and b are inconsistent');
end

R = 'hermitian';
if length(varargin) >= 1
  R = varargin{1};
end


% get dimensions etc.
dim = size(C,1);
K.s = dim;
if ~ (strcmp(R,'real') || strcmp(R,'re'))
  K.scomplex = [1];
end


%--- equality constraints ---
At = sparse([]);
for k=1:length(A)
    At = [At;vec(sparse(A{k})).'];
end


%--- objective function ---
C = vec(sparse(C));


%--- run SDP ---
[X,y,info] = sedumi(At, b, C, K);
X = mat(X);
obj = full(trace(mat(C)*X));