function dJ =...
    sssdpj(mode, n2, tn, dx, svind, sMind, X, Z, S, E, indf, iindf, precond, Tv, Ty, indsE)
%function dJ = ssdpj(mode, n2, tn, dx, svind, sMind, X, Z, S, E, indf, iindf, precond, Tv, Ty)
%sdpj: calculate J and J* depending on mode
%n2, tn: size of J
%svind: vector-matrix transformation (upper triangle)
%sMind: matrix-vector transformation
%X, Z: current data (X,Z)
%A = [S E], Q = [-S\E; I]
%indf: permutation vector with firs m values defining the well-conditioned block in A
%iindf: inverse permutation to get back to the original A
%precond: indicate which preconditioner is used
% 	precond = 0: no preconditioning
% 	precond = 1: diagonal preconditioner
%   precond >= 2: block diagonal preconditioner
%(Tv, Ty): parameters for preconditioner
% 	precond = 0: do not need  (Tv, Ty)
% 	precond = 1: (Tv, Ty) are diagonal elements of the preconditioner
%   precond >= 2: (Tv, Ty) are block diagonal matrices of the preconditioner
%mode = 1: 
%	dJ = J(dx)
%mode = 2:
% 	dJ = J*(dx)
%global svind; %assume X has a fixed size n by n
%x = xz(:,1);
%z = xz(:,2);
n = 2*tn - n2; %tn=n(n+1)/2 and n2=n^2
%Z = ssMat(z,n,sMind);
%X = ssMat(x,n,sMind);
m = size(S,1); %number of constraints
M = zeros(n,n); %matrix used to calculate Jacobian and its adjoint
T = zeros(n,n); %transpose matrix used in calculation
%calculate y = J(x)
	%indsE=find(sum(E)~=0);       % test by henry mar 12
	Ebar=E;    % test by henry mar 12
	%E=zeros(m,tn-m);   % testing
	%E(:,indsE)=Ebar;   % testing
if (mode==1)
	dv = dx(1:tn-m);
	dy = dx(tn-m+1:tn);
	if precond == 1
		dv = dv./Tv; %scale dv according to Tv
		dy = dy./Ty; %scale dy according to Ty
	else
		if precond >= 2
			dv = Tv\dv; 
			dy = Ty\dy;
		end
	end
	dX = [-S\(Ebar*dv(indsE)); dv];    % test by henry mar 12
	%dXbar = [-S\(Ebar*dv(indsE)); dv];    % test by henry mar 12
	%dX = [-S\(E*dv); dv];
	%if norm(dX-dXbar)>1e-15,   %test by henry
	%	'error'   %test by henry
	%	keyboard   %test by henry
	%end   %test by henry
	dX = dX(iindf); %shift the vector back to the original format
	nX = sssMat(dX,n,sMind);
	%dZ = [S'*dy; E'*dy];
	temp=zeros(tn-m,1);    % test by henry mar 12
	temp(indsE)= dy'*Ebar;    % test by henry mar 12
	dZ = [(dy'*S)'; temp];    % test by henry mar 12
	%dZbar = [(dy'*S)'; temp];    % test by henry mar 12
	%dZ = [(dy'*S)'; (dy'*E)'];
	dZ = dZ(iindf);
	nZ = sssMat(dZ,n,sMind);
	%M = Z*nX - nZ*X;
	M = Z*nX;
	dJ = M(:);
	M = nZ*X;
	dJ = dJ - M(:); %vectorize the matrix M
end
%calculate y = J^*(x)
if (mode==2)
	%M = zeros(n,n);
	M(:) = dx;
	ZM = Z*M;
	T = ZM';
	ZM = ZM + T;
	ZM = ZM/2;
	zm = sssvec(ZM,svind);
	zm = zm(indf); %shift the vector to the new format
	%XM = X*M';
	XM = M*X;
	T = XM';
	XM = XM + T;
	XM = XM/2;
	xm = sssvec(XM,svind);
	xm = xm(indf);
	if precond == 0
		%dJ = [-E'*((S')\zm(1:m))+zm(m+1:end); -(S*xm(1:m)+E*xm(m+1:end))];
		dJ = [-(((zm(1:m))'/S)*E)'+zm(m+1:end); -(S*xm(1:m)+E*xm(m+1:end))];
	else
		if precond == 1
			%scale results according to (Tv, Ty)
			%dJ = [(-E'*((S')\zm(1:m))+zm(m+1:end))./Tv; -(S*xm(1:m)+E*xm(m+1:end))./Ty]; 
			tempz=zeros(tn-m,1);   % test by henry
			tempz(indsE)= -(((zm(1:m))'/S)*Ebar); %test by henry
			tempx=xm(m+1:end); %testing
			tempv=tempx(indsE); %testing
			tempvv=Ebar*tempv; %testing
			dJ = [(tempz+zm(m+1:end))./Tv; -(S*xm(1:m)+tempvv)./Ty]; %testing 
			%dJbar = [(tempz+zm(m+1:end))./Tv; -(S*xm(1:m)+tempvv)./Ty]; %testing 
			%dJ = [(-(((zm(1:m))'/S)*E)'+zm(m+1:end))./Tv; -(S*xm(1:m)+E*xm(m+1:end))./Ty]; 
			%if norm(dJbar-dJ)>1e-14, %testing
			%	'errordJ' %testing
			%	keyboard %testing
			%end %testing
		else
			%%% still to do for new Ebar????
			if precond >= 2
				%dJ = [((Tv')\(-E'*((S')\zm(1:m))+zm(m+1:end))); -((Ty')\(S*xm(1:m)+E*xm(m+1:end)))];
				dJ = [((-((zm(1:m))'/S)*E+(zm(m+1:end))')/Tv)'; -((S*xm(1:m)+E*xm(m+1:end))'/Ty)'];
			end
		end
	end
end

