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