function  [X, y, Z, muu, itersinfo, meritfn] = srssdp(C, A, b, toler, precond, crosstoler, lsqr, itoler, X_0, y_0)
%function  [X, y, Z, muu, itersinfo, meritfn] = srssdp(C, A, b, toler, precond, crosstoler, lsqr, itoler, X_0, y_0)
%% 	INPUT:
%
%%   C  -  given real symmetric n by n matrix
%%   A  -  given m real symmetric n by n matrices in the form of Matlab cell structure
%	 b  -  given real m by 1 vector  
%
%%   toler - relative tolerance for stopping criteria of the pd-ip algorithm
%    precond - indicator for type of preconditioning 
%                                   0 for no preconditioning (default)
%                                   1 currently for opt diag precond 
%									2 currently for Q-less QR block diagonal precond
%    crosstoler - (optional input) tolerance for start of 
%                         			crossover to sigmaa=0 and no pos. check. 
%                       			default is .1 (see below)
%	 lsqr - indicator for least-square solver
%									0 for lsmr (default)
%									1 for lsqr
%	 itoler - parameter to select initial solution for X and Z ((1/toler)^{itoler}I)
%									default is itoler = 0 (I is the initial solution)
%	 X_0 - optional initial primal feasible solution
%	 y_0 - optional initial dual feasible solution
%
%% 	OUTPUT:
%%	 
%%	 X - optimum matrix
%%   y, Z - optimum dual vector,matrix
%%   muu - duality gap at end (trace XZ)
%%   itersinfo - 2 dim vector with
%%           	(number of iterations used, crossover iteration no.)
%%   meritfn - if included as desired input (optional), this is a 
%%             	vector with the value of ||F||| at each iteration
%%             	and a plot at completion.
%%
% 	PARAMETERS:
%
%      iterbnd   - upper bound on number of iterations
%
% 	Copyright (c) 2010 by
%		Xuan Vinh Doan, Serge Kruk, and Henry Wolkowicz
%
% 	DEVELOPMENTS:
%
% 	 December 16, 2009:
%          Start to modify and develop the algorithm.
%	 January 27, 2010:
%		   Include the option for infeasible starting point
%	 January 30, 2010:
%		   Change how to perform ssMat with sMind indices (inverse of svind indices)
%		   Change bases for both row and null space of A: A = [I AQ] and Q = [-AQ; I]
% 		by finding a well-conditioned m by m block in the original A
%		   Pass both permutation indf and the inverse permutation iindf (for the well-conditioned
%		block in A) in ssdpj.m to calculate the Jacobian and its adjoint
%		   Keep X and Z as current available and pass them to ssdpj.m
%	 February 05, 2010:
%		   Change to sparse version by adding prefix s- for all files
%		   Keep bases as A = [S E] and Q = [-S\E; I] under a permutation
%		   Pass S and E to lsqr instead of AQ
%	 February 08, 2010:
%		   Find the permutation both in row and column of A so that A = [S E], S is approximately upper triangle
%		   Change A, b, and y_0 according to permutation in row
%	 February 11, 2010:
%		   Try to find sparse representation of A using msparse.m (with 1-norm optimization), the result is good
%		   but 1-norm optimization problems take too long to compute
%    March 08, 2010:
%		   Change the way how to find permutations for row and column in A (mtriangle.m) by taking into consideration
%		   the condition numbers of resulting diagonal matrices
%	 March 10, 2010:
%		   Develop diagonal preconditioner. Need to update ssdplsqr.m and sssdpj.m 
%    March 12, 2010:
%		   Develop block diagonal preconditioner. sssdpx.m and sssdpz.m calculate X(.) and Z(.). Q-less QR is used
%	 March 18, 2010:
% 		   Change crossover condition to include measure of infeasibility. Add display information
% 	 March 22, 2010:
% 		   Develop block diagonal preconditioner using incomplete Cholesky. Update sssdpj.m
%	 March 28, 2010:
%		   Use ssdpxi.m and ssdpzi.m to calculate matrix for J, X, and Z transformations to use in preconditioners
%	 June 07, 2010:
%		   Change from lsqr to lsmr. Use ssdplsmr.m code
%	 August 24, 2010:
%		   Update calculation of preconditioners, remove all matrix transposes in sssdpj.m
%	 August 25, 2010:
%		   Add DIMACS error measures at the end
%	 August 26, 2010:
%		   Remove full-rank checking for matrix mA (to test with larger matrices)
%	 August 31, 2010:
%		   Keep geometric scale for step sizes as in general version
%	 October 26, 2010:
%		   Change time calculation from cputime to tic and toc 
%	 October 29, 2010:
%		   Add option to use either lsmr or lsqr
%	 	   Add option to use (1/toler)^{itoler}I for both Z and X as initial solution
%	 November 25, 2010:
%		   Complete version 0.1
%

