function [X, y, Z, C, A, b, X_0, y_0] = testsrssdp(n, m, toler, precond);
% generate C, A, b randomly (along with initial solutions X_0 and y_0) and solve the following problem
% min <C,X> subject to A(X) = b, X >= 0
% n: size of matrix X
% m: number of constraints
% toler: relative tolerance
% A: cell of m symmetric matrices
% C: symmetric objective matrix
% b: vector m by 1 right-hand side
% X_0, y_0: initial solutions

rand('state',sum(100*clock));
%construct random matrix A
%densitymA = .5/(n+m);
densitymA = .25/n;
n2 = n^2;
tn = n*(n+1)/2;
r = 0;
iter = 1;
while (r < m & iter<10)
	mA = tril(sprandn(m,m,densitymA))+speye(m); %first component is a lower triangle matrix
	mA = [sprandn(m,(tn-m),densitymA) mA];
	%check rank of mA
	r = sprank(mA);
	if r < m
		disp('The constraints are not linearly independent. Generate new ones ...');
		iter = iter + 1;
	end 
	if (iter==10)
		disp('Please try again with smaller m');
		return;
	end
end
%spy(mA);
disp(['[number of zero columns in mA, cols mA is ',...
                         num2str([sum(sum(mA)==0) size(mA,2)])]);
svind = find(triu(ones(n),1)); %linear indices of the strictly upper triangle
sMind = zeros(n2,1); %matrix indices with respect to linear indices of the diagonal and the strictly upper triangle
sMind(svind) = 1:tn-n; %strictly upper triangle
sMind(1:n+1:n2) = tn-n+1:tn;%the diagonal
%calculate the linear indices of the strictly lower triangle
indM = zeros(n); %matrix indices
indM(:) = 1:n2;
indM = indM';
lsvind = indM(svind); %strictly lower triangular 
sMind(lsvind) = 1:tn-n; %strictly lower triangular
%construct random matrices A
A = cell(m,1);
for i=1:m,
    A{i}=sssMat(mA(i,:),n,sMind);  
end
%construct random righ-hand side and intial solution X_0>=0
%B = -2*rand(n) + 1;
[q,r] = qr(randn(n));
X_0 = q*diag(1+10*rand(n,1))*q';
X_0 = (X_0+X_0')/2;
%X_0 = B*B';
b = zeros(m,1);
for i=1:m
	b(i) = trace(A{i}*X_0);
end
%construct random objective and initial solution y_0
y_0 = -2*rand(m,1) + 1;
[q,r] = qr(randn(n));
Z_0 = q*diag(1+10*rand(n,1))*q';
Z_0 = (Z_0+Z_0')/2;
%B = -2*rand(n) + 1;
%Z_0 = B*B';
C = Z_0;
for i=1:m
	C = C + y_0(i)*A{i};
end
%solve the problem using rsdp with all default parameters
profile on
tic
[X, y, Z, muu, itersinfo, meritfn] = srssdp(C, A, b, toler, precond);
toc
profile report
