Mercurial > octave
annotate scripts/sparse/spaugment.m @ 27919:1891570abac8
update Octave Project Developers copyright for the new year
In files that have the "Octave Project Developers" copyright notice,
update for 2020.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Mon, 06 Jan 2020 22:29:51 -0500 |
parents | b442ec6dda5c |
children | bd51beb6205e |
rev | line source |
---|---|
27919
1891570abac8
update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
27918
diff
changeset
|
1 ## Copyright (C) 2008-2020 The Octave Project Developers |
27918
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
2 ## |
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
3 ## See the file COPYRIGHT.md in the top-level directory of this distribution |
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
4 ## or <https://octave.org/COPYRIGHT.html/>. |
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
5 ## |
7681 | 6 ## |
7 ## This file is part of Octave. | |
8 ## | |
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
9 ## Octave is free software: you can redistribute it and/or modify it |
7681 | 10 ## under the terms of the GNU General Public License as published by |
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
11 ## the Free Software Foundation, either version 3 of the License, or |
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
12 ## (at your option) any later version. |
7681 | 13 ## |
14 ## Octave is distributed in the hope that it will be useful, but | |
15 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
16 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
17 ## GNU General Public License for more details. |
7681 | 18 ## |
19 ## You should have received a copy of the GNU General Public License | |
20 ## along with Octave; see the file COPYING. If not, see | |
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
21 ## <https://www.gnu.org/licenses/>. |
7681 | 22 |
23 ## -*- texinfo -*- | |
20852
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20711
diff
changeset
|
24 ## @deftypefn {} {@var{s} =} spaugment (@var{A}, @var{c}) |
18809
53af80da6781
doc: Update documentation of sparse functions including seealso links.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
25 ## Create the augmented matrix of @var{A}. |
53af80da6781
doc: Update documentation of sparse functions including seealso links.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
26 ## |
53af80da6781
doc: Update documentation of sparse functions including seealso links.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
27 ## This is given by |
7681 | 28 ## |
29 ## @example | |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
30 ## @group |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
31 ## [@var{c} * eye(@var{m}, @var{m}), @var{A}; |
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
32 ## @var{A}', zeros(@var{n}, @var{n})] |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
33 ## @end group |
7681 | 34 ## @end example |
35 ## | |
36 ## @noindent | |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
37 ## This is related to the least squares solution of |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
38 ## @code{@var{A} \ @var{b}}, by |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
39 ## |
7681 | 40 ## @example |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
41 ## @group |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
42 ## @var{s} * [ @var{r} / @var{c}; x] = [ @var{b}, zeros(@var{n}, columns(@var{b})) ] |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
43 ## @end group |
7681 | 44 ## @end example |
45 ## | |
46 ## @noindent | |
47 ## where @var{r} is the residual error | |
48 ## | |
49 ## @example | |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
50 ## @var{r} = @var{b} - @var{A} * @var{x} |
7681 | 51 ## @end example |
52 ## | |
20164
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
53 ## As the matrix @var{s} is symmetric indefinite it can be factorized with |
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
54 ## @code{lu}, and the minimum norm solution can therefore be found without the |
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
55 ## need for a @code{qr} factorization. As the residual error will be |
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
56 ## @code{zeros (@var{m}, @var{m})} for underdetermined problems, and example |
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
57 ## can be |
7681 | 58 ## |
59 ## @example | |
60 ## @group | |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
61 ## m = 11; n = 10; mn = max (m, n); |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
62 ## A = spdiags ([ones(mn,1), 10*ones(mn,1), -ones(mn,1)], |
8516 | 63 ## [-1, 0, 1], m, n); |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
64 ## x0 = A \ ones (m,1); |
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
65 ## s = spaugment (A); |
7681 | 66 ## [L, U, P, Q] = lu (s); |
67 ## x1 = Q * (U \ (L \ (P * [ones(m,1); zeros(n,1)]))); | |
68 ## x1 = x1(end - n + 1 : end); | |
69 ## @end group | |
70 ## @end example | |
71 ## | |
20164
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
72 ## To find the solution of an overdetermined problem needs an estimate of the |
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
73 ## residual error @var{r} and so it is more complex to formulate a minimum norm |
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
74 ## solution using the @code{spaugment} function. |
7681 | 75 ## |
20164
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
76 ## In general the left division operator is more stable and faster than using |
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
77 ## the @code{spaugment} function. |
18809
53af80da6781
doc: Update documentation of sparse functions including seealso links.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
78 ## @seealso{mldivide} |
7681 | 79 ## @end deftypefn |
80 | |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
81 function s = spaugment (A, c) |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
20852
diff
changeset
|
82 |
7681 | 83 if (nargin < 2) |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
84 if (issparse (A)) |
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
85 c = max (max (abs (A))) / 1000; |
7681 | 86 else |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
87 if (ndims (A) != 2) |
20711
7b608fadc663
Make error messages more specific about the variable and problem encountered.
Rik <rik@octave.org>
parents:
20164
diff
changeset
|
88 error ("spaugment: A must be a 2-D matrix"); |
7681 | 89 else |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
90 c = max (abs (A(:))) / 1000; |
7681 | 91 endif |
92 endif | |
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
93 elseif (! isscalar (c)) |
11472
1740012184f9
Use uppercase for variable names in error() strings to match Info documentation. Only m-files done.
Rik <octave@nomad.inbox5.com>
parents:
11471
diff
changeset
|
94 error ("spaugment: C must be a scalar"); |
7681 | 95 endif |
96 | |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
97 [m, n] = size (A); |
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
98 s = [ c * speye(m, m), A; A', sparse(n, n)]; |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
20852
diff
changeset
|
99 |
7681 | 100 endfunction |
101 | |
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
102 |
8871 | 103 %!testif HAVE_UMFPACK |
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
104 %! m = 11; n = 10; mn = max (m ,n); |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
105 %! A = spdiags ([ones(mn,1), 10*ones(mn,1), -ones(mn,1)],[-1,0,1], m, n); |
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
106 %! x0 = A \ ones (m,1); |
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
107 %! s = spaugment (A); |
7681 | 108 %! [L, U, P, Q] = lu (s); |
109 %! x1 = Q * (U \ (L \ (P * [ones(m,1); zeros(n,1)]))); | |
110 %! x1 = x1(end - n + 1 : end); | |
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
111 %! assert (x1, x0, 1e-6); |