ncvar.m

optimization toolbox

        % Find zero-variables
        constants = find(sys.lmi_variables==0);
        if ~isempty(constants);
            sys.lmi_variables(constants)=[];
            sys.basis(:,1) = sys.basis(:,1) + sum(sys.basis(:,1+constants),2);
            sys.basis(:,1+constants)=[];
        end
        if isempty(sys.lmi_variables)
            sys = full(reshape(sys.basis(:,1),sys.dim(1),sys.dim(2)));
        else
            sys = class(sys,'sdpvar');
        end
        return
    case 6 % Fast version for internal use
        sys.basis = varargin{5};
        sys.lmi_variables=varargin{4};
        sys.dim(1) = varargin{1};
        sys.dim(2) = varargin{2};
        sys.typeflag = varargin{6};
        sys.savedata = [];
        sys.extra = [];
        sys.extra.expanded = [];        
        sys.conicinfo = 0;
        % Find zero-variables
        constants = find(sys.lmi_variables==0);
        if ~isempty(constants);
            sys.lmi_variables(constants)=[];
            sys.basis(:,1) = sys.basis(:,1) + sum(sys.basis(:,1+constants),2);
            sys.basis(:,1+constants)=[];
        end
        if isempty(sys.lmi_variables)
            sys = full(reshape(sys.basis(:,1),sys.dim(1),sys.dim(2)));
        else
            sys = class(sys,'sdpvar');
        end
        return
    case 7 % Fast version for internal use
        sys.basis = varargin{5};
        sys.lmi_variables=varargin{4};
        sys.dim(1) = varargin{1};
        sys.dim(2) = varargin{2};
        sys.typeflag = varargin{6};
        sys.savedata = [];
        sys.extra = varargin{7};
        sys.extra.expanded = [];        
        sys.conicinfo = varargin{7};
        % Find zero-variables
        constants = find(sys.lmi_variables==0);
        if ~isempty(constants);
            sys.lmi_variables(constants)=[];
            sys.basis(:,1) = sys.basis(:,1) + sum(sys.basis(:,1+constants),2);
            sys.basis(:,1+constants)=[];
        end
        if isempty(sys.lmi_variables)
            sys = full(reshape(sys.basis(:,1),sys.dim(1),sys.dim(2)));
        else
            sys = class(sys,'sdpvar');
        end
        return

    otherwise
        error('Wrong number of arguments in sdpvar creation');
end

if isempty(n) | isempty(m)
    error('Size must be integer valued')
end;
if ~((abs((n-ceil(n)))+ abs((m-ceil(m))))==0)
    error('Size must be integer valued')
end

nonCommutingTable = yalmip('nonCommutingTable');
[monomtable,variabletype] = yalmip('monomtable');
lmi_variables = (1:nvar)+max(size(nonCommutingTable,1),size(monomtable,1));

