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