home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPARFUN.DI$ / SPRANDN.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  4.2 KB  |  128 lines
  1. function R = sprandn(arg1,n,density,rc)
  2. %SPRAND    Random sparse matrices.
  3. %    R = SPRANDN(S) has the same sparsity structure as S, but normally
  4. %        distributed random entries.
  5. %    R = SPRANDN(m,n,density) is a random, m-by-n, sparse matrix with 
  6. %        approximately density*m*n normally distributed nonzero entries.
  7. %    R = SPRANDN(m,n,density,rc) also has reciprocal condition number
  8. %        approximately equal to rc.  R is constructed from a sum of
  9. %        matrices of rank one. 
  10. %
  11. %        If rc is a vector of length lr <= min(m,n), then R has 
  12. %        rc as its first lr singular values, all others are zero.
  13. %        In this case, R is generated by random plane rotations
  14. %        applied to a diagonal matrix with the given singular values
  15. %        It has a great deal of topological and algebraic structure.
  16. %
  17. %    See also SPRANDSYM.
  18.  
  19. %    Rob Schreiber, RIACS, and Cleve Moler, MathWorks, 9/10/90.
  20. %    Revised 1/28/91, 2/12/91, RS; 8/12/91, CBM.
  21. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  22.  
  23.  
  24. if nargin == 1
  25.    S = arg1;
  26.    [m,n] = size(S);
  27.    [i,j] = find(S);  i = i(:);  %workaround for bug in find
  28.    R = sparse(i,j,randn(length(i),1),m,n);
  29. elseif nargin == 2
  30.    error('Too many or not enough input arguments.')
  31. elseif nargin == 3
  32.    m = arg1;
  33.    nnzwanted = round(m * n * density);
  34.    rands = randn( nnzwanted, 1 );
  35.    i = fix( rand(nnzwanted, 1) * m ) + 1;
  36.    j = fix( rand(nnzwanted, 1) * n ) + 1;
  37.    R = sparse(i,j,rands,m,n);
  38. else    %    nargin == 4
  39.    m = arg1;
  40.    nnzwanted = round(density*m*n);
  41.    minm = min(m,n);
  42.    maxmn = max(m,n);
  43.    lr = length(rc);
  44.    if (lr > 1)  %  specified singular values
  45.       if (lr > minm) lr = minm; end;
  46.       %
  47.       %   To start, put a random nonzero in every column / row if
  48.       %   there is room enough
  49.       %
  50.       anz = zeros(maxmn, 1);
  51.       anz(1:lr) = rc(1:lr);
  52.       lr = minm;      %   if lr < minm then this adds some zero s.v. 
  53.       sqrt2o2 = sqrt(2)/2;
  54.       while (lr < min(maxmn, nnzwanted))
  55.          ndo = min([lr, nnzwanted-lr, maxmn-lr]);
  56.          anz(lr+1 : lr+ndo) = sqrt2o2 * anz(1:ndo);
  57.          anz(1 : ndo) = sqrt2o2 * anz(1:ndo);
  58.          lr = lr + ndo;
  59.       end
  60.       [ignore,p] = sort(rand(m,1));
  61.       [ignore,q] = sort(rand(n,1));
  62.       if (m < n),
  63.          p = p(1 + rem((0:n-1), m));
  64.       else
  65.          q = q(1 + rem((0:m-1), n));
  66.       end
  67.       R = sparse(p, q, anz, m, n);
  68.       %
  69.       %   random kogbetliantz rotations
  70.       %
  71.       while ( nnz(R) < .95*nnzwanted )
  72.           R = rkr(R);
  73.       end
  74.    else    %   condition number specified
  75.       if rc > 1, error('Reciprocal condition numver must be less than 1.'); end;
  76.       MEMFAC = 1.05;
  77.       NZF = 1.25;
  78.       i = zeros(fix(nnzwanted*MEMFAC),1); j = zeros(fix(nnzwanted*MEMFAC),1); 
  79.       anz = zeros(fix(nnzwanted*MEMFAC),1); nnza = 0;
  80.       %
  81.       %   To start, put a random nonzero in every column / row
  82.       %
  83.       how_big = 1;
  84.       [ignore,p] = sort(rand(m,1));
  85.       [ignore,q] = sort(rand(n,1));
  86.       ratio = rc ^ (1/(minm-1));
  87.       if (m < n),
  88.          for jj = 1:n
  89.             i(jj) = p(1 + rem(jj-1,m));
  90.             j(jj) = q(jj);
  91.             anz(jj) = how_big;
  92.             how_big = how_big * ratio;
  93.          end
  94.          nnza = n;
  95.       else
  96.          for ii = 1:m
  97.             j(ii) = q(1 + rem(ii-1,n));
  98.             i(ii) = p(ii);
  99.             anz(ii) = how_big;
  100.             how_big = how_big * ratio;
  101.          end
  102.          nnza = m;
  103.       end
  104.       %
  105.       %   Now add a sum of scaled outer products with random structure
  106.       %
  107.       how_big = .1;
  108.       for iouter = 1:minm,
  109.          nzpc = fix(1 + sqrt((m/n)*(NZF*nnzwanted - nnza)/(minm-iouter+1)));
  110.          nzpr = fix(1 + sqrt((n/m)*(NZF*nnzwanted - nnza)/(minm-iouter+1)));
  111.          ak = sprandn(m,1, nzpc/m);
  112.          bk = sprandn(n,1, nzpr/n);
  113.          [ik,jk,anzk] = find(how_big * (ak * bk'));
  114.          lk = length(ik);
  115.          i(nnza+1:nnza+lk) = ik; 
  116.          j(nnza+1:nnza+lk) = jk; 
  117.          anz(nnza+1:nnza+lk) = anzk;
  118.          nnza = nnza + lk;
  119.          how_big = how_big * ratio;
  120.          if (nnza > nnzwanted) break; end;
  121.       end
  122.       i = i(1:nnza);
  123.       j = j(1:nnza);
  124.       anz = anz(1:nnza);
  125.       R = sparse(i,j,anz);
  126.    end  %   rc gives condition number
  127. end     %   nargin == 4
  128.