for blk = 1:length(n)
    switch matrix_type

        case 'full'
            basis{blk} = [spalloc(n(blk)*m(blk),1,0) speye(n(blk)*m(blk))];%speye(nvar)];

        case 'full complex'
            basis = [spalloc(n*m,1,0) speye(nvar/2) speye(nvar/2)*sqrt(-1)];

        case 'symm'
            if 0
                basis = spalloc(n^2,1+nvar,n^2);
                l = 2;
                an_empty = spalloc(n,n,2);
                for i=1:n
                    temp = an_empty;
                    temp(i,i)=1;
                    basis(:,l)=temp(:);
                    l = l+1;
                    for j=i+1:n,
                        temp = an_empty;
                        temp(i,j)=1;
                        temp(j,i)=1;
                        basis(:,l)=temp(:);
                        l = l+1;
                    end
                end
            else
                % Hrm...fast but completely f*d up

                Y = reshape(1:n(blk)^2,n(blk),n(blk));
                Y = tril(Y);
                Y = (Y+Y')-diag(sparse(diag(Y)));
                [uu,oo,pp] = unique(Y(:));
                if 1
                    basis{blk} = sparse(1:n(blk)^2,pp+1,1);
                else                    
                    basis{blk} = lazybasis(n^2,1+(n*(n+1)/2),1:n(blk)^2,pp+1,ones(n(blk)^2,1));
                end
              
            end
            
        case 'symm complex'
            basis = spalloc(n^2,1+nvar,2);
            l = 2;
            an_empty = spalloc(n,n,2);
            for i=1:n
                temp = an_empty;
                temp(i,i)=1;
                basis(:,l)=temp(:);
                l = l+1;
                for j=i+1:n,
                    temp = an_empty;
                    temp(i,j)=1;
                    temp(j,i)=1;
                    basis(:,l)=temp(:);
                    l = l+1;
                end
            end
            for i=1:n
                temp = an_empty;
                temp(i,i)=sqrt(-1);
                basis(:,l)=temp(:);
                l = l+1;
                for j=i+1:n,
                    temp = an_empty;
                    temp(i,j)=sqrt(-1);
                    temp(j,i)=sqrt(-1);
                    basis(:,l)=temp(:);
                    l = l+1;
                end
            end

        case 'herm complex'
            basis = spalloc(n^2,1+nvar,2);
            l = 2;
            an_empty = spalloc(n,n,2);
            for i=1:n
                temp = an_empty;
                temp(i,i)=1;
                basis(:,l)=temp(:);
                l = l+1;
                for j=i+1:n,
                    temp = an_empty;
                    temp(i,j)=1;
                    temp(j,i)=1;
                    basis(:,l)=temp(:);
                    l = l+1;
                end
            end
            for i=1:n
                for j=i+1:n,
                    temp = an_empty;
                    temp(i,j)=sqrt(-1);
                    temp(j,i)=-sqrt(-1);
                    basis(:,l)=temp(:);
                    l = l+1;
                end
            end

        case 'skew'
            basis = spalloc(n^2,1+nvar,2);
            l = 2;
            an_empty = spalloc(n,n,2);
            for i=1:n
                for j=i+1:n,
                    temp = an_empty;
                    temp(i,j)=1;
                    temp(j,i)=-1;
                    basis(:,l)=temp(:);
                    l = l+1;
                end
            end

        case 'skew complex'
            basis = spalloc(n^2,1+nvar,2);
            l = 2;
            an_empty = spalloc(n,n,2);
            for i=1:n
                for j=i+1:n,
                    temp = an_empty;
                    temp(i,j)=1;
                    temp(j,i)=-1;
                    basis(:,l)=temp(:);
                    l = l+1;
                end
            end
            for i=1:n
                for j=i+1:n,
                    temp = an_empty;
                    temp(i,j)=sqrt(-1);
                    temp(j,i)=-sqrt(-1);
                    basis(:,l)=temp(:);
                    l = l+1;
                end
            end

        case 'toep'
            basis = spalloc(n^2,1+nvar,2);
            an_empty = spalloc(n,1,1);
            for i=1:n,
                v = an_empty;
                v(i)=1;
                temp = sparse(toeplitz(v));
                basis(:,i+1) = temp(:);
            end

            % Notice, complex Toeplitz not Hermitian
        case 'toep complex'
            basis = spalloc(n^2,1+nvar,2);
            an_empty = spalloc(n,1,1);
            for i=1:n,
                v = an_empty;
                v(i)=1;
                temp = sparse(toeplitz(v));
                basis(:,i+1) = temp(:);
            end
            for i=1:n,
                v = an_empty;
                v(i)=sqrt(-1);
                temp = sparse(toeplitz(v));
                basis(:,n+i+1) = temp(:);
            end

        case 'hank'
            basis = spalloc(n^2,1+nvar,2);
            an_empty = spalloc(n,1,1);
            for i=1:n,
                v = an_empty;
                v(i)=1;
                temp = sparse(hankel(v));
                basis(:,i+1) = temp(:);
            end

        case 'hank complex'
            basis = spalloc(n^2,1+nvar,2);
            an_empty = spalloc(n,1,1);
            for i=1:n,
                v = an_empty;
                v(i)=1;
                temp = sparse(hankel(v));
                basis(:,i+1) = temp(:);
            end
            for i=1:n,
                v = an_empty;
                v(i)=sqrt(-1);
                temp = sparse(hankel(v));
                basis(:,n+i+1) = temp(:);
            end

        case 'rhankel'
            basis = spalloc(n^2,1+nvar,2);
            an_empty = spalloc(2*n-1,1,1);
            for i=1:nvar,
                v = an_empty;
                v(i)=1;
                temp = sparse(hankel(v(1:n),[v(n);v(n+1:2*n-1)]));
                basis(:,i+1) = temp(:);
            end

        case 'rhankel complex'
            basis = spalloc(n^2,1+nvar,2);
            an_empty = spalloc(2*n-1,1,1);
            for i=1:nvar/2,
                v = an_empty;
                v(i)=1;
                temp = sparse(hankel(v(1:n),[v(n);v(n+1:2*n-1)]));
                basis(:,i+1) = temp(:);
            end
            for i=1:nvar/2,
                v = an_empty;
                v(i)=sqrt(-1);
                temp = sparse(hankel(v(1:n),[v(n);v(n+1:2*n-1)]));
                basis(:,nvar/2+i+1) = temp(:);
            end

        otherwise
            error('Bug! Report')
    end

end

% Update noncommuting table and monomtables
nonCommutingTable(lmi_variables,1) = nan;
nonCommutingTable(lmi_variables,2) = lmi_variables;
yalmip('nonCommutingTable',nonCommutingTable);
variabletype(lmi_variables) = 0;
monomtable(lmi_variables(end),lmi_variables(end)) = 0;
yalmip('setmonomtable',monomtable,variabletype);


% Create an object
if isa(basis,'cell')
    top = 1;
    for blk = 1:length(n)
        sys{blk}.basis=basis{blk};
        nn = size(sys{blk}.basis,2)-1;
        sys{blk}.lmi_variables = lmi_variables(top:top+nn-1);
        top = top + nn;
        sys{blk}.dim(1) = n(blk);
        sys{blk}.dim(2) = m(blk);
        sys{blk}.typeflag = 0;
        sys{blk}.savedata = [];
        sys{blk}.extra = [];
        sys{blk}.extra.expanded = [];        
        sys{blk}.conicinfo = conicinfo;
        sys{blk} = class(sys{blk},'ncvar');
    end
    if length(n)==1
        sys = sys{1};
    end
else
    sys.basis=basis;
    sys.lmi_variables = lmi_variables;
    sys.dim(1) = n;
    sys.dim(2) = m;
    sys.typeflag = 0;
    sys.savedata = [];
    sys.extra = [];
    sys.extra.expanded = [];    
    sys.conicinfo = conicinfo;
    sys = class(sys,'ncvar');
    if ~isreal(basis)
        % Add internal information about complex pairs
        complex_elements = find(any(imag(basis),2));
        complex_pairs = [];
        for i = 1:length(complex_elements)
            complex_pairs = [complex_pairs;lmi_variables(find(basis(complex_elements(i),:))-1)];
        end
        complex_pairs = uniquesafe(complex_pairs,'rows');
        yalmip('addcomplexpair',complex_pairs);
    end
end