Mercurial > octave-nkf

comparison scripts/sparse/private/__sprand__.m @ 20434:26bd6008fc9c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Rename __sprand_impl__.m to __sprand__.m * scripts/sparse/private/__sprand__.m: Renamed file. * scripts/sparse/private/__sprand_impl__.m: Removed file. * scripts/sparse/module.mk: Rename file in build system. * sprand.m, sprandn.m: Rename file in calling functions.
author Rik <rik@octave.org>
date Sun, 26 Jul 2015 08:50:04 -0700
parents scripts/sparse/private/__sprand_impl__.m@83792dd9bcc1
children
comparison
equal deleted inserted replaced
20432:1bc95d7148b7 20434:26bd6008fc9c
1 ## Copyright (C) 2004-2015 Paul Kienzle
2 ## Copyright (C) 2012 Jordi Gutiérrez Hermoso
3 ##
4 ## This file is part of Octave.
5 ##
6 ## Octave is free software; you can redistribute it and/or modify it
7 ## under the terms of the GNU General Public License as published by
8 ## the Free Software Foundation; either version 3 of the License, or (at
9 ## your option) any later version.
10 ##
11 ## Octave is distributed in the hope that it will be useful, but
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 ## General Public License for more details.
15 ##
16 ## You should have received a copy of the GNU General Public License
17 ## along with Octave; see the file COPYING. If not, see
18 ## <http://www.gnu.org/licenses/>.
19 ##
20 ## Original version by Paul Kienzle distributed as free software in the
21 ## public domain.
22
23 ## -*- texinfo -*-
24 ## @deftypefn {Function File} {} __sprand__ (@var{s}, @var{randfun})
25 ## @deftypefnx {Function File} {} __sprand__ (@var{m}, @var{n}, @var{d}, @var{fcnname}, @var{randfun})
26 ## @deftypefnx {Function File} {} __sprand__ (@var{m}, @var{n}, @var{d}, @var{rc}, @var{fcnname}, @var{randfun})
27 ## Undocumented internal function.
28 ## @end deftypefn
29
30 ## Actual implementation of sprand and sprandn happens here.
31
32 function S = __sprand__ (varargin)
33
34 if (nargin == 2)
35 [m, randfun] = deal (varargin{1:2});
36 [i, j] = find (m);
37 [nr, nc] = size (m);
38 S = sparse (i, j, randfun (size (i)), nr, nc);
39 else
40 if (nargin == 5)
41 [m, n, d, fcnname, randfun] = deal (varargin{:});
42 else
43 [m, n, d, rc, fcnname, randfun] = deal (varargin{:});
44 endif
45
46 if (! (isscalar (m) && m == fix (m) && m > 0))
47 error ("%s: M must be an integer greater than 0", fcnname);
48 endif
49 if (! (isscalar (n) && n == fix (n) && n > 0))
50 error ("%s: N must be an integer greater than 0", fcnname);
51 endif
52 if (d < 0 || d > 1)
53 error ("%s: density D must be between 0 and 1", fcnname);
54 endif
55
56 if (nargin == 5)
57 mn = m*n;
58 k = round (d*mn);
59 if (mn > sizemax ())
60 ## randperm will overflow, so use alternative methods
61
62 idx = unique (fix (rand (1.01*k, 1) * mn)) + 1;
63
64 ## idx contains random numbers in [1,mn]
65 ## Generate 1% more random values than necessary in order to reduce the
66 ## probability that there are less than k distinct values; maybe a
67 ## better strategy could be used but I don't think it's worth the price.
68
69 ## actual number of entries in S
70 k = min (length (idx), k);
71 j = floor ((idx(1:k) - 1) / m);
72 i = idx(1:k) - j * m;
73 j++;
74 else
75 idx = randperm (mn, k);
76 [i, j] = ind2sub ([m, n], idx);
77 endif
78
79 S = sparse (i, j, randfun (k, 1), m, n);
80
81 elseif (nargin == 6)
82 ## Create a matrix with specified reciprocal condition number.
83
84 if (! isscalar (rc) && ! isvector (rc))
85 error ("%s: RC must be a scalar or vector", fcnname);
86 endif
87
88 ## We want to reverse singular valued decomposition A=U*S*V'.
89 ## First, first S is constructed and then U = U1*U2*..Un and
90 ## V' = V1*V2*..Vn are seen as Jacobi rotation matrices with angles and
91 ## planes of rotation randomized. Repeatedly apply rotations until the
92 ## required density for A is achieved.
93
94 if (isscalar (rc))
95 if (rc < 0 || rc > 1)
96 error ("%s: reciprocal condition number RC must be between 0 and 1", fcnname);
97 endif
98 ## Reciprocal condition number is ratio of smallest SV to largest SV
99 ## Generate singular values randomly and sort them to build S
100 ## Random singular values in range [rc, 1].
101 v = rand (1, min (m,n)) * (1 - rc) + rc;
102 v(1) = 1;
103 v(end) = rc;
104 v = sort (v, "descend");
105 S = sparse (diag (v, m, n));
106 else
107 ## Only the min (m, n) greater singular values from rc vector are used.
108 if (length (rc) > min (m,n))
109 rc = rc(1:min(m, n));
110 endif
111 S = sparse (diag (sort (rc, "descend"), m, n));
112 endif
113
114 Uinit = speye (m);
115 Vinit = speye (n);
116 k = round (d*m*n);
117 while (nnz (S) < k)
118 if (m > 1)
119 ## Construct U randomized rotation matrix
120 rot_angleu = 2 * pi * rand ();
121 cu = cos (rot_angleu); su = sin (rot_angleu);
122 rndtmp = randperm (m, 2);
123 i = rndtmp(1); j = rndtmp(2);
124 U = Uinit;
125 U(i, i) = cu; U(i, j) = -su;
126 U(j, i) = su; U(j, j) = cu;
127 S = U * S;
128 endif
129 if (n > 1)
130 ## Construct V' randomized rotation matrix
131 rot_anglev = 2 * pi * rand ();
132 cv = cos (rot_anglev); sv = sin (rot_anglev);
133 rndtmp = randperm (n, 2);
134 i = rndtmp(1); j = rndtmp(2);
135 V = Vinit;
136 V(i, i) = cv; V(i, j) = sv;
137 V(j, i) = -sv; V(j, j) = cv;
138 S *= V;
139 endif
140 endwhile
141 endif
142 endif
143
144 endfunction
145