%%********************************************************************
%% infeaspt: generate an initial point for sdp.m
%%
%%  [X0,y0,Z0] = infeaspt(blk,At,C,b,options,scalefac);
%%
%%  options = 1  if want X0,Z0 to be scaled identity matrices
%%          = 2  if want X0,Z0 to be scalefac*(identity matrices).
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%********************************************************************

function [X0,y0,Z0] = infeaspt(blk,At,C,b,options,scalefac)
%%
if (nargin < 5); options = 1; end;
if (options == 1); scalefac = []; end;
if (options == 2) && (nargin < 6); scalefac = 1000; end;
if (scalefac <= 0); error('scalefac must a positive number'); end;
%%
if ~iscell(At); At = {At}; end;
if ~iscell(C);  C = {C}; end;
m = length(b);
if all(size(At) == [size(blk,1) m]);
    convertyes = zeros(size(blk,1),1);
    for p = 1:size(blk,1)
        if strcmp(blk{p,1},'s') && all(size(At{p,1}) == sum(blk{p,2}))
            convertyes(p) = 1;
        end
    end
    if any(convertyes)
        At = svec(blk,At,ones(size(blk,1),1));
    end
end;
%%
%%[blk,At,C,b] = validate(blk,At,C,b);
%%
X0 = cell(size(C)); Z0 = cell(size(C));
m = length(b);
for p = 1:size(blk,1);
    pblk = blk(p,:);
    blktmp = pblk{2};
    n = length(C{p});
    y0 = zeros(m,1);
    b2 = 1 + abs(b');
    if (options == 1);
        if strcmp(pblk{1},'s');
            normAni = [];
            X0{p} = sparse(n,n); Z0{p} = sparse(n,n);
            ss = [0, cumsum(blktmp)];
            tt = [0, cumsum(blktmp.*(blktmp+1)/2)];
            for i = 1:length(pblk{2})
                if ~isempty(At{p,1})
                    pos = tt(i)+1 : tt(i+1);
                    Ai = At{p,1}(pos,:);
                    normAni = 1+sqrt(sum(Ai.*Ai));
                end
                if (length(At(p,:)) >= 2) %% for low rank constraints
                    dd = At{p,3};
                    qq = [0, cumsum(pblk{3})]; normtmp = ones(1,length(pblk{3}));
                    idxD = [0; find(diff(dd(:,1))); size(dd,1)];
                    for k = 1:length(pblk{3})
                        idx = qq(k)+1 : qq(k+1);
                        idx2 = idxD(k)+1: idxD(k+1);
                        Ak = At{p,2}(:,idx);
                        ii = dd(idx2,2)-qq(k); %% undo cumulative indexing
                        jj = dd(idx2,3)-qq(k);
                        len = pblk{3}(k);
                        Dk = spconvert([ii,jj,dd(idx2,4); len,len,0]);
                        tmp = Ak'*Ak*Dk;
                        normtmp(1,k) = 1+sqrt(sum(sum(tmp.*tmp')));
                    end
                    normAni = [normAni, normtmp]; %#ok
                end
                pos = ss(i)+1 : ss(i+1);  ni = length(pos);
                tmp = C{p}(pos,pos);
                normCni = 1+sqrt(sum(sum(tmp.*tmp)));
                const  = 10; %%--- old: const = 1;
                constX = max([const,sqrt(ni),ni*(b2./normAni)]);
                constZ = max([const,sqrt(ni),normAni,normCni]);
                X0{p}(pos,pos) = constX*spdiags(1+1e-10*randmat(ni,1,0,'u'),0,ni,ni);
                Z0{p}(pos,pos) = constZ*spdiags(1+1e-10*randmat(ni,1,0,'u'),0,ni,ni);
            end
        elseif strcmp(pblk{1},'q');
            s = 1+[0, cumsum(blktmp)];
            len = length(blktmp);
            normC = 1+norm(C{p});
            normA = 1+sqrt(sum(At{p,1}.*At{p,1}));
            idenqX = zeros(sum(blktmp),1);
            idenqZ = zeros(sum(blktmp),1);
            idenqX(s(1:len)) = max([1,b2./normA])*sqrt(blktmp') ;
            idenqZ(s(1:len)) = max([sqrt(blktmp); max([normA,normC])*ones(1,len)])';
            idenqX(s(1:len)) = idenqX(s(1:len)).*(1+1e-10*randmat(len,1,0,'u'));
            idenqZ(s(1:len)) = idenqZ(s(1:len)).*(1+1e-10*randmat(len,1,0,'u'));
            X0{p} = idenqX;
            Z0{p} = idenqZ;
        elseif strcmp(pblk{1},'l');
            normC = 1+norm(C{p});
            normA = 1+sqrt(sum(At{p,1}.*At{p,1}));
            const = 10; %%--- old: const =1;
            constX = max([const,sqrt(n),sqrt(n)*b2./normA]);
            constZ = max([const,sqrt(n),normA,normC]);
            X0{p} = constX*(1+1e-10*randmat(n,1,0,'u'));
            Z0{p} = constZ*(1+1e-10*randmat(n,1,0,'u'));
        elseif strcmp(pblk{1},'u');
            X0{p} = sparse(n,1);
            Z0{p} = sparse(n,1);
        else
            error(' blk: some fields not specified correctly');
        end;
    elseif (options == 2);
        if strcmp(pblk{1},'s');
            n = sum(blktmp);
            X0{p} = scalefac*spdiags(1+1e-10*randmat(n,1,0,'u'),0,n,n);
            Z0{p} = scalefac*spdiags(1+1e-10*randmat(n,1,0,'u'),0,n,n);
        elseif strcmp(pblk{1},'q');
            s = 1+[0, cumsum(blktmp)];
            len = length(blktmp);
            idenq = zeros(sum(blktmp),1);
            idenq(s(1:len)) = 1+1e-10*randmat(len,1,0,'u');
            X0{p} = scalefac*idenq;
            Z0{p} = scalefac*idenq;
        elseif strcmp(pblk{1},'l');
            X0{p} = scalefac*(1+1e-10*randmat(n,1,0,'u'));
            Z0{p} = scalefac*(1+1e-10*randmat(n,1,0,'u'));
        elseif strcmp(pblk{1},'u');
            X0{p} = sparse(n,1);
            Z0{p} = sparse(n,1);
        else
            error(' blk: some fields not specified correctly');
        end
    end
end
%%********************************************************************
