# HG changeset patch # User Antonio Pino Robles # Date 1432764527 -7200 # Node ID 6f8c572f27fe9a2ec231dc452c6e7ff7ca2c8ef3 # Parent 7bd87990a8f40c6db26c079559d5db9dd1f2046b 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(). diff -r 7bd87990a8f4 -r 6f8c572f27fe scripts/special-matrix/gallery.m --- 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])