Mercurial > hg > octave-nkf

comparison scripts/sparse/sprandsym.m @ 6252:738c97e101eb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-01-24 19:58:46 by jwe]
author jwe
date Wed, 24 Jan 2007 19:58:46 +0000
parents 34f96dd5441b
children 2c85044aa63f
comparison
equal deleted inserted replaced
6251:27f85707054a 6252:738c97e101eb
30 ## generated wherever the matrix @var{S} is non-zero in its lower 30 ## generated wherever the matrix @var{S} is non-zero in its lower
31 ## triangular part. 31 ## triangular part.
32 ## @seealso{sprand, sprandn} 32 ## @seealso{sprand, sprandn}
33 ## @end deftypefn 33 ## @end deftypefn
34 34
35 function S = sprandsym(n,d) 35 function S = sprandsym (n, d)
36 if nargin == 1 36 if (nargin == 1)
37 [i,j,v,nr,nc] = spfind(tril(n)); 37 [i, j, v, nr, nc] = spfind (tril (n));
38 S = sparse(i,j,randn(size(v)),nr,nc); 38 S = sparse (i, j, randn (size (v)), nr, nc);
39 S = S + tril(S,-1)'; 39 S = S + tril (S, -1)';
40 elseif nargin == 2 40 elseif (nargin == 2)
41 m1 = floor(n/2); 41 m1 = floor (n/2);
42 n1 = m1 + 1; 42 n1 = m1 + rem (n, 2);
43 mn1 = m1*n1; 43 mn1 = m1*n1;
44 k1 = round(d*mn1); 44 k1 = round (d*mn1);
45 idx1=unique(fix(rand(min(k1*1.01,k1+10),1)*mn1))+1; 45 idx1 = unique (fix(rand(min(k1*1.01,k1+10),1)*mn1))+1;
46 # idx contains random numbers in [1,mn] 46 ## idx contains random numbers in [1,mn] generate 1% or 10 more
47 # generate 1% or 10 more random values than necessary 47 ## random values than necessary in order to reduce the probability
48 # in order to reduce the probability that there are less than k 48 ## that there are less than k distinct values; maybe a better
49 # distinct values; 49 ## strategy could be used but I don't think it's worth the price.
50 # maybe a better strategy could be used 50 k1 = min (length (idx1), k1); # actual number of entries in S
51 # but I don't think it's worth the price 51 j1 = floor ((idx1(1:k1)-1)/m1);
52 k1 = min(length(idx1),k1); # actual number of entries in S
53 j1 = floor((idx1(1:k1)-1)/m1);
54 i1 = idx1(1:k1) - j1*m1; 52 i1 = idx1(1:k1) - j1*m1;
55 53
56 n2 = ceil(n/2); 54 n2 = ceil (n/2);
57 nn2 = n2*n2; 55 nn2 = n2*n2;
58 k2 = round(d*nn2); 56 k2 = round (d*nn2);
59 idx2=unique(fix(rand(min(k2*1.01,k1+10),1)*nn2))+1; 57 idx2 = unique (fix (rand (min (k2*1.01, k1+10), 1)*nn2)) + 1;
60 k2 = min(length(idx2),k2); 58 k2 = min (length (idx2), k2);
61 j2 = floor((idx2(1:k2)-1)/n2); 59 j2 = floor ((idx2(1:k2)-1)/n2);
62 i2 = idx2(1:k2) - j2*n2; 60 i2 = idx2(1:k2) - j2*n2;
63 61
64 if isempty(i1) && isempty(i2) 62 if (isempty (i1) && isempty (i2))
65 S = sparse(n,n); 63 S = sparse (n, n);
66 else 64 else
67 S1 = sparse(i1,j1+1,randn(k1,1),m1,n1); 65 S1 = sparse (i1, j1+1, randn (k1, 1), m1, n1);
68 S = [tril(S1), sparse(m1,m1); ... 66 S = [tril(S1), sparse(m1,m1); ...
69 sparse(i2,j2+1,randn(k2,1),n2,n2), triu(S1,1)']; 67 sparse(i2,j2+1,randn(k2,1),n2,n2), triu(S1,1)'];
70 S = S + tril(S,-1)'; 68 S = S + tril (S, -1)';
71 endif 69 endif
72 else 70 else
73 print_usage (); 71 print_usage ();
74 endif 72 endif
75 endfunction 73 endfunction