Mercurial > hg > octave-nkf
annotate scripts/linear-algebra/krylov.m @ 7795:df9519e9990c
Handle single precision eps values
| author | David Bateman <dbateman@free.fr> |
|---|---|
| date | Mon, 12 May 2008 22:57:11 +0200 |
| parents | ffca9912dc82 |
| children | bc982528de11 |
| rev | line source |
|---|---|
| 7017 | 1 ## Copyright (C) 1993, 1998, 1999, 2000, 2002, 2003, 2005, 2006, 2007 |
| 2 ## 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 | |
| 25 ## [v a*v a^2*v ... a^(k+1)*v] | |
| 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) |
| 3240 | 71 pflg = 0; # default permutation flag |
| 72 endif | |
| 3457 | 73 |
| 3240 | 74 if(nargin < 4) |
| 75 eps1 = defeps; # default tolerance parameter | |
| 76 endif | |
| 3273 | 77 |
| 3457 | 78 if (isempty (eps1)) |
| 79 eps1 = defeps; | |
| 80 endif | |
| 81 | |
| 4030 | 82 na = issquare (A); |
| 3457 | 83 if (! na) |
| 84 error ("A(%d x %d) must be square", rows (A), columns (A)); | |
| 3240 | 85 endif |
| 3225 | 86 |
| 7201 | 87 [m, kb] = size (V); |
| 3457 | 88 if (m != na) |
| 89 error("A(%d x %d), V(%d x %d): argument dimensions do not match", | |
| 90 na, na, m, kb) | |
| 3240 | 91 endif |
| 3233 | 92 |
| 4030 | 93 if (! isscalar (k)) |
| 3457 | 94 error ("krylov: third argument must be a scalar integer"); |
| 95 endif | |
| 96 | |
| 97 Vnrm = norm (V, Inf); | |
| 3273 | 98 |
| 3457 | 99 ## check for trivial solution |
| 100 if (Vnrm == 0) | |
| 101 Uret = []; | |
| 4609 | 102 H = []; |
| 3457 | 103 nu = 0; |
| 104 return; | |
| 3211 | 105 endif |
| 106 | |
| 3240 | 107 # identify trivial null space |
| 3457 | 108 abm = max (abs ([A, V]')); |
| 109 zidx = find (abm == 0); | |
| 3240 | 110 |
| 111 # set up vector of pivot points | |
| 112 pivot_vec = 1:na; | |
| 3225 | 113 |
| 3240 | 114 iter = 0; |
| 3273 | 115 alpha = []; |
| 3240 | 116 nh = 0; |
| 3457 | 117 while (length(alpha) < na) && (columns(V) > 0) && (iter < k) |
| 3240 | 118 iter++; |
| 3211 | 119 |
| 3457 | 120 ## get orthogonal basis of V |
| 3240 | 121 jj = 1; |
| 3457 | 122 while (jj <= columns (V) && length (alpha) < na) |
| 123 ## index of next Householder reflection | |
| 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) |
| 131 ## insignificant column; delete | |
| 132 nv = columns (V); | |
| 133 if (jj != nv) | |
| 134 [V(:,jj), V(:,nv)] = swap (V(:,jj), V(:,nv)); | |
| 5775 | 135 ## FIXME -- H columns should be swapped too. Not done |
| 3457 | 136 ## since Block Hessenberg structure is lost anyway. |
| 3240 | 137 endif |
| 138 V = V(:,1:(nv-1)); | |
| 3457 | 139 ## one less reflection |
| 140 nu--; | |
| 3240 | 141 else |
| 3457 | 142 ## new householder reflection |
| 143 if (pflg) | |
| 144 ## locate max magnitude element in short_q | |
| 145 asq = abs (short_q); | |
| 146 maxv = max (asq); | |
| 6024 | 147 maxidx = find (asq == maxv, 1); |
| 148 pivot_idx = short_pv(maxidx); | |
| 3273 | 149 |
| 3457 | 150 ## see if need to change the pivot list |
| 151 if (pivot_idx != pivot_vec(nu)) | |
| 6024 | 152 swapidx = maxidx + (nu-1); |
| 3457 | 153 [pivot_vec(nu), pivot_vec(swapidx)] = ... |
| 154 swap (pivot_vec(nu), pivot_vec(swapidx)); | |
| 3240 | 155 endif |
| 156 endif | |
| 3273 | 157 |
| 3457 | 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 |
| 3240 | 166 # reduce V per the reflection |
| 167 V(idx,:) = V(idx,:) - av*hv*(hv' * V(idx,:)); | |
| 168 if(iter > 1) | |
| 5775 | 169 ## FIXME -- not done correctly for block case |
| 3240 | 170 H(nu,nu-1) = V(pivot_vec(nu),jj); |
| 171 endif | |
| 172 | |
| 3457 | 173 ## advance to next column of V |
| 174 jj++; | |
| 3240 | 175 endif |
| 176 endwhile | |
| 3211 | 177 |
| 3457 | 178 ## check for oversize V (due to full rank) |
| 179 if ((columns (V) > na) && (length (alpha) == na)) | |
| 180 ## trim to size | |
| 3273 | 181 V = V(:,1:na); |
| 3457 | 182 elseif (columns(V) > na) |
| 3273 | 183 krylov_V = V |
| 184 krylov_na = na | |
| 3457 | 185 krylov_length_alpha = length (alpha) |
| 186 error ("This case should never happen; submit a bug report"); | |
| 3273 | 187 endif |
| 188 | |
| 3457 | 189 if (columns (V) > 0) |
| 190 ## construct next Q and multiply | |
| 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 |
| 3457 | 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 | |
| 3457 | 205 ## multiply to get new vector; |
| 3240 | 206 V = A*Q; |
| 3457 | 207 ## project off of previous vectors |
| 208 nu = length (alpha); | |
| 209 for i = 1:nu | |
| 210 hv = U(:,i); | |
| 211 av = alpha(i); | |
| 3240 | 212 V = V - av*hv*(hv'*V); |
| 213 H(i,nu-columns(V)+(1:columns(V))) = V(pivot_vec(i),:); | |
| 7151 | 214 endfor |
| 3240 | 215 |
| 3211 | 216 endwhile |
| 217 | |
| 3457 | 218 ## Back out complete U matrix |
| 219 ## back out U matrix ; | |
| 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", | |
| 234 eps1); | |
| 3225 | 235 endif |
| 236 | |
| 3211 | 237 endfunction |