function [X,S]=TriangulationGsedumi(P,E,Eo)
    
    % Test a framework with configuration/position P, graph 
    % specified by edge set E, and an objective function specified by by 
    % edges in Eo. The distance of each edge in E will be set as a 
    % constraint; and the sum of edge distances in Eo will be maximized.
    % E could represent a triangulation graph, and Eo represents edge set 
    % of all virtual edges as in Shamsi, Taheri, Zhu and Ye 2011.
    %
    % coded on 12/29/2011.
    %
    % 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 which Euclidean distance is set as a constraint.
    % Eo: 2 x |Eo| matrix where each column specifies two point-idexes of 
    %    an edge to be maximized in the sum of all edges in Eo.
    %
    % Output:
    % X: d' x n configuration matrix satisfying the edge distance 
    % constraint, where the first point in X is the origin.
    % Thus, the second column of X is the positiion of the second point in
    % P, and so on.
    %
    % d'=d if and only if the framework is uniquely localizable with the
    % added objective function.
    %
    % S: a max-rank PSD (dual) stress matrix approximately orthogonal to X
    %
    % Solver: sedumi 1.1 by Jos Stern
    %
    % Check the dimension of the configuration and the number of edges
    %
    [d,n]=size(P);
    n=n-1;
    [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));
        if (i ~= 1), a(i-1)=1;end;
        if (j ~= 1), a(j-1)=-1;end;
        A=[A vec(a*a')];            
        dd= P(:,i)-P(:,j);
        rr= dd'*dd;        
        b=[b; rr];
    end;
    %
    % Set up the SDP objective using edges specified in Eo
    %
    [m,objn]=size(Eo);
    C=vec(sparse(n,n));
    for k =1:objn
        i = Eo(1,k);
        j = Eo(2,k);
        a=sparse(zeros(n,1));
        if (i ~= 1), a(i-1)=1;end;
        if (j ~= 1), a(j-1)=-1;end;
        C=C-vec(a*a');            
    end;
    % 
    % 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);
    %
    % Prepare the output solutions
    %
    Y = mat(Y);
    X = chol(Y);
    r = rank(Y,sqrt(toler));
    X =[zeros(r,1) X(1:r,:)];
    S = mat(C);
    for k =1:mtotal
        S=S-y(k)*mat(A(:,k));
    end;
    %
    % 
    % end of the function