mtimes.m

optimization toolbox

                                if before==0
                                    before = [];
                                else
                                    before = possibleOld(before);
                                end
                            else
                                before = possibleOld(findhash(possibleOldHash,current_hash));                                
%                                before = possibleOld(findhash(possibleOldHash,current_hash,length(possibleOldHash)));                                
                            end
                        else
                            before = findhash(new_mt_hash,current_hash); % first among new monomials
%                            before = findhash(new_mt_hash_transpose,current_hash,topp); % first among new monomials
                            if before % Try new if we failed finding any match
                                before=before+size(mt_hash,1);
                            else
                                before = possibleOld(findhash(possibleOldHash,current_hash));                               
%                                before = possibleOld(findhash(possibleOldHash,current_hash,length(possibleOldHash)));                               
                            end
                        end
                        if before
                            Z.lmi_variables(i) = before;
                        else
                            changed_mt=1;
                            isemptynew_mt_hash=0;
                            thesewhereactuallyused(acounter) = 1;
                            if ~bilinearproduct % We don't need to update until later since it not is used inside the loop
                                new_mt_hash = [new_mt_hash;current_hash];
 %                               new_mt_hash_transpose(topp) = current_hash;topp = topp + 1;
                            end
                            nvar = nvar + 1;
                            Z.lmi_variables(i) = nvar;                           
                        end
                        i = i+1;
                    end
                end
                if all(copytofrom)
                    Z.basis = [Z.basis allprodbase];
                else
                    Z.basis = [Z.basis allprodbase(:,find(copytofrom))];
                end
                if all(thesewhereactuallyused)
                    new_mt = [new_mt allmt_xplusy];
                else
                    new_mt = [new_mt allmt_xplusy(:,find(thesewhereactuallyused))];
                end
                if bilinearproduct
                    new_mt_hash = [new_mt_hash;allhash(find(thesewhereactuallyused))];
                end
            end