n = size(C,1); %matrix size
m = length(b); %number of constraints
n2 = n^2; %number of equations in G-N formulation
tn = n*(n+1)/2; %number of variables in G-N formulation
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';
rvind = indM(svind); %strictly lower triangular 
sMind(rvind) = 1:tn-n; %strictly lower triangular

%check the arguments
if ~exist('crosstoler','var')
   crosstoler = .1;    %default tolerance to start crossover
end
if ~exist('precond','var')
   precond = 0;        %default is no preconditioning
end
if ~exist('toler','var')
   toler = 1e-09;	   %default tolerance
end
if ~exist('lsqr','var')
	lsqr = 0;		%default solver is lsmr
end
if ~exist('itoler','var')
	itoler = 0;		%default initial solution is identity matrix
end
wantmeritfn = nargout>=6;

disp(['solve new problem with crosstoler dim-n, tn1=n(n+1)/2 tn1^2: ',...
   num2str(crosstoler),'  ',num2str(n),'  ',num2str(n*(n+1)/2),...
   ' ',num2str((n*(n+1)/2)^2)]);

%initial parameters
sigmaa = 1;     % centering parameter - changes adaptively
step = 1;       % initial stepsize - changes adaptively
				% initially was 2 - using 1 for cond. numb tests
tau1 = .05;     % initial backtrack value
				% using: step = min(.99997,1-min(tau1,abs(tau1*gapp)))*step;
iterbnd = 50;   % bound on number of major iterations
nnzchol = 0;    % nnz for cholinc

%%Initializations
mA = zeros(m,tn); %matrix form of linear operator A
for i=1:m
	mA(i,:) = sssvec(A{i},svind)';
end
%check the rank of A
% if (rank(mA)<m)
	% disp('A is not full-rank');
	% keyboard;
% end
%sparsify mA
mA = sparse(mA);
%find an orthonomal basis for A
%Q = null(mA);
fnum = 0.0;
indf = 1:tn;
%find a square submatrix with low condition number
disp('Find new bases for the row space and null space of A ...');
[mp, sp] = mtriangle(mA);
%change rows of A
mA = mA(mp,:);
mb = b(mp); % change b
my_0 = zeros(m);
if exist('y_0','var')
	my_0 = y_0(mp); % change y_0
end
%the inverse permutation of mp
imp = zeros(m,1);
imp(mp) = 1:m;

%change the columns
indf(1:m) = sp;
indf(m+1:end) = setdiff(1:tn,sp);
%A = [S E] partition
S = sparse(mA(:,indf(1:m)));
E = mA;
E(:,indf(1:m)) = []; %remove S
E = sparse(E);
% basis for the null space of A
Q = sparse([ -S\E ; speye(tn-m)]);
%the inverse permutation
iindf = zeros(tn,1);
iindf(indf) = 1:tn;
%shift the rows of Q
Q = Q(iindf,:);
%construct Q in matrix form
mQ = cell(tn-m,1);
for i=1:tn-m
	mQ{i} = sssMat(Q(:,i),n,sMind);
end
%sparsify Q
Q = sparse(Q);

%the constant value x_0 and c
%tmA = mA';
%x_0 = tmA*((mA*tmA)\mb);
x_0 = mA\mb;
c = sssvec(C,svind);
[Vc,lc] = eig(C); %eigenvalue decomposition of C

v = zeros(tn-m,1); %initial solution v
%there is no initial primal solution
if ~exist('X_0','var')
	x = sssvec(1/toler^itoler*eye(n),svind);
else
	%check semidefiniteness of primal solution
	[Rx,choltrueX] = chol(X_0);
	if choltrueX > 0 %switch to primal infeasibility
		disp('The initial solution X_0 is not positive semidefinite. Switch to primal infeasible approach');
		x = ssvec(1/toler^itoler*eye(n),svind);
	else
		x_0 = sssvec(X_0,svind);
		x = x_0;
	end
