function R = sprand(arg1,n,density,rc)
%SPRAND Sparse uniformly distributed random matrix.
%   R = SPRAND(S) has the same sparsity structure as S, but uniformly
%       distributed random entries.
%
%   R = SPRAND(m,n,density) is a random, m-by-n, sparse matrix with 
%       approximately density*m*n uniformly distributed nonzero entries.
%       SPRAND is designed to produce large matrices with small density
%       and will generate significantly fewer nonzeros than requested
%       if m*n is small or density is large.
%
%   R = SPRAND(m,n,density,rc) also has reciprocal condition number
%       approximately equal to rc.  R is constructed from a sum of
%       matrices of rank one. 
%
%       If rc is a vector of length lr <= min(m,n), then R has 
%       rc as its first lr singular values, all others are zero.
%       In this case, R is generated by random plane rotations
%       applied to a diagonal matrix with the given singular values
%       It has a great deal of topological and algebraic structure.
%
%   See also SPRANDN, SPRANDSYM.

%   Rob Schreiber, RIACS, and Cleve Moler, MathWorks, 9/10/90.
%   Revised 1/28/91, 2/12/91, RS; 8/12/91, CBM.
%   Copyright 1984-2003 The MathWorks, Inc. 
%   $Revision: 1.11.4.1 $  $Date: 2003/05/01 20:43:13 $

if nargin == 1
   S = arg1;
   [m,n] = size(S);
   [i,j] = find(S);
   R = sparse(i,j,rand(length(i),1),m,n);
elseif nargin == 2
   error('MATLAB:sprand:TwoInputs', 'Too many or not enough input arguments.')
elseif nargin == 3
   m = arg1;
   nnzwanted = round(m * n * min(density,1));
   i = fix( rand(nnzwanted, 1) * m ) + 1;
   j = fix( rand(nnzwanted, 1) * n ) + 1;
   ij = unique([i j],'rows');
   if ~isempty(ij)
      i = ij(:,1);
      j = ij(:,2);
   end
   rands = rand( length(i), 1 );
   R = sparse(i,j,rands,m,n);
else    %    nargin == 4
   m = arg1;
   nnzwanted = round(m * n * min(density,1));
   minm = min(m,n);
   maxmn = max(m,n);
   lr = length(rc);
   if (lr > 1)  %  specified singular values
      if (lr > minm), lr = minm; end;
      %
      %   To start, put a random nonzero in every column / row if
      %   there is room enough
      %
      anz = zeros(maxmn, 1);
      anz(1:lr) = rc(1:lr);
      lr = minm;      %   if lr < minm then this adds some zero s.v. 
      sqrt2o2 = sqrt(2)/2;
      while (lr < min(maxmn, nnzwanted))
         ndo = min([lr, nnzwanted-lr, maxmn-lr]);
         anz(lr+1 : lr+ndo) = sqrt2o2 * anz(1:ndo);
         anz(1 : ndo) = sqrt2o2 * anz(1:ndo);
         lr = lr + ndo;
      end
      [ignore,p] = sort(rand(m,1));
      [ignore,q] = sort(rand(n,1));
      if (m < n),
         p = p(1 + rem((0:n-1), m));
      else
         q = q(1 + rem((0:m-1), n));
      end
      R = sparse(p, q, anz, m, n);
      %
      %   random two-sided rotations
      %
      while ( nnz(R) < .95*nnzwanted )
          R = rjr(R,2);
      end
   else    %   condition number specified
      if rc > 1, error('MATLAB:sprand:RcondTooHigh',...
                       'Reciprocal condition number must be 1 or lower.'); 
      end
      MEMFAC = 1.05;
      NZF = 1.25;
      i = zeros(fix(nnzwanted*MEMFAC),1); j = zeros(fix(nnzwanted*MEMFAC),1); 
      anz = zeros(fix(nnzwanted*MEMFAC),1); nnza = 0;
      %
      %   To start, put a random nonzero in every column / row
      %
      how_big = 1;
      [ignore,p] = sort(rand(m,1));
      [ignore,q] = sort(rand(n,1));
      ratio = rc ^ (1/(minm-1));
      if (m < n),
         for jj = 1:n
            i(jj) = p(1 + rem(jj-1,m));
            j(jj) = q(jj);
            anz(jj) = how_big;
            how_big = how_big * ratio;
         end
         nnza = n;
      else
         for ii = 1:m
            j(ii) = q(1 + rem(ii-1,n));
            i(ii) = p(ii);
            anz(ii) = how_big;
            how_big = how_big * ratio;
         end
         nnza = m;
      end
      %
      %   Now add a sum of scaled outer products with random structure
      %
      how_big = .1;
      for iouter = 1:minm,
         nzpc = fix(1 + sqrt((m/n)*(NZF*nnzwanted - nnza)/(minm-iouter+1)));
         nzpr = fix(1 + sqrt((n/m)*(NZF*nnzwanted - nnza)/(minm-iouter+1)));
         ak = sprand(m,1, nzpc/m);
         bk = sprand(n,1, nzpr/n);
         [ik,jk,anzk] = find(how_big * (ak * bk'));
         lk = length(ik);
         i(nnza+1:nnza+lk) = ik; 
         j(nnza+1:nnza+lk) = jk; 
         anz(nnza+1:nnza+lk) = anzk;
         nnza = nnza + lk;
         how_big = how_big * ratio;
         if (nnza > nnzwanted), break; end;
      end
      i = i(1:nnza);
      j = j(1:nnza);
      anz = anz(1:nnza);
      R = sparse(i,j,anz);
   end  %   rc gives condition number
end     %   nargin == 4