subsref2 (A,i,j)

function C = subsref2 (A,i,j)
% SUBSREF2 often a faster way to compute C=A(i,j) when A is sparse.
% Linear indexing of a 2D array (A(i) when A is a matrix) is not supported.
%
% Usage:
%   C = subsref2 (A,i,j)        % same as C = A(i,j)
%   C = subsref2 (A,i)          % same as C = A(i,:)
%   C = subsref2 (A,{},j)       % same as C = A(:,j)
%   C = subsref2 (A,':',j)      % same as C = A(:,j)
%
% If the i or j inputs are not present, or if they are not numueric, then
% they are treated as a single ':', which means to return all the rows (in
% the case of i) or all the columns (in the case of j).  Best results are
% obtained in MATLAB R2008a (7.6) or later, which has a faster sparse-times-
% sparse matrix multiplication function than previous versions.  i and j can
% be column vectors or matrices, but best results are obtained when they are
% row vectors.
%
% Example:
%
%   A = sprand (10000, 20000, 0.01) ;
%   i = 4000:5000 ;
%   j = randperm (20000) ;   
%   tic ; B = A (i,j) ; toc
%   tic ; C = subsref2 (A,i,j) ; toc
%   B-C
%
% See also subsref, mtimes, sparse.

% Copyright 2009, Timothy A. Davis, University of Florida
% Future work:  the two calls to the "sparse" function can be replaced by
% a simple C mexFunction which would be much faster.

% check inputs
if (nargin < 1 || nargin > 3 || nargout > 1)
    error ('Usage: C = subsref2 (A,i,j)') ;
end
if (~isa (A, 'numeric') || ndims (A) > 2)
    error ('subsref2: A must be a 2D numeric matrix') ;
end

% determine if i and/or j are present
i_present = (nargin >= 2 && isa (i, 'numeric')) ;
j_present = (nargin == 3 && isa (j, 'numeric')) ;

[m n] = size (A) ;

% construct the permutation matrix I
if (i_present)
    if (issparse (i))
        error ('subsref2: i cannot be sparse') ;
    end
    if (size (i,1) ~= 1)
        % make sure i is a row vector
        i = i (:)' ;
    end
    i_len = length (i) ;
    I = sparse (1:i_len, i, 1, i_len, m) ;
end

% construct the permutation matrix J
if (j_present)
    if (issparse (j))
        error ('subsref2: j cannot be sparse') ;
    end
    if (size (j,1) ~= 1)
        % make sure j is a row vector
        j = j (:)' ;
    end
    j_len = length (j) ;
    J = sparse (j, 1:j_len, 1, n, j_len) ;
end

% perform the subsref using sparse matrix multiplication
if (i_present)
    if (j_present)
        % C = A (i,j)
        C = I * A * J ;
    else
        % C = A (i,:)
        C = I * A ;
    end
else
    if (j_present)
        % C = A (:,j)
        C = A * J ;
    else
        % C = A (:,:)
        C = A ;
    end
end

% convert to full if necessary
if (~issparse (A))
    C = full (C) ;
end