end
X = sssMat(x,n,sMind); %the original matrix variable
%there is no initial dual solution
if ~exist('y_0','var')
	y = zeros(m,1);
	Z = 1/toler^itoler*eye(n);
	z = sssvec(Z,svind);
else
	%check semidefiniteness of dual solution
	z = c - mA'*my_0;
	Z = sssMat(z,n,sMind);
	% Z = C;
	% for i=1:m
		% Z = Z - y_0(i)*A{i};
	% end
	% z = sssvec(Z,svind);
	[Rz,choltrueZ] = chol(Z);
	if choltrueZ > 0 %switch to dual infeasibility
		disp('The initial solution Z_0 is not positive semidefinite. Switch to dual infeasible approach');
		y = zeros(m,1);
		Z = 1/toler^itoler*eye(n);
		z = sssvec(Z,svind);
	else
		y = my_0;
	end
end
%calculate the residuals
rx = x - x_0 - Q*v;
rd = z - c + mA'*y;		
%calculate the optimality matrix ZX
ZX = Z*X;
gapp = trace(ZX);
muu=gapp/n; %initial parameter muu
ZXnorm = norm(ZX,'fro'); %for crossover
%calculate the current objective
optval = c'*x;

noiter = 0;
crossover = 0; %indicator for crossover algor.
itersinfo(2) = 0; %in case no crossover takes place

if wantmeritfn, 
	meritfn = []; 
end

stagnant = 0; %flag to check if the gap is not improved

%display header information
disp('noiter   -log10(relgapp)    step        sigmaa         optval         lsqrtime      cputimeprec    lsqr.iter  nnz-X  nnz-chol');
%Dy = [];
while ( ~crossover*(gapp/(abs(optval)+1) > toler)  & (noiter<iterbnd)  )...
         |(crossover*(ZXnorm)/(abs(optval)+1) > toler& (noiter<iterbnd));

	noiter=noiter+1;
	%compute RHS at current iterate and current muu estimate
	Fmu = ZX - sigmaa*muu*eye(n);
    Fmu = Fmu - Z*sssMat(rx,n,sMind);
    Fmu = Fmu - sssMat(rd,n,sMind)*X;
	
	%disp(['Residual norms ' num2str(norm(rx)) ' ' num2str(norm(rd))]);

	%find the search direction at the current iterate
	tfordirc = cputime;
	damp = 0;
	localsize = 0;
   %  change accuracy: adaptive depending on muu!!! and no iterations
   %%  not working properly???i.e. precond fails when adaptive - why???
   %%                 needs more accuracy for large dimensions??
   %%            ?????to add: nnzQ influence with dimensions????
   %atol   = 1.0e-6;
   %btol   = 1.0e-6;
   %%%?????need to include change in atol,btol after crossover?????
%   atol   = 1e-14;
%   btol   = max(1.0e-1*toler,...
%             (1.0e-6)*(min(muu,1))^(1/3)  ); 
%   btol = max(btol,(1e-14)/muu);
%%%%%%%%%%%%%%%%%%sept/04 working tolerances below%%%%%%%%%%%%%
    atol = max(1.0e-1*toler,...
			(1.0e-6)*(min(muu,1))^(1/3)  );  
    btol = max(1.0e-1*toler,...
			(1.0e-6)*(min(muu,1))^(1/3)  ); 
