function R = sprandsym(arg1, density, rcond, kind)
%SPRANDSYM Sparse random symmetric matrix.
%   R = SPRANDSYM(S) is a symmetric random matrix whose lower triangle
%       and diagonal have the same structure as S.  The elements are
%       normally distributed, with mean 0 and variance 1.
%
%   R = SPRANDSYM(n,density) is a symmetric random, n-by-n, sparse 
%       matrix with approximately density*n*n nonzeros; each entry is
%       the sum of one or more normally distributed random samples.
%
%   R = SPRANDSYM(n,density,rc) also has a reciprocal condition number
%       equal to rc.  The distribution of entries is
%       nonuniform; it is roughly symmetric about 0; all are in [-1,1].
%
%       If rc is a VECTOR of length n, then R has eigenvalues rc.
%       Thus, if rc is a positive (nonnegative) vector then R will 
%       be positive (nonnegative) definite. In either case, R is
%       generated by random Jacobi rotations applied to a diagonal
%       matrix with the given eigenvalues or condition number.  It has
%       a great deal of topological and algebraic structure.
%
%   R = SPRANDSYM(n, density, rc, kind) is positive definite.
%
%       If kind = 1, R is generated by random Jacobi rotation of
%          a positive definite diagonal matrix.
%          R has the desired condition number exactly.
%
%       If kind = 2, R is a shifted sum of outer products.
%          R has the desired condition number only
%          approximately, but has less structure.
%
%   R = SPRANDSYM(S,[],rc,3) has the same structure as the MATRIX S
%       and approximate condition number 1/rc.
%
%   See also SPRANDN, SPRAND.

%   Rob Schreiber, RIACS, and Cleve Moler, MathWorks, 2/5/91.
%   Copyright 1984-2004 The MathWorks, Inc. 
%   $Revision: 5.11.4.3 $  $Date: 2004/03/02 21:48:30 $

if nargin > 4
    help sprandsym
    error('Too many or not enough input arguments.')
elseif nargin == 1
    S = arg1;
    [m,n] = size(S);
    [ii,jj] = find(S);
    lower = find(ii > jj);  %  find lower triangle
    di = find( ii == jj );  %  find diagonal
    nd = length( di );
    nl = length( lower );
    randi = randn( nd, 1 );
    randl = randn( nl, 1 );
    i = [ii(lower); jj(lower); ii(di)];
    j = [jj(lower); ii(lower); ii(di)];
    r = [randl; randl; randi];
    R = sparse(i,j,r,m,m);
elseif nargin == 2,
    n = arg1;
    nl = round( (n*(n+1)/2) * density );
    rands = randn( nl, 1 );
    ii = fix( rand(nl, 1) * n ) + 1;
    jj = fix( rand(nl, 1) * n ) + 1;
    di = find( ii == jj );
    off = find( ii ~= jj );
    nd = length( di );
    no = length( off );
    randi = rands( 1:nd );
    rando = rands( nd+1 : nl );
    i = [ii(off); jj(off); ii(di)];
    j = [jj(off); ii(off); ii(di)];
    r = [rando; rando; randi];
    R = sparse(i,j,r,n,n);
