Mercurial > octave-nkf
annotate scripts/linear-algebra/krylov.m @ 20614:10ec79b47808
use new string_value method to handle value extraction errors
* __voronoi__.cc, chol.cc, colamd.cc, fftw.cc:
Use new string_value method.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Thu, 08 Oct 2015 18:15:56 -0400 |
| parents | 83792dd9bcc1 |
| children |
| rev | line source |
|---|---|
|
19731
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
19074
diff
changeset
|
1 ## Copyright (C) 1993-2015 Auburn University. All rights reserved. |
| 3427 | 2 ## |
| 3 ## This file is part of Octave. | |
| 4 ## | |
| 5 ## Octave is free software; you can redistribute it and/or modify it | |
| 7016 | 6 ## under the terms of the GNU General Public License as published by |
| 7 ## the Free Software Foundation; either version 3 of the License, or (at | |
| 8 ## your option) any later version. | |
| 3427 | 9 ## |
| 7016 | 10 ## Octave is distributed in the hope that it will be useful, but |
| 11 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
| 12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
| 13 ## General Public License for more details. | |
| 3427 | 14 ## |
| 15 ## You should have received a copy of the GNU General Public License | |
| 7016 | 16 ## along with Octave; see the file COPYING. If not, see |
| 17 ## <http://www.gnu.org/licenses/>. | |
| 3427 | 18 |
| 3439 | 19 ## -*- texinfo -*- |
|
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
|
20 ## @deftypefn {Function File} {[@var{u}, @var{h}, @var{nu}] =} krylov (@var{A}, @var{V}, @var{k}, @var{eps1}, @var{pflg}) |
| 5555 | 21 ## Construct an orthogonal basis @var{u} of block Krylov subspace |
| 22 ## | |
| 23 ## @example | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
24 ## [v a*v a^2*v @dots{} a^(k+1)*v] |
| 5555 | 25 ## @end example |
| 26 ## | |
| 27 ## @noindent | |
|
20195
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19731
diff
changeset
|
28 ## using Householder reflections to guard against loss of orthogonality. |
| 3427 | 29 ## |
|
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
|
30 ## If @var{V} is a vector, then @var{h} contains the Hessenberg matrix |
|
20195
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19731
diff
changeset
|
31 ## 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:
19731
diff
changeset
|
32 ## @code{rk = a*u(:,k)-u*h(:,k)}, and @nospell{@tcode{ek'}} is the vector |
| 7248 | 33 ## @code{[0, 0, @dots{}, 1]} of length @code{k}. Otherwise, @var{h} is |
| 34 ## meaningless. | |
| 35 ## | |
|
20195
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19731
diff
changeset
|
36 ## 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:
19731
diff
changeset
|
37 ## then @var{h} contains the Hessenberg matrix such that @code{a*u == u*h}. |
| 5555 | 38 ## |
|
20195
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19731
diff
changeset
|
39 ## 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:
19731
diff
changeset
|
40 ## (based on @var{eps1}). |
| 5555 | 41 ## |
|
20195
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19731
diff
changeset
|
42 ## 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:
19731
diff
changeset
|
43 ## contains the Hessenberg decomposition of @var{A}. |
| 5555 | 44 ## |
|
20195
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19731
diff
changeset
|
45 ## 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:
19731
diff
changeset
|
46 ## value is 1e-12. |
| 5555 | 47 ## |
|
20195
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19731
diff
changeset
|
48 ## 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:
19731
diff
changeset
|
49 ## improve numerical behavior. The default value is 0. |
| 3427 | 50 ## |
|
19074
0850b5212619
doc: Add @nospell macro around proper names in documentation.
Rik <rik@octave.org>
parents:
18891
diff
changeset
|
51 ## 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:
18891
diff
changeset
|
52 ## 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:
18891
diff
changeset
|
53 ## the 42nd IEEE Conference on Decision and Control, December 2003. |
| 3439 | 54 ## @end deftypefn |
| 55 | |
| 56 ## Author: A. Scottedward Hodel <a.s.hodel@eng.auburn.edu> | |
| 3240 | 57 |
| 4609 | 58 function [Uret, H, nu] = krylov (A, V, k, eps1, pflg); |
| 3240 | 59 |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
60 if (isa (A, "single") || isa (V, "single")) |
|
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
61 defeps = 1e-6; |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
62 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
63 defeps = 1e-12; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
64 endif |
| 3457 | 65 |
| 66 if (nargin < 3 || nargin > 5) | |
| 6046 | 67 print_usage (); |
| 3457 | 68 elseif (nargin < 5) |
| 8506 | 69 ## Default permutation flag. |
| 70 pflg = 0; | |
| 3240 | 71 endif |
| 3457 | 72 |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
73 if (nargin < 4) |
| 8506 | 74 ## Default tolerance parameter. |
| 75 eps1 = defeps; | |
| 3240 | 76 endif |
| 3273 | 77 |
| 3457 | 78 if (isempty (eps1)) |
| 79 eps1 = defeps; | |
| 80 endif | |
| 81 | |
|
9872
72d6e0de76c7
fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents:
9843
diff
changeset
|
82 if (! issquare (A) || isempty (A)) |
|
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
83 error ("krylov: A(%d x %d) must be a non-empty square matrix", rows (A), columns (A)); |
| 3240 | 84 endif |
|
9872
72d6e0de76c7
fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents:
9843
diff
changeset
|
85 na = rows (A); |
| 3225 | 86 |
| 7201 | 87 [m, kb] = size (V); |
| 3457 | 88 if (m != na) |
|
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
89 error ("krylov: A(%d x %d), V(%d x %d): argument dimensions do not match", |
|
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
90 na, na, m, kb); |
| 3240 | 91 endif |
| 3233 | 92 |
| 4030 | 93 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
|
94 error ("krylov: K must be a scalar integer"); |
| 3457 | 95 endif |
| 96 | |
| 97 Vnrm = norm (V, Inf); | |
| 3273 | 98 |
| 8506 | 99 ## check for trivial solution. |
| 3457 | 100 if (Vnrm == 0) |
| 101 Uret = []; | |
| 4609 | 102 H = []; |
| 3457 | 103 nu = 0; |
| 104 return; | |
| 3211 | 105 endif |
| 106 | |
| 8506 | 107 ## Identify trivial null space. |
| 3457 | 108 abm = max (abs ([A, V]')); |
| 109 zidx = find (abm == 0); | |
| 3240 | 110 |
| 8506 | 111 ## Set up vector of pivot points. |
| 3240 | 112 pivot_vec = 1:na; |
| 3225 | 113 |
| 3240 | 114 iter = 0; |
| 3273 | 115 alpha = []; |
| 3240 | 116 nh = 0; |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
117 while (length (alpha) < na) && (columns (V) > 0) && (iter < k) |
| 3240 | 118 iter++; |
| 3211 | 119 |
| 8506 | 120 ## Get orthogonal basis of V. |
| 3240 | 121 jj = 1; |
| 3457 | 122 while (jj <= columns (V) && length (alpha) < na) |
| 8506 | 123 ## 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
|
124 nu = length (alpha)+1; |
| 3273 | 125 |
| 3240 | 126 short_pv = pivot_vec(nu:na); |
| 127 q = V(:,jj); | |
| 128 short_q = q(short_pv); | |
| 3211 | 129 |
| 3457 | 130 if (norm (short_q) < eps1) |
| 10549 | 131 ## Insignificant column; delete. |
| 3457 | 132 nv = columns (V); |
| 133 if (jj != nv) | |
| 134 [V(:,jj), V(:,nv)] = swap (V(:,jj), V(:,nv)); | |
|
18891
7bbe3658c5ef
maint: Use "FIXME:" coding convention in m-files.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
135 ## FIXME: H columns should be swapped too. |
|
7bbe3658c5ef
maint: Use "FIXME:" coding convention in m-files.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
136 ## Not done since Block Hessenberg structure is lost anyway. |
| 3240 | 137 endif |
| 138 V = V(:,1:(nv-1)); | |
| 10549 | 139 ## One less reflection. |
| 3457 | 140 nu--; |
| 3240 | 141 else |
| 10549 | 142 ## New householder reflection. |
| 3457 | 143 if (pflg) |
| 8506 | 144 ## Locate max magnitude element in short_q. |
| 3457 | 145 asq = abs (short_q); |
| 146 maxv = max (asq); | |
| 6024 | 147 maxidx = find (asq == maxv, 1); |
| 148 pivot_idx = short_pv(maxidx); | |
| 3273 | 149 |
| 10549 | 150 ## See if need to change the pivot list. |
| 3457 | 151 if (pivot_idx != pivot_vec(nu)) |
| 6024 | 152 swapidx = maxidx + (nu-1); |
| 3457 | 153 [pivot_vec(nu), pivot_vec(swapidx)] = ... |
| 10549 | 154 swap (pivot_vec(nu), pivot_vec(swapidx)); |
| 3240 | 155 endif |
| 156 endif | |
| 3273 | 157 |
| 10549 | 158 ## Isolate portion of vector for reflection. |
| 3240 | 159 idx = pivot_vec(nu:na); |
| 160 jdx = pivot_vec(1:nu); | |
| 161 | |
| 3457 | 162 [hv, av, z] = housh (q(idx), 1, 0); |
| 3240 | 163 alpha(nu) = av; |
| 164 U(idx,nu) = hv; | |
| 3211 | 165 |
| 8506 | 166 ## Reduce V per the reflection. |
| 3240 | 167 V(idx,:) = V(idx,:) - av*hv*(hv' * V(idx,:)); |
| 168 if(iter > 1) | |
|
18891
7bbe3658c5ef
maint: Use "FIXME:" coding convention in m-files.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
169 ## FIXME: not done correctly for block case. |
| 3240 | 170 H(nu,nu-1) = V(pivot_vec(nu),jj); |
| 171 endif | |
| 172 | |
| 8506 | 173 ## Advance to next column of V. |
| 3457 | 174 jj++; |
| 3240 | 175 endif |
| 176 endwhile | |
| 3211 | 177 |
| 8506 | 178 ## Check for oversize V (due to full rank). |
| 3457 | 179 if ((columns (V) > na) && (length (alpha) == na)) |
| 8506 | 180 ## Trim to size. |
| 3273 | 181 V = V(:,1:na); |
| 3457 | 182 elseif (columns(V) > na) |
|
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
183 krylov_V = V; |
|
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
184 krylov_na = na; |
|
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
185 krylov_length_alpha = length (alpha); |
|
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
186 error ("krylov: this case should never happen; submit a bug report"); |
| 3273 | 187 endif |
| 188 | |
| 3457 | 189 if (columns (V) > 0) |
| 8506 | 190 ## Construct next Q and multiply. |
| 3457 | 191 Q = zeros (size (V)); |
| 192 for kk = 1:columns (Q) | |
| 3240 | 193 Q(pivot_vec(nu-columns(Q)+kk),kk) = 1; |
| 194 endfor | |
| 3273 | 195 |
| 8506 | 196 ## Apply Householder reflections. |
| 3240 | 197 for ii = nu:-1:1 |
| 198 idx = pivot_vec(ii:na); | |
| 199 hv = U(idx,ii); | |
| 200 av = alpha(ii); | |
| 201 Q(idx,:) = Q(idx,:) - av*hv*(hv'*Q(idx,:)); | |
| 202 endfor | |
| 3211 | 203 endif |
| 204 | |
| 8506 | 205 ## Multiply to get new vector. |
| 3240 | 206 V = A*Q; |
| 8506 | 207 ## Project off of previous vectors. |
| 3457 | 208 nu = length (alpha); |
| 209 for i = 1:nu | |
| 210 hv = U(:,i); | |
| 211 av = alpha(i); | |
|
20266
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20195
diff
changeset
|
212 V -= av*hv*(hv'*V); |
| 3240 | 213 H(i,nu-columns(V)+(1:columns(V))) = V(pivot_vec(i),:); |
| 7151 | 214 endfor |
| 3240 | 215 |
| 3211 | 216 endwhile |
| 217 | |
| 8506 | 218 ## Back out complete U matrix. |
| 219 ## back out U matrix. | |
| 3457 | 220 j1 = columns (U); |
| 221 for i = j1:-1:1; | |
| 3240 | 222 idx = pivot_vec(i:na); |
| 223 hv = U(idx,i); | |
| 224 av = alpha(i); | |
| 3457 | 225 U(:,i) = zeros (na, 1); |
| 3240 | 226 U(idx(1),i) = 1; |
| 227 U(idx,i:j1) = U(idx,i:j1)-av*hv*(hv'*U(idx,i:j1)); | |
| 3211 | 228 endfor |
| 3273 | 229 |
| 3457 | 230 nu = length (alpha); |
| 3240 | 231 Uret = U; |
| 3457 | 232 if (max (max (abs (Uret(zidx,:)))) > 0) |
| 233 warning ("krylov: trivial null space corrupted; set pflg = 1 or eps1 > %e", | |
| 10549 | 234 eps1); |
| 3225 | 235 endif |
| 236 | |
| 3211 | 237 endfunction |
| 9843 | 238 |
| 239 | |
| 240 function [a1, b1] = swap (a, b) | |
| 241 | |
| 242 a1 = b; | |
| 243 b1 = a; | |
| 244 | |
| 245 endfunction | |
|
17338
1c89599167a6
maint: End m-files with 1 blank line.
Rik <rik@octave.org>
parents:
17268
diff
changeset
|
246 |