%%%%%%%%%%%%%%%%%%sept/04 working tolerances above%%%%%%%%%%%%%
   %atol   = max(1.0e-1*toler,...
   %          (1.0e-6)*sqrt(min(muu,1)));  %adaptive not working???why??
   %btol   = max(1.0e-1*toler,...
   %          (1.0e-6)*sqrt(min(muu,1)));
	conlim = 1.0e+10;
	itnlim = 10*tn;
	show   = 0;
   
	%PRECONDITIONING
    stimeprec = tic;
	Tv = [];
	Ty = [];
	%diagonal preconditioner
	if precond == 1
		Tv = zeros(tn-m,1);
		Ty = zeros(m,1);
		%calculate Jacobian J
		for i=1:tn-m
			%ei = zeros(tn-m,1);
			%ei(i) = 1;
			%Ji = sssdpj(1, n2, tn, ei, svind, sMind, X, Z, S, E, indf, iindf, 0, Tv, Ty);
			%Ji = sssdpzi(n, ei, sMind, Z, S, E(:,i), iindf); % J(ei) = Z(ei)
			Ji = Z*mQ{i};
			Tv(i) = norm(Ji,'fro');
			%keyboard;
		end
		for i=1:m
			%ei = zeros(tn,1);
			%ei(tn-m+i) = 1;
			%Ji = sssdpj(1, n2, tn, ei, svind, sMind, X, Z, S, E, indf, iindf, 0, Tv, Ty);
			%Ji = sssdpxi(n, sMind, X, S(i,:), E(i,:), iindf); % J(ei) = X(ei)
			Ji = A{mp(i)}*X;
			Ty(i) = norm(Ji,'fro');
			%keyboard;
		end
    else
		%block diagonal preconditioner
		if precond == 2
			%calculate transformation matrix for Z: R^{tn-m} -> M^n
			mZ = zeros(n2,tn-m);
			for i=1:tn-m
				ei = zeros(tn-m,1);
				ei(i) = 1;
				%Zi = sssdpz(n, ei, sMind, Z, S, E, iindf);
				%Zi = sssdpzi(n, ei, sMind, Z, S, E(:,i), iindf);
				Zi = Z*mQ{i};
				mZ(:,i) = Zi(:);
			end
			%Q-less QR decomposition
			Tv = triu(qr(mZ));
			Tv = Tv(1:tn-m,:);
			%cond(Tv)
			%calculate transformation matrix for X: R^m -> M^n
			mX = zeros(n2,m);
			for i=1:m
				%ei = zeros(m,1);
				%ei(i) = 1;
				%Xi = sssdpx(n, ei, sMind, X, S, E, iindf);
				%Xi = sssdpxi(n, sMind, X, S(i,:), E(i,:), iindf);
				Xi = -A{mp(i)}*X;
				mX(:,i) = Xi(:);
			end
			%Q-less QR decomposition
			Ty = triu(qr(mX));
			Ty = Ty(1:m,:);
			%cond(Ty)
		% else
			% if precond == 3
				% %calculate transformation matrix for Z: R^{tn-m} -> M^n
				% mZ = zeros(n2,tn-m);
				% for i=1:tn-m
					% ei = zeros(tn-m,1);
					% ei(i) = 1;
					% %Zi = sssdpz(n, ei, sMind, Z, S, E, iindf);
					% %Zi = sssdpzi(n, ei, sMind, Z, S, E(:,i), iindf);
					% Zi = Z*mQ{i};
					% mZ(:,i) = Zi(:);
				% end
				% %incomplete Cholesky decomposition
				% options = [];
				% options.droptol = toler;
				% options.michol = 0;
				% options.rdiag = 1;
				% tmZ = mZ';
				% ctmZ = tmZ*mZ;
				% ctmZ = sparse(ctmZ);
				% Tv = cholinc(ctmZ,options);
				% %Tv = cholinc(sparse(mZ'*mZ),options);
				% %Tv = chol(mZ'*mZ);
				% %condest(Tv)
				% %calculate transformation matrix for X: R^m -> M^n
				% mX = zeros(n2,m);
				% for i=1:m
					% %ei = zeros(m,1);
					% %ei(i) = 1;
					% %Xi = sssdpx(n, ei, sMind, X, S, E, iindf);
					% %Xi = sssdpxi(n, sMind, X, S(i,:), E(i,:), iindf);
					% Xi = -A{mp(i)}*X;
					% mX(:,i) = Xi(:);
				% end
				% %incomplete Cholesky decomposition
				% tmX = mX';
				% ctmX = tmX*mX;
				% ctmX = sparse(ctmX);
				% Ty = cholinc(ctmX,options);
				% %Ty = cholinc(sparse(mX'*mX),options);
				% %condest(Ty)
			% end
		end
	end
	timeprec=toc(stimeprec);       
   %%%%%%%%%%%%%%%%%%%
   % [ deltas, istop, lsqriter, rnorm, xnorm      ] =  ...
    % ssdplsqr( n2, tn, 'sssdpj', svind, sMind, X, Z, S, E, indf, iindf, 0, Tv, Ty, -Fmu(:), damp, atol, btol, conlim, itnlim, show );
	% deltav = deltas(1:tn-m);
	% deltay = deltas(tn-m+1:end);
	% deltav
	% deltay
	%assign data to calculate J and J* at each point (v,y)
    %xz = [x z]; 
	stlsqr = tic;
	%%[ deltas, istop, lsqriter, rnorm, xnorm, var ] =  ...  %  for var variable
	if (lsqr==1)
		[ deltas, istop, lsqriter, rnorm, xnorm      ] =  ...
		ssdplsqr( n2, tn, 'sssdpj', svind, sMind, X, Z, S, E, indf, iindf, precond, Tv, Ty, -Fmu(:), damp, atol, btol, conlim, itnlim, show );
	else
		% if (lsqr==0)
			[ deltas, istop, lsqriter, rnorm, normAr, normA, condA, xnorm      ] =  ...
			ssdplsmr( n2, tn, 'sssdpj', svind, sMind, X, Z, S, E, indf, iindf, precond, Tv, Ty, -Fmu(:), damp, atol, btol, conlim, itnlim, localsize, show );
		% else
			% %calculate transformation matrix for Z: R^{tn-m} -> M^n
			% mZ = zeros(n2,tn-m);
			% for i=1:tn-m
				% ei = zeros(tn-m,1);
				% ei(i) = 1;
				% %Zi = sssdpz(n, ei, sMind, Z, S, E, iindf);
				% %Zi = sssdpzi(n, ei, sMind, Z, S, E(:,i), iindf);
				% Zi = Z*mQ{i};
				% mZ(:,i) = Zi(:);
			% end
			% %calculate transformation matrix for X: R^m -> M^n
			% mX = zeros(n2,m);
			% for i=1:m
				% %ei = zeros(m,1);
				% %ei(i) = 1;
				% %Xi = sssdpx(n, ei, sMind, X, S, E, iindf);
				% %Xi = sssdpxi(n, sMind, X, S(i,:), E(i,:), iindf);
				% Xi = -A{mp(i)}*X;
				% mX(:,i) = Xi(:);
			% end
			% J = [mZ mX];
			% deltas = -J\Fmu(:);
			% lsqriter = 0;
		% end
	end
	%sdpj is used to calculate J(v,y) given A and Q and also J*(M)
	tlsqr = toc(stlsqr);
   %  disp(['solvlss-lsqr: itns rnorm, xnorm cputime:  ',num2str(itn),' ',...
   %      num2str(rnorm),' ',num2str(rnorm),' ',num2str(tlsqr)]);
   %deltax=deltas(1:tnm1);deltay=deltas(tnm1+1:length(deltas));
   %  disp(['solvlss-lsqr: atol btol :  ',num2str(atol),' ',num2str(btol)]);
   
   % if precond==2,
     % deltas(1:tnm1)=cholmat11\deltas(1:tnm1);
     % deltas(1:tnm1)=deltas(invpcols);
     % deltas(end-n+1:end)=deltas(end-n+1:end)./dataxy(end-n+1:end);
   % elseif precond==1,
   % deltas=deltas./dataxy;   % old diagonal preconditioner
   % elseif precond==0,   % no preconditioner
   % end
	%calculate the directions
	deltav = deltas(1:tn-m);
	deltay = deltas(tn-m+1:end);
	%disp(['Direction norms ' num2str(norm(deltav)) ' ' num2str(norm(deltay)) ', RHS norm ' num2str(norm(Fmu(:)))]);
	if precond == 1
		deltav = deltav./Tv; %scale dv back according to Tv;
		deltay = deltay./Ty; %scale dy back according to Ty;
	else
		if precond >= 2
			deltav = Tv\deltav;
			deltay = Ty\deltay;
		end
	end
	
	%disp(['Direction norms ' num2str(norm(deltav)) ' ' num2str(norm(deltay))]);
	%Dy = [Dy deltay];
	
    deltax = Q*deltav - rx;
	deltaX = sssMat(deltax,n,sMind);
    deltaz = -mA'*deltay - rd;
	deltaZ = sssMat(deltaz,n,sMind);
	%total time to find the directions
	dirtimec = cputime - tfordirc;
%%%%%%%%%%%%%%%%%%%%%%%
%%% start crossover - needs lots of polishing
	if ((norm(rx)+norm(rd)+gapp)/(abs(optval)+1) < crosstoler)  & ~crossover ,
		crossover = 1;
		step = 1; %step stays 1 after crossover starts
		sigmaa = 0; %stays 0 after crossover starts
		itersinfo(2) = noiter;
		disp(['*********CROSSover started at: ',num2str((norm(rx)+norm(rd)+gapp)/(abs(optval)+1))]);
   %% disp('noiter   -log10(relgapp)    step        sigmaa         optval         lsqrtime      cputimeprec    lsqr.iter  nnz-X  nnz-chol');
     %disp('noiter -log10(relnormZX)     step       sigmaa     optval          cputime.dir,');
	end
	
	if ~crossover,
%     disp(['starting the back tracking  ']);
		%vp = v + step*deltav;
		%yp = y + step*deltay;
        %xp = x + step*deltax;
        %zp = z + step*deltaz;
		Xp = X + step*deltaX;
		Zp = Z + step*deltaZ;
		%check the semidefiniteness of solution
        [Rx,choltrueX] = chol(Xp);
        [Rz,choltrueZ] = chol(Zp);	
		while choltrueX + choltrueZ > 0,
%%testing ?????
        %step=.9999*step;   % small backtrack to get bndry
                       % change to eig approach?????
%%testing ?????
			step = .8*step;  % original backtracking
			%vp = v + step*deltav;
			%yp = y + step*deltay;
            %xp = x + step*deltax;
            %zp = z + step*deltaz;
			Xp = X + step*deltaX;
			Zp = Z + step*deltaZ;
			[Rx,choltrueX] = chol(Xp);
			[Rz,choltrueZ] = chol(Zp);
		end
		% if choltrueX + choltrueZ > 0
			% stepL = 0;
			% stepH = step;
			% step = (stepL + stepH)/2;
			% while ((stepH-stepL)/2>1e-16)
				% Xp = X + step*deltaX;
				% Zp = Z + step*deltaZ;
				% [Rx,choltrueX] = chol(Xp);
				% [Rz,choltrueZ] = chol(Zp);
				% if choltrueX + choltrueZ > 0
					% stepH = step;
				% else
					% stepL = step;
				% end
				% step = (stepL + stepH)/2;
			% end
			% step = stepL;
			% if (step==0)
				% disp(['to USER from srssdp.m: ERROR; current step size is too small']);
				% return;
			% end
		% end
  %  step=.999997*step;     % original backtrack strategy
%%testing ?????
     %step=min(1.07,1.32-min(tau1,abs(tau1*gapp)   )    )*step;
%%testing ?????
		step = min(.99997,1-min(tau1,abs(tau1*gapp)))*step; %proper way
  %  step=.97*step;
	end
%%% end crossover modif - below for sigmaa as well.
	%update all the variables
	v = v + step*deltav;
	y = y + step*deltay;
    x = x + step*deltax;
    z = z + step*deltaz;
	X = X + step*deltaX;
	Z = Z + step*deltaZ;
	%keyboard;
	% v = vp;
	% y = yp;
    % x = xp;
    % z = zp;
	% X = Xp;
	% Z = Zp;
    rx = x - x_0 - Q*v;
    rd = z - c + mA'*y;
	%update the optimality matrix ZX
    ZX = Z*X;
	ZXnorm = norm(ZX,'fro'); %%for crossover check
	gapp = trace(ZX);
	muu = gapp/(2*n);
	%calculate the current objective
    optval = c'*x;
%  reler is relative error: duality-gapp/n1/(optval+1)
	if crossover,
		muu = ZXnorm/(abs(optval)+1);  % need positive for lsqr tolerance atol,btol

		fprintf(' %3g      %6.2e         %6.2e    %8.2e     %5.4e    %2.4e     %2.4e   %5g   %5g   %5g\n', ...
			noiter,-log10(muu), ...
			1,0,optval,tlsqr,timeprec,lsqriter,nnz(X),nnzchol);
	else
		fprintf(' %3g      %6.2e         %6.2e    %8.2e     %5.4e    %2.4e     %2.4e   %5g   %5g   %5g\n', ...
		    noiter,-log10(full(gapp)/(optval+1)), ...
			step,sigmaa,optval,tlsqr,timeprec,lsqriter,nnz(X),nnzchol);


%disp([sprintf(' %3g      %6.2e         %4g    %10g     %5.4e      %2.4e  %2.4e %2.4e', ...
%       noiter,-log10(ZXnorm/(abs(optval)+1)), ...
%        1,0,optval,dirtimec)]);
%  else
%disp([sprintf(' %3g      %6.2e       %7.4e    %2.1e     %2.4e      %2.4e  %2.4e %2.4e   %2.4', ...
%       noiter,-log10(gapp/(optval+1)), ...
%        step,sigmaa,optval,dirtimec,min(min(eig(Z),eig(X))))]);
% update muu and sigmaa depending on stepsize
		if step < .001,
			muu = muu * 2.1;
			sigmaa = 1;
			step = .9;  %muu increased
		elseif step < .01, 
			muu = muu * 2;
			sigmaa = 1;
			step = 1;  % muu unchanged
		elseif step < .25, 
			muu = muu * 1.5; 
			sigmaa = 1;
			step = 1; 
%%   testing for influence of condition number on lsqr
%%% steplength restricted to one when using lsqr - seems to help
%%% the next 3 lines are changed to get steplength <=1
%% This appears to be TRUE in case of crossover soon after,
%%   i.e. stay well centered till crossover.
%     elseif step < .95, muu = muu * 1.2; sigmaa=1-.3*min(step,1);step=1.5;
%     elseif step < 1.05, muu = muu *.5; sigmaa=1-.3*min(step,1); step=2.0;
%     else  sigmaa=1-.3*step;step=2;muu=muu*.8;
		elseif step < .95, 
			muu = muu * 1.2; 
			sigmaa = 1-.3*min(step,1);
			step = 1.0;
		elseif step < 1.05, 
			muu = muu *.5; 
			sigmaa = 1-.3*min(step,1); 
			step = 1.0;
		else  
			sigmaa = 1-.3*step;
			step = 1;
			muu = muu*.8;
		end
	end
%%% end crossover modif - see above as well
	if wantmeritfn, 
		meritfn = [meritfn ZXnorm];
	end
%%% check to see whether the lsrqr is stagnant
	if (precond==2 & ~crossover & lsqriter==itnlim) %& gapp/(abs(optval)+1) < 1e-10)
		disp(['to USER from srssdp.m: Possible maximum tolerance reached - alogrithms ends']);
		stagnant = 1;
		break;
	end
end          % end of main while loop
%save dtmp.mat Dy
%calculate the final solution
%X = sssMat(x,n,sMind);
%Z = sssMat(z,n,sMind);
%y = (T')\y; %change back to the original basis of the row space of A
y = y(imp); %change back to the original row setting
%  if wantmeritfn,
%    plot(-log10(meritfn/(abs(optval)+1)));
%    title('plot of -log10(norm(ZX)/(abs(optval)+1))')
%  end
if noiter < iterbnd,
	disp('to USER from srssdp.m: Tolerance achieved - algorithm ends')
else
	disp(['to USER from srssdp.m: ERROR;  exceeded iteration bound:  ',...
		num2str(noiter)]);
end
disp(['for crossover; final rel. norm(ZX) is: ',num2str(ZXnorm/(abs(optval)+1))]); 
disp(['for crossover; rel. min of eigs is: ',...
	num2str( min(  min(eig(X)),min(eig(Z)) )/(abs(optval)+1)  ) ]); 
itersinfo(1) = noiter;
%DIMACS errors
pfeas = b;
for i=1:m
	pfeas(i) = pfeas(i) - sum(sum(A{i}.*X));
end
err1 = norm(pfeas)/(1+max(abs(b)));
err2 = max([0 -min(eig(X))/(1+max(abs(b)))]);
dfeas = Z - C;
for i=1:m
	dfeas = dfeas + y(i)*A{i};
end
err3 = norm(dfeas,'fro')/(1+max(abs(C(:))));
err4 = max([0 -min(eig(Z))/(1+max(abs(C(:))))]);
err5 = (sum(sum(C.*X))-b'*y)/(1+abs(sum(sum(C.*X)))+abs(b'*y));
err6 = sum(sum(Z.*X))/(1+abs(sum(sum(C.*X)))+abs(b'*y));
disp(['DIMACS error measures: ' num2str(err1) ' ' num2str(err2) ' ' num2str(err3) ' ' num2str(err4) ' ' num2str(err5) ' ' num2str(err6)]);
