function [Y,S]=URFrameworkTest(P,E)
    
    % Test if a bar-framework with configuration/position P and graph 
    % specified by edge set E is universal rigid (UR); 
    % coded on 2/6/2011. Updated 2/10/2011 to output the configuration
    % matrix
    
    % Input:
    % P: d x n matrix where jth column is the position of jth point in R^d
    % E: 2 x |E| matrix where each column specifies two point-idexes of 
    %    an edge 
    
    % Output:
    % Y: a max-rank configuration matrix satisfying the edge distance 
    % constraint approximately
    %
    % rank(Y)=d+1 if and only if the framework is UR, if the graph contains
    % a (d+1) clique and the (d+1) points are in general position; see
    % So and Ye 2005.
    %
    % S: a max-rank PSD (dual) stress matrix approximately orthogonal to Y
    %
    % rank(S)=n-d-1 implies that the framework is UR, if all points are in 
    % general position; see Alfakih and Ye 2010, and Connelly 1999 and
    % Gortler/Thurston 2009 for points in generic position.
    %
    
    % Check the dimension of the configuration and the number of edges
    %
    [d,n]=size(P);
    [m,mtotal]=size(E);
    toler = 1.e-8;
    %
    % Compute edge distances specified by E and set up the SDP inputs
    %
    b=[];
    A=[];
    for k =1:mtotal
        i = E(1,k);
        j = E(2,k);
        a=sparse(zeros(n,1));
        a(i)=1;
        a(j)=-1;
        A=[A vec(a*a')];            
        dd= P(:,i)-P(:,j);
        rr= dd'*dd;        
        b=[b; rr];
    end;
    C=vec(sparse(n,n));
    % 
    % Set up algorithm parameter to call the SDP solver sedumi
    %
    K.s=[n];
    pars.eps=toler;
    pars.stepdif=0;
    [Y,y,info] = sedumi(A,b,C,K,pars);
    Y = mat(Y);
    S = zeros(n,n);
    for k =1:mtotal
        S=S-y(k)*mat(A(:,k));
    end;
    %
    % Prepare the outputs
    %
    YR = chol(Y);
    rankY = rank(Y,sqrt(toler))
    r = rankY;
    Y =YR(1:r,:);
    r = - min([min(eig(S)),0]);
    S = (r+toler^2)*eye(n,n)+S;
    SR = chol(S);
    rankS = rank(S,sqrt(toler))
    r = rankS;
    SR=SR(1:r,:);
    S = SR'*SR;
    if (r>0) & (r==n-d-1),
        YRP = SR\(-SR*Y');
        Y=Y+YRP';
    end
    % end of the function