mtimes.m

matlab波形优化算法经常要用到的matlab toolbox工具箱:yalmip

                            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
            
            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)
            Z = fix_variable_order(Z);
%             
%             % 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);
%         
%          if exist('ZZ')
%                 full(Z.basis)
%                 full(ZZ.basis)
%             end

    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;
        if (n_X==1) & is(X,'lpcone') & (n_Y == m_Y) & (size(X.basis,1)==size(X.basis,2-1)) & isequal(X.basis*[0 1:size(X.basis,2)-1]',(1:size(X.basis,2)-1)')
            % special case to speed up x'*Q, Q square. typically
            % encountered in large-scale QPs
            Z.basis = [X.basis(:,1) Y.'];
        else
            Z.basis = kron(Y.',speye(n_X))*X.basis;
        end
        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



function Z = fix_variable_order(Z)
% 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

function Z = super_fast_quadratic_multiplication(X,Y,mt,oldvariabletype,mt_hash,hash);

Q = X.basis(:,2:end);
R = Y.basis(:,2:end);
Q = (Q.')*R;
Q = Q + Q.' - diag(diag(Q));
indicies = find(tril(ones(length(Q))));
Q = Q(indicies);
Q = Q(:).';
% Q = Q + Q' - diag(diag(Q));
%  Q = Q(find(tril(ones(length(Q)))));
%  Q = Q(:)';

n = length(X.lmi_variables);
V = mt(X.lmi_variables,:);
if 1
    m1 = kron((1:n)',ones(n,1));
    m2 = kron(ones(n,1),(1:n)');
    r = reshape(1:n^2,n,n);
    r = r(find(tril(r)));
    m1 = m1(r);
    m2 = m2(r);
    VV = V';
    VV = VV(:,m1) + VV(:,m2);
    new_mt = VV';
else
    new_mt = kron(V,ones(n,1)) + kron(ones(n,1),V);
    r = reshape(1:n^2,n,n);
    new_mt = new_mt(r(find(tril(r))),:);
end
%newvariabletype = 3 - sum(new_mt | new_mt,2);
newvariabletype = max((new_mt),[],2);
yalmip('setmonomtable',[mt;new_mt],[oldvariabletype newvariabletype'],[mt_hash;new_mt*hash],hash);
Z = X;
varbase = [(X.basis(:,1).')*Y.basis(:,2:end)+(Y.basis(:,1).')*X.basis(:,2:end) Q];
Z.basis = [(X.basis(:,1).')*Y.basis(:,1) varbase(find(varbase))];
Z.lmi_variables = [X.lmi_variables size(mt,1) + (1:length(Q))];
Z.lmi_variables = Z.lmi_variables(find(varbase));
Z.dim(1) = 1;
Z.dim(2) = 1;
Z.conicinfo = [0 0];