Mercurial > octave-antonio

comparison scripts/special-matrix/gallery.m @ 20169:6f8c572f27fe draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
20168:7bd87990a8f4 20169:6f8c572f27fe
155 ## [@var{imin}, @var{imax}]. 155 ## [@var{imin}, @var{imax}].
156 ## 156 ##
157 ## The second input is a matrix of dimensions describing the size of the output. 157 ## The second input is a matrix of dimensions describing the size of the output.
158 ## The dimensions can also be input as comma-separated arguments. 158 ## The dimensions can also be input as comma-separated arguments.
159 ## 159 ##
160 ## The input @var{j} is an integer index in the range [0, 2^32-1]. The 160 ## The input @var{j} is an integer index in the range [0, 2^32-1]. The values
161 ## values of the output matrix are always exactly the same 161 ## of the output matrix are always exactly the same (reproducibility) for a
162 ## (reproducibility) for a given size input and @var{j} index. 162 ## given size input and @var{j} index.
163 ## 163 ##
164 ## The final optional argument determines the class of the resulting matrix. 164 ## The final optional argument determines the class of the resulting matrix.
165 ## Possible values for @var{class}: @qcode{"uint8"}, @qcode{"uint16"}, 165 ## Possible values for @var{class}: @qcode{"uint8"}, @qcode{"uint16"},
166 ## @qcode{"uint32"}, @qcode{"int8"}, @qcode{"int16"}, int32", @qcode{"single"}, 166 ## @qcode{"uint32"}, @qcode{"int8"}, @qcode{"int16"}, int32", @qcode{"single"},
167 ## @qcode{"double"}. The default is @qcode{"double"}. 167 ## @qcode{"double"}. The default is @qcode{"double"}.
179 ## 179 ##
180 ## @end deftypefn 180 ## @end deftypefn
181 ## 181 ##
182 ## @deftypefn {Function File} {@var{a} =} gallery ("ipjfact", @var{n}) 182 ## @deftypefn {Function File} {@var{a} =} gallery ("ipjfact", @var{n})
183 ## @deftypefnx {Function File} {@var{a} =} gallery ("ipjfact", @var{n}, @var{k}) 183 ## @deftypefnx {Function File} {@var{a} =} gallery ("ipjfact", @var{n}, @var{k})
184 ## Create an Hankel matrix with factorial elements. 184 ## Create a Hankel matrix with factorial elements.
185 ## 185 ##
186 ## @end deftypefn 186 ## @end deftypefn
187 ## 187 ##
188 ## @deftypefn {Function File} {@var{a} =} gallery ("jordbloc", @var{n}) 188 ## @deftypefn {Function File} {@var{a} =} gallery ("jordbloc", @var{n})
189 ## @deftypefnx {Function File} {@var{a} =} gallery ("jordbloc", @var{n}, @var{lambda}) 189 ## @deftypefnx {Function File} {@var{a} =} gallery ("jordbloc", @var{n}, @var{lambda})
255 ## (mean = 0, std = 1). 255 ## (mean = 0, std = 1).
256 ## 256 ##
257 ## The first input is a matrix of dimensions describing the size of the output. 257 ## The first input is a matrix of dimensions describing the size of the output.
258 ## The dimensions can also be input as comma-separated arguments. 258 ## The dimensions can also be input as comma-separated arguments.
259 ## 259 ##
260 ## The input @var{j} is an integer index in the range [0, 2^32-1]. The 260 ## The input @var{j} is an integer index in the range [0, 2^32-1]. The values
261 ## values of the output matrix are always exactly the same 261 ## of the output matrix are always exactly the same (reproducibility) for a
262 ## (reproducibility) for a given size input and @var{j} index. 262 ## given size input and @var{j} index.
263 ## 263 ##
264 ## The final optional argument determines the class of the resulting matrix. 264 ## The final optional argument determines the class of the resulting matrix.
265 ## Possible values for @var{class}: @qcode{"single"}, @qcode{"double"}. 265 ## Possible values for @var{class}: @qcode{"single"}, @qcode{"double"}.
266 ## The default is @qcode{"double"}. 266 ## The default is @qcode{"double"}.
267 ## 267 ##
378 ## (range [0,1]). 378 ## (range [0,1]).
379 ## 379 ##
380 ## The first input is a matrix of dimensions describing the size of the output. 380 ## The first input is a matrix of dimensions describing the size of the output.
381 ## The dimensions can also be input as comma-separated arguments. 381 ## The dimensions can also be input as comma-separated arguments.
382 ## 382 ##
383 ## The input @var{j} is an integer index in the range [0, 2^32-1]. The 383 ## The input @var{j} is an integer index in the range [0, 2^32-1]. The values
384 ## values of the output matrix are always exactly the same 384 ## of the output matrix are always exactly the same (reproducibility) for a
385 ## (reproducibility) for a given size input and @var{j} index. 385 ## given size input and @var{j} index.
386 ## 386 ##
387 ## The final optional argument determines the class of the resulting matrix. 387 ## The final optional argument determines the class of the resulting matrix.
388 ## Possible values for @var{class}: @qcode{"single"}, @qcode{"double"}. 388 ## Possible values for @var{class}: @qcode{"single"}, @qcode{"double"}.
389 ## The default is @qcode{"double"}. 389 ## The default is @qcode{"double"}.
390 ## 390 ##
1213 1213
1214 if (nargin < 1 || nargin > 3) 1214 if (nargin < 1 || nargin > 3)
1215 error ("gallery: 1 to 3 arguments are required for gearmat matrix."); 1215 error ("gallery: 1 to 3 arguments are required for gearmat matrix.");
1216 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) 1216 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
1217 error ("gallery: N must be an integer for gearmat matrix."); 1217 error ("gallery: N must be an integer for gearmat matrix.");
1218 elseif (! isnumeric (i) || ! isscalar (i) || i == 0 || abs (i) <= n) 1218 elseif (! isnumeric (i) || ! isscalar (i) || i == 0 || abs (i) > n)
1219 error ("gallery: I must be a nonzero scalar, and abs (I) <= N for gearmat matrix."); 1219 error ("gallery: I must be a nonzero scalar, and abs (I) <= N for gearmat matrix.");
1220 elseif (! isnumeric (j) || ! isscalar (j) || i == 0 || abs (j) <= n) 1220 elseif (! isnumeric (j) || ! isscalar (j) || i == 0 || abs (j) > n)
1221 error ("gallery: J must be a nonzero scalar, and abs (J) <= N for gearmat matrix."); 1221 error ("gallery: J must be a nonzero scalar, and abs (J) <= N for gearmat matrix.");
1222 endif 1222 endif
1223 1223
1224 A = diag (ones (n-1, 1), -1) + diag (ones (n-1, 1), 1); 1224 A = diag (ones (n-1, 1), -1) + diag (ones (n-1, 1), 1);
1225 A(1, abs (i)) = sign (i); 1225 A(1, abs (i)) = sign (i);
1269 error ("gallery: 1 to 2 arguments are required for hanowa matrix."); 1269 error ("gallery: 1 to 2 arguments are required for hanowa matrix.");
1270 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) 1270 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
1271 error ("gallery: N must be an integer for hanowa matrix."); 1271 error ("gallery: N must be an integer for hanowa matrix.");
1272 elseif (rem (n, 2) != 0) 1272 elseif (rem (n, 2) != 0)
1273 error ("gallery: N must be even for hanowa matrix."); 1273 error ("gallery: N must be even for hanowa matrix.");
1274 elseif (! isnumeric (lambda) || ! isscalar (lambda)) 1274 elseif (! isnumeric (d) || ! isscalar (d))
1275 error ("gallery: D must be a numeric scalar for hanowa matrix."); 1275 error ("gallery: D must be a numeric scalar for hanowa matrix.");
1276 endif 1276 endif
1277 1277
1278 m = n/2; 1278 m = n/2;
1279 A = [ d*eye(m) -diag(1:m) 1279 A = [ d*eye(m) -diag(1:m)
1600 ## W.F. Trench, Numerical solution of the eigenvalue problem 1600 ## W.F. Trench, Numerical solution of the eigenvalue problem
1601 ## for Hermitian Toeplitz matrices, SIAM J. Matrix Analysis and Appl., 1601 ## for Hermitian Toeplitz matrices, SIAM J. Matrix Analysis and Appl.,
1602 ## 10 (1989), pp. 135-146 (and see the references therein). 1602 ## 10 (1989), pp. 135-146 (and see the references therein).
1603 1603
1604 if (nargin < 1 || nargin > 2) 1604 if (nargin < 1 || nargin > 2)
1605 error ("gallery: 1 to 2 arguments are required for lauchli matrix."); 1605 error ("gallery: 1 to 2 arguments are required for kms matrix.");
1606 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) 1606 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
1607 error ("gallery: N must be an integer for lauchli matrix.") 1607 error ("gallery: N must be an integer for kms matrix.")
1608 elseif (! isscalar (mu)) 1608 elseif (! isscalar (rho))
1609 error ("gallery: MU must be a scalar for lauchli matrix.") 1609 error ("gallery: rho must be a scalar for kms matrix.")
1610 endif 1610 endif
1611 1611
1612 A = (1:n)'*ones (1,n); 1612 A = (1:n)'*ones (1,n);
1613 A = abs (A - A'); 1613 A = abs (A - A');
1614 A = rho .^ A; 1614 A = rho .^ A;
2013 2013
2014 if (nargin < 1 || nargin > 2) 2014 if (nargin < 1 || nargin > 2)
2015 error ("gallery: 1 to 2 arguments are required for pei matrix."); 2015 error ("gallery: 1 to 2 arguments are required for pei matrix.");
2016 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) 2016 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
2017 error ("gallery: N must be an integer for pei matrix."); 2017 error ("gallery: N must be an integer for pei matrix.");
2018 elseif (! isnumeric (w) || ! isscalar (w)) 2018 elseif (! isnumeric (alpha) || ! isscalar (alpha)) # change w to alpha
2019 error ("gallery: ALPHA must be a scalar for pei matrix."); 2019 error ("gallery: ALPHA must be a scalar for pei matrix.");
2020 endif 2020 endif
2021 2021
2022 P = alpha * eye (n) + ones (n); 2022 P = alpha * eye (n) + ones (n);
2023 endfunction 2023 endfunction
2854 A = A'; 2854 A = A';
2855 endif 2855 endif
2856 endfunction 2856 endfunction
2857 2857
2858 2858
2859 function B = qmult (A)
2860 ## QMULT Pre-multiply by random orthogonal matrix.
2861 ## QMULT(A) is Q*A where Q is a random real orthogonal matrix from
2862 ## the Haar distribution, of dimension the number of rows in A.
2863 ## Special case: if A is a scalar then QMULT(A) is the same as
2864 ## QMULT(EYE(A)).
2865 ##
2866 ##
2867 ## Called by RANDSVD.
2868 ##
2869 ## Reference:
2870 ## G.W. Stewart, The efficient generation of random
2871 ## orthogonal matrices with an application to condition estimators,
2872 ## SIAM J. Numer. Anal., 17 (1980), 403-409.
2873
2874 [n, m] = size (A);
2875
2876 # Handle scalar A.
2877 if (max (n,m) == 1)
2878 n = A;
2879 A = eye(n);
2880 endif
2881
2882 d = zeros (n);
2883
2884 for k = n-1:-1:1
2885 # Generate random Householder transformation.
2886 x = randn (n-k+1,1);
2887 s = norm (x);
2888 sgn = sign (x(1)) + (x(1) == 0); # Modification for sign(1)=1.
2889 s = sgn * s;
2890 d(k) = -sgn;
2891 x(1) = x(1) + s;
2892 beta = s * x(1);
2893 # Apply the transformation to A.
2894 y = x' * A(k:n,:);
2895 A(k:n,:) = A(k:n,:) - x * (y/beta);
2896 endfor
2897
2898 # Tidy up signs.
2899 for i=1:n-1
2900 A(i,:) = d(i) * A(i,:);
2901 endfor
2902 A(n,:) = A(n,:) * sign (randn);
2903 B = A;
2904 endfunction
2905
2906
2907
2859 ## BIST testing for just a few functions to verify that the main gallery 2908 ## BIST testing for just a few functions to verify that the main gallery
2860 ## dispatch function works. 2909 ## dispatch function works.
2861 %assert (gallery ("clement", 3), [0 1 0; 2 0 2; 0 1 0]) 2910 %assert (gallery ("clement", 3), [0 1 0; 2 0 2; 0 1 0])
2862 %assert (gallery ("invhess", 2), [1 -1; 1 2]) 2911 %assert (gallery ("invhess", 2), [1 -1; 1 2])
2863 2912