Mercurial > octave-nkf

annotate scripts/sparse/private/__sprand_impl__.m @ 17712:e3ac73be7894

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
__sprand_impl__.m : Fix off by one error in the column index for large matrices (bug #40325)
author David Bateman <dbateman@free.fr>
date Mon, 21 Oct 2013 20:28:24 +0200
parents 1c89599167a6
children d63878346099
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
14138
72c96de7a403 maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents: 13929
diff changeset
1 ## Copyright (C) 2004-2012 Paul Kienzle
72c96de7a403 maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents: 13929
diff changeset
2 ## Copyright (C) 2012 Jordi Gutiérrez Hermoso
13197
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
3 ##
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
4 ## This file is part of Octave.
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
5 ##
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
6 ## Octave is free software; you can redistribute it and/or modify it
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
7 ## under the terms of the GNU General Public License as published by
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
8 ## the Free Software Foundation; either version 3 of the License, or (at
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
9 ## your option) any later version.
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
10 ##
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
11 ## Octave is distributed in the hope that it will be useful, but
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
14 ## General Public License for more details.
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
15 ##
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
16 ## You should have received a copy of the GNU General Public License
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
17 ## along with Octave; see the file COPYING. If not, see
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
18 ## <http://www.gnu.org/licenses/>.
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
19 ##
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
20 ## Original version by Paul Kienzle distributed as free software in the
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
21 ## public domain.
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
22
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
23 ## -*- texinfo -*-
13929
9cae456085c2 Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents: 13202
diff changeset
24 ## @deftypefn {Function File} {} __sprand_impl__ (@var{s}, @var{randfun})
13202
92095e0ba372 Fix the copyright and docstring of __sprand_impl__
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13197
diff changeset
25 ## @deftypefnx {Function File} {} __sprand_impl__ (@var{m}, @var{n}, @var{d}, @var{funname}, @var{randfun})
13197
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
26 ## Undocumented internal function.
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
27 ## @end deftypefn
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
28
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
29 ## Actual implementation of sprand and sprandn happens here.
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
30
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
31 function S = __sprand_impl__ (varargin)
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
32
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
33 if (nargin == 2)
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
34 m = varargin{1};
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
35 randfun = varargin{2};
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
36 [i, j] = find (m);
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
37 [nr, nc] = size (m);
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
38 S = sparse (i, j, randfun (size (i)), nr, nc);
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
39 return;
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
40 endif
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
41
14868
5d3a684236b0 maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents: 14138
diff changeset
42 [m, n, d, funname, randfun] = deal (varargin{:});
13197
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
43
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
44 if (!(isscalar (m) && m == fix (m) && m > 0))
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
45 error ("%s: M must be an integer greater than 0", funname);
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
46 endif
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
47
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
48 if (!(isscalar (n) && n == fix (n) && n > 0))
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
49 error ("%s: N must be an integer greater than 0", funname);
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
50 endif
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
51
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
52 if (d < 0 || d > 1)
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
53 error ("%s: density D must be between 0 and 1", funname);
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
54 endif
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
55
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
56 mn = m*n;
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
57 k = round (d*mn);
16779
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
58 if (mn > sizemax ())
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
59 ## randperm will overflow, so use alternative methods
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
60
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
61 idx = unique (fix (rand (min (k*1.01, k+10), 1) * mn)) + 1;
13197
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
62
16779
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
63 ## idx contains random numbers in [1,mn]
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
64 ## generate 1% or 10 more random values than necessary in order to
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
65 ## reduce the probability that there are less than k distinct
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
66 ## values; maybe a better strategy could be used but I don't think
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
67 ## it's worth the price
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
68
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
69 ## actual number of entries in S
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
70 k = min (length (idx), k);
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
71 j = floor ((idx(1:k) - 1) / m);
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
72 i = idx(1:k) - j * m;
17712
e3ac73be7894 __sprand_impl__.m : Fix off by one error in the column index for large matrices (bug #40325)
David Bateman <dbateman@free.fr>
parents: 17338
diff changeset
73 j++;
16779
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
74 else
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
75 idx = randperm (mn, k);
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
76 [i, j] = ind2sub ([m, n], idx);
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
77 endif
8fce0ed4894a Specialize is_empty and numel methods for sparse matrices (debian bug #706376)
David Bateman <dbateman@free.fr>
parents: 14868
diff changeset
78
13197
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
79 S = sparse (i, j, randfun (k, 1), m, n);
6db186dfdeaa Refactor sprandn/sprand code, move common code to common function (bug #34352)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
80
14868
5d3a684236b0 maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents: 14138
diff changeset
81 endfunction
17338
1c89599167a6 maint: End m-files with 1 blank line.
Rik <rik@octave.org>
parents: 17176
diff changeset
82