24
File: octave.info, Node: Special Utility Matrices, Next: Famous Matrices, Prev: Rearranging Matrices, Up: Matrix Manipulation
26
16.3 Special Utility Matrices
27
=============================
29
-- Built-in Function: eye (N)
30
-- Built-in Function: eye (M, N)
31
-- Built-in Function: eye ([M N])
32
-- Built-in Function: eye (..., CLASS)
33
Return an identity matrix. If invoked with a single scalar
34
argument N, return a square NxN identity matrix. If supplied two
35
scalar arguments (M, N), 'eye' takes them to be the number of rows
36
and columns. If given a vector with two elements, 'eye' uses the
37
values of the elements as the number of rows and columns,
38
respectively. For example:
45
The following expressions all produce the same result:
51
eye (size ([1, 2; 3, 4])
53
The optional argument CLASS, allows 'eye' to return an array of the
56
val = zeros (n,m, "uint8")
58
Calling 'eye' with no arguments is equivalent to calling it with an
59
argument of 1. Any negative dimensions are treated as zero. These
60
odd definitions are for compatibility with MATLAB.
62
See also: *note speye: XREFspeye, *note ones: XREFones, *note
65
-- Built-in Function: ones (N)
66
-- Built-in Function: ones (M, N)
67
-- Built-in Function: ones (M, N, K, ...)
68
-- Built-in Function: ones ([M N ...])
69
-- Built-in Function: ones (..., CLASS)
70
Return a matrix or N-dimensional array whose elements are all 1.
71
If invoked with a single scalar integer argument N, return a square
72
NxN matrix. If invoked with two or more scalar integer arguments,
73
or a vector of integer values, return an array with the given
76
If you need to create a matrix whose values are all the same, you
77
should use an expression like
79
val_matrix = val * ones (m, n)
81
The optional argument CLASS specifies the class of the return array
82
and defaults to double. For example:
84
val = ones (m,n, "uint8")
86
See also: *note zeros: XREFzeros.
88
-- Built-in Function: zeros (N)
89
-- Built-in Function: zeros (M, N)
90
-- Built-in Function: zeros (M, N, K, ...)
91
-- Built-in Function: zeros ([M N ...])
92
-- Built-in Function: zeros (..., CLASS)
93
Return a matrix or N-dimensional array whose elements are all 0.
94
If invoked with a single scalar integer argument, return a square
95
NxN matrix. If invoked with two or more scalar integer arguments,
96
or a vector of integer values, return an array with the given
99
The optional argument CLASS specifies the class of the return array
100
and defaults to double. For example:
102
val = zeros (m,n, "uint8")
104
See also: *note ones: XREFones.
106
-- Function File: repmat (A, M)
107
-- Function File: repmat (A, M, N)
108
-- Function File: repmat (A, [M N])
109
-- Function File: repmat (A, [M N P ...])
110
Form a block matrix of size M by N, with a copy of matrix A as each
111
element. If N is not specified, form an M by M block matrix. For
112
copying along more than two dimensions, specify the number of times
113
to copy across each dimension M, N, P, ..., in a vector in the
116
See also: *note repelems: XREFrepelems.
118
-- Built-in Function: repelems (X, R)
119
Construct a vector of repeated elements from X. R is a 2xN integer
120
matrix specifying which elements to repeat and how often to repeat
123
Entries in the first row, R(1,j), select an element to repeat. The
124
corresponding entry in the second row, R(2,j), specifies the repeat
125
count. If X is a matrix then the columns of X are imagined to be
126
stacked on top of each other for purposes of the selection index.
127
A row vector is always returned.
129
Conceptually the result is calculated as follows:
132
for i = 1:columns (R)
133
y = [y, X(R(1,i)*ones(1, R(2,i)))];
136
See also: *note repmat: XREFrepmat, *note cat: XREFcat.
138
The functions 'linspace' and 'logspace' make it very easy to create
139
vectors with evenly or logarithmically spaced elements. *Note Ranges::.
141
-- Built-in Function: linspace (BASE, LIMIT)
142
-- Built-in Function: linspace (BASE, LIMIT, N)
143
Return a row vector with N linearly spaced elements between BASE
144
and LIMIT. If the number of elements is greater than one, then the
145
endpoints BASE and LIMIT are always included in the range. If BASE
146
is greater than LIMIT, the elements are stored in decreasing order.
147
If the number of points is not specified, a value of 100 is used.
149
The 'linspace' function always returns a row vector if both BASE
150
and LIMIT are scalars. If one, or both, of them are column
151
vectors, 'linspace' returns a matrix.
153
For compatibility with MATLAB, return the second argument (LIMIT)
154
if fewer than two values are requested.
156
See also: *note logspace: XREFlogspace.
158
-- Function File: logspace (A, B)
159
-- Function File: logspace (A, B, N)
160
-- Function File: logspace (A, pi, N)
161
Return a row vector with N elements logarithmically spaced from
162
10^A to 10^B. If N is unspecified it defaults to 50.
164
If B is equal to pi, the points are between 10^A and pi, _not_ 10^A
165
and 10^pi, in order to be compatible with the corresponding MATLAB
168
Also for compatibility with MATLAB, return the second argument B if
169
fewer than two values are requested.
171
See also: *note linspace: XREFlinspace.
173
-- Built-in Function: rand (N)
174
-- Built-in Function: rand (M, N, ...)
175
-- Built-in Function: rand ([M N ...])
176
-- Built-in Function: V = rand ("state")
177
-- Built-in Function: rand ("state", V)
178
-- Built-in Function: rand ("state", "reset")
179
-- Built-in Function: V = rand ("seed")
180
-- Built-in Function: rand ("seed", V)
181
-- Built-in Function: rand ("seed", "reset")
182
-- Built-in Function: rand (..., "single")
183
-- Built-in Function: rand (..., "double")
184
Return a matrix with random elements uniformly distributed on the
185
interval (0, 1). The arguments are handled the same as the
188
You can query the state of the random number generator using the
193
This returns a column vector V of length 625. Later, you can
194
restore the random number generator to the state V using the form
198
You may also initialize the state vector from an arbitrary vector
199
of length <= 625 for V. This new state will be a hash based on the
200
value of V, not V itself.
202
By default, the generator is initialized from '/dev/urandom' if it
203
is available, otherwise from CPU time, wall clock time, and the
204
current fraction of a second. Note that this differs from MATLAB,
205
which always initializes the state to the same state at startup.
206
To obtain behavior comparable to MATLAB, initialize with a
207
deterministic state vector in Octave's startup files (*note Startup
210
To compute the pseudo-random sequence, 'rand' uses the Mersenne
211
Twister with a period of 2^{19937}-1 (See M. Matsumoto and T.
212
Nishimura, 'Mersenne Twister: A 623-dimensionally equidistributed
213
uniform pseudorandom number generator', ACM Trans. on Modeling and
214
Computer Simulation Vol. 8, No. 1, pp. 3-30, January 1998,
215
<http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html>). Do
216
*not* use for cryptography without securely hashing several
217
returned values together, otherwise the generator state can be
218
learned after reading 624 consecutive values.
220
Older versions of Octave used a different random number generator.
221
The new generator is used by default as it is significantly faster
222
than the old generator, and produces random numbers with a
223
significantly longer cycle time. However, in some circumstances it
224
might be desirable to obtain the same random sequences as used by
225
the old generators. To do this the keyword "seed" is used to
226
specify that the old generators should be use, as in
230
which sets the seed of the generator to VAL. The seed of the
231
generator can be queried with
235
However, it should be noted that querying the seed will not cause
236
'rand' to use the old generators, only setting the seed will. To
237
cause 'rand' to once again use the new generators, the keyword
238
"state" should be used to reset the state of the 'rand'.
240
The state or seed of the generator can be reset to a new random
241
value using the "reset" keyword.
243
The class of the value returned can be controlled by a trailing
244
"double" or "single" argument. These are the only valid classes.
246
See also: *note randn: XREFrandn, *note rande: XREFrande, *note
247
randg: XREFrandg, *note randp: XREFrandp.
249
-- Function File: randi (IMAX)
250
-- Function File: randi (IMAX, N)
251
-- Function File: randi (IMAX, M, N, ...)
252
-- Function File: randi ([IMIN IMAX], ...)
253
-- Function File: randi (..., "CLASS")
254
Return random integers in the range 1:IMAX.
256
Additional arguments determine the shape of the return matrix.
257
When no arguments are specified a single random integer is
258
returned. If one argument N is specified then a square matrix (N x N)
259
is returned. Two or more arguments will return a multi-dimensional
260
matrix (M x N x ...).
262
The integer range may optionally be described by a two element
263
matrix with a lower and upper bound in which case the returned
264
integers will be on the interval [IMIN, IMAX].
266
The optional argument CLASS will return a matrix of the requested
267
type. The default is "double".
269
The following example returns 150 integers in the range 1-10.
271
ri = randi (10, 150, 1)
273
Implementation Note: 'randi' relies internally on 'rand' which uses
274
class "double" to represent numbers. This limits the maximum
275
integer (IMAX) and range (IMAX - IMIN) to the value returned by the
276
'bitmax' function. For IEEE floating point numbers this value is
279
See also: *note rand: XREFrand.
281
-- Built-in Function: randn (N)
282
-- Built-in Function: randn (M, N, ...)
283
-- Built-in Function: randn ([M N ...])
284
-- Built-in Function: V = randn ("state")
285
-- Built-in Function: randn ("state", V)
286
-- Built-in Function: randn ("state", "reset")
287
-- Built-in Function: V = randn ("seed")
288
-- Built-in Function: randn ("seed", V)
289
-- Built-in Function: randn ("seed", "reset")
290
-- Built-in Function: randn (..., "single")
291
-- Built-in Function: randn (..., "double")
292
Return a matrix with normally distributed random elements having
293
zero mean and variance one. The arguments are handled the same as
294
the arguments for 'rand'.
296
By default, 'randn' uses the Marsaglia and Tsang "Ziggurat
297
technique" to transform from a uniform to a normal distribution.
299
The class of the value returned can be controlled by a trailing
300
"double" or "single" argument. These are the only valid classes.
302
Reference: G. Marsaglia and W.W. Tsang, 'Ziggurat Method for
303
Generating Random Variables', J. Statistical Software, vol 5, 2000,
304
<http://www.jstatsoft.org/v05/i08/>)
306
See also: *note rand: XREFrand, *note rande: XREFrande, *note
307
randg: XREFrandg, *note randp: XREFrandp.
309
-- Built-in Function: rande (N)
310
-- Built-in Function: rande (M, N, ...)
311
-- Built-in Function: rande ([M N ...])
312
-- Built-in Function: V = rande ("state")
313
-- Built-in Function: rande ("state", V)
314
-- Built-in Function: rande ("state", "reset")
315
-- Built-in Function: V = rande ("seed")
316
-- Built-in Function: rande ("seed", V)
317
-- Built-in Function: rande ("seed", "reset")
318
-- Built-in Function: rande (..., "single")
319
-- Built-in Function: rande (..., "double")
320
Return a matrix with exponentially distributed random elements.
321
The arguments are handled the same as the arguments for 'rand'.
323
By default, 'randn' uses the Marsaglia and Tsang "Ziggurat
324
technique" to transform from a uniform to an exponential
327
The class of the value returned can be controlled by a trailing
328
"double" or "single" argument. These are the only valid classes.
330
Reference: G. Marsaglia and W.W. Tsang, 'Ziggurat Method for
331
Generating Random Variables', J. Statistical Software, vol 5, 2000,
332
<http://www.jstatsoft.org/v05/i08/>)
334
See also: *note rand: XREFrand, *note randn: XREFrandn, *note
335
randg: XREFrandg, *note randp: XREFrandp.
337
-- Built-in Function: randp (L, N)
338
-- Built-in Function: randp (L, M, N, ...)
339
-- Built-in Function: randp (L, [M N ...])
340
-- Built-in Function: V = randp ("state")
341
-- Built-in Function: randp ("state", V)
342
-- Built-in Function: randp ("state", "reset")
343
-- Built-in Function: V = randp ("seed")
344
-- Built-in Function: randp ("seed", V)
345
-- Built-in Function: randp ("seed", "reset")
346
-- Built-in Function: randp (..., "single")
347
-- Built-in Function: randp (..., "double")
348
Return a matrix with Poisson distributed random elements with mean
349
value parameter given by the first argument, L. The arguments are
350
handled the same as the arguments for 'rand', except for the
353
Five different algorithms are used depending on the range of L and
354
whether or not L is a scalar or a matrix.
356
For scalar L <= 12, use direct method.
357
W.H. Press, et al., 'Numerical Recipes in C', Cambridge
358
University Press, 1992.
360
For scalar L > 12, use rejection method.[1]
361
W.H. Press, et al., 'Numerical Recipes in C', Cambridge
362
University Press, 1992.
364
For matrix L <= 10, use inversion method.[2]
365
E. Stadlober, et al., WinRand source code, available via FTP.
367
For matrix L > 10, use patchwork rejection method.
368
E. Stadlober, et al., WinRand source code, available via FTP,
369
or H. Zechner, 'Efficient sampling from continuous and
370
discrete unimodal distributions', Doctoral Dissertation,
371
156pp., Technical University Graz, Austria, 1994.
373
For L > 1e8, use normal approximation.
374
L. Montanet, et al., 'Review of Particle Properties', Physical
375
Review D 50 p1284, 1994.
377
The class of the value returned can be controlled by a trailing
378
"double" or "single" argument. These are the only valid classes.
380
See also: *note rand: XREFrand, *note randn: XREFrandn, *note
381
rande: XREFrande, *note randg: XREFrandg.
383
-- Built-in Function: randg (N)
384
-- Built-in Function: randg (M, N, ...)
385
-- Built-in Function: randg ([M N ...])
386
-- Built-in Function: V = randg ("state")
387
-- Built-in Function: randg ("state", V)
388
-- Built-in Function: randg ("state", "reset")
389
-- Built-in Function: V = randg ("seed")
390
-- Built-in Function: randg ("seed", V)
391
-- Built-in Function: randg ("seed", "reset")
392
-- Built-in Function: randg (..., "single")
393
-- Built-in Function: randg (..., "double")
394
Return a matrix with 'gamma (A,1)' distributed random elements.
395
The arguments are handled the same as the arguments for 'rand',
396
except for the argument A.
398
This can be used to generate many distributions:
400
'gamma (a, b)' for 'a > -1', 'b > 0'
404
'beta (a, b)' for 'a > -1', 'b > -1'
407
r = r1 / (r1 + randg (b, 1))
413
'chisq (df)' for 'df > 0'
415
r = 2 * randg (df / 2)
417
't (df)' for '0 < df < inf' (use randn if df is infinite)
419
r = randn () / sqrt (2 * randg (df / 2) / df)
421
'F (n1, n2)' for '0 < n1', '0 < n2'
423
## r1 equals 1 if n1 is infinite
424
r1 = 2 * randg (n1 / 2) / n1
425
## r2 equals 1 if n2 is infinite
426
r2 = 2 * randg (n2 / 2) / n2
430
negative 'binomial (n, p)' for 'n > 0', '0 < p <= 1'
432
r = randp ((1 - p) / p * randg (n))
434
non-central 'chisq (df, L)', for 'df >= 0' and 'L > 0'
435
(use chisq if 'L = 0')
438
r(r > 0) = 2 * randg (r(r > 0))
439
r(df > 0) += 2 * randg (df(df > 0)/2)
441
'Dirichlet (a1, ... ak)'
443
r = (randg (a1), ..., randg (ak))
446
The class of the value returned can be controlled by a trailing
447
"double" or "single" argument. These are the only valid classes.
449
See also: *note rand: XREFrand, *note randn: XREFrandn, *note
450
rande: XREFrande, *note randp: XREFrandp.
452
The generators operate in the new or old style together, it is not
453
possible to mix the two. Initializing any generator with "state" or
454
"seed" causes the others to switch to the same style for future calls.
456
The state of each generator is independent and calls to different
457
generators can be interleaved without affecting the final result. For
460
rand ("state", [11, 22, 33]);
461
randn ("state", [44, 55, 66]);
467
rand ("state", [11, 22, 33]);
468
randn ("state", [44, 55, 66]);
476
produce equivalent results. When the generators are initialized in the
477
old style with "seed" only 'rand' and 'randn' are independent, because
478
the old 'rande', 'randg' and 'randp' generators make calls to 'rand' and
481
The generators are initialized with random states at start-up, so
482
that the sequences of random numbers are not the same each time you run
483
Octave.(1) If you really do need to reproduce a sequence of numbers
484
exactly, you can set the state or seed to a specific value.
486
If invoked without arguments, 'rand' and 'randn' return a single
487
element of a random sequence.
489
The original 'rand' and 'randn' functions use Fortran code from
490
RANLIB, a library of Fortran routines for random number generation,
491
compiled by Barry W. Brown and James Lovato of the Department of
492
Biomathematics at The University of Texas, M.D. Anderson Cancer Center,
495
-- Built-in Function: randperm (N)
496
-- Built-in Function: randperm (N, M)
497
Return a row vector containing a random permutation of '1:N'. If M
498
is supplied, return M unique entries, sampled without replacement
499
from '1:N'. The complexity is O(N) in memory and O(M) in time,
500
unless M < N/5, in which case O(M) memory is used as well. The
501
randomization is performed using rand(). All permutations are
504
See also: *note perms: XREFperms.
506
---------- Footnotes ----------
508
(1) The old versions of 'rand' and 'randn' obtain their initial seeds
509
from the system clock.
24
512
File: octave.info, Node: Famous Matrices, Prev: Special Utility Matrices, Up: Matrix Manipulation
26
514
16.4 Famous Matrices
7310
7806
(1) The CHOLMOD, UMFPACK and CXSPARSE packages were written by Tim
7311
7807
Davis and are available at <http://www.cise.ufl.edu/research/sparse/>
7314
File: octave.info, Node: Iterative Techniques, Next: Real Life Example, Prev: Sparse Linear Algebra, Up: Sparse Matrices
7316
22.3 Iterative Techniques Applied to Sparse Matrices
7317
====================================================
7319
The left division '\' and right division '/' operators, discussed in the
7320
previous section, use direct solvers to resolve a linear equation of the
7321
form 'X = A \ B' or 'X = B / A'. Octave equally includes a number of
7322
functions to solve sparse linear equations using iterative techniques.
7324
-- Function File: X = pcg (A, B, TOL, MAXIT, M1, M2, X0, ...)
7325
-- Function File: [X, FLAG, RELRES, ITER, RESVEC, EIGEST] = pcg (...)
7327
Solve the linear system of equations 'A * X = B' by means of the
7328
Preconditioned Conjugate Gradient iterative method. The input
7331
* A can be either a square (preferably sparse) matrix or a
7332
function handle, inline function or string containing the name
7333
of a function which computes 'A * X'. In principle, A should
7334
be symmetric and positive definite; if 'pcg' finds A not to be
7335
positive definite, a warning is printed and the FLAG output
7338
* B is the right-hand side vector.
7340
* TOL is the required relative tolerance for the residual error,
7341
'B - A * X'. The iteration stops if
7342
'norm (B - A * X)' <= TOL * norm (B). If TOL is omitted or
7343
empty then a tolerance of 1e-6 is used.
7345
* MAXIT is the maximum allowable number of iterations; if MAXIT
7346
is omitted or empty then a value of 20 is used.
7348
* M = M1 * M2 is the (left) preconditioning matrix, so that the
7349
iteration is (theoretically) equivalent to solving by 'pcg'
7350
'P * X = M \ B', with 'P = M \ A'. Note that a proper choice
7351
of the preconditioner may dramatically improve the overall
7352
performance of the method. Instead of matrices M1 and M2, the
7353
user may pass two functions which return the results of
7354
applying the inverse of M1 and M2 to a vector (usually this is
7355
the preferred way of using the preconditioner). If M1 is
7356
omitted or empty '[]' then no preconditioning is applied. If
7357
M2 is omitted, M = M1 will be used as a preconditioner.
7359
* X0 is the initial guess. If X0 is omitted or empty then the
7360
function sets X0 to a zero vector by default.
7362
The arguments which follow X0 are treated as parameters, and passed
7363
in a proper way to any of the functions (A or M) which are passed
7364
to 'pcg'. See the examples below for further details. The output
7367
* X is the computed approximation to the solution of
7370
* FLAG reports on the convergence. A value of 0 means the
7371
solution converged and the tolerance criterion given by TOL is
7372
satisfied. A value of 1 means that the MAXIT limit for the
7373
iteration count was reached. A value of 3 indicates that the
7374
(preconditioned) matrix was found not to be positive definite.
7376
* RELRES is the ratio of the final residual to its initial
7377
value, measured in the Euclidean norm.
7379
* ITER is the actual number of iterations performed.
7381
* RESVEC describes the convergence history of the method.
7382
'RESVEC(i,1)' is the Euclidean norm of the residual, and
7383
'RESVEC(i,2)' is the preconditioned residual norm, after the
7384
(I-1)-th iteration, 'I = 1, 2, ..., ITER+1'. The
7385
preconditioned residual norm is defined as 'norm (R) ^ 2 = R'
7386
* (M \ R)' where 'R = B - A * X', see also the description of
7387
M. If EIGEST is not required, only 'RESVEC(:,1)' is returned.
7389
* EIGEST returns the estimate for the smallest 'EIGEST(1)' and
7390
largest 'EIGEST(2)' eigenvalues of the preconditioned matrix
7391
'P = M \ A'. In particular, if no preconditioning is used,
7392
the estimates for the extreme eigenvalues of A are returned.
7393
'EIGEST(1)' is an overestimate and 'EIGEST(2)' is an
7394
underestimate, so that 'EIGEST(2) / EIGEST(1)' is a lower
7395
bound for 'cond (P, 2)', which nevertheless in the limit
7396
should theoretically be equal to the actual value of the
7397
condition number. The method which computes EIGEST works only
7398
for symmetric positive definite A and M, and the user is
7399
responsible for verifying this assumption.
7401
Let us consider a trivial problem with a diagonal matrix (we
7402
exploit the sparsity of A)
7405
A = diag (sparse (1:n));
7407
[l, u, p, q] = luinc (A, 1.e-3);
7409
EXAMPLE 1: Simplest use of 'pcg'
7413
EXAMPLE 2: 'pcg' with a function which computes 'A * X'
7415
function y = apply_a (x)
7419
x = pcg ("apply_a", b)
7421
EXAMPLE 3: 'pcg' with a preconditioner: L * U
7423
x = pcg (A, b, 1.e-6, 500, l*u)
7425
EXAMPLE 4: 'pcg' with a preconditioner: L * U. Faster than EXAMPLE
7426
3 since lower and upper triangular matrices are easier to invert
7428
x = pcg (A, b, 1.e-6, 500, l, u)
7430
EXAMPLE 5: Preconditioned iteration, with full diagnostics. The
7431
preconditioner (quite strange, because even the original matrix A
7432
is trivial) is defined as a function
7434
function y = apply_m (x)
7435
k = floor (length (x) - 2);
7437
y(1:k) = x(1:k) ./ [1:k]';
7440
[x, flag, relres, iter, resvec, eigest] = ...
7441
pcg (A, b, [], [], "apply_m");
7442
semilogy (1:iter+1, resvec);
7444
EXAMPLE 6: Finally, a preconditioner which depends on a parameter
7447
function y = apply_M (x, varargin)
7450
y(1:K) = x(1:K) ./ [1:K]';
7453
[x, flag, relres, iter, resvec, eigest] = ...
7454
pcg (A, b, [], [], "apply_m", [], [], 3)
7458
1. C.T. Kelley, 'Iterative Methods for Linear and Nonlinear
7459
Equations', SIAM, 1995. (the base PCG algorithm)
7461
2. Y. Saad, 'Iterative Methods for Sparse Linear Systems', PWS
7462
1996. (condition number estimate from PCG) Revised version of
7463
this book is available online at
7464
<http://www-users.cs.umn.edu/~saad/books.html>
7466
See also: *note sparse: XREFsparse, *note pcr: XREFpcr.
7468
-- Function File: X = pcr (A, B, TOL, MAXIT, M, X0, ...)
7469
-- Function File: [X, FLAG, RELRES, ITER, RESVEC] = pcr (...)
7471
Solve the linear system of equations 'A * X = B' by means of the
7472
Preconditioned Conjugate Residuals iterative method. The input
7475
* A can be either a square (preferably sparse) matrix or a
7476
function handle, inline function or string containing the name
7477
of a function which computes 'A * X'. In principle A should
7478
be symmetric and non-singular; if 'pcr' finds A to be
7479
numerically singular, you will get a warning message and the
7480
FLAG output parameter will be set.
7482
* B is the right hand side vector.
7484
* TOL is the required relative tolerance for the residual error,
7485
'B - A * X'. The iteration stops if 'norm (B - A * X) <= TOL
7486
* norm (B - A * X0)'. If TOL is empty or is omitted, the
7487
function sets 'TOL = 1e-6' by default.
7489
* MAXIT is the maximum allowable number of iterations; if '[]'
7490
is supplied for 'maxit', or 'pcr' has less arguments, a
7491
default value equal to 20 is used.
7493
* M is the (left) preconditioning matrix, so that the iteration
7494
is (theoretically) equivalent to solving by 'pcr' 'P * X = M \
7495
B', with 'P = M \ A'. Note that a proper choice of the
7496
preconditioner may dramatically improve the overall
7497
performance of the method. Instead of matrix M, the user may
7498
pass a function which returns the results of applying the
7499
inverse of M to a vector (usually this is the preferred way of
7500
using the preconditioner). If '[]' is supplied for M, or M is
7501
omitted, no preconditioning is applied.
7503
* X0 is the initial guess. If X0 is empty or omitted, the
7504
function sets X0 to a zero vector by default.
7506
The arguments which follow X0 are treated as parameters, and passed
7507
in a proper way to any of the functions (A or M) which are passed
7508
to 'pcr'. See the examples below for further details. The output
7511
* X is the computed approximation to the solution of 'A * X =
7514
* FLAG reports on the convergence. 'FLAG = 0' means the
7515
solution converged and the tolerance criterion given by TOL is
7516
satisfied. 'FLAG = 1' means that the MAXIT limit for the
7517
iteration count was reached. 'FLAG = 3' reports t 'pcr'
7518
breakdown, see [1] for details.
7520
* RELRES is the ratio of the final residual to its initial
7521
value, measured in the Euclidean norm.
7523
* ITER is the actual number of iterations performed.
7525
* RESVEC describes the convergence history of the method, so
7526
that 'RESVEC (i)' contains the Euclidean norms of the residual
7527
after the (I-1)-th iteration, 'I = 1,2, ..., ITER+1'.
7529
Let us consider a trivial problem with a diagonal matrix (we
7530
exploit the sparsity of A)
7533
A = sparse (diag (1:n));
7536
EXAMPLE 1: Simplest use of 'pcr'
7540
EXAMPLE 2: 'pcr' with a function which computes 'A * X'.
7542
function y = apply_a (x)
7546
x = pcr ("apply_a", b)
7548
EXAMPLE 3: Preconditioned iteration, with full diagnostics. The
7549
preconditioner (quite strange, because even the original matrix A
7550
is trivial) is defined as a function
7552
function y = apply_m (x)
7553
k = floor (length (x) - 2);
7555
y(1:k) = x(1:k) ./ [1:k]';
7558
[x, flag, relres, iter, resvec] = ...
7559
pcr (A, b, [], [], "apply_m")
7560
semilogy ([1:iter+1], resvec);
7562
EXAMPLE 4: Finally, a preconditioner which depends on a parameter
7565
function y = apply_m (x, varargin)
7568
y(1:k) = x(1:k) ./ [1:k]';
7571
[x, flag, relres, iter, resvec] = ...
7572
pcr (A, b, [], [], "apply_m"', [], 3)
7576
[1] W. Hackbusch, 'Iterative Solution of Large Sparse Systems of
7577
Equations', section 9.5.4; Springer, 1994
7579
See also: *note sparse: XREFsparse, *note pcg: XREFpcg.
7581
The speed with which an iterative solver converges to a solution can
7582
be accelerated with the use of a pre-conditioning matrix M. In this
7583
case the linear equation 'M^-1 * X = M^-1 * A \ B' is solved instead.
7584
Typical pre-conditioning matrices are partial factorizations of the
7587
-- Built-in Function: [L, U, P, Q] = luinc (A, '0')
7588
-- Built-in Function: [L, U, P, Q] = luinc (A, DROPTOL)
7589
-- Built-in Function: [L, U, P, Q] = luinc (A, OPTS)
7590
Produce the incomplete LU factorization of the sparse matrix A.
7591
Two types of incomplete factorization are possible, and the type is
7592
determined by the second argument to 'luinc'.
7594
Called with a second argument of '0', the zero-level incomplete
7595
LU factorization is produced. This creates a factorization of A
7596
where the position of the non-zero arguments correspond to the same
7597
positions as in the matrix A.
7599
Alternatively, the fill-in of the incomplete LU factorization can
7600
be controlled through the variable DROPTOL or the structure OPTS.
7601
The UMFPACK multifrontal factorization code by Tim A. Davis is used
7602
for the incomplete LU factorization, (availability
7603
<http://www.cise.ufl.edu/research/sparse/umfpack/>)
7605
DROPTOL determines the values below which the values in the LU
7606
factorization are dropped and replaced by zero. It must be a
7607
positive scalar, and any values in the factorization whose absolute
7608
value are less than this value are dropped, expect if leaving them
7609
increase the sparsity of the matrix. Setting DROPTOL to zero
7610
results in a complete LU factorization which is the default.
7612
OPTS is a structure containing one or more of the fields
7615
The drop tolerance as above. If OPTS only contains 'droptol'
7616
then this is equivalent to using the variable DROPTOL.
7619
A logical variable flagging whether to use the modified
7620
incomplete LU factorization. In the case that 'milu' is
7621
true, the dropped values are subtracted from the diagonal of
7622
the matrix U of the factorization. The default is 'false'.
7625
A logical variable that flags whether zero elements on the
7626
diagonal of U should be replaced with DROPTOL to attempt to
7627
avoid singular factors. The default is 'false'.
7630
Defines the pivot threshold in the interval [0,1]. Values
7631
outside that range are ignored.
7633
All other fields in OPTS are ignored. The outputs from 'luinc' are
7634
the same as for 'lu'.
7636
Given the string argument "vector", 'luinc' returns the values of P
7639
See also: *note sparse: XREFsparse, *note lu: XREFlu.
added
removed
doc/interpreter/octave.info-4