else % nargin > 2
    lrc = length(rcond);
    if (lrc == 1)
        if( rcond <= 0 || rcond > 1 )
            error('rc must be > 0 and <= 1.')
        end
    end
    if(nargin == 3 )
        n = arg1;
        nnzwanted = round(density*n*n);
        if (lrc == 1)   %   rcond is given
            ratio = ((-1)^nargin) * (rcond ^ (1/(n-1)));
            anz = (ratio .^ ( 0:(n-1)));
        else            % eigenvalues explicitly given
            anz = rcond;
        end
        R = sparse(1:n,1:n, anz);
        % for nargin == 3,  ratio < 0
        % so we start with an indefinite R;
        nnzr = n;
        %
        %   random jacobi rotations
        %
        while ( nnzr < .95*nnzwanted )
            R = rjr(R);
            nnzr = nnz(R);
        end
    else  % nargin == 4
        if(kind==3)
            %
            %  add a 2-by-2 rank-1 for each off diagonal
            %
            [ii,jj] = find(arg1); ii = ii(:);  % workaround
            [m, n] = size(arg1);
            lower = find(ii > jj);
            nl = length(lower);
            ratio = rcond ^ (1/(nl-1));
            ii = ii(lower);
            jj = jj(lower);
            i = zeros(n + 4*nl,1); 
            j = zeros(n + 4*nl,1); 
            anz = zeros(n + 4*nl,1); 
            [ignore,p] = sort(rand(n,1));
            i(1:n) = p;
            j(1:n) = p;
            en = sqrt((28/45) * (1 - ratio^(2*nl)) / (1 - ratio^2));
            anz(1:n) = rcond * en * ones(n,1);
            nnza = n;
            how_big = 1;
            for iouter = 1:nl
                rv = rand(2,1);
                op = how_big * rv * rv';
                op = op(:);
                ik = ii(iouter);
                jk = jj(iouter);
                i(nnza+1:nnza+4) = [ik; jk; ik; jk];
                j(nnza+1:nnza+4) = [ik; ik; jk; jk];
                anz(nnza+1:nnza+4) = op;
                nnza = nnza + 4;
                how_big = how_big * ratio;
            end
            R = sparse(i(1:nnza),j(1:nnza),anz(1:nnza));
            %  end of kind == 3
        else  %  kind == 1 or 2
            n = arg1;  
            if ~isscalar(n)
                error('First argument must be a scalar');
            end
            nnzwanted = round(density*n*n);
            ratio = ((-1)^nargin) * (rcond ^ (1/(n-1)));
            if (kind == 1)
                R = sparse(1:n,1:n, (ratio .^ ( 0:(n-1)))  );
                nnzr = n;
                %   random jacobi rotations
                while ( nnzr < .95*nnzwanted )
                    R = rjr(R);
                    nnzr = nnz(R);
                end
            elseif(kind==2)   %  positive definite via sum of outer products
                MEMFAC = 2;
                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
                %
                [ignore,p] = sort(rand(n,1));
                i(1:n) = p;
                j(1:n) = p;
                %
                % en <- Frobenius norm and max eval estimate for sprandsym
                %  Theory:  A is the sum of 2n terms with weight ratio^(-j)
                %  each has (density*n*n)/2n nonzeros with variance 1/3
                %  assuming all are separate nonzers, we calculate, by summing
                %  the series that the expected value of ||A||^2 is approximately
                %  en = (n*density/6) * (1-ratio^4n)/(1-ratio^2).  Further,
                %  we approximate max(eig(a)) by ||A||.  We note that
                %  the initial diagonal value is a very close approx to min(eig(a)).
                %  Thus, with the choice rcond * en we get the right cond(a) on
                %  average.
                %
                en = n * density / 6;
                en = en * (1 - ratio^(4*n)) / (1 - ratio^2);
                en = sqrt(en);
                %
                anz(1:n) = rcond * en * ones(n,1);
                nnza = n;   nnztrue = n;
                %
                %   Sum of scaled symmetric outer products with random structure
                %
                how_big = 1;
                for iouter = 1:2*n, 
                    nzpc = fix( 2 + sqrt((NZF*nnzwanted - nnza) / (2*n-iouter+1)));
                    bk = sprandn(n,1, nzpc/n);
                    [ik,jk,anzk] = find(how_big * (bk * 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;
                    nnztrue = nnztrue + (lk - nnz(bk));
                    how_big = how_big * ratio;
                    if (nnztrue > (nnzwanted*1.1)) break; end;
                end
                R = sparse(i(1:nnza),j(1:nnza),anz(1:nnza));
            end  %  kind == 2
        end   %  kind == 1 or 2
    end   %  nargin == 4
end      %  nargin > 2
if nargin>1,
    R = .5 * (R + R');    % make absolutely symmetric
end