changeset 20169:6f8c572f27fe draft

gallery.m: repair functions randsvd(), pei(), kms(), hanowa(), gearmat(). * scripts/special-matrix/gallery.m: add qmult() auxiliary function called by randsvd(), edit incorrect argument test for gearmat(), rename incorrect variable assignments in pei(), kms(), hanowa().
author Antonio Pino Robles <data.script93@gmail.com>
date Thu, 28 May 2015 00:08:47 +0200
parents 7bd87990a8f4
children af2b7695f1c4
files scripts/special-matrix/gallery.m
diffstat 1 files changed, 67 insertions(+), 18 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/special-matrix/gallery.m	Wed May 27 22:38:13 2015 +0200
+++ b/scripts/special-matrix/gallery.m	Thu May 28 00:08:47 2015 +0200
@@ -157,9 +157,9 @@
 ## The second input is a matrix of dimensions describing the size of the output.
 ## The dimensions can also be input as comma-separated arguments.
 ##
-## The input @var{j} is an integer index in the range [0, 2^32-1].  The
-## values of the output matrix are always exactly the same
-## (reproducibility) for a given size input and @var{j} index.
+## The input @var{j} is an integer index in the range [0, 2^32-1].  The values
+## of the output matrix are always exactly the same (reproducibility) for a
+## given size input and @var{j} index.
 ##
 ## The final optional argument determines the class of the resulting matrix.
 ## Possible values for @var{class}: @qcode{"uint8"}, @qcode{"uint16"},
@@ -181,7 +181,7 @@
 ##
 ## @deftypefn  {Function File} {@var{a} =} gallery ("ipjfact", @var{n})
 ## @deftypefnx {Function File} {@var{a} =} gallery ("ipjfact", @var{n}, @var{k})
-## Create an Hankel matrix with factorial elements.
+## Create a Hankel matrix with factorial elements.
 ##
 ## @end deftypefn
 ##
@@ -257,9 +257,9 @@
 ## The first input is a matrix of dimensions describing the size of the output.
 ## The dimensions can also be input as comma-separated arguments.
 ##
-## The input @var{j} is an integer index in the range [0, 2^32-1].  The
-## values of the output matrix are always exactly the same
-## (reproducibility) for a given size input and @var{j} index.
+## The input @var{j} is an integer index in the range [0, 2^32-1].  The values
+## of the output matrix are always exactly the same (reproducibility) for a
+## given size input and @var{j} index.
 ##
 ## The final optional argument determines the class of the resulting matrix.
 ## Possible values for @var{class}: @qcode{"single"}, @qcode{"double"}.
@@ -380,9 +380,9 @@
 ## The first input is a matrix of dimensions describing the size of the output.
 ## The dimensions can also be input as comma-separated arguments.
 ##
-## The input @var{j} is an integer index in the range [0, 2^32-1].  The
-## values of the output matrix are always exactly the same
-## (reproducibility) for a given size input and @var{j} index.
+## The input @var{j} is an integer index in the range [0, 2^32-1].  The values
+## of the output matrix are always exactly the same (reproducibility) for a
+## given size input and @var{j} index.
 ##
 ## The final optional argument determines the class of the resulting matrix.
 ## Possible values for @var{class}: @qcode{"single"}, @qcode{"double"}.
@@ -1215,9 +1215,9 @@
     error ("gallery: 1 to 3 arguments are required for gearmat matrix.");
   elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
     error ("gallery: N must be an integer for gearmat matrix.");
-  elseif (! isnumeric (i) || ! isscalar (i) || i == 0 || abs (i) <= n)
+  elseif (! isnumeric (i) || ! isscalar (i) || i == 0 || abs (i) > n)
     error ("gallery: I must be a nonzero scalar, and abs (I) <= N for gearmat matrix.");
-  elseif (! isnumeric (j) || ! isscalar (j) || i == 0 || abs (j) <= n)
+  elseif (! isnumeric (j) || ! isscalar (j) || i == 0 || abs (j) > n)
     error ("gallery: J must be a nonzero scalar, and abs (J) <= N for gearmat matrix.");
   endif
 
@@ -1271,7 +1271,7 @@
     error ("gallery: N must be an integer for hanowa matrix.");
   elseif (rem (n, 2) != 0)
     error ("gallery: N must be even for hanowa matrix.");
-  elseif (! isnumeric (lambda) || ! isscalar (lambda))
+  elseif (! isnumeric (d) || ! isscalar (d))
     error ("gallery: D must be a numeric scalar for hanowa matrix.");
   endif
 
@@ -1602,11 +1602,11 @@
   ##    10 (1989), pp. 135-146 (and see the references therein).
 
   if (nargin < 1 || nargin > 2)
-    error ("gallery: 1 to 2 arguments are required for lauchli matrix.");
+    error ("gallery: 1 to 2 arguments are required for kms matrix.");
   elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
-    error ("gallery: N must be an integer for lauchli matrix.")
-  elseif (! isscalar (mu))
-    error ("gallery: MU must be a scalar for lauchli matrix.")
+    error ("gallery: N must be an integer for kms matrix.")
+  elseif (! isscalar (rho))
+    error ("gallery: rho must be a scalar for kms matrix.")
   endif
 
   A = (1:n)'*ones (1,n);
@@ -2015,7 +2015,7 @@
     error ("gallery: 1 to 2 arguments are required for pei matrix.");
   elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
     error ("gallery: N must be an integer for pei matrix.");
-  elseif (! isnumeric (w) || ! isscalar (w))
+  elseif (! isnumeric (alpha) || ! isscalar (alpha)) # change w to alpha 
     error ("gallery: ALPHA must be a scalar for pei matrix.");
   endif
 
@@ -2856,6 +2856,55 @@
 endfunction
 
 
+function B = qmult (A)
+  ## QMULT  Pre-multiply by random orthogonal matrix.
+  ##   QMULT(A) is Q*A where Q is a random real orthogonal matrix from
+  ##   the Haar distribution, of dimension the number of rows in A.
+  ##   Special case: if A is a scalar then QMULT(A) is the same as
+  ##   QMULT(EYE(A)).
+  ##
+  ##
+  ##   Called by RANDSVD.
+  ##       
+  ##   Reference:
+  ##   G.W. Stewart, The efficient generation of random
+  ##   orthogonal matrices with an application to condition estimators,
+  ##   SIAM J. Numer. Anal., 17 (1980), 403-409.
+  
+  [n, m] = size (A);
+  
+  #  Handle scalar A.
+  if (max (n,m) == 1)
+     n = A;
+     A = eye(n);
+  endif
+  
+  d = zeros (n);
+  
+  for k = n-1:-1:1
+    # Generate random Householder transformation.
+    x = randn (n-k+1,1);
+    s = norm (x);
+    sgn = sign (x(1)) + (x(1) == 0);    # Modification for sign(1)=1.
+    s = sgn * s;
+    d(k) = -sgn;
+    x(1) = x(1) + s;
+    beta = s * x(1);
+    # Apply the transformation to A.
+    y = x' * A(k:n,:);
+    A(k:n,:) = A(k:n,:) - x * (y/beta);
+  endfor
+  
+  # Tidy up signs.
+  for i=1:n-1
+      A(i,:) = d(i) * A(i,:);
+  endfor
+  A(n,:) = A(n,:) * sign (randn);
+  B = A;
+endfunction
+
+
+
 ## BIST testing for just a few functions to verify that the main gallery
 ## dispatch function works.
 %assert (gallery ("clement", 3), [0 1 0; 2 0 2; 0 1 0])