mtimes.m

optimization toolbox

function Z = mtimes(X,Y)
%MTIMES (overloaded)

% Author Johan L鰂berg
% $Id: mtimes.m,v 1.57 2006/10/19 12:42:36 joloef Exp $

% Check classes
X_is_spdvar = isa(X,'sdpvar');
Y_is_spdvar = isa(Y,'sdpvar');

% Convert block objects
if ~X_is_spdvar
    if isa(X,'blkvar')
        X = sdpvar(X);
        X_is_spdvar = isa(X,'sdpvar');
    end
end

if ~Y_is_spdvar
    if isa(Y,'blkvar')
        Y = sdpvar(Y);
        Y_is_spdvar = isa(Y,'sdpvar');
    end
end

% Lame special cases, make sure to reurn
% empty matrices in the sense that the
% used MATLAB version
if isempty(X)
    YY = full(reshape(Y.basis(:,1),Y.dim(1),Y.dim(2)));
    Z = X*YY;
    return
elseif isempty(Y)
    XX = full(reshape(X.basis(:,1),X.dim(1),X.dim(2)));
    Z = XX*Y;
    return
end


% Optimized calls in different order?
if X_is_spdvar & Y_is_spdvar
    manytimesfew = length(X.lmi_variables) > 5*length(Y.lmi_variables);
    if manytimesfew
        Z = (Y'*X')'; % Optimized for this order (few variables * many variables)
        return
    end
end