%             
%             try
%            new_mt_hash_transpose = new_mt_hash_transpose(1:topp-1);
%            norm(new_mt_hash_transpose-new_mt_hash);
%             catch
%                 1
%             end
            
            if ~isempty(new_mt)
                [i1,j1,k1] = find(mt);
                [ii1,jj1,kk1] = find(new_mt');
                uv = find(used_variables);uv=uv(:);
                mt = sparse([i1(:);ii1(:)+size(mt,1)],[j1(:);uv(jj1(:))],[k1(:);kk1(:)],size(mt,1)+size(new_mt,2),size(mt,2));
            end

            % We pre-allocated a sufficiently long, now pick the ones we
            % actually filled with values
            Z.lmi_variables = Z.lmi_variables(1:i-1);

            % Fucked up order (lmi_variables should be sorted)
            if any(diff(Z.lmi_variables)<0)
                [i,j]=sort(Z.lmi_variables);
                Z.basis = [Z.basis(:,1) Z.basis(:,j+1)];
                Z.lmi_variables = Z.lmi_variables(j);
            end

            [un_Z_vars2] = uniquestripped(Z.lmi_variables);
            if length(un_Z_vars2) < length(Z.lmi_variables)
                [un_Z_vars,hh,jj] = unique(Z.lmi_variables);
                if length(Z.lmi_variables) ~=length(un_Z_vars)
                    Z.basis = Z.basis*sparse([1 1+jj],[1 1+(1:length(jj))],ones(1,1+length(jj)))';
                    Z.lmi_variables = un_Z_vars;
                end
            end

            if changed_mt%~isequal(mt,oldmt)
                newmt = mt(size(oldmt,1)+1:end,:);
                nonlinear = ~(sum(newmt,2)==1 & sum(newmt~=0,2)==1);
                newvariabletype = spalloc(size(newmt,1),1,nnz(nonlinear))';
                nonlinearvariables = find(nonlinear);
                newvariabletype = sparse(nonlinearvariables,ones(length(nonlinearvariables),1),3,size(newmt,1),1)';
                if ~isempty(nonlinear)
                    %mt = internal_sdpvarstate.monomtable;
                    %newvariabletype(nonlinear) = 3;
                    quadratic = sum(newmt,2)==2;
                    newvariabletype(quadratic) = 2;
                    bilinear = max(newmt,[],2)<=1;
                    newvariabletype(bilinear & quadratic) = 1;
                    sigmonial = any(0>newmt,2) | any(newmt-fix(newmt),2);
                    newvariabletype(sigmonial) = 4;
                end
                yalmip('setmonomtable',mt,[oldvariabletype newvariabletype],[mt_hash;new_mt_hash],oldhash);
            end

            if ~(x_isscalar | y_isscalar)
                Z.dim(1) = X.dim(1);
                Z.dim(2) = Y.dim(2);
            else
                Z.dim(1) = max(X.dim(1),Y.dim(1));
                Z.dim(2) = max(X.dim(2),Y.dim(2));
            end
        catch
            error(lasterr)
        end
        % Reset info about conic terms
        Z.conicinfo = [0 0];
        Z = clean(Z);

    case 2

        n_X = X.dim(1);
        m_X = X.dim(2);
        [n_Y,m_Y] = size(Y);

        x_isscalar =  (n_X*m_X==1);
        y_isscalar =  (n_Y*m_Y==1);

        if ~x_isscalar
            if ((m_X~= n_Y & ~y_isscalar))
                error('Inner matrix dimensions must agree.')
            end
        end

        n = n_X;
        m = m_Y;
        Z = X;

        if x_isscalar
            if y_isscalar
                if Y==0
                    Z = 0;
                    return
                else
                    Z.basis = Z.basis*Y;
                    % Reset info about conic terms
                    Z.conicinfo = [0 0];
                    return                                
                end                    
            else
                Z.dim(1) = n_Y;
                Z.dim(2) = m_Y;
                Z.basis = kron(Z.basis,Y(:));
                Z.conicinfo = [0 0];
                Z = clean(Z);
                return
            end
        elseif y_isscalar
            Z.dim(1) = n_X;
            Z.dim(2) = m_X;
            Z.basis = Z.basis*Y;
            Z.conicinfo = [0 0];
            Z = clean(Z);
            return
        end

        Z.dim(1) = n;
        Z.dim(2) = m;
        Z.basis = kron(Y.',speye(n_X))*X.basis;
        Z.conicinfo = [0 0];
        Z = clean(Z);

    case 1

        n_Y = Y.dim(1);
        m_Y = Y.dim(2);
        [n_X,m_X] = size(X);

        x_isscalar =  (n_X*m_X==1);
        y_isscalar =  (n_Y*m_Y==1);

        if ~x_isscalar
            if ((m_X~= n_Y & ~y_isscalar))
                error('Inner matrix dimensions must agree.')
            end
        end

        n = n_X;
        m = m_Y;
        Z = Y;

        % Special cases
        if x_isscalar
            if y_isscalar
                if X==0
                    Z = 0;
                    return
                else
                    Z.basis = Z.basis*X;
                    Z.conicinfo = [0 0];
                    return
                end
            else
                Z.dim(1) = n_Y;
                Z.dim(2) = m_Y;
                try
                    Z.basis = sparse(X)*Y.basis;
                catch
                    % This works better when low on memory in some cases
                    [i,j,k] = find(Y.basis);
                    Z.basis = sparse(i,j,X*k,size(Y.basis,1),size(Y.basis,2));
                end
                Z.conicinfo = [0 0];
                Z = clean(Z);
                return
            end
        elseif y_isscalar
            Z.dim(1) = n_X;
            Z.dim(2) = m_X;
            Z.basis = X(:)*Y.basis;
            Z = clean(Z);
            return
        end

        if m_Y==1
            Z.basis = X*Y.basis;
        else
            try               
                speyemy = speye(m_Y);
                kronX = kron(speyemy,X);
                Z.basis = kronX*Y.basis;
            catch
                for i = 1:size(Y.basis,2);
                    dummy = X*reshape(Y.basis(:,i),Y.dim(1),Y.dim(2));
                    Z.basis(:,i) = dummy(:);
                end
            end
        end
        Z.dim(1) = n;
        Z.dim(2) = m;
        Z.conicinfo = [0 0];
        Z = clean(Z);

    otherwise
        error('Logical error in mtimes. Report bug')
end



function Z=clean(X)
temp = any(X.basis,1);
temp = temp(2:end);
index = find(temp);
if ~isempty(index)
    Z = X;
    if length(index)~=length(Z.lmi_variables)
        Z.basis = Z.basis(:,[1 1+index]);
        Z.lmi_variables = Z.lmi_variables(index);
    end
else
    Z = full(reshape(X.basis(:,1),X.dim(1),X.dim(2)));
end