function [W, S, obj] = SymTriOrthNMF(A, r, max_iter, check_step, W0, S0)
% the symmetric orthogonal 3-factor NMF
% He Zhang, Sep. 2012
% Ref: C. Ding et al.: Orthogonal nonnegative matrix t-factorizations for clustering, 
% SIGKDD, ACM, 2006

n = size(A,1);
rand('seed', 0);
if ~exist('W0', 'var') || isempty(W0)
    W0 = rand(n,r);
end    
if ~exist('S0', 'var') || isempty(S0)
    S0 = W0'*A*W0;
end
W = W0;
S = S0;
for iter=1:max_iter
    AWS = A*W*S;
    W = W .* AWS ./ (W*(W'*AWS)+eps);
    W = max(W,eps);
    WW = W'*W;
    S = S .* (W'*A*W) ./ (WW*S*WW+eps);
    S = max(S,eps);
    if mod(iter, check_step)==0
        obj = norm(A,'fro')^2 -2*trace(W'*A*W*S) + norm(W'*W*S, 'fro')^2;
        fprintf('SymTriOrthNMF: iter=%d, obj=%6.5f\n', iter, obj);
    end
end

if nargout>2
    obj = norm(A,'fro')^2 -2*trace(W'*A*W*S) + norm(W'*W*S, 'fro')^2;
end
