Mercurial > octave-nkf

comparison scripts/special-matrix/gallery.m @ 20331:0b9d23557506

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
gallery: fix randsvd by adding missing dependency qmult(). * scripts/special-matrix/gallery.m (randsvd) was copied from the Test Matrix toolbox by Nicholas J. Higham. It made use of qmult() which was also part of that toolbox but was left behind. This qmult() implementation is also recovered from the Test Matrix toolbox. Note that Octave itself used to have qmult() which is part of the legacy quaternion package (now also removed from the new quaternion package). See cset 21904fe299c8 for when qmult() was removed from Octave. Also fix the default value for KL and KU so that a 2 element vector can be used as N.
author Carnë Draug <carandraug@octave.org>
date Fri, 03 Jul 2015 16:18:33 +0100
parents 557979395ca9
children 26fc9bbb8762
comparison
equal deleted inserted replaced
20330:557979395ca9 20331:0b9d23557506
2148 error ("gallery: unknown K '%d' for smoke matrix.", k); 2148 error ("gallery: unknown K '%d' for smoke matrix.", k);
2149 endswitch 2149 endswitch
2150 2150
2151 endfunction 2151 endfunction
2152 2152
2153 function A = randsvd (n, kappa = sqrt (1/eps), mode = 3, kl = n-1, ku = kl) 2153 function A = randsvd (n, kappa = sqrt (1/eps), mode = 3, kl = max (n) -1,
2154 ku = kl)
2154 ## RANDSVD Random matrix with pre-assigned singular values. 2155 ## RANDSVD Random matrix with pre-assigned singular values.
2155 ## RANDSVD(N, KAPPA, MODE, KL, KU) is a (banded) random matrix of order N 2156 ## RANDSVD(N, KAPPA, MODE, KL, KU) is a (banded) random matrix of order N
2156 ## with COND(A) = KAPPA and singular values from the distribution MODE. 2157 ## with COND(A) = KAPPA and singular values from the distribution MODE.
2157 ## N may be a 2-vector, in which case the matrix is N(1)-by-N(2). 2158 ## N may be a 2-vector, in which case the matrix is N(1)-by-N(2).
2158 ## Available types: 2159 ## Available types:
2849 if (flip) 2850 if (flip)
2850 A = A'; 2851 A = A';
2851 endif 2852 endif
2852 endfunction 2853 endfunction
2853 2854
2855 ## NOTE: qmult is part of the Test Matrix Toolbox and is used by randsvd()
2856 function B = qmult (A)
2857 ## QMULT Pre-multiply by random orthogonal matrix.
2858 ## QMULT(A) is Q*A where Q is a random real orthogonal matrix from
2859 ## the Haar distribution, of dimension the number of rows in A.
2860 ## Special case: if A is a scalar then QMULT(A) is the same as
2861 ## QMULT(EYE(A)).
2862 ##
2863 ## Called by RANDSVD.
2864 ##
2865 ## Reference:
2866 ## G.W. Stewart, The efficient generation of random
2867 ## orthogonal matrices with an application to condition estimators,
2868 ## SIAM J. Numer. Anal., 17 (1980), 403-409.
2869
2870 [n, m] = size (A);
2871
2872 ## Handle scalar A
2873 if (isscalar (A))
2874 n = A;
2875 A = eye (n);
2876 endif
2877
2878 d = zeros (n);
2879
2880 for k = n-1:-1:1
2881 ## Generate random Householder transformation.
2882 x = randn (n-k+1, 1);
2883 s = norm (x);
2884 sgn = sign (x(1)) + (x(1) == 0); # Modification for sign(1)=1.
2885 s = sgn*s;
2886 d(k) = -sgn;
2887 x(1) = x(1) + s;
2888 beta = s*x(1);
2889
2890 ## Apply the transformation to A.
2891 y = x'*A(k:n,:);
2892 A(k:n,:) = A(k:n,:) - x*(y/beta);
2893 endfor
2894
2895 ## Tidy up signs
2896 for i = 1:n-1
2897 A(i,:) = d(i)*A(i,:);
2898 endfor
2899 A(n,:) = A(n,:) * sign (randn);
2900 B = A;
2901 endfunction
2854 2902
2855 ## BIST testing for just a few functions to verify that the main gallery 2903 ## BIST testing for just a few functions to verify that the main gallery
2856 ## dispatch function works. 2904 ## dispatch function works.
2857 %assert (gallery ("clement", 3), [0 1 0; 2 0 2; 0 1 0]) 2905 %assert (gallery ("clement", 3), [0 1 0; 2 0 2; 0 1 0])
2858 %assert (gallery ("invhess", 2), [1 -1; 1 2]) 2906 %assert (gallery ("invhess", 2), [1 -1; 1 2])