% Different code for
% 1 : SDPVAR * DOUBLE
% 2 : DOUBLE * SDPVAR
% 3 : SDPVAR * SDPVAR
switch 2*X_is_spdvar+Y_is_spdvar
    case 3
        try
            x_isscalar =  (X.dim(1)*X.dim(2)==1);
            y_isscalar =  (Y.dim(1)*Y.dim(2)==1);

            % Optimized unique
            all_lmi_variables = uniquestripped([X.lmi_variables Y.lmi_variables]);

            % Create clean SDPVAR object
            Z = X;Z.dim(1) = 1;Z.dim(2) = 1;Z.lmi_variables = all_lmi_variables;Z.basis = [];

            % Awkward code due to bug in ML6.5
            Xbase = reshape(X.basis(:,1),X.dim(1),X.dim(2));
            Ybase = reshape(Y.basis(:,1),Y.dim(1),Y.dim(2));
            if x_isscalar
                Xbase = sparse(full(Xbase));
            end
            if y_isscalar
                Ybase = sparse(full(Ybase));
            end

            index_X = double(ismembc(all_lmi_variables,X.lmi_variables));
            index_Y = double(ismembc(all_lmi_variables,Y.lmi_variables));
            iX=find(index_X);
            iY=find(index_Y);
            index_X(iX)=1:length(iX);index_X=index_X(:);
            index_Y(iY)=1:length(iY);index_Y=index_Y(:);

            ny = Y.dim(1);
            my = Y.dim(2);
            nx = X.dim(1);
            mx = X.dim(2);

            % Pre-allocate sufficiently long
            Z.lmi_variables = [Z.lmi_variables zeros(1,length(X.lmi_variables)*length(Y.lmi_variables))];

            % Pre-calc identity (used a lot
            speyemy = sparse(1:my,1:my,1,my,my);

            % Linear terms
            inner_vector_product = (X.dim(1)==1 & Y.dim(2)==1 & (X.dim(2) == Y.dim(1)));
            if inner_vector_product
                base1=Xbase*Y.basis;base1=base1(2:end);
                base2=Ybase.'*X.basis;base2=base2(2:end);
                [i1,j1,k1]=find(base1);
                [i2,j2,k2]=find(base2);
                base1 = sparse(i1,iY(j1),k1,1,length(all_lmi_variables));
                base2 = sparse(i2,iX(j2),k2,1,length(all_lmi_variables));
                Z.basis = [Xbase*Ybase base1+base2];
            else
                base0 = Xbase*Ybase;
                if x_isscalar
                    base1 = Xbase*Y.basis(:,2:end);
                    base2 = Ybase(:)*X.basis(:,2:end);
                elseif y_isscalar
                    base1 = Xbase(:)*Y.basis;base1=base1(:,2:end);
                    base2 = X.basis*Ybase;base2=base2(:,2:end);
                else
                    base1 = kron(speyemy,Xbase)*Y.basis(:,2:end);
                    base2 = kron(Ybase.',speye(nx))*X.basis(:,2:end);
                end
                [i1,j1,k1]=find(base1);
                [i2,j2,k2]=find(base2);
                base1 = sparse(i1,iY(j1),k1,size(base0(:),1),length(all_lmi_variables));
                base2 = sparse(i2,iX(j2),k2,size(base0(:),1),length(all_lmi_variables));
                Z.basis = [base0(:) base1+base2];
            end

            % Loop start for nonlinear terms
            i = length(all_lmi_variables)+1;

            [mt,oldvariabletype,mt_hash,hash] = yalmip('monomtable');

            % Check if problem is bilinear. We can exploit this later
            % to improve performance significantly...
            bilinearproduct = 0;
            candofastlocation  = 0;
            if all(oldvariabletype(X.lmi_variables)==0) & all(oldvariabletype(Y.lmi_variables)==0)
               % if isempty(intersect(X.lmi_variables,Y.lmi_variables))
                if ~any(ismembc(X.lmi_variables,Y.lmi_variables))
                    bilinearproduct = 1;
                    try
                        dummy = ismembc2(1,1); % Not available in all versions (needed in ismember)
                        candofastlocation = 1;
                    catch
                    end
                end
            end

            oldmt = mt;
            local_mt = mt(all_lmi_variables,:);
            used_variables = any(local_mt,1);
            local_mt = local_mt(:,used_variables);

            possibleOld = find(any(mt(:,used_variables),2));
            if all(oldvariabletype <=3)
                % All monomials have non-negative integer powers
                % no chance of x^2*x^-1, hence all products
                % are nonlinear
                possibleOld = possibleOld(find(oldvariabletype(possibleOld)));
                if size(possibleOld,1)==0
                    possibleOld = [];
                end
            end

            if bilinearproduct & ~isempty(possibleOld)
                if length(X.lmi_variables)<=length(Y.lmi_variables)
                    temp = mt(:,X.lmi_variables);
                    temp = temp(possibleOld,:);
                    possibleOld=possibleOld(find(any(temp,2)));
                else
                    temp = mt(:,Y.lmi_variables);
                    temp = temp(possibleOld,:);
                    possibleOld=possibleOld(find(any(temp,2)));
                end
            end

            theyvars = find(index_Y);
            thexvars = find(index_X);

            possibleOldHash = mt_hash(possibleOld);
            oldhash = hash;
            hash = hash(used_variables);
            new_mt_hash = [];
            new_mt_hash_transpose = [];
            new_mt = [];
            changed_mt = 0;
            local_mt = local_mt';
            nvar = size(mt,1);

            for ix = thexvars(:)'

                mt_x = local_mt(:,ix);

                testthese = theyvars(:)';

                % Compute [vec(Xbasis*Ybasis1) vec(Xbasis*Ybasis2) ...]
                % in one shot using vectorization and Kronecker tricks
                % Optimized and treat special case scalar*matrix etc
                if x_isscalar
                    Xibase = X.basis(:,1+index_X(ix));
                    allprodbase = Xibase * Y.basis(:,1+index_Y(testthese));
                elseif y_isscalar
                    Xibase = X.basis(:,1+index_X(ix));
                    allprodbase =  Xibase * Y.basis(:,1+index_Y(testthese));
                elseif inner_vector_product
                    Xibase = X.basis(:,1+index_X(ix)).';
                    allprodbase = Xibase*Y.basis(:,1+index_Y(testthese));
                else
                    Xibase = reshape(X.basis(:,1+index_X(ix)),nx,mx);                  
                    temp = kron(speyemy,Xibase);
                    allprodbase = temp * Y.basis(:,1+index_Y(testthese));
                end

                % Keep non-zero matrices
                nonzeromatrices = find(sum(abs(allprodbase),1)>1e-12);
                testthese = testthese(nonzeromatrices);
                allprodbase = allprodbase(:,nonzeromatrices);

                % Some data for vectorization
                nyvars = length(testthese);
                if prod(size(mt_x))==1 % Bug in Solaris and Linux, ML 6.X
                    allmt_xplusy = local_mt(:,testthese) + sparse(repmat(full(mt_x),1,nyvars));
                else
                    allmt_xplusy = local_mt(:,testthese) + repmat(mt_x,1,nyvars);
                end
                allhash = allmt_xplusy'*hash;
                thesewhereactuallyused = zeros(1,nyvars);
                copytofrom  = ones(1,nyvars);
                acounter =  0;


                % Special case : x*inv(x) and similiar...
                sum_to_constant = abs(allhash)<eps;
                add_these = find(sum_to_constant);
                if ~isempty(add_these)
                    prodbase = allprodbase(:,add_these);
                    Z.basis(:,1) = Z.basis(:,1) + sum(prodbase,2);
                    copytofrom(add_these) = 0;
                end
                indicies = find(~sum_to_constant);
                indicies = indicies(:)';

                allbefore_in_old = 1;
                if bilinearproduct & candofastlocation
                    [dummy,allbefore_in_old] = ismember(allhash,possibleOldHash);
                end


                if bilinearproduct & candofastlocation & (nnz(allbefore_in_old)==0)
                    % All nonlinear variables are new, so we can create them at once
                    changed_mt=1;
                    thesewhereactuallyused = thesewhereactuallyused+1;
                    Z.lmi_variables(i:(i+length(indicies)-1)) = (nvar+1):(nvar+length(indicies));
                    nvar = nvar + length(indicies);
                    i = i + length(indicies);
                else
                    isemptynew_mt_hash = isempty(new_mt_hash);                    
                   
                    for acounter = indicies

                        current_hash = allhash(acounter);

                        % Ok, braze your self for some horrible special case
                        % treatment etc...
                        if isemptynew_mt_hash | bilinearproduct % only search among old monomials
                            if bilinearproduct & candofastlocation
                                before = allbefore_in_old(acounter);