home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPARFUN.DI$ / SPRANDSY.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  7.1 KB  |  209 lines
  1. function R = sprandsym(arg1, density, rcond, kind)
  2. %SPRANDSYM Random sparse symmetric matrices
  3. %
  4. %    R = SPRANDSYM(S) is a symmetric random matrix whose lower triangle
  5. %        and diagonal have the same structure as S.  The elements are
  6. %        normally distributed, with mean 0 and variance 1.
  7. %
  8. %    R = SPRANDSYM(n,density) is a symmetric random, n-by-n, sparse 
  9. %        matrix with approximately density*n*n nonzeros; each entry is
  10. %        the sum of one or more normally distributed random samples.
  11. %
  12. %    R = SPRANDSYM(n,density,rcond) also has a reciprocal condition number
  13. %        equal to rcond.  The distribution of entries is
  14. %        nonuniform; it is roughly symmetric about 0; all are in [-1,1].
  15. %
  16. %        If rcond is a VECTOR of length n, then R has eigenvalues rcond.
  17. %        Thus, if rcond is a positive (nonnegative) vector then R will 
  18. %        be positive (nonnegative) definite
  19. %
  20. %        In either case, R is generated by random Jacobi rotations
  21. %        applied to a diagonal matrix with the given eigenvalues or
  22. %        condition number.  It has a great deal of topological and
  23. %        algebraic structure.
  24. %
  25. %    R = SPRANDSYM(n, density, rcond, kind) is positive definite.
  26. %
  27. %        If kind = 1, R is generated by random Jacobi rotation of
  28. %           a positive definite diagonal matrix.
  29. %           R has the desired condition number exactly.
  30. %
  31. %        If kind = 2, R is a shifted sum of outer products.
  32. %           R has the desired condition number only
  33. %           approximately, but has less structure.
  34. %
  35. %        If R = SPRANDSYM(S, density, rcond, 3), 
  36. %           then R has the same structure as the MATRIX S and
  37. %           approximate condition number 1/rcond.  density is ignored.
  38. %
  39. %    See also SPRANDN.
  40.  
  41. %    Rob Schreiber, RIACS, and Cleve Moler, MathWorks, 2/5/91.
  42. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  43.  
  44. if nargin > 4
  45.   help sprandsym
  46.   error('Too many or not enough input arguments.')
  47. elseif nargin == 1
  48.    S = arg1;
  49.    [m,n] = size(S);
  50.    [ii,jj] = find(S);  ii = ii(:);  %workaround find bug
  51.    lower = find(ii > jj);  %  find lower triangle
  52.    di = find( ii == jj );  %  find diagonal
  53.    nd = length( di );
  54.    nl = length( lower );
  55.    randi = randn( nd, 1 );
  56.    randl = randn( nl, 1 );
  57.    i = [ii(lower); jj(lower); ii(di)];
  58.    j = [jj(lower); ii(lower); ii(di)];
  59.    r = [randl; randl; randi];
  60.    R = sparse(i,j,r,m,m);
  61. elseif nargin == 2,
  62.    n = arg1;
  63.    nl = round( (n*(n+1)/2) * density );
  64.    rands = randn( nl, 1 );
  65.    ii = fix( rand(nl, 1) * n ) + 1;
  66.    jj = fix( rand(nl, 1) * n ) + 1;
  67.    di = find( ii == jj );
  68.    off = find( ii ~= jj );
  69.    nd = length( di );
  70.    no = length( off );
  71.    randi = rands( 1:nd );
  72.    rando = rands( nd+1 : nl );
  73.    i = [ii(off); jj(off); ii(di)];
  74.    j = [jj(off); ii(off); ii(di)];
  75.    r = [rando; rando; randi];
  76.    R = sparse(i,j,r,n,n);
  77. else % nargin > 2
  78.    lrc = length(rcond);
  79.    if (lrc == 1)
  80.       if( (1/rcond) < 1 )
  81.          error('Reciprocal condition numbver must be less than 1.')
  82.       end
  83.    end
  84.    if(nargin == 3 )
  85.       n = arg1;
  86.       nnzwanted = round(density*n*n);
  87.       if (lrc == 1)   %   rcond is given
  88.          ratio = ((-1)^nargin) * (rcond ^ (1/(n-1)));
  89.          anz = (ratio .^ ( 0:(n-1)));
  90.       else            % eigenvalues explicitly given
  91.          anz = rcond;
  92.       end
  93.       R = sparse(1:n,1:n, anz);
  94.       % for nargin == 3,  ratio < 0
  95.       % so we start with an indefinite R;
  96.       nnzr = n;
  97.       %
  98.       %   random jacobi rotations
  99.       %
  100.       while ( nnzr < .95*nnzwanted )
  101.           R = rjr(R);
  102.           nnzr = nnz(R);
  103.       end
  104.    else  % nargin == 4
  105.     if(kind==3)
  106. %
  107. %  add a 2-by-2 rank-1 for each off diagonal
  108. %
  109.       [ii,jj] = find(arg1); ii = ii(:);  % workaround
  110.       [m, n] = size(arg1);
  111.       lower = find(ii > jj);
  112.       nl = length(lower);
  113.       ratio = rcond ^ (1/(nl-1));
  114.       ii = ii(lower);
  115.       jj = jj(lower);
  116.       i = zeros(n + 4*nl,1); 
  117.       j = zeros(n + 4*nl,1); 
  118.       anz = zeros(n + 4*nl,1); 
  119.       [ignore,p] = sort(rand(n,1));
  120.       i(1:n) = p;
  121.       j(1:n) = p;
  122.       en = sqrt((28/45) * (1 - ratio^(2*nl)) / (1 - ratio^2));
  123.       anz(1:n) = rcond * en * ones(n,1);
  124.       nnza = n;
  125.       how_big = 1;
  126.       for iouter = 1:nl
  127.          rv = rand(2,1);
  128.          op = how_big * rv * rv';
  129.          op = op(:);
  130.          ik = ii(iouter);
  131.          jk = jj(iouter);
  132.          i(nnza+1:nnza+4) = [ik; jk; ik; jk];
  133.          j(nnza+1:nnza+4) = [ik; ik; jk; jk];
  134.          anz(nnza+1:nnza+4) = op;
  135.          nnza = nnza + 4;
  136.          how_big = how_big * ratio;
  137.       end
  138.       R = sparse(i(1:nnza),j(1:nnza),anz(1:nnza));
  139.     %  end of kind == 3
  140.     else  %  kind == 1 or 2
  141.       n = arg1;  [mn,nn]=size(n);
  142.       if (mn > 1 | nn > 1)
  143.          error('First argument must be a scalar');
  144.       end
  145.       nnzwanted = round(density*n*n);
  146.       ratio = ((-1)^nargin) * (rcond ^ (1/(n-1)));
  147.       if (kind == 1)
  148.         R = sparse(1:n,1:n, (ratio .^ ( 0:(n-1)))  );
  149.         nnzr = n;
  150.         %   random jacobi rotations
  151.         while ( nnzr < .95*nnzwanted )
  152.             R = rjr(R);
  153.             nnzr = nnz(R);
  154.         end
  155.       elseif(kind==2)   %  positive definite via sum of outer products
  156.         MEMFAC = 2;
  157.         NZF = 1.25;
  158.         i = zeros(fix(nnzwanted*MEMFAC),1); 
  159.         j = zeros(fix(nnzwanted*MEMFAC),1); 
  160.         anz = zeros(fix(nnzwanted*MEMFAC),1); 
  161.         nnza = 0;
  162.         %
  163.         %   To start, put a random nonzero in every column / row
  164.         %
  165.         [ignore,p] = sort(rand(n,1));
  166.         i(1:n) = p;
  167.         j(1:n) = p;
  168.         %
  169.         % en <- Frobenius norm and max eval estimate for sprandsym
  170.         %  Theory:  A is the sum of 2n terms with weight ratio^(-j)
  171.         %  each has (density*n*n)/2n nonzeros with variance 1/3
  172.         %  assuming all are separate nonzers, we calculate, by summing
  173.         %  the series that the expected value of ||A||^2 is approximately
  174.         %  en = (n*density/6) * (1-ratio^4n)/(1-ratio^2).  Further,
  175.         %  we approximate max(eig(a)) by ||A||.  We note that
  176.         %  the initial diagonal value is a very close approx to min(eig(a)).
  177.         %  Thus, with the choice rcond * en we get the right cond(a) on
  178.         %  average.
  179.         %
  180.         en = n * density / 6;
  181.         en = en * (1 - ratio^(4*n)) / (1 - ratio^2);
  182.         en = sqrt(en);
  183.         %
  184.         anz(1:n) = rcond * en * ones(n,1);
  185.         nnza = n;   nnztrue = n;
  186.         %
  187.         %   Sum of scaled symmetric outer products with random structure
  188.         %
  189.         how_big = 1;
  190.         for iouter = 1:2*n, 
  191.            nzpc = fix( 2 + sqrt((NZF*nnzwanted - nnza) / (2*n-iouter+1)));
  192.            bk = sprandn(n,1, nzpc/n);
  193.            [ik,jk,anzk] = find(how_big * (bk * bk'));
  194.            lk = length(ik);
  195.            i(nnza+1:nnza+lk) = ik; 
  196.            j(nnza+1:nnza+lk) = jk; 
  197.            anz(nnza+1:nnza+lk) = anzk;
  198.            nnza = nnza + lk;
  199.            nnztrue = nnztrue + (lk - nnz(bk));
  200.            how_big = how_big * ratio;
  201.            if (nnztrue > (nnzwanted*1.1)) break; end;
  202.         end
  203.         R = sparse(i(1:nnza),j(1:nnza),anz(1:nnza));
  204.       end  %  kind == 2
  205.     end   %  kind == 1 or 2
  206.    end   %  nargin == 4
  207.    R = .5 * (R + R');    % make absolutely symmetric
  208. end      %  nargin > 2
  209.