Mercurial > octave-nkf

comparison scripts/special-matrix/gallery.m @ 19630:0e1f5a750d00

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Periodic merge of gui-release to default.
author John W. Eaton <jwe@octave.org>
date Tue, 20 Jan 2015 10:24:46 -0500
parents cdfc8bc9ab62 446c46af4b42
children 4197fc428c7d
comparison
equal deleted inserted replaced
19626:37d37297acf8 19630:0e1f5a750d00
1294 ## Must have |s| = norm(x), v = x+s*e_1, and 1294 ## Must have |s| = norm(x), v = x+s*e_1, and
1295 ## x'y = x'Hx =(x'Hx)' real => arg(s) = arg(x(1)). 1295 ## x'y = x'Hx =(x'Hx)' real => arg(s) = arg(x(1)).
1296 ## So take s = sign(x(1))*norm(x) (which avoids cancellation). 1296 ## So take s = sign(x(1))*norm(x) (which avoids cancellation).
1297 ## v'v = (x(1)+s)^2 + x(2)^2 + ... + x(n)^2 1297 ## v'v = (x(1)+s)^2 + x(2)^2 + ... + x(n)^2
1298 ## = 2*norm(x)*(norm(x) + |x(1)|). 1298 ## = 2*norm(x)*(norm(x) + |x(1)|).
1299 ## 1299 ##
1300 ## References: 1300 ## References:
1301 ## G.H. Golub and C.F. Van Loan, Matrix Computations, second edition, 1301 ## G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
1302 ## Johns Hopkins University Press, Baltimore, Maryland, 1989. 1302 ## Johns Hopkins University Press, Baltimore, Maryland, 1989.
1303 ## G.W. Stewart, Introduction to Matrix Computations, Academic Press, 1303 ## G.W. Stewart, Introduction to Matrix Computations, Academic Press,
1304 ## New York, 1973, 1304 ## New York, 1973,
1343 "for integerdata matrix."]); 1343 "for integerdata matrix."]);
1344 endif 1344 endif
1345 jidx = varargin{end-1}; 1345 jidx = varargin{end-1};
1346 svec = [varargin{1:end-1}]; 1346 svec = [varargin{1:end-1}];
1347 varargin(end-1) = []; 1347 varargin(end-1) = [];
1348 else 1348 else
1349 error (["gallery: J must be an integer in the range [0, 2^32-1] " ... 1349 error (["gallery: J must be an integer in the range [0, 2^32-1] " ...
1350 "for integerdata matrix"]); 1350 "for integerdata matrix"]);
1351 endif 1351 endif
1352 1352
1353 if (! (isnumeric (jidx) && isscalar (jidx) 1353 if (! (isnumeric (jidx) && isscalar (jidx)
1357 "for integerdata matrix"]); 1357 "for integerdata matrix"]);
1358 endif 1358 endif
1359 1359
1360 ## Save and restore random state. Initialization done so that reproducible 1360 ## Save and restore random state. Initialization done so that reproducible
1361 ## data is available from gallery depending on the jidx and size vector. 1361 ## data is available from gallery depending on the jidx and size vector.
1362 randstate = rand ("state"); 1362 randstate = rand ("state");
1363 unwind_protect 1363 unwind_protect
1364 rand ("state", svec); 1364 rand ("state", svec);
1365 A = randi (varargin{:}); 1365 A = randi (varargin{:});
1366 unwind_protect_cleanup 1366 unwind_protect_cleanup
1367 rand ("state", randstate); 1367 rand ("state", randstate);
1461 ## The determinant and inverse are known explicitly. 1461 ## The determinant and inverse are known explicitly.
1462 ## If a second output argument is present, d = DET(A) is returned: 1462 ## If a second output argument is present, d = DET(A) is returned:
1463 ## [A, d] = IPJFACT(N, K); 1463 ## [A, d] = IPJFACT(N, K);
1464 ## 1464 ##
1465 ## Suggested by P. R. Graves-Morris. 1465 ## Suggested by P. R. Graves-Morris.
1466 ## 1466 ##
1467 ## Reference: 1467 ## Reference:
1468 ## M.J.C. Gover, The explicit inverse of factorial Hankel matrices, 1468 ## M.J.C. Gover, The explicit inverse of factorial Hankel matrices,
1469 ## Dept. of Mathematics, University of Bradford, 1993. 1469 ## Dept. of Mathematics, University of Bradford, 1993.
1470 1470
1471 if (nargin < 1 || nargin > 2) 1471 if (nargin < 1 || nargin > 2)
1862 "for normaldata matrix."]); 1862 "for normaldata matrix."]);
1863 endif 1863 endif
1864 jidx = varargin{end-1}; 1864 jidx = varargin{end-1};
1865 svec = [varargin{1:end-1}]; 1865 svec = [varargin{1:end-1}];
1866 varargin(end-1) = []; 1866 varargin(end-1) = [];
1867 else 1867 else
1868 error (["gallery: J must be an integer in the range [0, 2^32-1] " ... 1868 error (["gallery: J must be an integer in the range [0, 2^32-1] " ...
1869 "for normaldata matrix"]); 1869 "for normaldata matrix"]);
1870 endif 1870 endif
1871 1871
1872 if (! (isnumeric (jidx) && isscalar (jidx) 1872 if (! (isnumeric (jidx) && isscalar (jidx)
1876 "for normaldata matrix"]); 1876 "for normaldata matrix"]);
1877 endif 1877 endif
1878 1878
1879 ## Save and restore random state. Initialization done so that reproducible 1879 ## Save and restore random state. Initialization done so that reproducible
1880 ## data is available from gallery depending on the jidx and size vector. 1880 ## data is available from gallery depending on the jidx and size vector.
1881 randstate = randn ("state"); 1881 randstate = randn ("state");
1882 unwind_protect 1882 unwind_protect
1883 randn ("state", svec); 1883 randn ("state", svec);
1884 A = randn (varargin{:}); 1884 A = randn (varargin{:});
1885 unwind_protect_cleanup 1885 unwind_protect_cleanup
1886 randn ("state", randstate); 1886 randn ("state", randstate);
1981 ## 1981 ##
1982 ## At the Second SIAM Conference on Linear Algebra, Raleigh, N.C., 1982 ## At the Second SIAM Conference on Linear Algebra, Raleigh, N.C.,
1983 ## 1985, Cleve Moler noted that most of the singular values of 1983 ## 1985, Cleve Moler noted that most of the singular values of
1984 ## PARTER(N) are very close to PI. An explanation of the phenomenon 1984 ## PARTER(N) are very close to PI. An explanation of the phenomenon
1985 ## was given by Parter; see also the paper by Tyrtyshnikov. 1985 ## was given by Parter; see also the paper by Tyrtyshnikov.
1986 ## 1986 ##
1987 ## References: 1987 ## References:
1988 ## The MathWorks Newsletter, Volume 1, Issue 1, March 1986, page 2. 1988 ## The MathWorks Newsletter, Volume 1, Issue 1, March 1986, page 2.
1989 ## S.V. Parter, On the distribution of the singular values of Toeplitz 1989 ## S.V. Parter, On the distribution of the singular values of Toeplitz
1990 ## matrices, Linear Algebra and Appl., 80 (1986), pp. 115-130. 1990 ## matrices, Linear Algebra and Appl., 80 (1986), pp. 115-130.
1991 ## E.E. Tyrtyshnikov, Cauchy-Toeplitz matrices and some applications, 1991 ## E.E. Tyrtyshnikov, Cauchy-Toeplitz matrices and some applications,
2087 ## unitary Hessenberg matrix with positive subdiagonal elements in terms 2087 ## unitary Hessenberg matrix with positive subdiagonal elements in terms
2088 ## of 2N-1 real parameters (the Schur parametrization). 2088 ## of 2N-1 real parameters (the Schur parametrization).
2089 ## This M-file handles the real case only and is intended simply as a 2089 ## This M-file handles the real case only and is intended simply as a
2090 ## convenient way to generate random or non-random orthogonal Hessenberg 2090 ## convenient way to generate random or non-random orthogonal Hessenberg
2091 ## matrices. 2091 ## matrices.
2092 ## 2092 ##
2093 ## Reference: 2093 ## Reference:
2094 ## W.B. Gragg, The QR algorithm for unitary Hessenberg matrices, 2094 ## W.B. Gragg, The QR algorithm for unitary Hessenberg matrices,
2095 ## J. Comp. Appl. Math., 16 (1986), pp. 1-8. 2095 ## J. Comp. Appl. Math., 16 (1986), pp. 1-8.
2096 2096
2097 if (nargin != 1) 2097 if (nargin != 1)
2583 "for uniformdata matrix."]); 2583 "for uniformdata matrix."]);
2584 endif 2584 endif
2585 jidx = varargin{end-1}; 2585 jidx = varargin{end-1};
2586 svec = [varargin{1:end-1}]; 2586 svec = [varargin{1:end-1}];
2587 varargin(end-1) = []; 2587 varargin(end-1) = [];
2588 else 2588 else
2589 error (["gallery: J must be an integer in the range [0, 2^32-1] " ... 2589 error (["gallery: J must be an integer in the range [0, 2^32-1] " ...
2590 "for uniformdata matrix"]); 2590 "for uniformdata matrix"]);
2591 endif 2591 endif
2592 2592
2593 if (! (isnumeric (jidx) && isscalar (jidx) 2593 if (! (isnumeric (jidx) && isscalar (jidx)
2597 "for uniformdata matrix"]); 2597 "for uniformdata matrix"]);
2598 endif 2598 endif
2599 2599
2600 ## Save and restore random state. Initialization done so that reproducible 2600 ## Save and restore random state. Initialization done so that reproducible
2601 ## data is available from gallery depending on the jidx and size vector. 2601 ## data is available from gallery depending on the jidx and size vector.
2602 randstate = rand ("state"); 2602 randstate = rand ("state");
2603 unwind_protect 2603 unwind_protect
2604 rand ("state", svec); 2604 rand ("state", svec);
2605 A = rand (varargin{:}); 2605 A = rand (varargin{:});
2606 unwind_protect_cleanup 2606 unwind_protect_cleanup
2607 rand ("state", randstate); 2607 rand ("state", randstate);
2713 2713
2714 rho = 100 * rand (nx, ny); 2714 rho = 100 * rand (nx, ny);
2715 2715
2716 for j = 1:ny 2716 for j = 1:ny
2717 for i = 1:nx 2717 for i = 1:nx
2718 ## 2718 ##
2719 ## For the element (I,J), determine the indices of the 8 nodes. 2719 ## For the element (I,J), determine the indices of the 8 nodes.
2720 ## 2720 ##
2721 nn(1) = 3*j*nx + 2*i + 2*j + 1; 2721 nn(1) = 3*j*nx + 2*i + 2*j + 1;
2722 nn(2) = nn(1) - 1; 2722 nn(2) = nn(1) - 1;
2723 nn(3) = nn(2) - 1; 2723 nn(3) = nn(2) - 1;
2724 nn(4) = (3*j - 1) * nx + 2*j + i - 1; 2724 nn(4) = (3*j - 1) * nx + 2*j + i - 1;
2725 nn(5) = 3 * (j-1) * nx + 2*i + 2*j - 3; 2725 nn(5) = 3 * (j-1) * nx + 2*i + 2*j - 3;
2815 ## 2815 ##
2816 ## Called by RANDSVD. 2816 ## Called by RANDSVD.
2817 ## This is a `standard' reduction. Cf. reduction to bidiagonal form 2817 ## This is a `standard' reduction. Cf. reduction to bidiagonal form
2818 ## prior to computing the SVD. This code is a little wasteful in that 2818 ## prior to computing the SVD. This code is a little wasteful in that
2819 ## it computes certain elements which are immediately set to zero! 2819 ## it computes certain elements which are immediately set to zero!
2820 ## 2820 ##
2821 ## Reference: 2821 ## Reference:
2822 ## G.H. Golub and C.F. Van Loan, Matrix Computations, second edition, 2822 ## G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
2823 ## Johns Hopkins University Press, Baltimore, Maryland, 1989. 2823 ## Johns Hopkins University Press, Baltimore, Maryland, 1989.
2824 ## Section 5.4.3. 2824 ## Section 5.4.3.
2825 2825