Mercurial > octave-nkf
annotate scripts/linear-algebra/krylov.m @ 12312:b10ea6efdc58 release-3-4-x ss-3-3-91
version is now 3.3.91
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Mon, 31 Jan 2011 08:36:58 -0500 |
| parents | b0084095098e |
| children | 990762e784fe |
| rev | line source |
|---|---|
| 11523 | 1 ## Copyright (C) 1993-2011 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 | |
| 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 |
| 7248 | 31 ## such that @code{a*u == u*h+rk*ek'}, in which @code{rk = |
| 32 ## a*u(:,k)-u*h(:,k)}, and @code{ek'} is the vector | |
| 33 ## @code{[0, 0, @dots{}, 1]} of length @code{k}. Otherwise, @var{h} is | |
| 34 ## meaningless. | |
| 35 ## | |
|
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
|
36 ## If @var{V} is a vector and @var{k} is greater than |
| 7248 | 37 ## @code{length(A)-1}, then @var{h} contains the Hessenberg matrix such |
| 38 ## that @code{a*u == u*h}. | |
| 5555 | 39 ## |
|
10791
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
40 ## The value of @var{nu} is the dimension of the span of the Krylov |
| 5555 | 41 ## subspace (based on @var{eps1}). |
| 42 ## | |
| 43 ## If @var{b} is a vector and @var{k} is greater than @var{m-1}, then | |
|
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
|
44 ## @var{h} contains the Hessenberg decomposition of @var{A}. |
| 5555 | 45 ## |
| 46 ## The optional parameter @var{eps1} is the threshold for zero. The | |
| 47 ## default value is 1e-12. | |
| 48 ## | |
| 49 ## If the optional parameter @var{pflg} is nonzero, row pivoting is used | |
| 50 ## to improve numerical behavior. The default value is 0. | |
| 3427 | 51 ## |
|
10791
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
52 ## Reference: A. Hodel, P. Misra, @cite{Partial Pivoting in the Computation of |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
53 ## Krylov Subspaces of Large Sparse Systems}, Proceedings of the 42nd IEEE |
|
10791
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
54 ## Conference on Decision and Control, December 2003. |
| 3439 | 55 ## @end deftypefn |
| 56 | |
| 57 ## Author: A. Scottedward Hodel <a.s.hodel@eng.auburn.edu> | |
| 3240 | 58 |
| 4609 | 59 function [Uret, H, nu] = krylov (A, V, k, eps1, pflg); |
| 3240 | 60 |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
61 if (isa (A, "single") || isa (V, "single")) |
|
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
62 defeps = 1e-6; |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
63 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
64 defeps = 1e-12; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7248
diff
changeset
|
65 endif |
| 3457 | 66 |
| 67 if (nargin < 3 || nargin > 5) | |
| 6046 | 68 print_usage (); |
| 3457 | 69 elseif (nargin < 5) |
| 8506 | 70 ## Default permutation flag. |
| 71 pflg = 0; | |
| 3240 | 72 endif |
| 3457 | 73 |
| 3240 | 74 if(nargin < 4) |
| 8506 | 75 ## Default tolerance parameter. |
| 76 eps1 = defeps; | |
| 3240 | 77 endif |
| 3273 | 78 |
| 3457 | 79 if (isempty (eps1)) |
| 80 eps1 = defeps; | |
| 81 endif | |
| 82 | |
|
9872
72d6e0de76c7
fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents:
9843
diff
changeset
|
83 if (! issquare (A) || isempty (A)) |
|
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
84 error ("krylov: A(%d x %d) must be a non-empty square matrix", rows (A), columns (A)); |
| 3240 | 85 endif |
|
9872
72d6e0de76c7
fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents:
9843
diff
changeset
|
86 na = rows (A); |
| 3225 | 87 |
| 7201 | 88 [m, kb] = size (V); |
| 3457 | 89 if (m != na) |
|
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
90 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
|
91 na, na, m, kb); |
| 3240 | 92 endif |
| 3233 | 93 |
| 4030 | 94 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
|
95 error ("krylov: K must be a scalar integer"); |
| 3457 | 96 endif |
| 97 | |
| 98 Vnrm = norm (V, Inf); | |
| 3273 | 99 |
| 8506 | 100 ## check for trivial solution. |
| 3457 | 101 if (Vnrm == 0) |
| 102 Uret = []; | |
| 4609 | 103 H = []; |
| 3457 | 104 nu = 0; |
| 105 return; | |
| 3211 | 106 endif |
| 107 | |
| 8506 | 108 ## Identify trivial null space. |
| 3457 | 109 abm = max (abs ([A, V]')); |
| 110 zidx = find (abm == 0); | |
| 3240 | 111 |
| 8506 | 112 ## Set up vector of pivot points. |
| 3240 | 113 pivot_vec = 1:na; |
| 3225 | 114 |
| 3240 | 115 iter = 0; |
| 3273 | 116 alpha = []; |
| 3240 | 117 nh = 0; |
| 3457 | 118 while (length(alpha) < na) && (columns(V) > 0) && (iter < k) |
| 3240 | 119 iter++; |
| 3211 | 120 |
| 8506 | 121 ## Get orthogonal basis of V. |
| 3240 | 122 jj = 1; |
| 3457 | 123 while (jj <= columns (V) && length (alpha) < na) |
| 8506 | 124 ## Index of next Householder reflection. |
| 3457 | 125 nu = length(alpha)+1; |
| 3273 | 126 |
| 3240 | 127 short_pv = pivot_vec(nu:na); |
| 128 q = V(:,jj); | |
| 129 short_q = q(short_pv); | |
| 3211 | 130 |
| 3457 | 131 if (norm (short_q) < eps1) |
| 10549 | 132 ## Insignificant column; delete. |
| 3457 | 133 nv = columns (V); |
| 134 if (jj != nv) | |
| 135 [V(:,jj), V(:,nv)] = swap (V(:,jj), V(:,nv)); | |
| 10549 | 136 ## FIXME -- H columns should be swapped too. Not done |
| 137 ## since Block Hessenberg structure is lost anyway. | |
| 3240 | 138 endif |
| 139 V = V(:,1:(nv-1)); | |
| 10549 | 140 ## One less reflection. |
| 3457 | 141 nu--; |
| 3240 | 142 else |
| 10549 | 143 ## New householder reflection. |
| 3457 | 144 if (pflg) |
| 8506 | 145 ## Locate max magnitude element in short_q. |
| 3457 | 146 asq = abs (short_q); |
| 147 maxv = max (asq); | |
| 6024 | 148 maxidx = find (asq == maxv, 1); |
| 149 pivot_idx = short_pv(maxidx); | |
| 3273 | 150 |
| 10549 | 151 ## See if need to change the pivot list. |
| 3457 | 152 if (pivot_idx != pivot_vec(nu)) |
| 6024 | 153 swapidx = maxidx + (nu-1); |
| 3457 | 154 [pivot_vec(nu), pivot_vec(swapidx)] = ... |
| 10549 | 155 swap (pivot_vec(nu), pivot_vec(swapidx)); |
| 3240 | 156 endif |
| 157 endif | |
| 3273 | 158 |
| 10549 | 159 ## Isolate portion of vector for reflection. |
| 3240 | 160 idx = pivot_vec(nu:na); |
| 161 jdx = pivot_vec(1:nu); | |
| 162 | |
| 3457 | 163 [hv, av, z] = housh (q(idx), 1, 0); |
| 3240 | 164 alpha(nu) = av; |
| 165 U(idx,nu) = hv; | |
| 3211 | 166 |
| 8506 | 167 ## Reduce V per the reflection. |
| 3240 | 168 V(idx,:) = V(idx,:) - av*hv*(hv' * V(idx,:)); |
| 169 if(iter > 1) | |
| 10549 | 170 ## FIXME -- not done correctly for block case. |
| 3240 | 171 H(nu,nu-1) = V(pivot_vec(nu),jj); |
| 172 endif | |
| 173 | |
| 8506 | 174 ## Advance to next column of V. |
| 3457 | 175 jj++; |
| 3240 | 176 endif |
| 177 endwhile | |
| 3211 | 178 |
| 8506 | 179 ## Check for oversize V (due to full rank). |
| 3457 | 180 if ((columns (V) > na) && (length (alpha) == na)) |
| 8506 | 181 ## Trim to size. |
| 3273 | 182 V = V(:,1:na); |
| 3457 | 183 elseif (columns(V) > na) |
|
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
184 krylov_V = V; |
|
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
185 krylov_na = na; |
|
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
186 krylov_length_alpha = length (alpha); |
|
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
187 error ("krylov: this case should never happen; submit a bug report"); |
| 3273 | 188 endif |
| 189 | |
| 3457 | 190 if (columns (V) > 0) |
| 8506 | 191 ## Construct next Q and multiply. |
| 3457 | 192 Q = zeros (size (V)); |
| 193 for kk = 1:columns (Q) | |
| 3240 | 194 Q(pivot_vec(nu-columns(Q)+kk),kk) = 1; |
| 195 endfor | |
| 3273 | 196 |
| 8506 | 197 ## Apply Householder reflections. |
| 3240 | 198 for ii = nu:-1:1 |
| 199 idx = pivot_vec(ii:na); | |
| 200 hv = U(idx,ii); | |
| 201 av = alpha(ii); | |
| 202 Q(idx,:) = Q(idx,:) - av*hv*(hv'*Q(idx,:)); | |
| 203 endfor | |
| 3211 | 204 endif |
| 205 | |
| 8506 | 206 ## Multiply to get new vector. |
| 3240 | 207 V = A*Q; |
| 8506 | 208 ## Project off of previous vectors. |
| 3457 | 209 nu = length (alpha); |
| 210 for i = 1:nu | |
| 211 hv = U(:,i); | |
| 212 av = alpha(i); | |
| 3240 | 213 V = V - av*hv*(hv'*V); |
| 214 H(i,nu-columns(V)+(1:columns(V))) = V(pivot_vec(i),:); | |
| 7151 | 215 endfor |
| 3240 | 216 |
| 3211 | 217 endwhile |
| 218 | |
| 8506 | 219 ## Back out complete U matrix. |
| 220 ## back out U matrix. | |
| 3457 | 221 j1 = columns (U); |
| 222 for i = j1:-1:1; | |
| 3240 | 223 idx = pivot_vec(i:na); |
| 224 hv = U(idx,i); | |
| 225 av = alpha(i); | |
| 3457 | 226 U(:,i) = zeros (na, 1); |
| 3240 | 227 U(idx(1),i) = 1; |
| 228 U(idx,i:j1) = U(idx,i:j1)-av*hv*(hv'*U(idx,i:j1)); | |
| 3211 | 229 endfor |
| 3273 | 230 |
| 3457 | 231 nu = length (alpha); |
| 3240 | 232 Uret = U; |
| 3457 | 233 if (max (max (abs (Uret(zidx,:)))) > 0) |
| 234 warning ("krylov: trivial null space corrupted; set pflg = 1 or eps1 > %e", | |
| 10549 | 235 eps1); |
| 3225 | 236 endif |
| 237 | |
| 3211 | 238 endfunction |
| 9843 | 239 |
| 240 | |
| 241 function [a1, b1] = swap (a, b) | |
| 242 | |
| 243 a1 = b; | |
| 244 b1 = a; | |
| 245 | |
| 246 endfunction |