Mercurial > octave-antonio
annotate scripts/special-matrix/gallery.m @ 16816:12005245b645
doc: Periodic grammarcheck of documentation.
* doc/interpreter/basics.txi, doc/interpreter/container.txi,
doc/interpreter/contrib.txi, doc/interpreter/diagperm.txi,
doc/interpreter/errors.txi, doc/interpreter/install.txi,
doc/interpreter/sparse.txi, libinterp/corefcn/ellipj.cc,
libinterp/corefcn/mappers.cc, libinterp/corefcn/regexp.cc,
libinterp/corefcn/tril.cc, libinterp/dldfcn/__init_fltk__.cc,
libinterp/dldfcn/fftw.cc, libinterp/interpfcn/dirfns.cc,
libinterp/interpfcn/input.cc, libinterp/octave-value/ov-fcn-inline.cc,
libinterp/octave-value/ov-struct.cc, scripts/@ftp/cd.m,
scripts/general/interp1.m, scripts/general/num2str.m, scripts/image/ind2rgb.m,
scripts/image/rgb2ind.m, scripts/io/importdata.m, scripts/io/textread.m,
scripts/java/javamem.m, scripts/linear-algebra/condest.m,
scripts/linear-algebra/onenormest.m, scripts/miscellaneous/error_ids.m,
scripts/miscellaneous/getfield.m, scripts/miscellaneous/setfield.m,
scripts/plot/area.m, scripts/plot/pcolor.m, scripts/plot/stairs.m,
scripts/set/powerset.m, scripts/sparse/bicg.m, scripts/sparse/bicgstab.m,
scripts/sparse/cgs.m, scripts/specfun/ellipke.m,
scripts/special-matrix/gallery.m, scripts/strings/strjoin.m,
scripts/strings/strsplit.m, scripts/testfun/__have_feature__.m,
scripts/testfun/__printf_assert__.m, scripts/testfun/__prog_output_assert__.m,
scripts/testfun/__run_test_suite__.m: grammarcheck documentation.
author | Rik <rik@octave.org> |
---|---|
date | Sat, 22 Jun 2013 19:47:32 -0700 |
parents | 7268845c0a1e |
children | e39f00a32dc7 |
rev | line source |
---|---|
16634 | 1 ## Copyright (C) 1989-1995 Nicholas .J. Higham |
2 ## Copyright (C) 2013 Carnë Draug | |
3 ## | |
4 ## This file is part of Octave. | |
5 ## | |
6 ## Octave is free software; you can redistribute it and/or modify it | |
7 ## under the terms of the GNU General Public License as published by | |
8 ## the Free Software Foundation; either version 3 of the License, or (at | |
9 ## your option) any later version. | |
10 ## | |
11 ## Octave is distributed in the hope that it will be useful, but | |
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
14 ## General Public License for more details. | |
15 ## | |
16 ## You should have received a copy of the GNU General Public License | |
17 ## along with Octave; see the file COPYING. If not, see | |
18 ## <http://www.gnu.org/licenses/>. | |
19 | |
20 ## -*- texinfo -*- | |
21 ## @deftypefn {Function File} {} gallery (@var{name}) | |
22 ## @deftypefnx {Function File} {} gallery (@var{name}, @var{args}) | |
23 ## Create interesting matrices for testing. | |
24 ## | |
25 ## | |
26 ## @end deftypefn | |
27 ## | |
28 ## @deftypefn {Function File} {@var{c} =} gallery ("cauchy", @var{x}) | |
29 ## @deftypefnx {Function File} {@var{c} =} gallery ("cauchy", @var{x}, @var{y}) | |
30 ## Create a Cauchy matrix. | |
31 ## | |
32 ## @end deftypefn | |
33 ## | |
34 ## @deftypefn {Function File} {@var{c} =} gallery ("chebspec", @var{n}) | |
35 ## @deftypefnx {Function File} {@var{c} =} gallery ("chebspec", @var{n}, @var{k}) | |
36 ## Create a Chebyshev spectral differentiation matrix. | |
37 ## | |
38 ## @end deftypefn | |
39 ## | |
40 ## @deftypefn {Function File} {@var{c} =} gallery ("chebvand", @var{p}) | |
41 ## @deftypefnx {Function File} {@var{c} =} gallery ("chebvand", @var{m}, @var{p}) | |
42 ## Create a Vandermonde-like matrix for the Chebyshev polynomials. | |
43 ## | |
44 ## @end deftypefn | |
45 ## | |
46 ## @deftypefn {Function File} {@var{a} =} gallery ("chow", @var{n}) | |
47 ## @deftypefnx {Function File} {@var{a} =} gallery ("chow", @var{n}, @var{alpha}) | |
48 ## @deftypefnx {Function File} {@var{a} =} gallery ("chow", @var{n}, @var{alpha}, @var{delta}) | |
49 ## Create a Chow matrix -- a singular Toeplitz lower Hessenberg matrix. | |
50 ## | |
51 ## @end deftypefn | |
52 ## | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
53 ## @deftypefn {Function File} {@var{c} =} gallery ("circul", @var{v}) |
16634 | 54 ## Create a circulant matrix. |
55 ## | |
56 ## @end deftypefn | |
57 ## | |
58 ## @deftypefn {Function File} {@var{a} =} gallery ("clement", @var{n}) | |
59 ## @deftypefnx {Function File} {@var{a} =} gallery ("clement", @var{n}, @var{k}) | |
60 ## Create a tridiagonal matrix with zero diagonal entries. | |
61 ## | |
62 ## @end deftypefn | |
63 ## | |
64 ## @deftypefn {Function File} {@var{c} =} gallery ("compar", @var{a}) | |
65 ## @deftypefnx {Function File} {@var{c} =} gallery ("compar", @var{a}, @var{k}) | |
66 ## Create a comparison matrix. | |
67 ## | |
68 ## @end deftypefn | |
69 ## | |
70 ## @deftypefn {Function File} {@var{a} =} gallery ("condex", @var{n}) | |
71 ## @deftypefnx {Function File} {@var{a} =} gallery ("condex", @var{n}, @var{k}) | |
72 ## @deftypefnx {Function File} {@var{a} =} gallery ("condex", @var{n}, @var{k}, @var{theta}) | |
73 ## Create a `counterexample' matrix to a condition estimator. | |
74 ## | |
75 ## @end deftypefn | |
76 ## | |
77 ## @deftypefn {Function File} {@var{a} =} gallery ("cycol", [@var{m} @var{n}]) | |
78 ## @deftypefnx {Function File} {@var{a} =} gallery ("cycol", @var{n}) | |
79 ## @deftypefnx {Function File} {@var{a} =} gallery (@dots{}, @var{k}) | |
80 ## Create a matrix whose columns repeat cyclically. | |
81 ## | |
82 ## @end deftypefn | |
83 ## | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
84 ## @deftypefn {Function File} {[@var{c}, @var{d}, @var{e}] =} gallery ("dorr", @var{n}) |
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
85 ## @deftypefnx {Function File} {[@var{c}, @var{d}, @var{e}] =} gallery ("dorr", @var{n}, @var{theta}) |
16634 | 86 ## @deftypefnx {Function File} {@var{a} =} gallery ("dorr", @dots{}) |
87 ## Create a diagonally dominant, ill conditioned, tridiagonal matrix. | |
88 ## | |
89 ## @end deftypefn | |
90 ## | |
91 ## @deftypefn {Function File} {@var{a} =} gallery ("dramadah", @var{n}) | |
92 ## @deftypefnx {Function File} {@var{a} =} gallery ("dramadah", @var{n}, @var{k}) | |
93 ## Create a (0, 1) matrix whose inverse has large integer entries. | |
94 ## | |
95 ## @end deftypefn | |
96 ## | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
97 ## @deftypefn {Function File} {@var{a} =} gallery ("fiedler", @var{c}) |
16634 | 98 ## Create a symmetric Fiedler matrix. |
99 ## | |
100 ## @end deftypefn | |
101 ## | |
102 ## @deftypefn {Function File} {@var{a} =} gallery ("forsythe", @var{n}) | |
103 ## @deftypefnx {Function File} {@var{a} =} gallery ("forsythe", @var{n}, @var{alpha}) | |
104 ## @deftypefnx {Function File} {@var{a} =} gallery ("forsythe", @var{n}, @var{alpha}, @var{lambda}) | |
105 ## Create a Forsythe matrix (a perturbed Jordan block). | |
106 ## | |
107 ## @end deftypefn | |
108 ## | |
109 ## @deftypefn {Function File} {@var{f} =} gallery ("frank", @var{n}) | |
110 ## @deftypefnx {Function File} {@var{f} =} gallery ("frank", @var{n}, @var{k}) | |
111 ## Create a Frank matrix (ill conditioned eigenvalues). | |
112 ## | |
113 ## @end deftypefn | |
114 ## | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
115 ## @deftypefn {Function File} {@var{c} =} gallery ("gcdmat", @var{n}) |
16634 | 116 ## Create a greatest common divisor matrix. |
117 ## | |
118 ## @var{c} is an @var{n}-by-@var{n} matrix whose values correspond to the | |
119 ## greatest common divisor of its coordinate values, i.e., @var{c}(i,j) | |
120 ## correspond @code{gcd (i, j)}. | |
121 ## @end deftypefn | |
122 ## | |
123 ## @deftypefn {Function File} {@var{a} =} gallery ("gearmat", @var{n}) | |
124 ## @deftypefnx {Function File} {@var{a} =} gallery ("gearmat", @var{n}, @var{i}) | |
125 ## @deftypefnx {Function File} {@var{a} =} gallery ("gearmat", @var{n}, @var{i}, @var{j}) | |
126 ## Create a Gear matrix. | |
127 ## | |
128 ## @end deftypefn | |
129 ## | |
130 ## @deftypefn {Function File} {@var{g} =} gallery ("grcar", @var{n}) | |
131 ## @deftypefnx {Function File} {@var{g} =} gallery ("grcar", @var{n}, @var{k}) | |
132 ## Create a Toeplitz matrix with sensitive eigenvalues. | |
133 ## | |
134 ## @end deftypefn | |
135 ## | |
136 ## @deftypefn {Function File} {@var{a} =} gallery ("hanowa", @var{n}) | |
137 ## @deftypefnx {Function File} {@var{a} =} gallery ("hanowa", @var{n}, @var{d}) | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
138 ## Create a matrix whose eigenvalues lie on a vertical line in the complex |
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
139 ## plane. |
16634 | 140 ## |
141 ## @end deftypefn | |
142 ## | |
143 ## @deftypefn {Function File} {@var{v} =} gallery ("house", @var{x}) | |
144 ## @deftypefnx {Function File} {[@var{v}, @var{beta}] =} gallery ("house", @var{x}) | |
145 ## Create a householder matrix. | |
146 ## | |
147 ## @end deftypefn | |
148 ## | |
149 ## @deftypefn {Function File} {@var{a} =} gallery ("invhess", @var{x}) | |
150 ## @deftypefnx {Function File} {@var{a} =} gallery ("invhess", @var{x}, @var{y}) | |
151 ## Create the inverse of an upper Hessenberg matrix. | |
152 ## | |
153 ## @end deftypefn | |
154 ## | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
155 ## @deftypefn {Function File} {@var{a} =} gallery ("invol", @var{n}) |
16634 | 156 ## Create an involutory matrix. |
157 ## | |
158 ## @end deftypefn | |
159 ## | |
160 ## @deftypefn {Function File} {@var{a} =} gallery ("ipjfact", @var{n}) | |
161 ## @deftypefnx {Function File} {@var{a} =} gallery ("ipjfact", @var{n}, @var{k}) | |
162 ## Create an Hankel matrix with factorial elements. | |
163 ## | |
164 ## @end deftypefn | |
165 ## | |
166 ## @deftypefn {Function File} {@var{a} =} gallery ("jordbloc", @var{n}) | |
167 ## @deftypefnx {Function File} {@var{a} =} gallery ("jordbloc", @var{n}, @var{lambda}) | |
168 ## Create a Jordan block. | |
169 ## | |
170 ## @end deftypefn | |
171 ## | |
172 ## @deftypefn {Function File} {@var{u} =} gallery ("kahan", @var{n}) | |
173 ## @deftypefnx {Function File} {@var{u} =} gallery ("kahan", @var{n}, @var{theta}) | |
174 ## @deftypefnx {Function File} {@var{u} =} gallery ("kahan", @var{n}, @var{theta}, @var{pert}) | |
175 ## Create a Kahan matrix (upper trapezoidal). | |
176 ## | |
177 ## @end deftypefn | |
178 ## | |
179 ## @deftypefn {Function File} {@var{a} =} gallery ("kms", @var{n}) | |
180 ## @deftypefnx {Function File} {@var{a} =} gallery ("kms", @var{n}, @var{rho}) | |
181 ## Create a Kac-Murdock-Szego Toeplitz matrix. | |
182 ## | |
183 ## @end deftypefn | |
184 ## | |
185 ## @deftypefn {Function File} {@var{b} =} gallery ("krylov", @var{a}) | |
186 ## @deftypefnx {Function File} {@var{b} =} gallery ("krylov", @var{a}, @var{x}) | |
187 ## @deftypefnx {Function File} {@var{b} =} gallery ("krylov", @var{a}, @var{x}, @var{j}) | |
188 ## Create a Krylov matrix. | |
189 ## | |
190 ## @end deftypefn | |
191 ## | |
192 ## @deftypefn {Function File} {@var{a} =} gallery ("lauchli", @var{n}) | |
193 ## @deftypefnx {Function File} {@var{a} =} gallery ("lauchli", @var{n}, @var{mu}) | |
194 ## Create a Lauchli matrix (rectangular). | |
195 ## | |
196 ## @end deftypefn | |
197 ## | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
198 ## @deftypefn {Function File} {@var{a} =} gallery ("lehmer", @var{n}) |
16634 | 199 ## Create a Lehmer matrix (symmetric positive definite). |
200 ## | |
201 ## @end deftypefn | |
202 ## | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
203 ## @deftypefn {Function File} {@var{t} =} gallery ("lesp", @var{n}) |
16634 | 204 ## Create a tridiagonal matrix with real, sensitive eigenvalues. |
205 ## | |
206 ## @end deftypefn | |
207 ## | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
208 ## @deftypefn {Function File} {@var{a} =} gallery ("lotkin", @var{n}) |
16634 | 209 ## Create a Lotkin matrix. |
210 ## | |
211 ## @end deftypefn | |
212 ## | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
213 ## @deftypefn {Function File} {@var{a} =} gallery ("minij", @var{n}) |
16634 | 214 ## Create a symmetric positive definite matrix MIN(i,j). |
215 ## | |
216 ## @end deftypefn | |
217 ## | |
218 ## @deftypefn {Function File} {@var{a} =} gallery ("moler", @var{n}) | |
219 ## @deftypefnx {Function File} {@var{a} =} gallery ("moler", @var{n}, @var{alpha}) | |
220 ## Create a Moler matrix (symmetric positive definite). | |
221 ## | |
222 ## @end deftypefn | |
223 ## | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
224 ## @deftypefn {Function File} {[@var{a}, @var{t}] =} gallery ("neumann", @var{n}) |
16634 | 225 ## Create a singular matrix from the discrete Neumann problem (sparse). |
226 ## | |
227 ## @end deftypefn | |
228 ## | |
229 ## @deftypefn {Function File} {@var{q} =} gallery ("orthog", @var{n}) | |
230 ## @deftypefnx {Function File} {@var{q} =} gallery ("orthog", @var{n}, @var{k}) | |
231 ## Create orthogonal and nearly orthogonal matrices. | |
232 ## | |
233 ## @end deftypefn | |
234 ## | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
235 ## @deftypefn {Function File} {@var{a} =} gallery ("parter", @var{n}) |
16634 | 236 ## Create a Parter matrix (a Toeplitz matrix with singular values near pi). |
237 ## | |
238 ## @end deftypefn | |
239 ## | |
240 ## @deftypefn {Function File} {@var{p} =} gallery ("pei", @var{n}) | |
241 ## @deftypefnx {Function File} {@var{p} =} gallery ("pei", @var{n}, @var{alpha}) | |
242 ## Create a Pei matrix. | |
243 ## | |
244 ## @end deftypefn | |
245 ## | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
246 ## @deftypefn {Function File} {@var{a} =} gallery ("Poisson", @var{n}) |
16634 | 247 ## Create a block tridiagonal matrix from Poisson's equation (sparse). |
248 ## | |
249 ## @end deftypefn | |
250 ## | |
251 ## @deftypefn {Function File} {@var{a} =} gallery ("prolate", @var{n}) | |
252 ## @deftypefnx {Function File} {@var{a} =} gallery ("prolate", @var{n}, @var{w}) | |
253 ## Create a prolate matrix (symmetric, ill-conditioned Toeplitz matrix). | |
254 ## | |
255 ## @end deftypefn | |
256 ## | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
257 ## @deftypefn {Function File} {@var{h} =} gallery ("randhess", @var{x}) |
16634 | 258 ## Create a random, orthogonal upper Hessenberg matrix. |
259 ## | |
260 ## @end deftypefn | |
261 ## | |
262 ## @deftypefn {Function File} {@var{a} =} gallery ("rando", @var{n}) | |
263 ## @deftypefnx {Function File} {@var{a} =} gallery ("rando", @var{n}, @var{k}) | |
264 ## Create a random matrix with elements -1, 0 or 1. | |
265 ## | |
266 ## @end deftypefn | |
267 ## | |
268 ## @deftypefn {Function File} {@var{a} =} gallery ("randsvd", @var{n}) | |
269 ## @deftypefnx {Function File} {@var{a} =} gallery ("randsvd", @var{n}, @var{kappa}) | |
270 ## @deftypefnx {Function File} {@var{a} =} gallery ("randsvd", @var{n}, @var{kappa}, @var{mode}) | |
271 ## @deftypefnx {Function File} {@var{a} =} gallery ("randsvd", @var{n}, @var{kappa}, @var{mode}, @var{kl}) | |
272 ## @deftypefnx {Function File} {@var{a} =} gallery ("randsvd", @var{n}, @var{kappa}, @var{mode}, @var{kl}, @var{ku}) | |
273 ## Create a random matrix with pre-assigned singular values. | |
274 ## | |
275 ## @end deftypefn | |
276 ## | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
277 ## @deftypefn {Function File} {@var{a} =} gallery ("redheff", @var{n}) |
16634 | 278 ## Create a zero and ones matrix of Redheffer associated with the Riemann |
279 ## hypothesis. | |
280 ## | |
281 ## @end deftypefn | |
282 ## | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
283 ## @deftypefn {Function File} {@var{a} =} gallery ("riemann", @var{n}) |
16634 | 284 ## Create a matrix associated with the Riemann hypothesis. |
285 ## | |
286 ## @end deftypefn | |
287 ## | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
288 ## @deftypefn {Function File} {@var{a} =} gallery ("ris", @var{n}) |
16634 | 289 ## Create a symmetric Hankel matrix. |
290 ## | |
291 ## @end deftypefn | |
292 ## | |
293 ## @deftypefn {Function File} {@var{a} =} gallery ("smoke", @var{n}) | |
294 ## @deftypefnx {Function File} {@var{a} =} gallery ("smoke", @var{n}, @var{k}) | |
295 ## Create a complex matrix, with a `smoke ring' pseudospectrum. | |
296 ## | |
297 ## @end deftypefn | |
298 ## | |
299 ## @deftypefn {Function File} {@var{t} =} gallery ("toeppd", @var{n}) | |
300 ## @deftypefnx {Function File} {@var{t} =} gallery ("toeppd", @var{n}, @var{m}) | |
301 ## @deftypefnx {Function File} {@var{t} =} gallery ("toeppd", @var{n}, @var{m}, @var{w}) | |
302 ## @deftypefnx {Function File} {@var{t} =} gallery ("toeppd", @var{n}, @var{m}, @var{w}, @var{theta}) | |
303 ## Create a symmetric positive definite Toeplitz matrix. | |
304 ## | |
305 ## @end deftypefn | |
306 ## | |
307 ## @deftypefn {Function File} {@var{p} =} gallery ("toeppen", @var{n}) | |
308 ## @deftypefnx {Function File} {@var{p} =} gallery ("toeppen", @var{n}, @var{a}) | |
309 ## @deftypefnx {Function File} {@var{p} =} gallery ("toeppen", @var{n}, @var{a}, @var{b}) | |
310 ## @deftypefnx {Function File} {@var{p} =} gallery ("toeppen", @var{n}, @var{a}, @var{b}, @var{c}) | |
311 ## @deftypefnx {Function File} {@var{p} =} gallery ("toeppen", @var{n}, @var{a}, @var{b}, @var{c}, @var{d}) | |
312 ## @deftypefnx {Function File} {@var{p} =} gallery ("toeppen", @var{n}, @var{a}, @var{b}, @var{c}, @var{d}, @var{e}) | |
313 ## Create a pentadiagonal Toeplitz matrix (sparse). | |
314 ## | |
315 ## @end deftypefn | |
316 ## | |
317 ## @deftypefn {Function File} {@var{a} =} gallery ("tridiag", @var{x}, @var{y}, @var{z}) | |
318 ## @deftypefnx {Function File} {@var{a} =} gallery ("tridiag", @var{n}) | |
319 ## @deftypefnx {Function File} {@var{a} =} gallery ("tridiag", @var{n}, @var{c}, @var{d}, @var{e}) | |
320 ## Create a tridiagonal matrix (sparse). | |
321 ## | |
322 ## @end deftypefn | |
323 ## | |
324 ## @deftypefn {Function File} {@var{t} =} gallery ("triw", @var{n}) | |
325 ## @deftypefnx {Function File} {@var{t} =} gallery ("triw", @var{n}, @var{alpha}) | |
326 ## @deftypefnx {Function File} {@var{t} =} gallery ("triw", @var{n}, @var{alpha}, @var{k}) | |
327 ## Create an upper triangular matrix discussed by Kahan, Golub and Wilkinson. | |
328 ## | |
329 ## @end deftypefn | |
330 ## | |
331 ## @deftypefn {Function File} {@var{a} =} gallery ("wathen", @var{nx}, @var{ny}) | |
332 ## @deftypefnx {Function File} {@var{a} =} gallery ("wathen", @var{nx}, @var{ny}, @var{k}) | |
333 ## Create the Wathen matrix. | |
334 ## | |
335 ## @end deftypefn | |
336 ## | |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
16766
diff
changeset
|
337 ## @deftypefn {Function File} {[@var{a}, @var{b}] =} gallery ("wilk", @var{n}) |
16634 | 338 ## Create various specific matrices devised/discussed by Wilkinson. |
339 ## | |
340 ## @end deftypefn | |
341 | |
342 ## Code for most of the individual matrices (except binomial, gcdmat, | |
343 ## integerdata, leslie, normaldata, randcolu, randcorr, randjorth, sampling, | |
344 ## uniformdata) by Nicholas .J. Higham <Nicholas.J.Higham@manchester.ac.uk> | |
345 ## Adapted for Octave and into single gallery function by Carnë Draug | |
346 | |
347 function [varargout] = gallery (name, varargin) | |
348 | |
349 if (nargin < 1) | |
350 print_usage (); | |
351 elseif (! ischar (name)) | |
352 error ("gallery: NAME must be a string."); | |
353 endif | |
354 | |
355 ## NOTE: there isn't a lot of input check in the individual functions | |
356 ## that actually build the functions. This is by design. The original | |
357 ## code by Higham did not perform it and was propagated to Matlab, so | |
358 ## for compatibility, we also don't make it. For example, arguments | |
359 ## that behave as switches, and in theory accepting a value of 0 or 1, | |
360 ## will use a value of 0, for any value other than 1 (only check made | |
361 ## is if the value is equal to 1). It will often also accept string | |
362 ## values instead of numeric. Only input check added was where it | |
363 ## would be causing an error anyway. | |
364 | |
365 ## we will always want to return at least 1 output | |
366 n_out = nargout; | |
367 if (n_out == 0) | |
368 n_out = 1; | |
369 endif | |
370 | |
371 switch (tolower (name)) | |
372 case "binomial" | |
373 error ("gallery: matrix %s not implemented.", name); | |
374 case "cauchy" , [varargout{1:n_out}] = cauchy (varargin{:}); | |
375 case "chebspec" , [varargout{1:n_out}] = chebspec (varargin{:}); | |
376 case "chebvand" , [varargout{1:n_out}] = chebvand (varargin{:}); | |
377 case "chow" , [varargout{1:n_out}] = chow (varargin{:}); | |
378 case "circul" , [varargout{1:n_out}] = circul (varargin{:}); | |
379 case "clement" , [varargout{1:n_out}] = clement (varargin{:}); | |
380 case "compar" , [varargout{1:n_out}] = compar (varargin{:}); | |
381 case "condex" , [varargout{1:n_out}] = condex (varargin{:}); | |
382 case "cycol" , [varargout{1:n_out}] = cycol (varargin{:}); | |
383 case "dorr" , [varargout{1:n_out}] = dorr (varargin{:}); | |
384 case "dramadah" , [varargout{1:n_out}] = dramadah (varargin{:}); | |
385 case "fiedler" , [varargout{1:n_out}] = fiedler (varargin{:}); | |
386 case "forsythe" , [varargout{1:n_out}] = forsythe (varargin{:}); | |
387 case "frank" , [varargout{1:n_out}] = frank (varargin{:}); | |
388 case "gearmat" , [varargout{1:n_out}] = gearmat (varargin{:}); | |
389 case "gcdmat" , [varargout{1:n_out}] = gcdmat (varargin{:}); | |
390 case "grcar" , [varargout{1:n_out}] = grcar (varargin{:}); | |
391 case "hanowa" , [varargout{1:n_out}] = hanowa (varargin{:}); | |
392 case "house" , [varargout{1:n_out}] = house (varargin{:}); | |
393 case "integerdata" | |
394 error ("gallery: matrix %s not implemented.", name); | |
395 case "invhess" , [varargout{1:n_out}] = invhess (varargin{:}); | |
396 case "invol" , [varargout{1:n_out}] = invol (varargin{:}); | |
397 case "ipjfact" , [varargout{1:n_out}] = ipjfact (varargin{:}); | |
398 case "jordbloc" , [varargout{1:n_out}] = jordbloc (varargin{:}); | |
399 case "kahan" , [varargout{1:n_out}] = kahan (varargin{:}); | |
400 case "kms" , [varargout{1:n_out}] = kms (varargin{:}); | |
401 case "krylov" , [varargout{1:n_out}] = krylov (varargin{:}); | |
402 case "lauchli" , [varargout{1:n_out}] = lauchli (varargin{:}); | |
403 case "lehmer" , [varargout{1:n_out}] = lehmer (varargin{:}); | |
404 case "leslie" | |
405 error ("gallery: matrix %s not implemented.", name); | |
406 case "lesp" , [varargout{1:n_out}] = lesp (varargin{:}); | |
407 case "lotkin" , [varargout{1:n_out}] = lotkin (varargin{:}); | |
408 case "minij" , [varargout{1:n_out}] = minij (varargin{:}); | |
409 case "moler" , [varargout{1:n_out}] = moler (varargin{:}); | |
410 case "neumann" , [varargout{1:n_out}] = neumann (varargin{:}); | |
411 case "normaldata" | |
412 error ("gallery: matrix %s not implemented.", name); | |
413 case "orthog" , [varargout{1:n_out}] = orthog (varargin{:}); | |
414 case "parter" , [varargout{1:n_out}] = parter (varargin{:}); | |
415 case "pei" , [varargout{1:n_out}] = pei (varargin{:}); | |
416 case "poisson" , [varargout{1:n_out}] = poisson (varargin{:}); | |
417 case "prolate" , [varargout{1:n_out}] = prolate (varargin{:}); | |
418 case "randcolu" | |
419 error ("gallery: matrix %s not implemented.", name); | |
420 case "randcorr" | |
421 error ("gallery: matrix %s not implemented.", name); | |
422 case "randhess" , [varargout{1:n_out}] = randhess (varargin{:}); | |
423 case "randjorth" | |
424 error ("gallery: matrix %s not implemented.", name); | |
425 case "rando" , [varargout{1:n_out}] = rando (varargin{:}); | |
426 case "randsvd" , [varargout{1:n_out}] = randsvd (varargin{:}); | |
427 case "redheff" , [varargout{1:n_out}] = redheff (varargin{:}); | |
428 case "riemann" , [varargout{1:n_out}] = riemann (varargin{:}); | |
429 case "ris" , [varargout{1:n_out}] = ris (varargin{:}); | |
430 case "sampling" | |
431 error ("gallery: matrix %s not implemented.", name); | |
432 case "smoke" , [varargout{1:n_out}] = smoke (varargin{:}); | |
433 case "toeppd" , [varargout{1:n_out}] = toeppd (varargin{:}); | |
434 case "toeppen" , [varargout{1:n_out}] = toeppen (varargin{:}); | |
435 case "tridiag" , [varargout{1:n_out}] = tridiag (varargin{:}); | |
436 case "triw" , [varargout{1:n_out}] = triw (varargin{:}); | |
437 case "uniformdata" | |
438 error ("gallery: matrix %s not implemented.", name); | |
439 case "wathen" , [varargout{1:n_out}] = wathen (varargin{:}); | |
440 case "wilk" , [varargout{1:n_out}] = wilk (varargin{:}); | |
441 otherwise | |
442 error ("gallery: unknown matrix with NAME %s", name); | |
443 endswitch | |
444 | |
445 endfunction | |
446 | |
447 function C = cauchy (x, y) | |
448 ##CAUCHY Cauchy matrix. | |
449 ## C = CAUCHY(X, Y), where X, Y are N-vectors, is the N-by-N matrix | |
450 ## with C(i,j) = 1/(X(i)+Y(j)). By default, Y = X. | |
451 ## Special case: if X is a scalar CAUCHY(X) is the same as CAUCHY(1:X). | |
452 ## Explicit formulas are known for DET(C) (which is nonzero if X and Y | |
453 ## both have distinct elements) and the elements of INV(C). | |
454 ## C is totally positive if 0 < X(1) < ... < X(N) and | |
455 ## 0 < Y(1) < ... < Y(N). | |
456 ## | |
457 ## References: | |
458 ## N.J. Higham, Accuracy and Stability of Numerical Algorithms, | |
459 ## Society for Industrial and Applied Mathematics, Philadelphia, PA, | |
460 ## USA, 1996; sec. 26.1. | |
461 ## D.E. Knuth, The Art of Computer Programming, Volume 1, | |
462 ## Fundamental Algorithms, second edition, Addison-Wesley, Reading, | |
463 ## Massachusetts, 1973, p. 36. | |
464 ## E.E. Tyrtyshnikov, Cauchy-Toeplitz matrices and some applications, | |
465 ## Linear Algebra and Appl., 149 (1991), pp. 1-18. | |
466 ## O. Taussky and M. Marcus, Eigenvalues of finite matrices, in | |
467 ## Survey of Numerical Analysis, J. Todd, ed., McGraw-Hill, New York, | |
468 ## pp. 279-313, 1962. (States the totally positive property on p. 295.) | |
469 | |
470 if (nargin < 1 || nargin > 2) | |
471 error ("gallery: 1 or 2 arguments are required for cauchy matrix."); | |
472 elseif (! isnumeric (x)) | |
473 error ("gallery: X must be numeric for cauchy matrix."); | |
474 elseif (nargin == 2 && ! isnumeric (y)) | |
475 error ("gallery: Y must be numeric for cauchy matrix."); | |
476 endif | |
477 | |
478 n = numel (x); | |
479 if (isscalar (x) && fix (x) == x) | |
480 n = x; | |
481 x = 1:n; | |
482 elseif (n > 1 && isvector (x)) | |
483 ## do nothing | |
484 else | |
485 error ("gallery: X be an integer or a vector for cauchy matrix."); | |
486 endif | |
487 | |
488 if (nargin == 1) | |
489 y = x; | |
490 endif | |
491 | |
492 ## Ensure x and y are column vectors | |
493 x = x(:); | |
494 y = y(:); | |
495 if (numel (x) != numel (y)) | |
496 error ("gallery: X and Y must be vectors of same length for cauchy matrix."); | |
497 endif | |
498 | |
499 C = x * ones (1, n) + ones (n, 1) * y.'; | |
500 C = ones (n) ./ C; | |
501 endfunction | |
502 | |
503 function C = chebspec (n, k = 0) | |
504 ## CHEBSPEC Chebyshev spectral differentiation matrix. | |
505 ## C = CHEBSPEC(N, K) is a Chebyshev spectral differentiation | |
506 ## matrix of order N. K = 0 (the default) or 1. | |
507 ## For K = 0 (`no boundary conditions'), C is nilpotent, with | |
508 ## C^N = 0 and it has the null vector ONES(N,1). | |
509 ## C is similar to a Jordan block of size N with eigenvalue zero. | |
510 ## For K = 1, C is nonsingular and well-conditioned, and its eigenvalues | |
511 ## have negative real parts. | |
512 ## For both K, the computed eigenvector matrix X from EIG is | |
513 ## ill-conditioned (MESH(REAL(X)) is interesting). | |
514 ## | |
515 ## References: | |
516 ## C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zang, Spectral | |
517 ## Methods in Fluid Dynamics, Springer-Verlag, Berlin, 1988; p. 69. | |
518 ## L.N. Trefethen and M.R. Trummer, An instability phenomenon in | |
519 ## spectral methods, SIAM J. Numer. Anal., 24 (1987), pp. 1008-1023. | |
520 ## D. Funaro, Computing the inverse of the Chebyshev collocation | |
521 ## derivative, SIAM J. Sci. Stat. Comput., 9 (1988), pp. 1050-1057. | |
522 | |
523 if (nargin < 1 || nargin > 2) | |
524 error ("gallery: 1 to 2 arguments are required for chebspec matrix."); | |
525 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
526 error ("gallery: N must be an integer for chebspec matrix."); | |
527 elseif (! isnumeric (k) || ! isscalar (k)) | |
528 error ("gallery: K must be a scalar for chebspec matrix."); | |
529 endif | |
530 | |
531 ## k = 1 case obtained from k = 0 case with one bigger n. | |
532 switch (k) | |
533 case (0), # do nothing | |
534 case (1), n = n + 1; | |
535 otherwise | |
16766
7268845c0a1e
avoid backquote in error messages, some uses in doc strings
John W. Eaton <jwe@octave.org>
parents:
16734
diff
changeset
|
536 error ("gallery: unknown K '%d' for chebspec matrix.", k); |
16634 | 537 endswitch |
538 | |
539 n = n-1; | |
540 C = zeros (n+1); | |
541 | |
542 one = ones (n+1, 1); | |
543 x = cos ((0:n)' * (pi/n)); | |
544 d = ones (n+1, 1); | |
545 d(1) = 2; | |
546 d(n+1) = 2; | |
547 | |
548 ## eye(size(C)) on next line avoids div by zero. | |
549 C = (d * (one./d)') ./ (x*one'-one*x' + eye (size (C))); | |
550 | |
551 ## Now fix diagonal and signs. | |
552 C(1,1) = (2*n^2+1)/6; | |
553 for i = 2:n+1 | |
554 if (rem (i, 2) == 0) | |
555 C(:,i) = -C(:,i); | |
556 C(i,:) = -C(i,:); | |
557 endif | |
558 if (i < n+1) | |
559 C(i,i) = -x(i)/(2*(1-x(i)^2)); | |
560 else | |
561 C(n+1,n+1) = -C(1,1); | |
562 endif | |
563 endfor | |
564 | |
565 if (k == 1) | |
566 C = C(2:n+1,2:n+1); | |
567 endif | |
568 endfunction | |
569 | |
570 function C = chebvand (m, p) | |
571 ## CHEBVAND Vandermonde-like matrix for the Chebyshev polynomials. | |
572 ## C = CHEBVAND(P), where P is a vector, produces the (primal) | |
573 ## Chebyshev Vandermonde matrix based on the points P, | |
574 ## i.e., C(i,j) = T_{i-1}(P(j)), where T_{i-1} is the Chebyshev | |
575 ## polynomial of degree i-1. | |
576 ## CHEBVAND(M,P) is a rectangular version of CHEBVAND(P) with M rows. | |
577 ## Special case: If P is a scalar then P equally spaced points on | |
578 ## [0,1] are used. | |
579 ## | |
580 ## Reference: | |
581 ## N.J. Higham, Stability analysis of algorithms for solving confluent | |
582 ## Vandermonde-like systems, SIAM J. Matrix Anal. Appl., 11 (1990), | |
583 ## pp. 23-41. | |
584 | |
585 if (nargin < 1 || nargin > 2) | |
586 error ("gallery: 1 or 2 arguments are required for chebvand matrix."); | |
587 endif | |
588 | |
589 ## because the order of the arguments changes if nargin is 1 or 2 ... | |
590 | |
591 if (nargin == 1) | |
592 p = m; | |
593 endif | |
594 | |
595 n = numel (p); | |
596 if (! isnumeric (p)) | |
597 error ("gallery: P must be numeric for chebvand matrix."); | |
598 elseif (isscalar (p) && fix (p) == p) | |
599 n = p; | |
600 p = linspace (0, 1, n); | |
601 elseif (n > 1 && isvector (p)) | |
602 ## do nothing | |
603 endif | |
604 p = p(:).'; # Ensure p is a row vector. | |
605 | |
606 if (nargin == 1) | |
607 m = n; | |
608 elseif (! isnumeric (m) || ! isscalar (m)) | |
609 error ("gallery: M must be a scalar for chebvand matrix."); | |
610 endif | |
611 | |
612 C = ones (m, n); | |
613 if (m != 1) | |
614 C(2,:) = p; | |
615 ## Use Chebyshev polynomial recurrence. | |
616 for i = 3:m | |
617 C(i,:) = 2.*p.*C(i-1,:) - C(i-2,:); | |
618 endfor | |
619 endif | |
620 endfunction | |
621 | |
622 function A = chow (n, alpha = 1, delta = 0) | |
623 ## CHOW Chow matrix - a singular Toeplitz lower Hessenberg matrix. | |
624 ## A = CHOW(N, ALPHA, DELTA) is a Toeplitz lower Hessenberg matrix | |
625 ## A = H(ALPHA) + DELTA*EYE, where H(i,j) = ALPHA^(i-j+1). | |
626 ## H(ALPHA) has p = FLOOR(N/2) zero eigenvalues, the rest being | |
627 ## 4*ALPHA*COS( k*PI/(N+2) )^2, k=1:N-p. | |
628 ## Defaults: ALPHA = 1, DELTA = 0. | |
629 ## | |
630 ## References: | |
631 ## T.S. Chow, A class of Hessenberg matrices with known | |
632 ## eigenvalues and inverses, SIAM Review, 11 (1969), pp. 391-395. | |
633 ## G. Fairweather, On the eigenvalues and eigenvectors of a class of | |
634 ## Hessenberg matrices, SIAM Review, 13 (1971), pp. 220-221. | |
635 | |
636 if (nargin < 1 || nargin > 3) | |
637 error ("gallery: 1 to 3 arguments are required for chow matrix."); | |
638 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
639 error ("gallery: N must be an integer for chow matrix."); | |
640 elseif (! isnumeric (alpha) || ! isscalar (alpha)) | |
641 error ("gallery: ALPHA must be a scalar for chow matrix."); | |
642 elseif (! isnumeric (delta) || ! isscalar (delta)) | |
643 error ("gallery: DELTA must be a scalar for chow matrix."); | |
644 endif | |
645 | |
646 A = toeplitz (alpha.^(1:n), [alpha 1 zeros(1, n-2)]) + delta * eye (n); | |
647 endfunction | |
648 | |
649 function C = circul (v) | |
650 ## CIRCUL Circulant matrix. | |
651 ## C = CIRCUL(V) is the circulant matrix whose first row is V. | |
652 ## (A circulant matrix has the property that each row is obtained | |
653 ## from the previous one by cyclically permuting the entries one step | |
654 ## forward; it is a special Toeplitz matrix in which the diagonals | |
655 ## `wrap round'.) | |
656 ## Special case: if V is a scalar then C = CIRCUL(1:V). | |
657 ## The eigensystem of C (N-by-N) is known explicitly. If t is an Nth | |
658 ## root of unity, then the inner product of V with W = [1 t t^2 ... t^N] | |
659 ## is an eigenvalue of C, and W(N:-1:1) is an eigenvector of C. | |
660 ## | |
661 ## Reference: | |
662 ## P.J. Davis, Circulant Matrices, John Wiley, 1977. | |
663 | |
664 if (nargin != 1) | |
665 error ("gallery: 1 argument is required for circul matrix."); | |
666 elseif (! isnumeric (v)) | |
667 error ("gallery: V must be numeric for circul matrix."); | |
668 endif | |
669 | |
16734
67b67fc0969a
gallery: fix bug from typo in variable name
Carnë Draug <carandraug@octave.org>
parents:
16634
diff
changeset
|
670 n = numel (v); |
67b67fc0969a
gallery: fix bug from typo in variable name
Carnë Draug <carandraug@octave.org>
parents:
16634
diff
changeset
|
671 if (isscalar (v) && fix (v) == v) |
16634 | 672 n = v; |
673 v = 1:n; | |
16734
67b67fc0969a
gallery: fix bug from typo in variable name
Carnë Draug <carandraug@octave.org>
parents:
16634
diff
changeset
|
674 elseif (n > 1 && isvector (v)) |
16634 | 675 ## do nothing |
676 else | |
677 error ("gallery: X must be a scalar or a vector for circul matrix."); | |
678 endif | |
679 | |
680 v = v(:).'; # Make sure v is a row vector | |
681 C = toeplitz ([v(1) v(n:-1:2)], v); | |
682 endfunction | |
683 | |
684 function A = clement (n, k = 0) | |
685 ## CLEMENT Clement matrix - tridiagonal with zero diagonal entries. | |
686 ## CLEMENT(N, K) is a tridiagonal matrix with zero diagonal entries | |
687 ## and known eigenvalues. It is singular if N is odd. About 64 | |
688 ## percent of the entries of the inverse are zero. The eigenvalues | |
689 ## are plus and minus the numbers N-1, N-3, N-5, ..., (1 or 0). | |
690 ## For K = 0 (the default) the matrix is unsymmetric, while for | |
691 ## K = 1 it is symmetric. | |
692 ## CLEMENT(N, 1) is diagonally similar to CLEMENT(N). | |
693 ## | |
694 ## Similar properties hold for TRIDIAG(X,Y,Z) where Y = ZEROS(N,1). | |
695 ## The eigenvalues still come in plus/minus pairs but they are not | |
696 ## known explicitly. | |
697 ## | |
698 ## References: | |
699 ## P.A. Clement, A class of triple-diagonal matrices for test | |
700 ## purposes, SIAM Review, 1 (1959), pp. 50-52. | |
701 ## A. Edelman and E. Kostlan, The road from Kac's matrix to Kac's | |
702 ## random polynomials. In John~G. Lewis, editor, Proceedings of | |
703 ## the Fifth SIAM Conference on Applied Linear Algebra Society | |
704 ## for Industrial and Applied Mathematics, Philadelphia, 1994, | |
705 ## pp. 503-507. | |
706 ## O. Taussky and J. Todd, Another look at a matrix of Mark Kac, | |
707 ## Linear Algebra and Appl., 150 (1991), pp. 341-360. | |
708 | |
709 if (nargin < 1 || nargin > 2) | |
710 error ("gallery: 1 or 2 arguments are required for clement matrix."); | |
711 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
712 error ("gallery: N must be an integer for clement matrix."); | |
713 elseif (! isnumeric (k) || ! isscalar (k)) | |
714 error ("gallery: K must be a numeric scalar for clement matrix."); | |
715 endif | |
716 | |
717 n = n-1; | |
718 x = n:-1:1; | |
719 z = 1:n; | |
720 | |
721 if (k == 0) | |
722 A = diag (x, -1) + diag (z, 1); | |
723 elseif (k == 1) | |
724 y = sqrt (x.*z); | |
725 A = diag (y, -1) + diag (y, 1); | |
726 else | |
727 error ("gallery: K must have a value of 0 or 1 for clement matrix."); | |
728 endif | |
729 endfunction | |
730 | |
731 function C = compar (A, k = 0) | |
732 ## COMP Comparison matrices. | |
733 ## COMP(A) is DIAG(B) - TRIL(B,-1) - TRIU(B,1), where B = ABS(A). | |
734 ## COMP(A, 1) is A with each diagonal element replaced by its | |
735 ## absolute value, and each off-diagonal element replaced by minus | |
736 ## the absolute value of the largest element in absolute value in | |
737 ## its row. However, if A is triangular COMP(A, 1) is too. | |
738 ## COMP(A, 0) is the same as COMP(A). | |
739 ## COMP(A) is often denoted by M(A) in the literature. | |
740 ## | |
741 ## Reference (e.g.): | |
742 ## N.J. Higham, A survey of condition number estimation for | |
743 ## triangular matrices, SIAM Review, 29 (1987), pp. 575-596. | |
744 | |
745 if (nargin < 1 || nargin > 2) | |
746 error ("gallery: 1 or 2 arguments are required for compar matrix."); | |
747 elseif (! isnumeric (A) || ndims (A) != 2) | |
748 error ("gallery: A must be a 2D matrix for compar matrix."); | |
749 elseif (! isnumeric (k) || ! isscalar (k)) | |
750 error ("gallery: K must be a numeric scalar for compar matrix."); | |
751 endif | |
752 | |
753 [m, n] = size (A); | |
754 p = min (m, n); | |
755 | |
756 if (k == 0) | |
757 ## This code uses less temporary storage than | |
758 ## the `high level' definition above. | |
759 C = -abs (A); | |
760 for j = 1:p | |
761 C(j,j) = abs (A(j,j)); | |
762 endfor | |
763 | |
764 elseif (k == 1) | |
765 C = A'; | |
766 for j = 1:p | |
767 C(k,k) = 0; | |
768 endfor | |
769 mx = max (abs (C)); | |
770 C = -mx'*ones (1, n); | |
771 for j = 1:p | |
772 C(j,j) = abs (A(j,j)); | |
773 endfor | |
774 if (all (A == tril (A))), C = tril (C); endif | |
775 if (all (A == triu (A))), C = triu (C); endif | |
776 | |
777 else | |
778 error ("gallery: K must have a value of 0 or 1 for compar matrix."); | |
779 endif | |
780 | |
781 endfunction | |
782 | |
783 function A = condex (n, k = 4, theta = 100) | |
784 ## CONDEX `Counterexamples' to matrix condition number estimators. | |
785 ## CONDEX(N, K, THETA) is a `counterexample' matrix to a condition | |
786 ## estimator. It has order N and scalar parameter THETA (default 100). | |
787 ## If N is not equal to the `natural' size of the matrix then | |
788 ## the matrix is padded out with an identity matrix to order N. | |
789 ## The matrix, its natural size, and the estimator to which it applies | |
790 ## are specified by K (default K = 4) as follows: | |
791 ## K = 1: 4-by-4, LINPACK (RCOND) | |
792 ## K = 2: 3-by-3, LINPACK (RCOND) | |
793 ## K = 3: arbitrary, LINPACK (RCOND) (independent of THETA) | |
794 ## K = 4: N >= 4, SONEST (Higham 1988) | |
795 ## (Note that in practice the K = 4 matrix is not usually a | |
796 ## counterexample because of the rounding errors in forming it.) | |
797 ## | |
798 ## References: | |
799 ## A.K. Cline and R.K. Rew, A set of counter-examples to three | |
800 ## condition number estimators, SIAM J. Sci. Stat. Comput., | |
801 ## 4 (1983), pp. 602-611. | |
802 ## N.J. Higham, FORTRAN codes for estimating the one-norm of a real or | |
803 ## complex matrix, with applications to condition estimation | |
804 ## (Algorithm 674), ACM Trans. Math. Soft., 14 (1988), pp. 381-396. | |
805 | |
806 if (nargin < 1 || nargin > 3) | |
807 error ("gallery: 1 to 3 arguments are required for condex matrix."); | |
808 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
809 error ("gallery: N must be an integer for condex matrix."); | |
810 elseif (! isnumeric (k) || ! isscalar (k)) | |
811 error ("gallery: K must be a numeric scalar for condex matrix."); | |
812 elseif (! isnumeric (theta) || ! isscalar (theta)) | |
813 error ("gallery: THETA must be a numeric scalar for condex matrix."); | |
814 endif | |
815 | |
816 if (k == 1) # Cline and Rew (1983), Example B. | |
817 A = [1 -1 -2*theta 0 | |
818 0 1 theta -theta | |
819 0 1 1+theta -(theta+1) | |
820 0 0 0 theta]; | |
821 | |
822 elseif (k == 2) # Cline and Rew (1983), Example C. | |
823 A = [1 1-2/theta^2 -2 | |
824 0 1/theta -1/theta | |
825 0 0 1]; | |
826 | |
827 elseif (k == 3) # Cline and Rew (1983), Example D. | |
828 A = gallery ("triw", n, -1)'; | |
829 A(n,n) = -1; | |
830 | |
831 elseif (k == 4) # Higham (1988), p. 390. | |
832 x = ones (n, 3); # First col is e | |
833 x(2:n,2) = zeros (n-1, 1); # Second col is e(1) | |
834 | |
835 ## Third col is special vector b in SONEST | |
836 x(:, 3) = (-1).^[0:n-1]' .* ( 1 + [0:n-1]'/(n-1) ); | |
837 | |
838 Q = orth (x); # Q*Q' is now the orthogonal projector onto span(e(1),e,b)). | |
839 P = eye (n) - Q*Q'; | |
840 A = eye (n) + theta*P; | |
841 | |
842 else | |
16766
7268845c0a1e
avoid backquote in error messages, some uses in doc strings
John W. Eaton <jwe@octave.org>
parents:
16734
diff
changeset
|
843 error ("gallery: unknown estimator K '%d' for condex matrix.", k); |
16634 | 844 endif |
845 | |
846 ## Pad out with identity as necessary. | |
847 m = columns (A); | |
848 if (m < n) | |
849 for i = n:-1:m+1 | |
850 A(i,i) = 1; | |
851 endfor | |
852 endif | |
853 endfunction | |
854 | |
855 function A = cycol (n, k) | |
856 ## CYCOL Matrix whose columns repeat cyclically. | |
857 ## A = CYCOL([M N], K) is an M-by-N matrix of the form A = B(1:M,1:N) | |
858 ## where B = [C C C...] and C = RANDN(M, K). Thus A's columns repeat | |
859 ## cyclically, and A has rank at most K. K need not divide N. | |
860 ## K defaults to ROUND(N/4). | |
861 ## CYCOL(N, K), where N is a scalar, is the same as CYCOL([N N], K). | |
862 ## | |
863 ## This type of matrix can lead to underflow problems for Gaussian | |
864 ## elimination: see NA Digest Volume 89, Issue 3 (January 22, 1989). | |
865 | |
866 if (nargin < 1 || nargin > 2) | |
867 error ("gallery: 1 or 2 arguments are required for cycol matrix."); | |
868 elseif (! isnumeric (n) || all (numel (n) != [1 2]) || fix (n) != n) | |
869 error ("gallery: N must be a 1 or 2 element integer for cycol matrix."); | |
870 elseif (! isnumeric (k) || ! isscalar (k)) | |
871 error ("gallery: K must be a scalar for cycol matrix."); | |
872 endif | |
873 | |
874 ## Parameter n specifies dimension: m-by-n | |
875 m = n(1); | |
876 n = n(end); | |
877 | |
878 if (nargin < 2) | |
879 k = max (round (n/4), 1); | |
880 endif | |
881 | |
882 A = randn (m, k); | |
883 for i = 2:ceil (n/k) | |
884 A = [A A(:,1:k)]; | |
885 endfor | |
886 A = A(:,1:n); | |
887 endfunction | |
888 | |
889 function [c, d, e] = dorr (n, theta = 0.01) | |
890 ## DORR Dorr matrix - diagonally dominant, ill conditioned, tridiagonal. | |
891 ## [C, D, E] = DORR(N, THETA) returns the vectors defining a row diagonally | |
892 ## dominant, tridiagonal M-matrix that is ill conditioned for small | |
893 ## values of the parameter THETA >= 0. | |
894 ## If only one output parameter is supplied then | |
895 ## C = FULL(TRIDIAG(C,D,E)), i.e., the matrix iself is returned. | |
896 ## The columns of INV(C) vary greatly in norm. THETA defaults to 0.01. | |
897 ## The amount of diagonal dominance is given by (ignoring rounding errors): | |
898 ## COMP(C)*ONES(N,1) = THETA*(N+1)^2 * [1 0 0 ... 0 1]'. | |
899 ## | |
900 ## Reference: | |
901 ## F.W. Dorr, An example of ill-conditioning in the numerical | |
902 ## solution of singular perturbation problems, Math. Comp., 25 (1971), | |
903 ## pp. 271-283. | |
904 | |
905 if (nargin < 1 || nargin > 2) | |
906 error ("gallery: 1 or 2 arguments are required for dorr matrix."); | |
907 elseif (! isscalar (n) || ! isnumeric (n) || fix (n) != n) | |
908 error ("gallery: N must be an integer for dorr matrix."); | |
909 elseif (! isscalar (theta) || ! isnumeric (theta)) | |
910 error ("gallery: THETA must be a numeric scalar for dorr matrix."); | |
911 endif | |
912 | |
913 c = zeros (n, 1); | |
914 e = c; | |
915 d = c; | |
916 ## All length n for convenience. Make c, e of length n-1 later. | |
917 | |
918 h = 1/(n+1); | |
919 m = floor ((n+1)/2); | |
920 term = theta/h^2; | |
921 | |
922 i = (1:m)'; | |
923 c(i) = -term * ones (m, 1); | |
924 e(i) = c(i) - (0.5-i*h)/h; | |
925 d(i) = -(c(i) + e(i)); | |
926 | |
927 i = (m+1:n)'; | |
928 e(i) = -term * ones (n-m, 1); | |
929 c(i) = e(i) + (0.5-i*h)/h; | |
930 d(i) = -(c(i) + e(i)); | |
931 | |
932 c = c(2:n); | |
933 e = e(1:n-1); | |
934 | |
935 if (nargout <= 1) | |
936 c = tridiag (c, d, e); | |
937 endif | |
938 endfunction | |
939 | |
940 function A = dramadah (n, k = 1) | |
941 ## DRAMADAH A (0,1) matrix whose inverse has large integer entries. | |
942 ## An anti-Hadamard matrix A is a matrix with elements 0 or 1 for | |
943 ## which MU(A) := NORM(INV(A),'FRO') is maximal. | |
944 ## A = DRAMADAH(N, K) is an N-by-N (0,1) matrix for which MU(A) is | |
945 ## relatively large, although not necessarily maximal. | |
946 ## Available types (the default is K = 1): | |
947 ## K = 1: A is Toeplitz, with ABS(DET(A)) = 1, and MU(A) > c(1.75)^N, | |
948 ## where c is a constant. | |
949 ## K = 2: A is upper triangular and Toeplitz. | |
950 ## The inverses of both types have integer entries. | |
951 ## | |
952 ## Another interesting (0,1) matrix: | |
953 ## K = 3: A has maximal determinant among (0,1) lower Hessenberg | |
954 ## matrices: det(A) = the n'th Fibonacci number. A is Toeplitz. | |
955 ## The eigenvalues have an interesting distribution in the complex | |
956 ## plane. | |
957 ## | |
958 ## References: | |
959 ## R.L. Graham and N.J.A. Sloane, Anti-Hadamard matrices, | |
960 ## Linear Algebra and Appl., 62 (1984), pp. 113-137. | |
961 ## L. Ching, The maximum determinant of an nxn lower Hessenberg | |
962 ## (0,1) matrix, Linear Algebra and Appl., 183 (1993), pp. 147-153. | |
963 | |
964 if (nargin < 1 || nargin > 2) | |
965 error ("gallery: 1 to 2 arguments are required for dramadah matrix."); | |
966 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
967 error ("gallery: N must be an integer for dramadah matrix."); | |
968 elseif (! isnumeric (k) || ! isscalar (k)) | |
969 error ("gallery: K must be a numeric scalar for dramadah matrix."); | |
970 endif | |
971 | |
972 switch (k) | |
973 case (1) # Toeplitz | |
974 c = ones (n, 1); | |
975 for i = 2:4:n | |
976 m = min (1, n-i); | |
977 c(i:i+m) = zeros (m+1, 1); | |
978 endfor | |
979 r = zeros (n, 1); | |
980 r(1:4) = [1 1 0 1]; | |
981 if (n < 4) | |
982 r = r(1:n); | |
983 endif | |
984 A = toeplitz (c, r); | |
985 | |
986 case (2) # Upper triangular and Toeplitz | |
987 c = zeros (n, 1); | |
988 c(1) = 1; | |
989 r = ones (n, 1); | |
990 for i= 3:2:n | |
991 r(i) = 0; | |
992 endfor | |
993 A = toeplitz (c, r); | |
994 | |
995 case (3) # Lower Hessenberg | |
996 c = ones (n, 1); | |
997 for i= 2:2:n | |
998 c(i) = 0; | |
999 endfor | |
1000 A = toeplitz (c, [1 1 zeros(1,n-2)]); | |
1001 | |
1002 otherwise | |
16766
7268845c0a1e
avoid backquote in error messages, some uses in doc strings
John W. Eaton <jwe@octave.org>
parents:
16734
diff
changeset
|
1003 error ("gallery: unknown K '%d' for dramadah matrix.", k); |
16634 | 1004 endswitch |
1005 endfunction | |
1006 | |
1007 function A = fiedler (c) | |
1008 ## FIEDLER Fiedler matrix - symmetric. | |
1009 ## A = FIEDLER(C), where C is an n-vector, is the n-by-n symmetric | |
1010 ## matrix with elements ABS(C(i)-C(j)). | |
1011 ## Special case: if C is a scalar, then A = FIEDLER(1:C) | |
1012 ## (i.e. A(i,j) = ABS(i-j)). | |
1013 ## Properties: | |
1014 ## FIEDLER(N) has a dominant positive eigenvalue and all the other | |
1015 ## eigenvalues are negative (Szego, 1936). | |
1016 ## Explicit formulas for INV(A) and DET(A) are given by Todd (1977) | |
1017 ## and attributed to Fiedler. These indicate that INV(A) is | |
1018 ## tridiagonal except for nonzero (1,n) and (n,1) elements. | |
1019 ## [I think these formulas are valid only if the elements of | |
1020 ## C are in increasing or decreasing order---NJH.] | |
1021 ## | |
1022 ## References: | |
1023 ## G. Szego, Solution to problem 3705, Amer. Math. Monthly, | |
1024 ## 43 (1936), pp. 246-259. | |
1025 ## J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra, | |
1026 ## Birkhauser, Basel, and Academic Press, New York, 1977, p. 159. | |
1027 | |
1028 if (nargin != 1) | |
1029 error ("gallery: 1 argument is required for fiedler matrix."); | |
1030 elseif (! isnumeric (c)) | |
1031 error ("gallery: C must be numeric for fiedler matrix."); | |
1032 endif | |
1033 | |
1034 n = numel (c); | |
1035 if (isscalar (c) && fix (c) == c) | |
1036 n = c; | |
1037 c = 1:n; | |
1038 elseif (n > 1 && isvector (c)) | |
1039 ## do nothing | |
1040 else | |
1041 error ("gallery: C must be an integer or a vector for fiedler matrix."); | |
1042 endif | |
1043 c = c(:).'; # Ensure c is a row vector. | |
1044 | |
1045 A = ones (n, 1) * c; | |
1046 A = abs (A - A.'); # NB. array transpose. | |
1047 endfunction | |
1048 | |
1049 function A = forsythe (n, alpha = sqrt (eps), lambda = 0) | |
1050 ## FORSYTHE Forsythe matrix - a perturbed Jordan block. | |
1051 ## FORSYTHE(N, ALPHA, LAMBDA) is the N-by-N matrix equal to | |
1052 ## JORDBLOC(N, LAMBDA) except it has an ALPHA in the (N,1) position. | |
1053 ## It has the characteristic polynomial | |
1054 ## DET(A-t*EYE) = (LAMBDA-t)^N - (-1)^N ALPHA. | |
1055 ## ALPHA defaults to SQRT(EPS) and LAMBDA to 0. | |
1056 | |
1057 if (nargin < 1 || nargin > 3) | |
1058 error ("gallery: 1 to 3 arguments are required for forsythe matrix."); | |
1059 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1060 error ("gallery: N must be an integer for forsythe matrix."); | |
1061 elseif (! isnumeric (alpha) || ! isscalar (alpha)) | |
1062 error ("gallery: ALPHA must be a numeric scalar for forsythe matrix."); | |
1063 elseif (! isnumeric (lambda) || ! isscalar (lambda)) | |
1064 error ("gallery: LAMBDA must be a numeric scalar for forsythe matrix."); | |
1065 endif | |
1066 | |
1067 A = jordbloc (n, lambda); | |
1068 A(n,1) = alpha; | |
1069 endfunction | |
1070 | |
1071 function F = frank (n, k = 0) | |
1072 ## FRANK Frank matrix---ill conditioned eigenvalues. | |
1073 ## F = FRANK(N, K) is the Frank matrix of order N. It is upper | |
1074 ## Hessenberg with determinant 1. K = 0 is the default; if K = 1 the | |
1075 ## elements are reflected about the anti-diagonal (1,N)--(N,1). | |
1076 ## F has all positive eigenvalues and they occur in reciprocal pairs | |
1077 ## (so that 1 is an eigenvalue if N is odd). | |
1078 ## The eigenvalues of F may be obtained in terms of the zeros of the | |
1079 ## Hermite polynomials. | |
1080 ## The FLOOR(N/2) smallest eigenvalues of F are ill conditioned, | |
1081 ## the more so for bigger N. | |
1082 ## | |
1083 ## DET(FRANK(N)') comes out far from 1 for large N---see Frank (1958) | |
1084 ## and Wilkinson (1960) for discussions. | |
1085 ## | |
1086 ## This version incorporates improvements suggested by W. Kahan. | |
1087 ## | |
1088 ## References: | |
1089 ## W.L. Frank, Computing eigenvalues of complex matrices by determinant | |
1090 ## evaluation and by methods of Danilewski and Wielandt, J. Soc. | |
1091 ## Indust. Appl. Math., 6 (1958), pp. 378-392 (see pp. 385, 388). | |
1092 ## G.H. Golub and J.H. Wilkinson, Ill-conditioned eigensystems and the | |
1093 ## computation of the Jordan canonical form, SIAM Review, 18 (1976), | |
1094 ## pp. 578-619 (Section 13). | |
1095 ## H. Rutishauser, On test matrices, Programmation en Mathematiques | |
1096 ## Numeriques, Editions Centre Nat. Recherche Sci., Paris, 165, | |
1097 ## 1966, pp. 349-365. Section 9. | |
1098 ## J.H. Wilkinson, Error analysis of floating-point computation, | |
1099 ## Numer. Math., 2 (1960), pp. 319-340 (Section 8). | |
1100 ## J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University | |
1101 ## Press, 1965 (pp. 92-93). | |
1102 ## The next two references give details of the eigensystem, as does | |
1103 ## Rutishauser (see above). | |
1104 ## P.J. Eberlein, A note on the matrices denoted by B_n, SIAM J. Appl. | |
1105 ## Math., 20 (1971), pp. 87-92. | |
1106 ## J.M. Varah, A generalization of the Frank matrix, SIAM J. Sci. Stat. | |
1107 ## Comput., 7 (1986), pp. 835-839. | |
1108 | |
1109 if (nargin < 1 || nargin > 2) | |
1110 error ("gallery: 1 to 2 arguments are required for frank matrix."); | |
1111 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1112 error ("gallery: N must be an integer for frank matrix."); | |
1113 elseif (! isnumeric (k) || ! isscalar (k)) | |
1114 error ("gallery: K must be a numeric scalar for frank matrix."); | |
1115 endif | |
1116 | |
1117 p = n:-1:1; | |
1118 F = triu (p(ones (n, 1), :) - diag (ones (n-1, 1), -1), -1); | |
1119 | |
1120 switch (k) | |
1121 case (0), # do nothing | |
1122 case (1), F = F(p,p)'; | |
1123 otherwise | |
1124 error ("gallery: K must have a value of 0 or 1 for frank matrix."); | |
1125 endswitch | |
1126 endfunction | |
1127 | |
1128 function c = gcdmat (n) | |
1129 if (nargin != 1) | |
1130 error ("gallery: 1 argument is required for gcdmat matrix."); | |
1131 elseif (! isscalar (n) || ! isnumeric (n) || fix (n) != n) | |
1132 error ("gallery: N must be an integer for gcdmat matrix."); | |
1133 endif | |
1134 c = gcd (repmat ((1:n)', [1 n]), repmat (1:n, [n 1])); | |
1135 endfunction | |
1136 | |
1137 function A = gearmat (n, i = n, j = -n) | |
1138 ## NOTE: this function was named gearm in the original Test Matrix Toolbox | |
1139 ## GEARMAT Gear matrix. | |
1140 ## A = GEARMAT(N,I,J) is the N-by-N matrix with ones on the sub- and | |
1141 ## super-diagonals, SIGN(I) in the (1,ABS(I)) position, SIGN(J) | |
1142 ## in the (N,N+1-ABS(J)) position, and zeros everywhere else. | |
1143 ## Defaults: I = N, j = -N. | |
1144 ## All eigenvalues are of the form 2*COS(a) and the eigenvectors | |
1145 ## are of the form [SIN(w+a), SIN(w+2a), ..., SIN(w+Na)]. | |
1146 ## The values of a and w are given in the reference below. | |
1147 ## A can have double and triple eigenvalues and can be defective. | |
1148 ## GEARMAT(N) is singular. | |
1149 ## | |
1150 ## (GEAR is a Simulink function, hence GEARMAT for Gear matrix.) | |
1151 ## Reference: | |
1152 ## C.W. Gear, A simple set of test matrices for eigenvalue programs, | |
1153 ## Math. Comp., 23 (1969), pp. 119-125. | |
1154 | |
1155 if (nargin < 1 || nargin > 3) | |
1156 error ("gallery: 1 to 3 arguments are required for gearmat matrix."); | |
1157 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1158 error ("gallery: N must be an integer for gearmat matrix."); | |
1159 elseif (! isnumeric (i) || ! isscalar (i) || i == 0 || abs (i) <= n) | |
1160 error ("gallery: I must be a non-zero scalar, and abs (I) <= N for gearmat matrix."); | |
1161 elseif (! isnumeric (j) || ! isscalar (j) || i == 0 || abs (j) <= n) | |
1162 error ("gallery: J must be a non-zero scalar, and abs (J) <= N for gearmat matrix."); | |
1163 endif | |
1164 | |
1165 A = diag (ones (n-1, 1), -1) + diag (ones (n-1, 1), 1); | |
1166 A(1, abs (i)) = sign (i); | |
1167 A(n, n+1 - abs (j)) = sign (j); | |
1168 endfunction | |
1169 | |
1170 function G = grcar (n, k = 3) | |
1171 ## GRCAR Grcar matrix - a Toeplitz matrix with sensitive eigenvalues. | |
1172 ## GRCAR(N, K) is an N-by-N matrix with -1s on the | |
1173 ## subdiagonal, 1s on the diagonal, and K superdiagonals of 1s. | |
1174 ## The default is K = 3. The eigenvalues of this matrix form an | |
1175 ## interesting pattern in the complex plane (try PS(GRCAR(32))). | |
1176 ## | |
1177 ## References: | |
1178 ## J.F. Grcar, Operator coefficient methods for linear equations, | |
1179 ## Report SAND89-8691, Sandia National Laboratories, Albuquerque, | |
1180 ## New Mexico, 1989 (Appendix 2). | |
1181 ## N.M. Nachtigal, L. Reichel and L.N. Trefethen, A hybrid GMRES | |
1182 ## algorithm for nonsymmetric linear systems, SIAM J. Matrix Anal. | |
1183 ## Appl., 13 (1992), pp. 796-825. | |
1184 | |
1185 if (nargin < 1 || nargin > 2) | |
1186 error ("gallery: 1 to 2 arguments are required for grcar matrix."); | |
1187 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1188 error ("gallery: N must be an integer for grcar matrix."); | |
1189 elseif (! isnumeric (k) || ! isscalar (k)) | |
1190 error ("gallery: K must be a numeric scalar for grcar matrix."); | |
1191 endif | |
1192 | |
1193 G = tril (triu (ones (n)), k) - diag (ones (n-1, 1), -1); | |
1194 endfunction | |
1195 | |
1196 function A = hanowa (n, d = -1) | |
1197 ## HANOWA A matrix whose eigenvalues lie on a vertical line in the complex plane. | |
1198 ## HANOWA(N, d) is the N-by-N block 2x2 matrix (thus N = 2M must be even) | |
1199 ## [d*EYE(M) -DIAG(1:M) | |
1200 ## DIAG(1:M) d*EYE(M)] | |
1201 ## It has complex eigenvalues lambda(k) = d +/- k*i (1 <= k <= M). | |
1202 ## Parameter d defaults to -1. | |
1203 ## | |
1204 ## Reference: | |
1205 ## E. Hairer, S.P. Norsett and G. Wanner, Solving Ordinary | |
1206 ## Differential Equations I: Nonstiff Problems, Springer-Verlag, | |
1207 ## Berlin, 1987. (pp. 86-87) | |
1208 | |
1209 if (nargin < 1 || nargin > 2) | |
1210 error ("gallery: 1 to 2 arguments are required for hanowa matrix."); | |
1211 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1212 error ("gallery: N must be an integer for hanowa matrix."); | |
1213 elseif (rem (n, 2) != 0) | |
1214 error ("gallery: N must be even for hanowa matrix."); | |
1215 elseif (! isnumeric (lambda) || ! isscalar (lambda)) | |
1216 error ("gallery: D must be a numeric scalar for hanowa matrix."); | |
1217 endif | |
1218 | |
1219 m = n/2; | |
1220 A = [ d*eye(m) -diag(1:m) | |
1221 diag(1:m) d*eye(m) ]; | |
1222 endfunction | |
1223 | |
1224 function [v, beta] = house (x) | |
1225 ## HOUSE Householder matrix. | |
1226 ## If [v, beta] = HOUSE(x) then H = EYE - beta*v*v' is a Householder | |
1227 ## matrix such that Hx = -sign(x(1))*norm(x)*e_1. | |
1228 ## NB: If x = 0 then v = 0, beta = 1 is returned. | |
1229 ## x can be real or complex. | |
1230 ## sign(x) := exp(i*arg(x)) ( = x./abs(x) when x ~= 0). | |
1231 ## | |
1232 ## Theory: (textbook references Golub & Van Loan 1989, 38-43; | |
1233 ## Stewart 1973, 231-234, 262; Wilkinson 1965, 48-50). | |
1234 ## Hx = y: (I - beta*v*v')x = -s*e_1. | |
1235 ## Must have |s| = norm(x), v = x+s*e_1, and | |
1236 ## x'y = x'Hx =(x'Hx)' real => arg(s) = arg(x(1)). | |
1237 ## So take s = sign(x(1))*norm(x) (which avoids cancellation). | |
1238 ## v'v = (x(1)+s)^2 + x(2)^2 + ... + x(n)^2 | |
1239 ## = 2*norm(x)*(norm(x) + |x(1)|). | |
1240 ## | |
1241 ## References: | |
1242 ## G.H. Golub and C.F. Van Loan, Matrix Computations, second edition, | |
1243 ## Johns Hopkins University Press, Baltimore, Maryland, 1989. | |
1244 ## G.W. Stewart, Introduction to Matrix Computations, Academic Press, | |
1245 ## New York, 1973, | |
1246 ## J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University | |
1247 ## Press, 1965. | |
1248 | |
1249 if (nargin != 1) | |
1250 error ("gallery: 1 argument is required for house matrix."); | |
1251 elseif (! isnumeric (x) || ! isvector (x) || numel (x) <= 1) | |
1252 error ("gallery: X must be a vector for house matrix."); | |
1253 endif | |
1254 | |
1255 ## must be a column vector | |
1256 x = x(:); | |
1257 | |
1258 s = norm (x) * (sign (x(1)) + (x(1) == 0)); # Modification for sign (0) == 1. | |
1259 v = x; | |
1260 if (s == 0) | |
1261 ## Quit if x is the zero vector. | |
1262 beta = 1; | |
1263 else | |
1264 v(1) = v(1) + s; | |
1265 beta = 1/(s'*v(1)); # NB the conjugated s. | |
1266 ## beta = 1/(abs (s) * (abs (s) +abs(x(1)) would guarantee beta real. | |
1267 ## But beta as above can be non-real (due to rounding) only when x is complex. | |
1268 endif | |
1269 endfunction | |
1270 | |
1271 function A = invhess (x, y) | |
1272 ## INVHESS Inverse of an upper Hessenberg matrix. | |
1273 ## INVHESS(X, Y), where X is an N-vector and Y an N-1 vector, | |
1274 ## is the matrix whose lower triangle agrees with that of | |
1275 ## ONES(N,1)*X' and whose strict upper triangle agrees with | |
1276 ## that of [1 Y]*ONES(1,N). | |
1277 ## The matrix is nonsingular if X(1) ~= 0 and X(i+1) ~= Y(i) | |
1278 ## for all i, and its inverse is an upper Hessenberg matrix. | |
1279 ## If Y is omitted it defaults to -X(1:N-1). | |
1280 ## Special case: if X is a scalar INVHESS(X) is the same as | |
1281 ## INVHESS(1:X). | |
1282 ## | |
1283 ## References: | |
1284 ## F.N. Valvi and V.S. Geroyannis, Analytic inverses and | |
1285 ## determinants for a class of matrices, IMA Journal of Numerical | |
1286 ## Analysis, 7 (1987), pp. 123-128. | |
1287 ## W.-L. Cao and W.J. Stewart, A note on inverses of Hessenberg-like | |
1288 ## matrices, Linear Algebra and Appl., 76 (1986), pp. 233-240. | |
1289 ## Y. Ikebe, On inverses of Hessenberg matrices, Linear Algebra and | |
1290 ## Appl., 24 (1979), pp. 93-97. | |
1291 ## P. Rozsa, On the inverse of band matrices, Integral Equations and | |
1292 ## Operator Theory, 10 (1987), pp. 82-95. | |
1293 | |
1294 if (nargin < 1 || nargin > 2) | |
1295 error ("gallery: 1 to 2 arguments are required for invhess matrix."); | |
1296 elseif (! isnumeric (x)) | |
1297 error ("gallery: X must be numeric for invhess matrix."); | |
1298 endif | |
1299 | |
1300 if (isscalar (x) && fix (x) == x) | |
1301 n = x; | |
1302 x = 1:n; | |
1303 elseif (! isscalar (x) && isvector (x)) | |
1304 n = numel (n); | |
1305 else | |
1306 error ("gallery: X must be an integer scalar, or a vector for invhess matrix."); | |
1307 endif | |
1308 | |
1309 if (nargin < 2) | |
1310 y = -x(1:end-1); | |
1311 elseif (! isvector (y) || numel (y) != numel (x) -1) | |
1312 error ("gallery: Y must be a vector of length -1 than X for invhess matrix."); | |
1313 endif | |
1314 | |
1315 x = x(:); | |
1316 y = y(:); | |
1317 | |
1318 ## On next line, z = x'; A = z(ones(n,1),:) would be more efficient. | |
1319 A = ones (n, 1) * x'; | |
1320 for j = 2:n | |
1321 A(1:j-1,j) = y(1:j-1); | |
1322 endfor | |
1323 endfunction | |
1324 | |
1325 function A = invol (n) | |
1326 ## INVOL An involutory matrix. | |
1327 ## A = INVOL(N) is an N-by-N involutory (A*A = EYE(N)) and | |
1328 ## ill-conditioned matrix. | |
1329 ## It is a diagonally scaled version of HILB(N). | |
1330 ## NB: B = (EYE(N)-A)/2 and B = (EYE(N)+A)/2 are idempotent (B*B = B). | |
1331 ## | |
1332 ## Reference: | |
1333 ## A.S. Householder and J.A. Carpenter, The singular values | |
1334 ## of involutory and of idempotent matrices, Numer. Math. 5 (1963), | |
1335 ## pp. 234-237. | |
1336 | |
1337 if (nargin != 1) | |
1338 error ("gallery: 1 argument is required for invol matrix."); | |
1339 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1340 error ("gallery: N must be an integer for invol matrix."); | |
1341 endif | |
1342 | |
1343 A = hilb (n); | |
1344 | |
1345 d = -n; | |
1346 A(:, 1) = d * A(:, 1); | |
1347 | |
1348 for i = 1:n-1 | |
1349 d = -(n+i)*(n-i)*d/(i*i); | |
1350 A(i+1,:) = d * A(i+1,:); | |
1351 endfor | |
1352 endfunction | |
1353 | |
1354 function [A, detA] = ipjfact (n, k = 0) | |
1355 ## IPJFACT A Hankel matrix with factorial elements. | |
1356 ## A = IPJFACT(N, K) is the matrix with | |
1357 ## A(i,j) = (i+j)! (K = 0, default) | |
1358 ## A(i,j) = 1/(i+j)! (K = 1) | |
1359 ## Both are Hankel matrices. | |
1360 ## The determinant and inverse are known explicitly. | |
1361 ## If a second output argument is present, d = DET(A) is returned: | |
1362 ## [A, d] = IPJFACT(N, K); | |
1363 ## | |
1364 ## Suggested by P. R. Graves-Morris. | |
1365 ## | |
1366 ## Reference: | |
1367 ## M.J.C. Gover, The explicit inverse of factorial Hankel matrices, | |
1368 ## Dept. of Mathematics, University of Bradford, 1993. | |
1369 | |
1370 if (nargin < 1 || nargin > 2) | |
1371 error ("gallery: 1 to 2 arguments are required for ipjfact matrix."); | |
1372 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1373 error ("gallery: N must be an integer for ipjfact matrix."); | |
1374 elseif (! isnumeric (k) || ! isscalar (k)) | |
1375 error ("gallery: K must be a numeric scalar for ipjfact matrix."); | |
1376 endif | |
1377 | |
1378 c = cumprod (2:n+1); | |
1379 d = cumprod (n+1:2*n) * c(n-1); | |
1380 | |
1381 A = hankel (c, d); | |
1382 | |
1383 switch (k) | |
1384 case (0), # do nothing | |
1385 case (1), A = ones (n) ./ A; | |
1386 otherwise | |
1387 error ("gallery: K must have a value of 0 or 1 for ipjfact matrix."); | |
1388 endswitch | |
1389 | |
1390 if (nargout == 2) | |
1391 d = 1; | |
1392 | |
1393 if (k == 0) | |
1394 for i = 1:n-1 | |
1395 d = d * prod (1:i+1) * prod (1:n-i); | |
1396 endfor | |
1397 d = d * prod (1:n+1); | |
1398 | |
1399 elseif (k == 1) | |
1400 for i = 0:n-1 | |
1401 d = d * prod (1:i) / prod (1:n+1+i); | |
1402 endfor | |
1403 if (rem (n*(n-1)/2, 2)) | |
1404 d = -d; | |
1405 endif | |
1406 | |
1407 else | |
1408 error ("gallery: K must have a value of 0 or 1 for ipjfact matrix."); | |
1409 endif | |
1410 | |
1411 detA = d; | |
1412 endif | |
1413 endfunction | |
1414 | |
1415 function J = jordbloc (n, lambda = 1) | |
1416 ## JORDBLOC Jordan block. | |
1417 ## JORDBLOC(N, LAMBDA) is the N-by-N Jordan block with eigenvalue | |
1418 ## LAMBDA. LAMBDA = 1 is the default. | |
1419 | |
1420 if (nargin < 1 || nargin > 2) | |
1421 error ("gallery: 1 to 2 arguments are required for jordbloc matrix."); | |
1422 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1423 error ("gallery: N must be an integer for jordbloc matrix."); | |
1424 elseif (! isnumeric (lambda) || ! isscalar (lambda)) | |
1425 error ("gallery: LAMBDA must be a numeric scalar for jordbloc matrix."); | |
1426 endif | |
1427 | |
1428 J = lambda * eye (n) + diag (ones (n-1, 1), 1); | |
1429 endfunction | |
1430 | |
1431 function U = kahan (n, theta = 1.2, pert = 25) | |
1432 ## KAHAN Kahan matrix - upper trapezoidal. | |
1433 ## KAHAN(N, THETA) is an upper trapezoidal matrix | |
1434 ## that has some interesting properties regarding estimation of | |
1435 ## condition and rank. | |
1436 ## The matrix is N-by-N unless N is a 2-vector, in which case it | |
1437 ## is N(1)-by-N(2). | |
1438 ## The parameter THETA defaults to 1.2. | |
1439 ## The useful range of THETA is 0 < THETA < PI. | |
1440 ## | |
1441 ## To ensure that the QR factorization with column pivoting does not | |
1442 ## interchange columns in the presence of rounding errors, the diagonal | |
1443 ## is perturbed by PERT*EPS*diag( [N:-1:1] ). | |
1444 ## The default is PERT = 25, which ensures no interchanges for KAHAN(N) | |
1445 ## up to at least N = 90 in IEEE arithmetic. | |
1446 ## KAHAN(N, THETA, PERT) uses the given value of PERT. | |
1447 ## | |
1448 ## The inverse of KAHAN(N, THETA) is known explicitly: see | |
1449 ## Higham (1987, p. 588), for example. | |
1450 ## The diagonal perturbation was suggested by Christian Bischof. | |
1451 ## | |
1452 ## References: | |
1453 ## W. Kahan, Numerical linear algebra, Canadian Math. Bulletin, | |
1454 ## 9 (1966), pp. 757-801. | |
1455 ## N.J. Higham, A survey of condition number estimation for | |
1456 ## triangular matrices, SIAM Review, 29 (1987), pp. 575-596. | |
1457 | |
1458 if (nargin < 1 || nargin > 3) | |
1459 error ("gallery: 1 to 3 arguments are required for kahan matrix."); | |
1460 elseif (! isnumeric (n) || all (numel (n) != [1 2]) || fix (n) != n) | |
1461 error ("gallery: N must be a 1 or 2 element integer for kahan matrix."); | |
1462 elseif (! isnumeric (theta) || ! isscalar (theta)) | |
1463 error ("gallery: THETA must be a numeric scalar for kahan matrix."); | |
1464 elseif (! isnumeric (pert) || ! isscalar (pert)) | |
1465 error ("gallery: PERT must be a numeric scalar for kahan matrix."); | |
1466 endif | |
1467 | |
1468 ## Parameter n specifies dimension: r-by-n | |
1469 r = n(1); | |
1470 n = n(end); | |
1471 | |
1472 s = sin (theta); | |
1473 c = cos (theta); | |
1474 | |
1475 U = eye (n) - c * triu (ones (n), 1); | |
1476 U = diag (s.^[0:n-1]) * U + pert*eps* diag ([n:-1:1]); | |
1477 if (r > n) | |
1478 U(r,n) = 0; # Extend to an r-by-n matrix | |
1479 else | |
1480 U = U(1:r,:); # Reduce to an r-by-n matrix | |
1481 endif | |
1482 endfunction | |
1483 | |
1484 function A = kms (n, rho = 0.5) | |
1485 ## KMS Kac-Murdock-Szego Toeplitz matrix. | |
1486 ## A = KMS(N, RHO) is the N-by-N Kac-Murdock-Szego Toeplitz matrix with | |
1487 ## A(i,j) = RHO^(ABS((i-j))) (for real RHO). | |
1488 ## If RHO is complex, then the same formula holds except that elements | |
1489 ## below the diagonal are conjugated. | |
1490 ## RHO defaults to 0.5. | |
1491 ## Properties: | |
1492 ## A has an LDL' factorization with | |
1493 ## L = INV(TRIW(N,-RHO,1)'), | |
1494 ## D(i,i) = (1-ABS(RHO)^2)*EYE(N) except D(1,1) = 1. | |
1495 ## A is positive definite if and only if 0 < ABS(RHO) < 1. | |
1496 ## INV(A) is tridiagonal. | |
1497 ## | |
1498 ## Reference: | |
1499 ## W.F. Trench, Numerical solution of the eigenvalue problem | |
1500 ## for Hermitian Toeplitz matrices, SIAM J. Matrix Analysis and Appl., | |
1501 ## 10 (1989), pp. 135-146 (and see the references therein). | |
1502 | |
1503 if (nargin < 1 || nargin > 2) | |
1504 error ("gallery: 1 to 2 arguments are required for lauchli matrix."); | |
1505 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1506 error("gallery: N must be an integer for lauchli matrix.") | |
1507 elseif (! isscalar (mu)) | |
1508 error("gallery: MU must be a scalar for lauchli matrix.") | |
1509 endif | |
1510 | |
1511 A = (1:n)'*ones(1,n); | |
1512 A = abs(A - A'); | |
1513 A = rho .^ A; | |
1514 if imag(rho) | |
1515 A = conj(tril(A,-1)) + triu(A); | |
1516 endif | |
1517 endfunction | |
1518 | |
1519 function B = krylov (A, x, j) | |
1520 ## KRYLOV Krylov matrix. | |
1521 ## KRYLOV(A, x, j) is the Krylov matrix | |
1522 ## [x, Ax, A^2x, ..., A^(j-1)x], | |
1523 ## where A is an n-by-n matrix and x is an n-vector. | |
1524 ## Defaults: x = ONES(n,1), j = n. | |
1525 ## KRYLOV(n) is the same as KRYLOV(RANDN(n)). | |
1526 ## | |
1527 ## Reference: | |
1528 ## G.H. Golub and C.F. Van Loan, Matrix Computations, second edition, | |
1529 ## Johns Hopkins University Press, Baltimore, Maryland, 1989, p. 369. | |
1530 | |
1531 if (nargin < 1 || nargin > 3) | |
1532 error ("gallery: 1 to 3 arguments are required for krylov matrix."); | |
1533 elseif (! isnumeric (A) || ! issquare (A) || ndims (A) != 2) | |
1534 error ("gallery: A must be a square 2D matrix for krylov matrix."); | |
1535 endif | |
1536 | |
1537 n = length (A); | |
1538 if (isscalar (A)) | |
1539 n = A; | |
1540 A = randn (n); | |
1541 endif | |
1542 | |
1543 if (nargin < 2) | |
1544 x = ones (n, 1); | |
1545 elseif (! isvector (x) || numel (x) != n) | |
1546 error ("gallery: X must be a vector of length equal to A for krylov matrix."); | |
1547 endif | |
1548 | |
1549 if (nargin < 3) | |
1550 j = n; | |
1551 elseif (! isnumeric (j) || ! isscalar (j) || fix (j) != j) | |
1552 error ("gallery: J must be an integer for krylov matrix."); | |
1553 endif | |
1554 | |
1555 B = ones (n, j); | |
1556 B(:,1) = x(:); | |
1557 for i = 2:j | |
1558 B(:,i) = A*B(:,i-1); | |
1559 endfor | |
1560 endfunction | |
1561 | |
1562 function A = lauchli (n, mu = sqrt (eps)) | |
1563 ## LAUCHLI Lauchli matrix - rectangular. | |
1564 ## LAUCHLI(N, MU) is the (N+1)-by-N matrix [ONES(1,N); MU*EYE(N))]. | |
1565 ## It is a well-known example in least squares and other problems | |
1566 ## that indicates the dangers of forming A'*A. | |
1567 ## MU defaults to SQRT(EPS). | |
1568 ## | |
1569 ## Reference: | |
1570 ## P. Lauchli, Jordan-Elimination und Ausgleichung nach | |
1571 ## kleinsten Quadraten, Numer. Math, 3 (1961), pp. 226-240. | |
1572 | |
1573 if (nargin < 1 || nargin > 2) | |
1574 error ("gallery: 1 to 2 arguments are required for lauchli matrix."); | |
1575 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1576 error ("gallery: N must be an integer for lauchli matrix."); | |
1577 elseif (! isscalar (mu)) | |
1578 error ("gallery: MU must be a scalar for lauchli matrix."); | |
1579 endif | |
1580 | |
1581 A = [ones(1, n) | |
1582 mu*eye(n) ]; | |
1583 endfunction | |
1584 | |
1585 function A = lehmer (n) | |
1586 ## LEHMER Lehmer matrix - symmetric positive definite. | |
1587 ## A = LEHMER(N) is the symmetric positive definite N-by-N matrix with | |
1588 ## A(i,j) = i/j for j >= i. | |
1589 ## A is totally nonnegative. INV(A) is tridiagonal, and explicit | |
1590 ## formulas are known for its entries. | |
1591 ## N <= COND(A) <= 4*N*N. | |
1592 ## | |
1593 ## References: | |
1594 ## M. Newman and J. Todd, The evaluation of matrix inversion | |
1595 ## programs, J. Soc. Indust. Appl. Math., 6 (1958), pp. 466-476. | |
1596 ## Solutions to problem E710 (proposed by D.H. Lehmer): The inverse | |
1597 ## of a matrix, Amer. Math. Monthly, 53 (1946), pp. 534-535. | |
1598 ## J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra, | |
1599 ## Birkhauser, Basel, and Academic Press, New York, 1977, p. 154. | |
1600 | |
1601 if (nargin != 1) | |
1602 error ("gallery: 1 argument is required for lehmer matrix."); | |
1603 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1604 error ("gallery: N must be an integer for lehmer matrix."); | |
1605 endif | |
1606 | |
1607 A = ones (n, 1) * (1:n); | |
1608 A = A./A'; | |
1609 A = tril (A) + tril (A, -1)'; | |
1610 endfunction | |
1611 | |
1612 function T = lesp (n) | |
1613 ## LESP A tridiagonal matrix with real, sensitive eigenvalues. | |
1614 ## LESP(N) is an N-by-N matrix whose eigenvalues are real and smoothly | |
1615 ## distributed in the interval approximately [-2*N-3.5, -4.5]. | |
1616 ## The sensitivities of the eigenvalues increase exponentially as | |
1617 ## the eigenvalues grow more negative. | |
1618 ## The matrix is similar to the symmetric tridiagonal matrix with | |
1619 ## the same diagonal entries and with off-diagonal entries 1, | |
1620 ## via a similarity transformation with D = diag(1!,2!,...,N!). | |
1621 ## | |
1622 ## References: | |
1623 ## H.W.J. Lenferink and M.N. Spijker, On the use of stability regions in | |
1624 ## the numerical analysis of initial value problems, | |
1625 ## Math. Comp., 57 (1991), pp. 221-237. | |
1626 ## L.N. Trefethen, Pseudospectra of matrices, in Numerical Analysis 1991, | |
1627 ## Proceedings of the 14th Dundee Conference, | |
1628 ## D.F. Griffiths and G.A. Watson, eds, Pitman Research Notes in | |
1629 ## Mathematics, volume 260, Longman Scientific and Technical, Essex, | |
1630 ## UK, 1992, pp. 234-266. | |
1631 | |
1632 if (nargin != 1) | |
1633 error ("gallery: 1 argument is required for lesp matrix."); | |
1634 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1635 error ("gallery: N must be an integer for lesp matrix."); | |
1636 endif | |
1637 | |
1638 x = 2:n; | |
1639 T = full (tridiag (ones (size (x)) ./x, -(2*[x n+1]+1), x)); | |
1640 endfunction | |
1641 | |
1642 function A = lotkin (n) | |
1643 ## LOTKIN Lotkin matrix. | |
1644 ## A = LOTKIN(N) is the Hilbert matrix with its first row altered to | |
1645 ## all ones. A is unsymmetric, ill-conditioned, and has many negative | |
1646 ## eigenvalues of small magnitude. | |
1647 ## The inverse has integer entries and is known explicitly. | |
1648 ## | |
1649 ## Reference: | |
1650 ## M. Lotkin, A set of test matrices, MTAC, 9 (1955), pp. 153-161. | |
1651 | |
1652 if (nargin != 1) | |
1653 error ("gallery: 1 argument is required for lotkin matrix."); | |
1654 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1655 error ("gallery: N must be an integer for lotkin matrix."); | |
1656 endif | |
1657 | |
1658 A = hilb (n); | |
1659 A(1,:) = ones (1, n); | |
1660 endfunction | |
1661 | |
1662 function A = minij (n) | |
1663 ## MINIJ Symmetric positive definite matrix MIN(i,j). | |
1664 ## A = MINIJ(N) is the N-by-N symmetric positive definite matrix with | |
1665 ## A(i,j) = MIN(i,j). | |
1666 ## Properties, variations: | |
1667 ## INV(A) is tridiagonal: it is minus the second difference matrix | |
1668 ## except its (N,N) element is 1. | |
1669 ## 2*A-ONES(N) (Givens' matrix) has tridiagonal inverse and | |
1670 ## eigenvalues .5*sec^2([2r-1)PI/4N], r=1:N. | |
1671 ## (N+1)*ONES(N)-A also has a tridiagonal inverse. | |
1672 ## | |
1673 ## References: | |
1674 ## J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra, | |
1675 ## Birkhauser, Basel, and Academic Press, New York, 1977, p. 158. | |
1676 ## D.E. Rutherford, Some continuant determinants arising in physics and | |
1677 ## chemistry---II, Proc. Royal Soc. Edin., 63, A (1952), pp. 232-241. | |
1678 ## (For the eigenvalues of Givens' matrix.) | |
1679 | |
1680 if (nargin != 1) | |
1681 error ("gallery: 1 argument is required for minij matrix."); | |
1682 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1683 error ("gallery: N must be an integer for minij matrix."); | |
1684 endif | |
1685 | |
1686 A = min (ones (n, 1) * (1:n), (1:n)' * ones (1, n)); | |
1687 endfunction | |
1688 | |
1689 function A = moler (n, alpha = -1) | |
1690 ## MOLER Moler matrix - symmetric positive definite. | |
1691 ## A = MOLER(N, ALPHA) is the symmetric positive definite N-by-N matrix | |
1692 ## U'*U where U = TRIW(N, ALPHA). | |
1693 ## For ALPHA = -1 (the default) A(i,j) = MIN(i,j)-2, A(i,i) = i. | |
1694 ## A has one small eigenvalue. | |
1695 ## | |
1696 ## Nash (1990) attributes the ALPHA = -1 matrix to Moler. | |
1697 ## | |
1698 ## Reference: | |
1699 ## J.C. Nash, Compact Numerical Methods for Computers: Linear | |
1700 ## Algebra and Function Minimisation, second edition, Adam Hilger, | |
1701 ## Bristol, 1990 (Appendix 1). | |
1702 | |
1703 if (nargin < 1 || nargin > 2) | |
1704 error ("gallery: 1 to 2 arguments are required for moler matrix."); | |
1705 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1706 error ("gallery: N must be an integer for moler matrix."); | |
1707 elseif (! isscalar (alpha)) | |
1708 error ("gallery: ALPHA must be a scalar for moler matrix."); | |
1709 endif | |
1710 | |
1711 A = triw (n, alpha)' * triw (n, alpha); | |
1712 endfunction | |
1713 | |
1714 function [A, T] = neumann (n) | |
1715 ## NEUMANN Singular matrix from the discrete Neumann problem (sparse). | |
1716 ## NEUMANN(N) is the singular, row diagonally dominant matrix resulting | |
1717 ## from discretizing the Neumann problem with the usual five point | |
1718 ## operator on a regular mesh. | |
1719 ## It has a one-dimensional null space with null vector ONES(N,1). | |
1720 ## The dimension N should be a perfect square, or else a 2-vector, | |
1721 ## in which case the dimension of the matrix is N(1)*N(2). | |
1722 ## | |
1723 ## Reference: | |
1724 ## R.J. Plemmons, Regular splittings and the discrete Neumann | |
1725 ## problem, Numer. Math., 25 (1976), pp. 153-161. | |
1726 | |
1727 if (nargin != 1) | |
1728 error ("gallery: 1 argument is required for neumann matrix."); | |
1729 elseif (! isnumeric (n) || all (numel (n) != [1 2]) || fix (n) != n) | |
1730 error ("gallery: N must be a 1 or 2 element integer for neumann matrix."); | |
1731 endif | |
1732 | |
1733 if (isscalar (n)) | |
1734 m = sqrt (n); | |
1735 if (m^2 != n) | |
1736 error ("gallery: N must be a perfect square for neumann matrix."); | |
1737 endif | |
1738 n(1) = m; | |
1739 n(2) = m; | |
1740 endif | |
1741 | |
1742 T = tridiag (n(1), -1, 2, -1); | |
1743 T(1,2) = -2; | |
1744 T(n(1),n(1)-1) = -2; | |
1745 | |
1746 A = kron (T, eye (n(2))) + kron (eye (n(2)), T); | |
1747 endfunction | |
1748 | |
1749 function Q = orthog (n, k = 1) | |
1750 ## ORTHOG Orthogonal and nearly orthogonal matrices. | |
1751 ## Q = ORTHOG(N, K) selects the K'th type of matrix of order N. | |
1752 ## K > 0 for exactly orthogonal matrices, K < 0 for diagonal scalings of | |
1753 ## orthogonal matrices. | |
1754 ## Available types: (K = 1 is the default) | |
1755 ## K = 1: Q(i,j) = SQRT(2/(n+1)) * SIN( i*j*PI/(n+1) ) | |
1756 ## Symmetric eigenvector matrix for second difference matrix. | |
1757 ## K = 2: Q(i,j) = 2/SQRT(2*n+1)) * SIN( 2*i*j*PI/(2*n+1) ) | |
1758 ## Symmetric. | |
1759 ## K = 3: Q(r,s) = EXP(2*PI*i*(r-1)*(s-1)/n) / SQRT(n) (i=SQRT(-1)) | |
1760 ## Unitary, the Fourier matrix. Q^4 is the identity. | |
1761 ## This is essentially the same matrix as FFT(EYE(N))/SQRT(N)! | |
1762 ## K = 4: Helmert matrix: a permutation of a lower Hessenberg matrix, | |
1763 ## whose first row is ONES(1:N)/SQRT(N). | |
1764 ## K = 5: Q(i,j) = SIN( 2*PI*(i-1)*(j-1)/n ) + COS( 2*PI*(i-1)*(j-1)/n ). | |
1765 ## Symmetric matrix arising in the Hartley transform. | |
1766 ## K = -1: Q(i,j) = COS( (i-1)*(j-1)*PI/(n-1) ) | |
1767 ## Chebyshev Vandermonde-like matrix, based on extrema of T(n-1). | |
1768 ## K = -2: Q(i,j) = COS( (i-1)*(j-1/2)*PI/n) ) | |
1769 ## Chebyshev Vandermonde-like matrix, based on zeros of T(n). | |
1770 ## | |
1771 ## References: | |
1772 ## N.J. Higham and D.J. Higham, Large growth factors in Gaussian | |
1773 ## elimination with pivoting, SIAM J. Matrix Analysis and Appl., | |
1774 ## 10 (1989), pp. 155-164. | |
1775 ## P. Morton, On the eigenvectors of Schur's matrix, J. Number Theory, | |
1776 ## 12 (1980), pp. 122-127. (Re. ORTHOG(N, 3)) | |
1777 ## H.O. Lancaster, The Helmert Matrices, Amer. Math. Monthly, 72 (1965), | |
1778 ## pp. 4-12. | |
1779 ## D. Bini and P. Favati, On a matrix algebra related to the discrete | |
1780 ## Hartley transform, SIAM J. Matrix Anal. Appl., 14 (1993), | |
1781 ## pp. 500-507. | |
1782 | |
1783 if (nargin < 1 || nargin > 2) | |
1784 error ("gallery: 1 to 2 arguments are required for orthog matrix."); | |
1785 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1786 error ("gallery: N must be an integer for orthog matrix."); | |
1787 elseif (! isnumeric (k) || ! isscalar (k)) | |
1788 error ("gallery: K must be a numeric scalar for orthog matrix."); | |
1789 endif | |
1790 | |
1791 switch (k) | |
1792 case (1) | |
1793 ## E'vectors second difference matrix | |
1794 m = (1:n)'*(1:n) * (pi/(n+1)); | |
1795 Q = sin (m) * sqrt (2/(n+1)); | |
1796 | |
1797 case (2) | |
1798 m = (1:n)'*(1:n) * (2*pi/(2*n+1)); | |
1799 Q = sin (m) * (2/ sqrt (2*n+1)); | |
1800 | |
1801 case (3) | |
1802 ## Vandermonde based on roots of unity | |
1803 m = 0:n-1; | |
1804 Q = exp (m'*m*2*pi* sqrt (-1) / n) / sqrt (n); | |
1805 | |
1806 case (4) | |
1807 ## Helmert matrix | |
1808 Q = tril (ones (n)); | |
1809 Q(1,2:n) = ones (1, n-1); | |
1810 for i = 2:n | |
1811 Q(i,i) = -(i-1); | |
1812 end | |
1813 Q = diag (sqrt ([n 1:n-1] .* [1:n])) \ Q; | |
1814 | |
1815 case (5) | |
1816 ## Hartley matrix | |
1817 m = (0:n-1)'*(0:n-1) * (2*pi/n); | |
1818 Q = (cos (m) + sin (m)) / sqrt (n); | |
1819 | |
1820 case (-1) | |
1821 ## extrema of T(n-1) | |
1822 m = (0:n-1)'*(0:n-1) * (pi/(n-1)); | |
1823 Q = cos (m); | |
1824 | |
1825 case (-2) | |
1826 ## zeros of T(n) | |
1827 m = (0:n-1)'*(.5:n-.5) * (pi/n); | |
1828 Q = cos (m); | |
1829 | |
1830 otherwise | |
16766
7268845c0a1e
avoid backquote in error messages, some uses in doc strings
John W. Eaton <jwe@octave.org>
parents:
16734
diff
changeset
|
1831 error ("gallery: unknown K '%d' for orthog matrix.", k); |
16634 | 1832 endswitch |
1833 endfunction | |
1834 | |
1835 function A = parter (n) | |
1836 ## PARTER Parter matrix - a Toeplitz matrix with singular values near PI. | |
1837 ## PARTER(N) is the matrix with (i,j) element 1/(i-j+0.5). | |
1838 ## It is a Cauchy matrix and a Toeplitz matrix. | |
1839 ## | |
1840 ## At the Second SIAM Conference on Linear Algebra, Raleigh, N.C., | |
1841 ## 1985, Cleve Moler noted that most of the singular values of | |
1842 ## PARTER(N) are very close to PI. An explanation of the phenomenon | |
1843 ## was given by Parter; see also the paper by Tyrtyshnikov. | |
1844 ## | |
1845 ## References: | |
1846 ## The MathWorks Newsletter, Volume 1, Issue 1, March 1986, page 2. | |
1847 ## S.V. Parter, On the distribution of the singular values of Toeplitz | |
1848 ## matrices, Linear Algebra and Appl., 80 (1986), pp. 115-130. | |
1849 ## E.E. Tyrtyshnikov, Cauchy-Toeplitz matrices and some applications, | |
1850 ## Linear Algebra and Appl., 149 (1991), pp. 1-18. | |
1851 | |
1852 if (nargin != 1) | |
1853 error ("gallery: 1 argument is required for parter matrix."); | |
1854 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1855 error ("gallery: N must be an integer for parter matrix."); | |
1856 endif | |
1857 | |
1858 A = cauchy ((1:n) + 0.5, -(1:n)); | |
1859 endfunction | |
1860 | |
1861 function P = pei (n, alpha = 1) | |
1862 ## PEI Pei matrix. | |
1863 ## PEI(N, ALPHA), where ALPHA is a scalar, is the symmetric matrix | |
1864 ## ALPHA*EYE(N) + ONES(N). | |
1865 ## If ALPHA is omitted then ALPHA = 1 is used. | |
1866 ## The matrix is singular for ALPHA = 0, -N. | |
1867 ## | |
1868 ## Reference: | |
1869 ## M.L. Pei, A test matrix for inversion procedures, | |
1870 ## Comm. ACM, 5 (1962), p. 508. | |
1871 | |
1872 if (nargin < 1 || nargin > 2) | |
1873 error ("gallery: 1 to 2 arguments are required for pei matrix."); | |
1874 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1875 error ("gallery: N must be an integer for pei matrix."); | |
1876 elseif (! isnumeric (w) || ! isscalar (w)) | |
1877 error ("gallery: ALPHA must be a scalar for pei matrix."); | |
1878 endif | |
1879 | |
1880 P = alpha * eye (n) + ones (n); | |
1881 endfunction | |
1882 | |
1883 function A = poisson (n) | |
1884 ## POISSON Block tridiagonal matrix from Poisson's equation (sparse). | |
1885 ## POISSON(N) is the block tridiagonal matrix of order N^2 | |
1886 ## resulting from discretizing Poisson's equation with the | |
1887 ## 5-point operator on an N-by-N mesh. | |
1888 ## | |
1889 ## Reference: | |
1890 ## G.H. Golub and C.F. Van Loan, Matrix Computations, second edition, | |
1891 ## Johns Hopkins University Press, Baltimore, Maryland, 1989 | |
1892 ## (Section 4.5.4). | |
1893 | |
1894 if (nargin != 1) | |
1895 error ("gallery: 1 argument is required for poisson matrix."); | |
1896 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1897 error ("gallery: N must be an integer for poisson matrix."); | |
1898 endif | |
1899 | |
1900 S = tridiag (n, -1, 2, -1); | |
1901 I = speye (n); | |
1902 A = kron (I, S) + kron (S, I); | |
1903 endfunction | |
1904 | |
1905 function A = prolate (n, w = 0.25) | |
1906 ## PROLATE Prolate matrix - symmetric, ill-conditioned Toeplitz matrix. | |
1907 ## A = PROLATE(N, W) is the N-by-N prolate matrix with parameter W. | |
1908 ## It is a symmetric Toeplitz matrix. | |
1909 ## If 0 < W < 0.5 then | |
1910 ## - A is positive definite | |
1911 ## - the eigenvalues of A are distinct, lie in (0, 1), and | |
1912 ## tend to cluster around 0 and 1. | |
1913 ## W defaults to 0.25. | |
1914 ## | |
1915 ## Reference: | |
1916 ## J.M. Varah. The Prolate matrix. Linear Algebra and Appl., | |
1917 ## 187:269--278, 1993. | |
1918 | |
1919 if (nargin < 1 || nargin > 2) | |
1920 error ("gallery: 1 to 2 arguments are required for prolate matrix."); | |
1921 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
1922 error ("gallery: N must be an integer for prolate matrix."); | |
1923 elseif (! isnumeric (w) || ! isscalar (w)) | |
1924 error ("gallery: W must be a scalar for prolate matrix."); | |
1925 endif | |
1926 | |
1927 a = zeros (n, 1); | |
1928 a(1) = 2*w; | |
1929 a(2:n) = sin (2*pi*w*(1:n-1)) ./ (pi*(1:n-1)); | |
1930 | |
1931 A = toeplitz(a); | |
1932 endfunction | |
1933 | |
1934 function H = randhess (x) | |
1935 ## NOTE: this function was named ohess in the original Test Matrix Toolbox | |
1936 ## RANDHESS Random, orthogonal upper Hessenberg matrix. | |
1937 ## H = RANDHESS(N) is an N-by-N real, random, orthogonal | |
1938 ## upper Hessenberg matrix. | |
1939 ## Alternatively, H = RANDHESS(X), where X is an arbitrary real | |
1940 ## N-vector (N > 1) constructs H non-randomly using the elements | |
1941 ## of X as parameters. | |
1942 ## In both cases H is constructed via a product of N-1 Givens rotations. | |
1943 ## | |
1944 ## Note: See Gragg (1986) for how to represent an N-by-N (complex) | |
1945 ## unitary Hessenberg matrix with positive subdiagonal elements in terms | |
1946 ## of 2N-1 real parameters (the Schur parametrization). | |
1947 ## This M-file handles the real case only and is intended simply as a | |
1948 ## convenient way to generate random or non-random orthogonal Hessenberg | |
1949 ## matrices. | |
1950 ## | |
1951 ## Reference: | |
1952 ## W.B. Gragg, The QR algorithm for unitary Hessenberg matrices, | |
1953 ## J. Comp. Appl. Math., 16 (1986), pp. 1-8. | |
1954 | |
1955 if (nargin != 1) | |
1956 error ("gallery: 1 argument is required for randhess matrix."); | |
1957 elseif (! isnumeric (x) || ! isreal (x)) | |
1958 error ("gallery: N or X must be numeric real values for randhess matrix."); | |
1959 endif | |
1960 | |
1961 if (isscalar (x)) | |
1962 n = x; | |
1963 x = rand (n-1, 1) * 2*pi; | |
1964 H = eye (n); | |
1965 H(n,n) = sign (randn); | |
1966 elseif (isvector (x)) | |
1967 n = numel (x); | |
1968 H = eye (n); | |
1969 H(n,n) = sign (x(n)) + (x(n) == 0); # Second term ensures H(n,n) nonzero. | |
1970 else | |
1971 error ("gallery: N or X must be a scalar or a vector for randhess matrix."); | |
1972 endif | |
1973 | |
1974 for i = n:-1:2 | |
1975 ## Apply Givens rotation through angle x(i-1). | |
1976 theta = x(i-1); | |
1977 c = cos (theta); | |
1978 s = sin (theta); | |
1979 H([i-1 i], :) = [ c*H(i-1,:)+s*H(i,:) | |
1980 -s*H(i-1,:)+c*H(i,:) ]; | |
1981 endfor | |
1982 endfunction | |
1983 | |
1984 function A = rando (n, k = 1) | |
1985 ## RANDO Random matrix with elements -1, 0 or 1. | |
1986 ## A = RANDO(N, K) is a random N-by-N matrix with elements from | |
1987 ## one of the following discrete distributions (default K = 1): | |
1988 ## K = 1: A(i,j) = 0 or 1 with equal probability, | |
1989 ## K = 2: A(i,j) = -1 or 1 with equal probability, | |
1990 ## K = 3: A(i,j) = -1, 0 or 1 with equal probability. | |
1991 ## N may be a 2-vector, in which case the matrix is N(1)-by-N(2). | |
1992 | |
1993 if (nargin < 1 || nargin > 2) | |
1994 error ("gallery: 1 to 2 arguments are required for rando matrix."); | |
1995 elseif (! isnumeric (n) || all (numel (n) != [1 2]) || fix (n) != n) | |
1996 error ("gallery: N must be an integer for rando matrix."); | |
1997 elseif (! isnumeric (k) || ! isscalar (k)) | |
1998 error ("gallery: K must be a numeric scalar for smoke matrix."); | |
1999 endif | |
2000 | |
2001 ## Parameter n specifies dimension: m-by-n. | |
2002 m = n(1); | |
2003 n = n(end); | |
2004 | |
2005 switch (k) | |
2006 case (1), A = floor ( rand(m, n) + 0.5); # {0, 1} | |
2007 case (2), A = 2*floor ( rand(m, n) + 0.5) -1; # {-1, 1} | |
2008 case (3), A = round (3*rand(m, n) - 1.5); # {-1, 0, 1} | |
2009 otherwise | |
16766
7268845c0a1e
avoid backquote in error messages, some uses in doc strings
John W. Eaton <jwe@octave.org>
parents:
16734
diff
changeset
|
2010 error ("gallery: unknown K '%d' for smoke matrix.", k); |
16634 | 2011 endswitch |
2012 | |
2013 endfunction | |
2014 | |
2015 function A = randsvd (n, kappa = sqrt (1/eps), mode = 3, kl = n-1, ku = kl) | |
2016 ## RANDSVD Random matrix with pre-assigned singular values. | |
2017 ## RANDSVD(N, KAPPA, MODE, KL, KU) is a (banded) random matrix of order N | |
2018 ## with COND(A) = KAPPA and singular values from the distribution MODE. | |
2019 ## N may be a 2-vector, in which case the matrix is N(1)-by-N(2). | |
2020 ## Available types: | |
2021 ## MODE = 1: one large singular value, | |
2022 ## MODE = 2: one small singular value, | |
2023 ## MODE = 3: geometrically distributed singular values, | |
2024 ## MODE = 4: arithmetically distributed singular values, | |
2025 ## MODE = 5: random singular values with unif. dist. logarithm. | |
2026 ## If omitted, MODE defaults to 3, and KAPPA defaults to SQRT(1/EPS). | |
2027 ## If MODE < 0 then the effect is as for ABS(MODE) except that in the | |
2028 ## original matrix of singular values the order of the diagonal entries | |
2029 ## is reversed: small to large instead of large to small. | |
2030 ## KL and KU are the lower and upper bandwidths respectively; if they | |
2031 ## are omitted a full matrix is produced. | |
2032 ## If only KL is present, KU defaults to KL. | |
2033 ## Special case: if KAPPA < 0 then a random full symmetric positive | |
2034 ## definite matrix is produced with COND(A) = -KAPPA and | |
2035 ## eigenvalues distributed according to MODE. | |
2036 ## KL and KU, if present, are ignored. | |
2037 ## | |
2038 ## Reference: | |
2039 ## N.J. Higham, Accuracy and Stability of Numerical Algorithms, | |
2040 ## Society for Industrial and Applied Mathematics, Philadelphia, PA, | |
2041 ## USA, 1996; sec. 26.3. | |
2042 ## | |
2043 ## This routine is similar to the more comprehensive Fortran routine xLATMS | |
2044 ## in the following reference: | |
2045 ## J.W. Demmel and A. McKenney, A test matrix generation suite, | |
2046 ## LAPACK Working Note #9, Courant Institute of Mathematical Sciences, | |
2047 ## New York, 1989. | |
2048 | |
2049 if (nargin < 1 || nargin > 5) | |
2050 error ("gallery: 1 to 5 arguments are required for randsvd matrix."); | |
2051 elseif (! isnumeric (n) || all (numel (n) != [1 2]) || fix (n) != n) | |
2052 error ("gallery: N must be a 1 or 2 element integer vector for randsvd matrix."); | |
2053 elseif (! isnumeric (kappa) || ! isscalar (kappa)) | |
2054 error ("gallery: KAPPA must be a numeric scalar for randsvd matrix."); | |
2055 elseif (abs (kappa) < 1) | |
2056 error ("gallery: KAPPA must larger than or equal to 1 for randsvd matrix."); | |
2057 elseif (! isnumeric (mode) || ! isscalar (mode)) | |
2058 error ("gallery: MODE must be a numeric scalar for randsvd matrix."); | |
2059 elseif (! isnumeric (kl) || ! isscalar (kl)) | |
2060 error ("gallery: KL must be a numeric scalar for randsvd matrix."); | |
2061 elseif (! isnumeric (ku) || ! isscalar (ku)) | |
2062 error ("gallery: KU must be a numeric scalar for randsvd matrix."); | |
2063 endif | |
2064 | |
2065 posdef = 0; | |
2066 if (kappa < 0) | |
2067 posdef = 1; | |
2068 kappa = -kappa; | |
2069 endif | |
2070 | |
2071 ## Parameter n specifies dimension: m-by-n. | |
2072 m = n(1); | |
2073 n = n(end); | |
2074 p = min ([m n]); | |
2075 | |
2076 ## If A will be a vector | |
2077 if (p == 1) | |
2078 A = randn (m, n); | |
2079 A = A / norm (A); | |
2080 return | |
2081 end | |
2082 | |
2083 ## Set up vector sigma of singular values. | |
2084 switch (abs (mode)) | |
2085 case (1) | |
2086 sigma = ones (p, 1) ./ kappa; | |
2087 sigma(1) = 1; | |
2088 case (2) | |
2089 sigma = ones (p, 1); | |
2090 sigma(p) = 1 / kappa; | |
2091 case (3) | |
2092 factor = kappa^(-1/(p-1)); | |
2093 sigma = factor.^[0:p-1]; | |
2094 case (4) | |
2095 sigma = ones (p, 1) - (0:p-1)'/(p-1)*(1-1/kappa); | |
2096 case (5) | |
2097 ## In this case cond (A) <= kappa. | |
2098 rand ("uniform"); | |
2099 sigma = exp (-rand (p, 1) * log (kappa)); | |
2100 otherwise | |
16766
7268845c0a1e
avoid backquote in error messages, some uses in doc strings
John W. Eaton <jwe@octave.org>
parents:
16734
diff
changeset
|
2101 error ("gallery: unknown MODE '%d' for randsvd matrix.", mode); |
16634 | 2102 endswitch |
2103 | |
2104 ## Convert to diagonal matrix of singular values. | |
2105 if (mode < 0) | |
2106 sigma = sigma (p:-1:1); | |
2107 end | |
2108 sigma = diag (sigma); | |
2109 | |
2110 if (posdef) | |
2111 ## handle case where KAPPA was negative | |
2112 Q = qmult (p); | |
2113 A = Q' * sigma * Q; | |
2114 A = (A + A') / 2; # Ensure matrix is symmetric. | |
2115 return | |
2116 endif | |
2117 | |
2118 if (m != n) | |
2119 ## Expand to m-by-n diagonal matrix | |
2120 sigma(m, n) = 0; | |
2121 end | |
2122 | |
2123 if (kl == 0 && ku == 0) | |
2124 ## Diagonal matrix requested - nothing more to do. | |
2125 A = sigma; | |
2126 else | |
2127 ## A = U*sigma*V, where U, V are random orthogonal matrices from the | |
2128 ## Haar distribution. | |
2129 A = qmult (sigma'); | |
2130 A = qmult (A'); | |
2131 | |
2132 if (kl < n-1 || ku < n-1) | |
2133 ## Bandwidth reduction | |
2134 A = bandred (A, kl, ku); | |
2135 endif | |
2136 endif | |
2137 endfunction | |
2138 | |
2139 function A = redheff (n) | |
2140 ## REDHEFF A (0,1) matrix of Redheffer associated with the Riemann hypothesis. | |
2141 ## A = REDHEFF(N) is an N-by-N matrix of 0s and 1s defined by | |
2142 ## A(i,j) = 1 if j = 1 or if i divides j, | |
2143 ## A(i,j) = 0 otherwise. | |
2144 ## It has N - FLOOR(LOG2(N)) - 1 eigenvalues equal to 1, | |
2145 ## a real eigenvalue (the spectral radius) approximately SQRT(N), | |
2146 ## a negative eigenvalue approximately -SQRT(N), | |
2147 ## and the remaining eigenvalues are provably ``small''. | |
2148 ## Barrett and Jarvis (1992) conjecture that | |
2149 ## ``the small eigenvalues all lie inside the unit circle | |
2150 ## ABS(Z) = 1'', | |
2151 ## and a proof of this conjecture, together with a proof that some | |
2152 ## eigenvalue tends to zero as N tends to infinity, would yield | |
2153 ## a new proof of the prime number theorem. | |
2154 ## The Riemann hypothesis is true if and only if | |
2155 ## DET(A) = O( N^(1/2+epsilon) ) for every epsilon > 0 | |
2156 ## (`!' denotes factorial). | |
2157 ## See also RIEMANN. | |
2158 ## | |
2159 ## Reference: | |
2160 ## W.W. Barrett and T.J. Jarvis, | |
2161 ## Spectral Properties of a Matrix of Redheffer, | |
2162 ## Linear Algebra and Appl., 162 (1992), pp. 673-683. | |
2163 | |
2164 if (nargin != 1) | |
2165 error ("gallery: 1 argument is required for redheff matrix."); | |
2166 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
2167 error ("gallery: N must be an integer for redheff matrix."); | |
2168 endif | |
2169 | |
2170 i = (1:n)' * ones (1, n); | |
2171 A = ! rem (i', i); | |
2172 A(:,1) = ones (n, 1); | |
2173 endfunction | |
2174 | |
2175 function A = riemann (n) | |
2176 ## RIEMANN A matrix associated with the Riemann hypothesis. | |
2177 ## A = RIEMANN(N) is an N-by-N matrix for which the | |
2178 ## Riemann hypothesis is true if and only if | |
2179 ## DET(A) = O( N! N^(-1/2+epsilon) ) for every epsilon > 0 | |
2180 ## (`!' denotes factorial). | |
2181 ## A = B(2:N+1, 2:N+1), where | |
2182 ## B(i,j) = i-1 if i divides j and -1 otherwise. | |
2183 ## Properties include, with M = N+1: | |
2184 ## Each eigenvalue E(i) satisfies ABS(E(i)) <= M - 1/M. | |
2185 ## i <= E(i) <= i+1 with at most M-SQRT(M) exceptions. | |
2186 ## All integers in the interval (M/3, M/2] are eigenvalues. | |
2187 ## | |
2188 ## See also REDHEFF. | |
2189 ## | |
2190 ## Reference: | |
2191 ## F. Roesler, Riemann's hypothesis as an eigenvalue problem, | |
2192 ## Linear Algebra and Appl., 81 (1986), pp. 153-198. | |
2193 | |
2194 if (nargin != 1) | |
2195 error ("gallery: 1 argument is required for riemann matrix."); | |
2196 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
2197 error ("gallery: N must be an integer for riemann matrix."); | |
2198 endif | |
2199 | |
2200 n = n+1; | |
2201 i = (2:n)' * ones (1, n-1); | |
2202 j = i'; | |
2203 A = i .* (! rem (j, i)) - ones (n-1); | |
2204 endfunction | |
2205 | |
2206 function A = ris (n) | |
2207 ## NOTE: this function was named dingdong in the original Test Matrix Toolbox | |
2208 ## RIS Dingdong matrix - a symmetric Hankel matrix. | |
2209 ## A = RIS(N) is the symmetric N-by-N Hankel matrix with | |
2210 ## A(i,j) = 0.5/(N-i-j+1.5). | |
2211 ## The eigenvalues of A cluster around PI/2 and -PI/2. | |
2212 ## | |
2213 ## Invented by F.N. Ris. | |
2214 ## | |
2215 ## Reference: | |
2216 ## J.C. Nash, Compact Numerical Methods for Computers: Linear | |
2217 ## Algebra and Function Minimisation, second edition, Adam Hilger, | |
2218 ## Bristol, 1990 (Appendix 1). | |
2219 | |
2220 if (nargin != 1) | |
2221 error ("gallery: 1 argument is required for ris matrix."); | |
2222 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
2223 error ("gallery: N must be an integer for ris matrix."); | |
2224 endif | |
2225 | |
2226 p = -2*(1:n) + (n+1.5); | |
2227 A = cauchy (p); | |
2228 endfunction | |
2229 | |
2230 function A = smoke (n, k = 0) | |
2231 ## SMOKE Smoke matrix - complex, with a `smoke ring' pseudospectrum. | |
2232 ## SMOKE(N) is an N-by-N matrix with 1s on the | |
2233 ## superdiagonal, 1 in the (N,1) position, and powers of | |
2234 ## roots of unity along the diagonal. | |
2235 ## SMOKE(N, 1) is the same except for a zero (N,1) element. | |
2236 ## The eigenvalues of SMOKE(N, 1) are the N'th roots of unity; | |
2237 ## those of SMOKE(N) are the N'th roots of unity times 2^(1/N). | |
2238 ## | |
2239 ## Try PS(SMOKE(32)). For SMOKE(N, 1) the pseudospectrum looks | |
2240 ## like a sausage folded back on itself. | |
2241 ## GERSH(SMOKE(N, 1)) is interesting. | |
2242 ## | |
2243 ## Reference: | |
2244 ## L. Reichel and L.N. Trefethen, Eigenvalues and pseudo-eigenvalues of | |
2245 ## Toeplitz matrices, Linear Algebra and Appl., 162-164:153-185, 1992. | |
2246 | |
2247 if (nargin < 1 || nargin > 2) | |
2248 error ("gallery: 1 to 2 arguments are required for smoke matrix."); | |
2249 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
2250 error ("gallery: N must be an integer for smoke matrix."); | |
2251 elseif (! isnumeric (n) || ! isscalar (n)) | |
2252 error ("gallery: K must be a numeric scalar for smoke matrix."); | |
2253 endif | |
2254 | |
2255 w = exp(2*pi*i/n); | |
2256 A = diag( [w.^(1:n-1) 1] ) + diag(ones(n-1,1),1); | |
2257 | |
2258 switch (k) | |
2259 case (0), A(n,1) = 1; | |
2260 case (1), # do nothing | |
2261 otherwise, | |
2262 error ("gallery: K must have a value of 0 or 1 for smoke matrix."); | |
2263 endswitch | |
2264 endfunction | |
2265 | |
2266 function T = toeppd (n, m = n, w = rand (m,1), theta = rand (m,1)) | |
2267 ## NOTE: this function was named pdtoep in the original Test Matrix Toolbox | |
2268 ## TOEPPD Symmetric positive definite Toeplitz matrix. | |
2269 ## TOEPPD(N, M, W, THETA) is an N-by-N symmetric positive (semi-) | |
2270 ## definite (SPD) Toeplitz matrix, comprised of the sum of M rank 2 | |
2271 ## (or, for certain THETA, rank 1) SPD Toeplitz matrices. | |
2272 ## Specifically, | |
2273 ## T = W(1)*T(THETA(1)) + ... + W(M)*T(THETA(M)), | |
2274 ## where T(THETA(k)) has (i,j) element COS(2*PI*THETA(k)*(i-j)). | |
2275 ## Defaults: M = N, W = RAND(M,1), THETA = RAND(M,1). | |
2276 ## | |
2277 ## Reference: | |
2278 ## G. Cybenko and C.F. Van Loan, Computing the minimum eigenvalue of | |
2279 ## a symmetric positive definite Toeplitz matrix, SIAM J. Sci. Stat. | |
2280 ## Comput., 7 (1986), pp. 123-131. | |
2281 | |
2282 if (nargin < 1 || nargin > 4) | |
2283 error ("gallery: 1 to 4 arguments are required for toeppd matrix."); | |
2284 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
2285 error ("gallery: N must be a numeric integer for toeppd matrix."); | |
2286 elseif (! isnumeric (m) || ! isscalar (m) || fix (m) != m) | |
2287 error ("gallery: M must be a numeric integer for toeppd matrix."); | |
2288 elseif (numel (w) != m || numel (theta) != m) | |
2289 error ("gallery: W and THETA must be vectors of length M for toeppd matrix."); | |
2290 endif | |
2291 | |
2292 T = zeros (n); | |
2293 E = 2*pi * ((1:n)' * ones (1, n) - ones (n, 1) * (1:n)); | |
2294 | |
2295 for i = 1:m | |
2296 T = T + w(i) * cos (theta(i)*E); | |
2297 endfor | |
2298 endfunction | |
2299 | |
2300 function P = toeppen (n, a = 1, b = -10, c = 0, d = 10, e = 1) | |
2301 ## NOTE: this function was named pentoep in the original Test Matrix Toolbox | |
2302 ## TOEPPEN Pentadiagonal Toeplitz matrix (sparse). | |
2303 ## P = TOEPPEN(N, A, B, C, D, E) is the N-by-N pentadiagonal | |
2304 ## Toeplitz matrix with diagonals composed of the numbers | |
2305 ## A =: P(3,1), B =: P(2,1), C =: P(1,1), D =: P(1,2), E =: P(1,3). | |
2306 ## Default: (A,B,C,D,E) = (1,-10,0,10,1) (a matrix of Rutishauser). | |
2307 ## This matrix has eigenvalues lying approximately on | |
2308 ## the line segment 2*cos(2*t) + 20*i*sin(t). | |
2309 ## | |
2310 ## Interesting plots are | |
2311 ## PS(FULL(TOEPPEN(32,0,1,0,0,1/4))) - `triangle' | |
2312 ## PS(FULL(TOEPPEN(32,0,1/2,0,0,1))) - `propeller' | |
2313 ## PS(FULL(TOEPPEN(32,0,1/2,1,1,1))) - `fish' | |
2314 ## | |
2315 ## References: | |
2316 ## R.M. Beam and R.F. Warming, The asymptotic spectra of | |
2317 ## banded Toeplitz and quasi-Toeplitz matrices, SIAM J. Sci. | |
2318 ## Comput. 14 (4), 1993, pp. 971-1006. | |
2319 ## H. Rutishauser, On test matrices, Programmation en Mathematiques | |
2320 ## Numeriques, Editions Centre Nat. Recherche Sci., Paris, 165, | |
2321 ## 1966, pp. 349-365. | |
2322 | |
2323 if (nargin < 1 || nargin > 6) | |
2324 error ("gallery: 1 to 6 arguments are required for toeppen matrix."); | |
2325 elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n) | |
2326 error ("gallery: N must be a numeric integer for toeppen matrix."); | |
2327 elseif (any (cellfun (@(x) ! isnumeric (x) || ! isscalar (x), {a b c d e}))) | |
2328 error ("gallery: A, B, C, D and E must be numeric scalars for toeppen matrix."); | |
2329 endif | |
2330 | |
2331 P = spdiags ([a*ones(n,1) b*ones(n,1) c*ones(n,1) d*ones(n,1) e*ones(n,1)], | |
2332 -2:2, n, n); | |
2333 endfunction | |
2334 | |
2335 function T = tridiag (n, x = -1, y = 2, z = -1) | |
2336 ## TRIDIAG Tridiagonal matrix (sparse). | |
2337 ## TRIDIAG(X, Y, Z) is the tridiagonal matrix with subdiagonal X, | |
2338 ## diagonal Y, and superdiagonal Z. | |
2339 ## X and Z must be vectors of dimension one less than Y. | |
2340 ## Alternatively TRIDIAG(N, C, D, E), where C, D, and E are all | |
2341 ## scalars, yields the Toeplitz tridiagonal matrix of order N | |
2342 ## with subdiagonal elements C, diagonal elements D, and superdiagonal | |
2343 ## elements E. This matrix has eigenvalues (Todd 1977) | |
2344 ## D + 2*SQRT(C*E)*COS(k*PI/(N+1)), k=1:N. | |
2345 ## TRIDIAG(N) is the same as TRIDIAG(N,-1,2,-1), which is | |
2346 ## a symmetric positive definite M-matrix (the negative of the | |
2347 ## second difference matrix). | |
2348 ## | |
2349 ## References: | |
2350 ## J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra, | |
2351 ## Birkhauser, Basel, and Academic Press, New York, 1977, p. 155. | |
2352 ## D.E. Rutherford, Some continuant determinants arising in physics and | |
2353 ## chemistry---II, Proc. Royal Soc. Edin., 63, A (1952), pp. 232-241. | |
2354 | |
2355 if (nargin != 1 && nargin != 3 && nargin != 4) | |
2356 error ("gallery: 1, 3, or 4 arguments are required for tridiag matrix."); | |
2357 elseif (nargin == 3) | |
2358 z = y; | |
2359 y = x; | |
2360 x = n; | |
2361 endif | |
2362 | |
2363 ## Force column vectors | |
2364 x = x(:); | |
2365 y = y(:); | |
2366 z = z(:); | |
2367 | |
2368 if (isscalar (x) && isscalar (y) && isscalar (z)) | |
2369 x *= ones (n-1, 1); | |
2370 z *= ones (n-1, 1); | |
2371 y *= ones (n, 1); | |
2372 elseif (numel (y) != numel (x) + 1) | |
2373 error ("gallery: X must have one element less than Y for tridiag matrix."); | |
2374 elseif (numel (y) != numel (z) + 1) | |
2375 error ("gallery: Z must have one element less than Y for tridiag matrix."); | |
2376 endif | |
2377 | |
2378 ## T = diag (x, -1) + diag (y) + diag (z, 1); # For non-sparse matrix. | |
2379 n = numel (y); | |
2380 T = spdiags ([[x;0] y [0;z]], -1:1, n, n); | |
2381 endfunction | |
2382 | |
2383 function t = triw (n, alpha = -1, k = -1) | |
2384 ## TRIW Upper triangular matrix discussed by Wilkinson and others. | |
2385 ## TRIW(N, ALPHA, K) is the upper triangular matrix with ones on | |
2386 ## the diagonal and ALPHAs on the first K >= 0 superdiagonals. | |
2387 ## N may be a 2-vector, in which case the matrix is N(1)-by-N(2) and | |
2388 ## upper trapezoidal. | |
2389 ## Defaults: ALPHA = -1, | |
2390 ## K = N - 1 (full upper triangle). | |
2391 ## TRIW(N) is a matrix discussed by Kahan, Golub and Wilkinson. | |
2392 ## | |
2393 ## Ostrowski (1954) shows that | |
2394 ## COND(TRIW(N,2)) = COT(PI/(4*N))^2, | |
2395 ## and for large ABS(ALPHA), | |
2396 ## COND(TRIW(N,ALPHA)) is approximately ABS(ALPHA)^N*SIN(PI/(4*N-2)). | |
2397 ## | |
2398 ## Adding -2^(2-N) to the (N,1) element makes TRIW(N) singular, | |
2399 ## as does adding -2^(1-N) to all elements in the first column. | |
2400 ## | |
2401 ## References: | |
2402 ## G.H. Golub and J.H. Wilkinson, Ill-conditioned eigensystems and the | |
2403 ## computation of the Jordan canonical form, SIAM Review, | |
2404 ## 18(4), 1976, pp. 578-619. | |
2405 ## W. Kahan, Numerical linear algebra, Canadian Math. Bulletin, | |
2406 ## 9 (1966), pp. 757-801. | |
2407 ## A.M. Ostrowski, On the spectrum of a one-parametric family of | |
2408 ## matrices, J. Reine Angew. Math., 193 (3/4), 1954, pp. 143-160. | |
2409 ## J.H. Wilkinson, Singular-value decomposition---basic aspects, | |
2410 ## in D.A.H. Jacobs, ed., Numerical Software---Needs and Availability, | |
2411 ## Academic Press, London, 1978, pp. 109-135. | |
2412 | |
2413 if (nargin < 1 || nargin > 3) | |
2414 error ("gallery: 1 to 3 arguments are required for triw matrix."); | |
2415 elseif (! isnumeric (n) || all (numel (n) != [1 2])) | |
2416 error ("gallery: N must be a 1 or 2 elements vector for triw matrix."); | |
2417 elseif (! isscalar (alpha)) | |
2418 error ("gallery: ALPHA must be a scalar for triw matrix."); | |
2419 elseif (! isscalar (k) || ! isnumeric (k) || fix (k) != k) | |
2420 error ("gallery: K must be a numeric integer for triw matrix."); | |
2421 endif | |
2422 | |
2423 m = n(1); # Parameter n specifies dimension: m-by-n. | |
2424 n = n(end); | |
2425 | |
2426 t = tril (eye (m, n) + alpha * triu (ones (m, n), 1), k); | |
2427 endfunction | |
2428 | |
2429 function A = wathen (nx, ny, k = 0) | |
2430 ## # WATHEN returns the Wathen matrix. | |
2431 ## | |
2432 ## Discussion: | |
2433 ## | |
2434 ## The Wathen matrix is a finite element matrix which is sparse. | |
2435 ## | |
2436 ## The entries of the matrix depend in part on a physical quantity | |
2437 ## related to density. That density is here assigned random values between | |
2438 ## 0 and 100. | |
2439 ## | |
2440 ## A = WATHEN ( NX, NY ) is a sparse random N-by-N finite element matrix | |
2441 ## where N = 3*NX*NY + 2*NX + 2*NY + 1. | |
2442 ## | |
2443 ## A is the consistent mass matrix for a regular NX-by-NY | |
2444 ## grid of 8-node (serendipity) elements in 2 space dimensions. | |
2445 ## | |
2446 ## Here is an illustration for NX = 3, NX = 2: | |
2447 ## | |
2448 ## 23-24-25-26-27-28-29 | |
2449 ## | | | | | |
2450 ## 19 20 21 22 | |
2451 ## | | | | | |
2452 ## 12-13-14-15-16-17-18 | |
2453 ## | | | | | |
2454 ## 8 9 10 11 | |
2455 ## | | | | | |
2456 ## 1--2--3--4--5--6--7 | |
2457 ## | |
2458 ## For this example, the total number of nodes is, as expected, | |
2459 ## | |
2460 ## N = 3 * 3 * 2 + 2 * 2 + 2 * 3 + 1 = 29. | |
2461 ## | |
2462 ## A is symmetric positive definite for any (positive) values of | |
2463 ## the density, RHO(NX,NY), which is chosen randomly in this routine. | |
2464 ## | |
2465 ## In particular, if D = DIAG(DIAG(A)), then | |
2466 ## 0.25 <= EIG(INV(D)*A) <= 4.5 | |
2467 ## for any positive integers NX and NY and any densities RHO(NX,NY). | |
2468 ## | |
2469 ## A = WATHEN ( NX, NY, 1 ) returns the diagonally scaled matrix. | |
2470 ## | |
2471 ## Modified: | |
2472 ## | |
2473 ## 17 September 2007 | |
2474 ## | |
2475 ## Author: | |
2476 ## | |
2477 ## Nicholas Higham | |
2478 ## | |
2479 ## Reference: | |
2480 ## | |
2481 ## Nicholas Higham, | |
2482 ## Algorithm 694: A Collection of Test Matrices in MATLAB, | |
2483 ## ACM Transactions on Mathematical Software, | |
2484 ## Volume 17, Number 3, September 1991, pages 289-305. | |
2485 ## | |
2486 ## Andrew Wathen, | |
2487 ## Realistic eigenvalue bounds for the Galerkin mass matrix, | |
2488 ## IMA Journal of Numerical Analysis, | |
2489 ## Volume 7, 1987, pages 449-457. | |
2490 ## | |
2491 ## Parameters: | |
2492 ## | |
2493 ## Input, integer NX, NY, the number of elements in the X and Y directions | |
2494 ## of the finite element grid. NX and NY must each be at least 1. | |
2495 ## | |
2496 ## Optional input, integer K, is used to request that the diagonally scaled | |
2497 ## version of the matrix be returned. This happens if K is specified with | |
2498 ## the value 1. | |
2499 ## | |
2500 ## Output, sparse real A(N,N), the matrix. The dimension N is determined by | |
2501 ## NX and NY, as described above. A is stored in the MATLAB sparse matrix | |
2502 ## format. | |
2503 | |
2504 if (nargin < 2 || nargin > 3) | |
2505 error ("gallery: 2 or 3 arguments are required for wathen matrix."); | |
2506 elseif (! isnumeric (nx) || ! isscalar (nx) || nx < 1) | |
2507 error ("gallery: NX must be a positive scalar for wathen matrix."); | |
2508 elseif (! isnumeric (ny) || ! isscalar (ny) || ny < 1) | |
2509 error ("gallery: NY must be a positive scalar for wathen matrix."); | |
2510 elseif (! isscalar (k)) | |
2511 error ("gallery: K must be a scalar for wathen matrix."); | |
2512 endif | |
2513 | |
2514 e1 = [ 6 -6 2 -8 | |
2515 -6 32 -6 20 | |
2516 2 -6 6 -6 | |
2517 -8 20 -6 32 ]; | |
2518 | |
2519 e2 = [ 3 -8 2 -6 | |
2520 -8 16 -8 20 | |
2521 2 -8 3 -8 | |
2522 -6 20 -8 16 ]; | |
2523 | |
2524 e = [ e1 e2 | |
2525 e2' e1] / 45; | |
2526 | |
2527 n = 3*nx*ny + 2*nx + 2*ny + 1; | |
2528 | |
2529 A = sparse (n, n); | |
2530 | |
2531 rho = 100 * rand (nx, ny); | |
2532 | |
2533 for j = 1:ny | |
2534 for i = 1:nx | |
2535 ## | |
2536 ## For the element (I,J), determine the indices of the 8 nodes. | |
2537 ## | |
2538 nn(1) = 3*j*nx + 2*i + 2*j + 1; | |
2539 nn(2) = nn(1) - 1; | |
2540 nn(3) = nn(2) - 1; | |
2541 nn(4) = (3*j - 1) * nx + 2*j + i - 1; | |
2542 nn(5) = 3 * (j-1) * nx + 2*i + 2*j - 3; | |
2543 nn(6) = nn(5) + 1; | |
2544 nn(7) = nn(6) + 1; | |
2545 nn(8) = nn(4) + 1; | |
2546 | |
2547 em = e * rho(i,j); | |
2548 | |
2549 for krow = 1:8 | |
2550 for kcol = 1:8 | |
2551 A(nn(krow),nn(kcol)) = A(nn(krow),nn(kcol)) + em(krow,kcol); | |
2552 endfor | |
2553 endfor | |
2554 | |
2555 endfor | |
2556 endfor | |
2557 | |
2558 ## If requested, return A with diagonal scaling. | |
2559 if (k) | |
2560 A = diag (diag (A)) \ A; | |
2561 endif | |
2562 endfunction | |
2563 | |
2564 function [A, b] = wilk (n) | |
2565 ## WILK Various specific matrices devised/discussed by Wilkinson. | |
2566 ## [A, b] = WILK(N) is the matrix or system of order N. | |
2567 ## N = 3: upper triangular system Ux=b illustrating inaccurate solution. | |
2568 ## N = 4: lower triangular system Lx=b, ill-conditioned. | |
2569 ## N = 5: HILB(6)(1:5,2:6)*1.8144. Symmetric positive definite. | |
2570 ## N = 21: W21+, tridiagonal. Eigenvalue problem. | |
2571 ## | |
2572 ## References: | |
2573 ## J.H. Wilkinson, Error analysis of direct methods of matrix inversion, | |
2574 ## J. Assoc. Comput. Mach., 8 (1961), pp. 281-330. | |
2575 ## J.H. Wilkinson, Rounding Errors in Algebraic Processes, Notes on Applied | |
2576 ## Science No. 32, Her Majesty's Stationery Office, London, 1963. | |
2577 ## J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University | |
2578 ## Press, 1965. | |
2579 | |
2580 if (nargin != 1) | |
2581 error ("gallery: 1 argument is required for wilk matrix."); | |
2582 elseif (! isnumeric (n) || ! isscalar (n)) | |
2583 error ("gallery: N must be a numeric scalar for wilk matrix."); | |
2584 endif | |
2585 | |
2586 if (n == 3) | |
2587 ## Wilkinson (1961) p.323. | |
2588 A = [ 1e-10 0.9 -0.4 | |
2589 0 0.9 -0.4 | |
2590 0 0 1e-10 ]; | |
2591 | |
2592 b = [ 0 | |
2593 0 | |
2594 1]; | |
2595 | |
2596 elseif (n == 4) | |
2597 ## Wilkinson (1963) p.105. | |
2598 A = [0.9143e-4 0 0 0 | |
2599 0.8762 0.7156e-4 0 0 | |
2600 0.7943 0.8143 0.9504e-4 0 | |
2601 0.8017 0.6123 0.7165 0.7123e-4]; | |
2602 | |
2603 b = [0.6524 | |
2604 0.3127 | |
2605 0.4186 | |
2606 0.7853]; | |
2607 | |
2608 elseif (n == 5) | |
2609 ## Wilkinson (1965), p.234. | |
2610 A = hilb (6); | |
2611 A = A(1:5, 2:6) * 1.8144; | |
2612 | |
2613 elseif (n == 21) | |
2614 ## Wilkinson (1965), p.308. | |
2615 E = diag (ones (n-1, 1), 1); | |
2616 m = (n-1)/2; | |
2617 A = diag (abs (-m:m)) + E + E'; | |
2618 | |
2619 else | |
16766
7268845c0a1e
avoid backquote in error messages, some uses in doc strings
John W. Eaton <jwe@octave.org>
parents:
16734
diff
changeset
|
2620 error ("gallery: unknown N '%d' for wilk matrix.", n); |
16634 | 2621 endif |
2622 endfunction | |
2623 | |
2624 ## NOTE: bandred is part of the Test Matrix Toolbox and is used by randsvd() | |
2625 function A = bandred (A, kl, ku) | |
2626 ## BANDRED Band reduction by two-sided unitary transformations. | |
2627 ## B = BANDRED(A, KL, KU) is a matrix unitarily equivalent to A | |
2628 ## with lower bandwidth KL and upper bandwidth KU | |
2629 ## (i.e. B(i,j) = 0 if i > j+KL or j > i+KU). | |
2630 ## The reduction is performed using Householder transformations. | |
2631 ## If KU is omitted it defaults to KL. | |
2632 ## | |
2633 ## Called by RANDSVD. | |
2634 ## This is a `standard' reduction. Cf. reduction to bidiagonal form | |
2635 ## prior to computing the SVD. This code is a little wasteful in that | |
2636 ## it computes certain elements which are immediately set to zero! | |
2637 ## | |
2638 ## Reference: | |
2639 ## G.H. Golub and C.F. Van Loan, Matrix Computations, second edition, | |
2640 ## Johns Hopkins University Press, Baltimore, Maryland, 1989. | |
2641 ## Section 5.4.3. | |
2642 | |
2643 ## Check for special case where order of left/right transformations matters. | |
2644 ## Easiest approach is to work on the transpose, flipping back at the end. | |
2645 flip = false; | |
2646 if (ku == 0) | |
2647 flip = true; | |
2648 A = A'; | |
2649 [ku, kl] = deal (kl, ku); | |
2650 endif | |
2651 | |
2652 [m, n] = size (A); | |
2653 | |
2654 for j = 1:min (min (m, n), max (m-kl-1, n-ku-1)) | |
2655 if (j+kl+1 <= m) | |
2656 [v, beta] = house (A(j+kl:m,j)); | |
2657 temp = A(j+kl:m,j:n); | |
2658 A(j+kl:m,j:n) = temp - beta*v*(v'*temp); | |
2659 A(j+kl+1:m,j) = zeros (m-j-kl, 1); | |
2660 endif | |
2661 | |
2662 if (j+ku+1 <= n) | |
2663 [v, beta] = house (A(j,j+ku:n)'); | |
2664 temp = A(j:m,j+ku:n); | |
2665 A(j:m,j+ku:n) = temp - beta*(temp*v)*v'; | |
2666 A(j,j+ku+1:n) = zeros (1, n-j-ku); | |
2667 endif | |
2668 endfor | |
2669 | |
2670 if (flip) | |
2671 A = A'; | |
2672 endif | |
2673 endfunction |