Issue when passing sparse arrays to function: "Sparse single array arithmetic operations are not supported"

Hello,

I am running a code where I define some sparse arrays that are pages of a cell structure, which are subsequently fed to a function in a for loop. I got the error "Sparse single array arithmetic operations are not supported" when running Matlab on cluster. However, I do not use any (at least explicitly) single point variables. Furthermore, this problem does not happen if I run the code on my local workstation instead of the cluster.

The error I got is:

'Error using SSPOD6_partial

Sparse single array arithmetic operations are not supported.

Error in run_SPOD_airfoil_case_4A_6 (line 377)

SSPOD6_partial(squeeze(UVW_pert_fft_span_only(:,t,bb)),k,window,W,sqrtW,winWeight,...'

Here is a snippet of the code giving the problem:

            % definition of the sparse arrays
            t_idx       = zeros(nBlks,1,Nb);                                              % sliding, relative time index for each block
            for blk_i = 1:nBlks
                t_idx(blk_i,:,:) =  t_idx(blk_i,:,:) - (blk_i-1)*dn + 1;
            end
            X_sum        = cell(nBlks,Nb);
            [X_sum{:,:}] = deal(sparse(zeros(nvars*numelx,nFreq)));

            U_hat        = zeros(nvars*numelx,nFreq,k,Nb);

            S_hat        = cell(Nb);
            [S_hat{:}]   = deal(sparse(zeros(k,nFreq)));

            mu           = cell(Nb);
            [mu{:}]      = deal(sparse(zeros(nvars*numelx,1)));

            mse_prev    = nan(1e3,k,nFreq,Nb);
            proj_prev   = nan(nFreq,1e3,k,Nb);
            blk_i       = ones(1,Nb);
            t_i         = zeros(1,Nb);            

           % location of issue
           for bb = 1:Nb
		        for t = 1:n4
		            [S_hat{bb},t_i(bb),blk_i(bb),mu{bb},X_sum(:,bb),...
		             U_hat(:,:,:,bb),t_idx(:,:,bb),mse_prev(:,:,:,bb),proj_prev(:,:,:,bb)] =...
		            SSPOD6_partial(squeeze(UVW_pert_fft_span_only(:,t,bb)),k,window,W,sqrtW,winWeight,...
		            nDFT,nFreq,dn,nBlks,t_idx(:,:,bb),X_sum(:,bb),U_hat(:,:,:,bb),S_hat{bb},...
		            mu{bb},Fourier,mse_prev(:,:,:,bb),proj_prev(:,:,:,bb),blk_i(bb),t_i(bb),update_avg);
		        end
            end

The function SSPOD6_partial is defined as

function [S_hat,t_i,blk_i,mu,X_sum,U_hat,t_idx,mse_prev,proj_prev] =...
    SSPOD6_partial(getdata,k,window,W,sqrtW,winWeight,...
    nDFT,nFreq,dn,nBlks,t_idx,X_sum,U_hat,S_hat,...
    mu,Fourier,mse_prev,proj_prev,blk_i,t_i,update_avg)

z           = zeros(1,k);    

X_SPOD_prev      = bsxfun(@times,U_hat,1./sqrtW);                                % rescale such that <U_i,U_j>_E = U_i^H*W*U_j = delta_ij

t_i = t_i + 1;
% disp(['Time ' num2str(t_i) ' / block ' num2str(blk_i)])

% Get new snapshot and abort if data stream runs dry
x_new   = getdata;
x_new   = x_new(:);

% Update sample mean
if update_avg
    mu_old  = mu;
    mu      = ((t_i-1)*mu_old + x_new)/t_i;
end

% Update incomplete Fourier sums, eqn (17)
update  = false;
for blk_j = 1:nBlks
    if t_idx(blk_j)>0
        X_sum{blk_j}      = X_sum{blk_j} + window(t_idx(blk_j))*Fourier(t_idx(blk_j),:).*x_new;
    end
    % check if sum is completed, and if so, initiate update
    if t_idx(blk_j)==nDFT      
        update              = true;
        X_hat               = X_sum{blk_j}; 
        X_sum{blk_j}        = 0;                                              
        t_idx(blk_j)        = min(t_idx) - dn;
    else
        t_idx(blk_j)        = t_idx(blk_j) + 1;
    end
end

% Update basis if a Fourier sum is completed 
if update
    blk_i               = blk_i + 1;
    
    % subtract mean contribution to Fourier sum
    if update_avg
        for row_idx=1:nDFT
            X_hat               = X_hat - window(row_idx)*Fourier(row_idx,:).*mu;
        end
    end
    
    % correct for windowing function and apply 1/nDFT factor
    X_hat               = winWeight/nDFT*X_hat;
         
    if blk_i==1
        % initialize basis with first vector
%         disp('--> Initializing left singular vectors')    
        U_hat(:,:,1)    = bsxfun(@times,X_hat,sqrtW);
        S_hat(1,:)      = sum(abs(U_hat(:,:,1).^2));
    else
        % update basis
%         disp('--> Updating left singular vectors')
        S_hat_prev  = S_hat;
        for fi = 1:nFreq
            
            x       = X_hat(:,fi).*sqrtW;                               % new data (weighted)
            U       = squeeze(U_hat(:,fi,:));                           % old basis
            S       = diag(squeeze(S_hat(:,fi)));                       % old singular values
            
            Ux      = U'*x;                                             % product U^H*x needed in eqns. (27,32)            
            u_p     = x-U*Ux;                                           % orthogonal complement to U, eqn. (27)
            abs_up  = sqrt(u_p'*u_p);                                   % norm of orthogonal complement
            u_new   = u_p/abs_up;                                       % normalized orthogonal complement
            
            % build K matric and compute its SVD, eqn. (32) 
            K       = sqrt(blk_i/(blk_i+1)^2)*[sqrt(blk_i+1)*S Ux; z abs_up];                
            [Up,Sp,~] = svd(full(K),'econ');

            Up = sparse(Up);
            Sp = sparse(Sp);
            
            % update U as in eqn. (33)
            U       = ([U u_new]*Up);                                   % for simplicity, we could not rotate here and instead update U<-[U p] and Up<-[Up 0;0 1]*Up and rotate later; see Brand (LAA ,2006, section 4.1)
            
            % best rank-k approximation, eqn. (37)
            U_hat(:,fi,:)   = U(:,1:k);
            S_hat(:,fi)     = diag(Sp(1:k,1:k));
            
        end
        X_hat(:,:)      = 0;                                            % reset Fourier sum
    end
    
    X_SPOD      = bsxfun(@times,U_hat,1./sqrtW);
    
    % Convergence
    for fi = 1:nFreq
        proj_fi                     = bsxfun(@times,squeeze(X_SPOD_prev(:,fi,:)),W)'*squeeze(X_SPOD(:,fi,:));
        [proj_prev(fi,blk_i,:),~]   = max(abs(proj_fi));
    end
    mse_prev(blk_i,:,:)             = (abs(S_hat_prev.^2 - S_hat.^2).^2)./S_hat_prev.^2;
    
end

Thank you for your answer.

[Walter Roberson]

It is odd that it is not giving you a line number and line of code inside SSPOD6_partial even though you have the source code. That hints that there might be a different SSPOD6_partial being executed. Try using

which -all SSPOD6_partial

to explore what versions are known to MATLAB.

[Thales Coelho Leite Fava]

I've run

which -all SSPOD6_partial

in the running folder. There is only one version known to Matlab, which corresponds to the file

/xxx/matlab/SSPOD6_partial.m

as expected by me. I wonder if the problem occurs before SSPOD6_partial is actually executed, i.e., in the phase where the input arguments are passed (therefore the absence of a line number for the error inside SSPOD6_partial).

[Walter Roberson]

In order for SSPOD6_partial to be mentioned in the error message, some routine named SSPOD6_partial needs to have started running (including possibly SSPOD6_partial.p or SSPOD6_partial.slx or SSPOD6_partial.mexa64 ); or else SSPOD6_partial needs to have been invoked with too many parameters; or else the algorithm for resolving class methods needs to decide that there is no SSPOD6_partial applicable for a class of the dominant parameter of the function call; or else (not sure on this one) SSPOD6_partial needs to have an arguments block [I have not experimented with that possibility myself]

If the problem were occuring while preparing the input arguments for the call to SSPOD6_partial then the error would be reported against the calling routine. For example if the problem was being triggered in the section of the line that says squeeze(UVW_pert_fft_span_only(:,t,bb)) then the error would be reported against squeeze or against run_SPOD_airfoil_case_4A_6

[Christine Tobler]

The error message here would be thrown if a double-precision sparse matrix is combined with a single-precision dense array, for example in a call to + or times:

x = sparse([1 0; 2 3]);
y = single([1 2; 3 4]);
x+y

Error using +
Sparse single array arithmetic operations are not supported.

This is because following the datatype propagation rules in MATLAB, the output would need to be a single-precision sparse matrix, which is not supported.

I could see this error happening in your SSPOD6_partial function if some of the inputs are single-precision. It's still quite mysterious why you would have different behavior when running on the cluster. I'm not an expert on running MATLAB on cluster, but is it possible that something is configured in your cluster setup to convert inputs to single-precision for performance?

One idea would be to call class on all input variables to SSPOD6_partial to check if any of them are single-precision when running on the cluster.