# HG changeset patch
# User Eduardo Ramos (edu159) <eduradical951@gmail.com>
# Date 1395491021 -3600
# Node ID 35a5e7740a6dcd6f707551d27fc2763a7e960924
# Parent  5032ac119d52001a26e3bdf28abaf49b1d9a7f63
Added implementation for 4th argument of sprand/sprandn (bug #41839).

* __sprand_impl__.m: Implementation done here

* sprand.m: Added documentation and some tests.

* sprandn.m: Added documentation and some tests.

diff -r 5032ac119d52 -r 35a5e7740a6d scripts/sparse/private/__sprand_impl__.m
--- a/scripts/sparse/private/__sprand_impl__.m	Sun Mar 23 12:41:15 2014 +0100
+++ b/scripts/sparse/private/__sprand_impl__.m	Sat Mar 22 13:23:41 2014 +0100
@@ -23,6 +23,7 @@
 ## -*- texinfo -*-
 ## @deftypefn  {Function File} {} __sprand_impl__ (@var{s}, @var{randfun})
 ## @deftypefnx {Function File} {} __sprand_impl__ (@var{m}, @var{n}, @var{d}, @var{funname}, @var{randfun})
+## @deftypefnx {Function File} {} __sprand_impl__ (@var{m}, @var{n}, @var{d}, @var{rc}, @var{funname}, @var{randfun})
 ## Undocumented internal function.
 ## @end deftypefn
 
@@ -37,46 +38,104 @@
     [nr, nc] = size (m);
     S = sparse (i, j, randfun (size (i)), nr, nc);
     return;
-  endif
-
-  [m, n, d, funname, randfun] = deal (varargin{:});
+  else
+    if (nargin == 5)
+      [m, n, d, funname, randfun] = deal(varargin{:});
+    else 
+      [m, n, d, rc, funname, randfun] = deal(varargin{:});
+    endif
 
-  if (!(isscalar (m) && m == fix (m) && m > 0))
-    error ("%s: M must be an integer greater than 0", funname);
-  endif
+    if (!(isscalar (m) && m == fix (m) && m > 0))
+      error ("%s: M must be an integer greater than 0", funname);
+    endif
+
+    if (!(isscalar (n) && n == fix (n) && n > 0))
+      error ("%s: N must be an integer greater than 0", funname);
+    endif
+
+    if (d < 0 || d > 1)
+      error ("%s: density D must be between 0 and 1", funname);
+    endif
+
 
-  if (!(isscalar (n) && n == fix (n) && n > 0))
-    error ("%s: N must be an integer greater than 0", funname);
-  endif
+    if (nargin == 5)
+      mn = m*n;
+      k = round (d*mn);
+      if (mn > sizemax ())
+        ## randperm will overflow, so use alternative methods
+
+        idx = unique (fix (rand (min (k*1.01, k+10), 1) * mn)) + 1;
 
-  if (d < 0 || d > 1)
-    error ("%s: density D must be between 0 and 1", funname);
-  endif
+        ## idx contains random numbers in [1,mn]
+        ## generate 1% or 10 more random values than necessary in order to
+        ## reduce the probability that there are less than k distinct
+        ## values; maybe a better strategy could be used but I don't think
+        ## it's worth the price
+        
+        ## actual number of entries in S
+        k = min (length (idx), k);
+        j = floor ((idx(1:k) - 1) / m);
+        i = idx(1:k) - j * m;
+        j++;
+      else
+        idx = randperm (mn, k);
+        [i, j] = ind2sub ([m, n], idx);
+      endif
 
-  mn = m*n;
-  k = round (d*mn);
-  if (mn > sizemax ())
-    ## randperm will overflow, so use alternative methods
+      S = sparse (i, j, randfun (k, 1), m, n);
+    elseif (nargin == 6)
+      ## We assume that we want to reverse A=U*S*V' so firstly S is constructed
+      ## and then U = U1*U2*..Un and V' = V1*V2*..Vn  are seen as Jacobi rotation matrices with angles and
+      ## planes of rotation randomized. In the nth step the density required for A is achieved.
 
-    idx = unique (fix (rand (min (k*1.01, k+10), 1) * mn)) + 1;
-
-    ## idx contains random numbers in [1,mn]
-    ## generate 1% or 10 more random values than necessary in order to
-    ## reduce the probability that there are less than k distinct
-    ## values; maybe a better strategy could be used but I don't think
-    ## it's worth the price
-    
-    ## actual number of entries in S
-    k = min (length (idx), k);
-    j = floor ((idx(1:k) - 1) / m);
-    i = idx(1:k) - j * m;
-    j++;
-  else
-    idx = randperm (mn, k);
-    [i, j] = ind2sub ([m, n], idx);
+      mynnz = round (m * n * d);
+      if (!isscalar(rc)) 
+        ## Only the min(m, n) greater singular values from rc vector are used. Needed to be compliant.
+        if (length (rc) > min (m,n))
+          rc = rc(1:min(m, n));
+        endif
+        S = sparse (diag (sort (rc, 'descend'), m, n));
+      else
+        if(rc < 0 || rc > 1)
+          error ("%s: reciprocal condition number rc must be between 0 and 1", funname);
+        endif
+        ## Generate the singular values randomly and sort them to build S
+        for (i = 1:min(m, n))
+          ## Randon singular values between 1 and rc.
+          v(i) = rand () * (1 - rc) + rc; 
+        endfor
+        v(1) = 1;
+        v(end) = rc;
+        v = sort (v, 'descend');
+        S = sparse (diag (v, m, n));
+      endif
+      while (nnz(S) < mynnz) 
+        [mm, nn] = size(S);
+        rot_angleu = 2 * randfun () * pi;
+        rot_anglev = 2 * randfun () * pi;
+        cu = cos (rot_angleu); cv = cos (rot_anglev); 
+        su = sin (rot_angleu); sv = sin (rot_anglev);
+        ## Rotation related with U
+        i = fix (rand () * m) + 1;
+        do
+          ## If j==i rotation matrix would be no longer that kind
+          j = fix (rand () * m) + 1;   
+        until (j != i)
+        U = sparse (eye (m,m));
+        U(i, i) = cu; U(i, j) = -su;
+        U(j, i) = su; U(j, j) = cu;
+        S = U * S;
+        ## Rotation related with V'
+        i = fix (rand () * nn) + 1;
+        do
+          j = fix (rand () * nn) + 1; 
+        until(j != i)
+        V = sparse (eye (n, n));
+        V(i, i) = cv; V(i, j) = sv;
+        V(j, i) = -sv; V(j, j) = cv;
+        S = S * V;
+      endwhile
+    endif
   endif
-
-  S = sparse (i, j, randfun (k, 1), m, n);
-
 endfunction
 
diff -r 5032ac119d52 -r 35a5e7740a6d scripts/sparse/sprand.m
--- a/scripts/sparse/sprand.m	Sun Mar 23 12:41:15 2014 +0100
+++ b/scripts/sparse/sprand.m	Sat Mar 22 13:23:41 2014 +0100
@@ -21,6 +21,7 @@
 
 ## -*- texinfo -*-
 ## @deftypefn  {Function File} {} sprand (@var{m}, @var{n}, @var{d})
+## @deftypefnx  {Function File} {} sprand (@var{m}, @var{n}, @var{d}, @var{rc})
 ## @deftypefnx {Function File} {} sprand (@var{s})
 ## Generate a random sparse matrix.  The size of the matrix will be
 ## @var{m} by @var{n}, with a density of values given by @var{d}.
@@ -29,6 +30,11 @@
 ##
 ## If called with a single matrix argument, a random sparse matrix is
 ## generated wherever the matrix @var{S} is non-zero.
+##
+## If called with the rc parameter as a scalar, a random sparse matrix with a
+## reciprocal condition number rc is generated. If rc is a vector, then it 
+## contains the first singular values of the matrix generated (length(rc) <= min(m, n)).
+##
 ## @seealso{sprandn, sprandsym}
 ## @end deftypefn
 
@@ -42,12 +48,14 @@
 ## David Bateman
 ##      2004-10-20      Texinfo help and copyright message
 
-function S = sprand (m, n, d)
+function S = sprand (m, n, d, rc)
 
   if (nargin == 1 )
     S = __sprand_impl__ (m, @rand);
   elseif ( nargin == 3)
     S = __sprand_impl__ (m, n, d, "sprand", @rand);
+  elseif (nargin == 4)
+    S = __sprand_impl__ (m, n, d, rc, "sprand", @rand);
   else
     print_usage ();
   endif
@@ -68,6 +76,19 @@
 %! assert (sort (j), [2 3 3]');
 %! assert (all (v > 0 & v < 1));
 
+%% Test 4-input calling form
+%!test
+%! d = rand();
+%! s1 = sprand (100, 100, d, 0.4);
+%! rc = [5, 4, 3, 2, 1, 0.1];
+%! s2 = sprand (100, 100, d, rc);
+%! s3 = sprand (6, 4, d, rc);
+%! assert (svd (s2)'(1:length (rc)), rc, sqrt (eps)); 
+%! assert (1/cond (s1), 0.4, sqrt (eps));
+%! assert (nnz (s1) / (100*100), d, 0.02);
+%! assert (nnz (s2) / (100*100), d, 0.02); 
+%! assert (svd (s3)', [5 4 3 2], sqrt (eps));
+
 %% Test input validation
 %!error sprand ()
 %!error sprand (1, 2)
@@ -80,6 +101,8 @@
 %!error sprand (3, 0, 0.5)
 %!error sprand (3, 3, -1)
 %!error sprand (3, 3, 2)
+%!error sprand (2, 2, 0.2, -1)
+%!error sprand (2, 2, 0.2, 2)
 
 %% Test very large, very low density matrix doesn't fail 
 %!test
diff -r 5032ac119d52 -r 35a5e7740a6d scripts/sparse/sprandn.m
--- a/scripts/sparse/sprandn.m	Sun Mar 23 12:41:15 2014 +0100
+++ b/scripts/sparse/sprandn.m	Sat Mar 22 13:23:41 2014 +0100
@@ -21,6 +21,7 @@
 
 ## -*- texinfo -*-
 ## @deftypefn  {Function File} {} sprandn (@var{m}, @var{n}, @var{d})
+## @deftypefnx {Function File} {} sprandn (@var{m}, @var{n}, @var{d}, @var{rc})
 ## @deftypefnx {Function File} {} sprandn (@var{s})
 ## Generate a random sparse matrix.  The size of the matrix will be
 ## @var{m} by @var{n}, with a density of values given by @var{d}.
@@ -29,17 +30,25 @@
 ##
 ## If called with a single matrix argument, a random sparse matrix is
 ## generated wherever the matrix @var{S} is non-zero.
+## 
+## If called with the rc parameter, then a random sparse matrix with a
+## reciprocal condition number rc is generated if rc is scalar. If rc 
+## is a vector, then it contains the first singular values of the matrix
+## generated (length(rc) <= min(m, n)).
+## 
 ## @seealso{sprand, sprandsym}
 ## @end deftypefn
 
 ## Author: Paul Kienzle <pkienzle@users.sf.net>
 
-function S = sprandn (m, n, d)
+function S = sprandn (m, n, d, rc)
 
   if (nargin == 1 )
     S = __sprand_impl__ (m, @randn);
   elseif ( nargin == 3)
     S = __sprand_impl__ (m, n, d, "sprandn", @randn);
+  elseif (nargin == 4)
+    S = __sprand_impl__ (m, n, d, rc, "sprandn", @randn);
   else
     print_usage ();
   endif
@@ -59,6 +68,19 @@
 %! assert (sort (i), [1 2 3]');
 %! assert (sort (j), [2 3 3]');
 
+%% Test 4-input calling form
+%!test
+%! d = rand();
+%! s1 = sprandn (100, 100, d, 0.4);
+%! rc = [5, 4, 3, 2, 1, 0.1];
+%! s2 = sprandn (100, 100, d, rc);
+%! s3 = sprandn (6, 4, d, rc);
+%! assert (svd (s2)'(1:length (rc)), rc, sqrt (eps)); 
+%! assert (1/cond (s1), 0.4, sqrt (eps));
+%! assert (nnz (s1) / (100*100), d, 0.02); 
+%! assert (nnz (s2) / (100*100), d, 0.02); 
+%! assert (svd (s3)', [5 4 3 2], sqrt (eps));
+
 %% Test input validation
 %!error sprandn ()
 %!error sprandn (1, 2)
@@ -71,6 +93,8 @@
 %!error sprandn (3, 0, 0.5)
 %!error sprandn (3, 3, -1)
 %!error sprandn (3, 3, 2)
+%!error sprandn (2, 2, 0.2, -1)
+%!error sprandn (2, 2, 0.2, 2)
 
 %% Test very large, very low density matrix doesn't fail 
 %!test