Mercurial > octave
annotate scripts/linear-algebra/krylov.m @ 30564:796f54d4ddbf stable
update Octave Project Developers copyright for the new year
In files that have the "Octave Project Developers" copyright notice,
update for 2021.
In all .txi and .texi files except gpl.txi and gpl.texi in the
doc/liboctave and doc/interpreter directories, change the copyright
to "Octave Project Developers", the same as used for other source
files. Update copyright notices for 2022 (not done since 2019). For
gpl.txi and gpl.texi, change the copyright notice to be "Free Software
Foundation, Inc." and leave the date at 2007 only because this file
only contains the text of the GPL, not anything created by the Octave
Project Developers.
Add Paul Thomas to contributors.in.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Tue, 28 Dec 2021 18:22:40 -0500 |
| parents | 7854d5752dd2 |
| children | 597f3ee61a48 |
| rev | line source |
|---|---|
|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
1 ######################################################################## |
|
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
2 ## |
|
30564
796f54d4ddbf
update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
29359
diff
changeset
|
3 ## Copyright (C) 1993-2022 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
|
4 ## |
|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
5 ## See the file COPYRIGHT.md in the top-level directory of this |
|
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
6 ## distribution or <https://octave.org/copyright/>. |
| 3427 | 7 ## |
| 8 ## This file is part of Octave. | |
| 9 ## | |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
10 ## Octave is free software: you can redistribute it and/or modify it |
| 7016 | 11 ## 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
|
12 ## the Free Software Foundation, either version 3 of the License, or |
|
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22737
diff
changeset
|
13 ## (at your option) any later version. |
| 3427 | 14 ## |
| 7016 | 15 ## Octave is distributed in the hope that it will be useful, but |
| 16 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
|
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22737
diff
changeset
|
17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22737
diff
changeset
|
18 ## GNU General Public License for more details. |
| 3427 | 19 ## |
| 20 ## You should have received a copy of the GNU General Public License | |
| 7016 | 21 ## 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
|
22 ## <https://www.gnu.org/licenses/>. |
|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
23 ## |
|
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
24 ######################################################################## |
| 3427 | 25 |
| 3439 | 26 ## -*- texinfo -*- |
|
20852
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20735
diff
changeset
|
27 ## @deftypefn {} {[@var{u}, @var{h}, @var{nu}] =} krylov (@var{A}, @var{V}, @var{k}, @var{eps1}, @var{pflg}) |
|
25106
d7ad543255c5
doc: Shorten very long first sentences of docstrings (bug #53388).
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
28 ## Construct an orthogonal basis @var{u} of a block Krylov subspace. |
|
d7ad543255c5
doc: Shorten very long first sentences of docstrings (bug #53388).
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
29 ## |
|
d7ad543255c5
doc: Shorten very long first sentences of docstrings (bug #53388).
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
30 ## The block Krylov subspace has the following form: |
| 5555 | 31 ## |
| 32 ## @example | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
33 ## [v a*v a^2*v @dots{} a^(k+1)*v] |
| 5555 | 34 ## @end example |
| 35 ## | |
| 36 ## @noindent | |
|
25106
d7ad543255c5
doc: Shorten very long first sentences of docstrings (bug #53388).
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
37 ## The construction is made with Householder reflections to guard against loss |
|
d7ad543255c5
doc: Shorten very long first sentences of docstrings (bug #53388).
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
38 ## of orthogonality. |
| 3427 | 39 ## |
|
11470
eb9e0b597d61
Use common names for variables in documentation and code for a few more m-script files.
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
40 ## If @var{V} is a vector, then @var{h} contains the Hessenberg matrix |
|
20160
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
41 ## such that @nospell{@tcode{a*u == u*h+rk*ek'}}, in which |
|
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
42 ## @code{rk = a*u(:,k)-u*h(:,k)}, and @nospell{@tcode{ek'}} is the vector |
|
25472
9771111f04f4
doc: Use @var rather than @code to mark inputs to functions in docstrings.
Rik <rik@octave.org>
parents:
25106
diff
changeset
|
43 ## @code{[0, 0, @dots{}, 1]} of length @var{k}. Otherwise, @var{h} is |
| 7248 | 44 ## meaningless. |
| 45 ## | |
|
20160
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
46 ## If @var{V} is a vector and @var{k} is greater than @code{length (A) - 1}, |
|
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
47 ## then @var{h} contains the Hessenberg matrix such that @code{a*u == u*h}. |
| 5555 | 48 ## |
|
20160
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
49 ## The value of @var{nu} is the dimension of the span of the Krylov subspace |
|
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
50 ## (based on @var{eps1}). |
| 5555 | 51 ## |
|
20160
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
52 ## If @var{b} is a vector and @var{k} is greater than @var{m-1}, then @var{h} |
|
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
53 ## contains the Hessenberg decomposition of @var{A}. |
| 5555 | 54 ## |
|
20160
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
55 ## The optional parameter @var{eps1} is the threshold for zero. The default |
|
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
56 ## value is 1e-12. |
| 5555 | 57 ## |
|
20160
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
58 ## If the optional parameter @var{pflg} is nonzero, row pivoting is used to |
|
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
59 ## improve numerical behavior. The default value is 0. |
| 3427 | 60 ## |
|
19040
0850b5212619
doc: Add @nospell macro around proper names in documentation.
Rik <rik@octave.org>
parents:
18857
diff
changeset
|
61 ## Reference: @nospell{A. Hodel, P. Misra}, @cite{Partial Pivoting in the |
|
0850b5212619
doc: Add @nospell macro around proper names in documentation.
Rik <rik@octave.org>
parents:
18857
diff
changeset
|
62 ## Computation of Krylov Subspaces of Large Sparse Systems}, Proceedings of |
|
0850b5212619
doc: Add @nospell macro around proper names in documentation.
Rik <rik@octave.org>
parents:
18857
diff
changeset
|
63 ## the 42nd IEEE Conference on Decision and Control, December 2003. |
| 3439 | 64 ## @end deftypefn |
| 65 | |
|
22767
212333a97d8d
maint: Remove dangling ';' from m-file function declarations.
Rik <rik@octave.org>
parents:
22755
diff
changeset
|
66 function [Uret, H, nu] = krylov (A, V, k, eps1, pflg) |
| 3240 | 67 |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
68 if (isa (A, "single") || isa (V, "single")) |
|
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
69 defeps = 1e-6; |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
70 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
71 defeps = 1e-12; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
72 endif |
| 3457 | 73 |
|
28789
28de41192f3c
Eliminate unneeded verification of nargin, nargout in m-files.
Rik <rik@octave.org>
parents:
27978
diff
changeset
|
74 if (nargin < 3) |
| 6046 | 75 print_usage (); |
| 3457 | 76 elseif (nargin < 5) |
| 8506 | 77 ## Default permutation flag. |
| 78 pflg = 0; | |
| 3240 | 79 endif |
| 3457 | 80 |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
81 if (nargin < 4) |
| 8506 | 82 ## Default tolerance parameter. |
| 83 eps1 = defeps; | |
| 3240 | 84 endif |
| 3273 | 85 |
| 3457 | 86 if (isempty (eps1)) |
| 87 eps1 = defeps; | |
| 88 endif | |
| 89 | |
|
9872
72d6e0de76c7
fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents:
9843
diff
changeset
|
90 if (! issquare (A) || isempty (A)) |
|
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
91 error ("krylov: A(%d x %d) must be a non-empty square matrix", rows (A), columns (A)); |
| 3240 | 92 endif |
|
9872
72d6e0de76c7
fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents:
9843
diff
changeset
|
93 na = rows (A); |
| 3225 | 94 |
| 7201 | 95 [m, kb] = size (V); |
| 3457 | 96 if (m != na) |
|
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
97 error ("krylov: A(%d x %d), V(%d x %d): argument dimensions do not match", |
|
21031
66a08c3cafe3
maint: Follow Octave coding conventions in m-files.
Rik <rik@octave.org>
parents:
20852
diff
changeset
|
98 na, na, m, kb); |
| 3240 | 99 endif |
| 3233 | 100 |
| 4030 | 101 if (! isscalar (k)) |
|
11472
1740012184f9
Use uppercase for variable names in error() strings to match Info documentation. Only m-files done.
Rik <octave@nomad.inbox5.com>
parents:
11470
diff
changeset
|
102 error ("krylov: K must be a scalar integer"); |
| 3457 | 103 endif |
| 104 | |
| 105 Vnrm = norm (V, Inf); | |
| 3273 | 106 |
| 8506 | 107 ## check for trivial solution. |
| 3457 | 108 if (Vnrm == 0) |
| 109 Uret = []; | |
| 4609 | 110 H = []; |
| 3457 | 111 nu = 0; |
| 112 return; | |
| 3211 | 113 endif |
| 114 | |
| 8506 | 115 ## Identify trivial null space. |
| 3457 | 116 abm = max (abs ([A, V]')); |
| 117 zidx = find (abm == 0); | |
| 3240 | 118 |
| 8506 | 119 ## Set up vector of pivot points. |
| 3240 | 120 pivot_vec = 1:na; |
| 3225 | 121 |
| 3240 | 122 iter = 0; |
| 3273 | 123 alpha = []; |
| 3240 | 124 nh = 0; |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
125 while (length (alpha) < na) && (columns (V) > 0) && (iter < k) |
|
20735
418ae0cb752f
Replace ++,-- with in-place operators for performance.
Rik <rik@octave.org>
parents:
20231
diff
changeset
|
126 iter += 1; |
| 3211 | 127 |
| 8506 | 128 ## Get orthogonal basis of V. |
| 3240 | 129 jj = 1; |
| 3457 | 130 while (jj <= columns (V) && length (alpha) < na) |
| 8506 | 131 ## Index of next Householder reflection. |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
132 nu = length (alpha)+1; |
| 3273 | 133 |
| 3240 | 134 short_pv = pivot_vec(nu:na); |
| 135 q = V(:,jj); | |
| 136 short_q = q(short_pv); | |
| 3211 | 137 |
| 3457 | 138 if (norm (short_q) < eps1) |
| 10549 | 139 ## Insignificant column; delete. |
| 3457 | 140 nv = columns (V); |
| 141 if (jj != nv) | |
| 142 [V(:,jj), V(:,nv)] = swap (V(:,jj), V(:,nv)); | |
|
18857
7bbe3658c5ef
maint: Use "FIXME:" coding convention in m-files.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
143 ## FIXME: H columns should be swapped too. |
|
7bbe3658c5ef
maint: Use "FIXME:" coding convention in m-files.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
144 ## Not done since Block Hessenberg structure is lost anyway. |
| 3240 | 145 endif |
| 146 V = V(:,1:(nv-1)); | |
| 10549 | 147 ## One less reflection. |
|
20735
418ae0cb752f
Replace ++,-- with in-place operators for performance.
Rik <rik@octave.org>
parents:
20231
diff
changeset
|
148 nu -= 1; |
| 3240 | 149 else |
| 10549 | 150 ## New householder reflection. |
| 3457 | 151 if (pflg) |
| 8506 | 152 ## Locate max magnitude element in short_q. |
| 3457 | 153 asq = abs (short_q); |
| 154 maxv = max (asq); | |
| 6024 | 155 maxidx = find (asq == maxv, 1); |
| 156 pivot_idx = short_pv(maxidx); | |
| 3273 | 157 |
| 10549 | 158 ## See if need to change the pivot list. |
| 3457 | 159 if (pivot_idx != pivot_vec(nu)) |
| 6024 | 160 swapidx = maxidx + (nu-1); |
| 3457 | 161 [pivot_vec(nu), pivot_vec(swapidx)] = ... |
| 10549 | 162 swap (pivot_vec(nu), pivot_vec(swapidx)); |
| 3240 | 163 endif |
| 164 endif | |
| 3273 | 165 |
| 10549 | 166 ## Isolate portion of vector for reflection. |
| 3240 | 167 idx = pivot_vec(nu:na); |
| 168 jdx = pivot_vec(1:nu); | |
| 169 | |
| 3457 | 170 [hv, av, z] = housh (q(idx), 1, 0); |
| 3240 | 171 alpha(nu) = av; |
| 172 U(idx,nu) = hv; | |
| 3211 | 173 |
| 8506 | 174 ## Reduce V per the reflection. |
| 3240 | 175 V(idx,:) = V(idx,:) - av*hv*(hv' * V(idx,:)); |
|
22737
7abc25e6206a
maint: Clean up code base to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
176 if (iter > 1) |
|
18857
7bbe3658c5ef
maint: Use "FIXME:" coding convention in m-files.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
177 ## FIXME: not done correctly for block case. |
| 3240 | 178 H(nu,nu-1) = V(pivot_vec(nu),jj); |
| 179 endif | |
| 180 | |
| 8506 | 181 ## Advance to next column of V. |
|
20735
418ae0cb752f
Replace ++,-- with in-place operators for performance.
Rik <rik@octave.org>
parents:
20231
diff
changeset
|
182 jj += 1; |
| 3240 | 183 endif |
| 184 endwhile | |
| 3211 | 185 |
| 8506 | 186 ## Check for oversize V (due to full rank). |
| 3457 | 187 if ((columns (V) > na) && (length (alpha) == na)) |
| 8506 | 188 ## Trim to size. |
| 3273 | 189 V = V(:,1:na); |
|
28912
0de38a6ef693
maint: Use Octave convention of space after function name in scripts dir.
Rik <rik@octave.org>
parents:
28789
diff
changeset
|
190 elseif (columns (V) > na) |
|
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
191 krylov_V = V; |
|
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
192 krylov_na = na; |
|
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
193 krylov_length_alpha = length (alpha); |
|
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
194 error ("krylov: this case should never happen; submit a bug report"); |
| 3273 | 195 endif |
| 196 | |
| 3457 | 197 if (columns (V) > 0) |
| 8506 | 198 ## Construct next Q and multiply. |
| 3457 | 199 Q = zeros (size (V)); |
| 200 for kk = 1:columns (Q) | |
|
28912
0de38a6ef693
maint: Use Octave convention of space after function name in scripts dir.
Rik <rik@octave.org>
parents:
28789
diff
changeset
|
201 Q(pivot_vec(nu-columns (Q)+kk),kk) = 1; |
| 3240 | 202 endfor |
| 3273 | 203 |
| 8506 | 204 ## Apply Householder reflections. |
| 3240 | 205 for ii = nu:-1:1 |
| 206 idx = pivot_vec(ii:na); | |
| 207 hv = U(idx,ii); | |
| 208 av = alpha(ii); | |
| 209 Q(idx,:) = Q(idx,:) - av*hv*(hv'*Q(idx,:)); | |
| 210 endfor | |
| 3211 | 211 endif |
| 212 | |
| 8506 | 213 ## Multiply to get new vector. |
| 3240 | 214 V = A*Q; |
| 8506 | 215 ## Project off of previous vectors. |
| 3457 | 216 nu = length (alpha); |
| 217 for i = 1:nu | |
| 218 hv = U(:,i); | |
| 219 av = alpha(i); | |
|
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20160
diff
changeset
|
220 V -= av*hv*(hv'*V); |
|
28912
0de38a6ef693
maint: Use Octave convention of space after function name in scripts dir.
Rik <rik@octave.org>
parents:
28789
diff
changeset
|
221 H(i,nu-columns (V)+(1:columns (V))) = V(pivot_vec(i),:); |
| 7151 | 222 endfor |
| 3240 | 223 |
| 3211 | 224 endwhile |
| 225 | |
| 8506 | 226 ## Back out complete U matrix. |
| 227 ## back out U matrix. | |
| 3457 | 228 j1 = columns (U); |
|
26268
6dd232798997
maint: Remove useless ';' from end of for, if, while, etc. statements.
Rik <rik@octave.org>
parents:
25472
diff
changeset
|
229 for i = j1:-1:1 |
| 3240 | 230 idx = pivot_vec(i:na); |
| 231 hv = U(idx,i); | |
| 232 av = alpha(i); | |
| 3457 | 233 U(:,i) = zeros (na, 1); |
| 3240 | 234 U(idx(1),i) = 1; |
| 235 U(idx,i:j1) = U(idx,i:j1)-av*hv*(hv'*U(idx,i:j1)); | |
| 3211 | 236 endfor |
| 3273 | 237 |
| 3457 | 238 nu = length (alpha); |
| 3240 | 239 Uret = U; |
| 3457 | 240 if (max (max (abs (Uret(zidx,:)))) > 0) |
| 241 warning ("krylov: trivial null space corrupted; set pflg = 1 or eps1 > %e", | |
| 10549 | 242 eps1); |
| 3225 | 243 endif |
| 244 | |
| 3211 | 245 endfunction |
| 9843 | 246 |
| 247 | |
| 248 function [a1, b1] = swap (a, b) | |
| 249 | |
| 250 a1 = b; | |
| 251 b1 = a; | |
| 252 | |
| 253 endfunction |