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