function w = wcond(A,toler,kind)
%function w = wcond(A,toler,kind)
%wcond: calculate the w-condition number of A = X'X pos def, 
%       sqrt(trace(X'X)/n/(det(X'X))^{1/n})
%INPUT:
%  A,X    n x n pos. def.; X is m x n full column rank
%  toler  tolerance to find w if det(A) jumps from 0 to Inf
%  kind   =1 if A pos. def. is the input, 
%		  =2 if X rectangular is input
%default is assuming pos. def. A is input
if nargin==3,
	if kind==1,
		disp(['starting wcond assuming given input is A pos. def.'])
	elseif kind==2,
		disp(['starting wcond assuming given input is X rect.'])
		[m,n] = size(A);
		if (m>=n)
			A = A'*A;
		else
			A = A*A';
		end
	else
		disp('ERROR: input kind is not 1 or 2')
		w = nan;
		return
	end
elseif nargin==2
    disp(['starting wcond assuming given input is A pos. def.'])
else
    disp(['starting wcond assuming given input is A pos. def.'])
	toler=log(n)*1e-10;  % default tolerance
	disp(['default tolerance for omega is  ',num2str(toler)]);
end
%check the input
[m,n]=size(A);
if m ~= n, 
	disp('ERROR: A is not square. Please check input.');
	w = nan;
	return
end
if normest(A-A')>1e-10,
	disp('ERROR: A is not symmetric. Please check input.');
	w = nan;
	return
end
%limit number of iterations to check w = Inf
maxlim = ceil((log2(realmax)-log2(realmin))); 
w = trace(A)/n;
%lu decomposition of A
[~,u] = lu(full(A));
%use the formulation det(A)=abs(prod(diag(u)))
A = diag(u);
dA = prod(A);
%if det(A) is finite, calculate w directly
if (isfinite(dA) & isfinite(log(dA)))
	w = w/abs(dA)^(1/n);
	w = sqrt(w);
	return;
end
%calculate det(A)^{1/n} using binary search
%find upper and lower bound for scalar a
if isinf(dA)
	aL = 0;
	aH = 1;
else
%if det(A) = 0, calculate aL: det(aL*A) = 0 while det(2*aL*A) ~= 0
	aL = 1;
	aH = 2;
	dA = prod(aH*A);
	index = 1;
	while (isfinite(dA) & isinf(log(dA)) & index <= maxlim)
		aL = aH;
		aH = 2*aH;
		dA = prod(aH*A);
		index = index + 1;
	end
	%if det(aH*A) = 0
	if (isfinite(dA) & isinf(log(dA)))
		disp('ERROR: A is not pos. def. Please check input.');
		w = Inf;
		return;
	end
end
%binary search to find the scalar a such that det(aA) is finite
a = aH;
%if det(aH*A) is not finite
while ((isinf(dA) | isinf(log(dA))) & (aH-aL)/2>toler)
	a = (aL + aH)/2;
	dA = prod(a*A);
	if (isfinite(dA) & isfinite(log(dA)))
		w = a*w/abs(dA)^(1/n);
		w = sqrt(w);
		return;
	end
	%update aL and aH
	if (isinf(dA))
		aH = a;
	else
		aL = a;
	end
end
%if det(a*A) and log(det(a*A)) are finite
if ((aH-aL)/2 > toler)
	w = a*w/abs(dA)^(1/n);
	w = sqrt(w);
else
%1/det(A)^{1/n} is between aL and aH
	w = w*((aH+aL)/2);
	w = sqrt(w);
end