--- /dev/null
+## Copyright (C) 1993-2012 Auburn University. All rights reserved.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{u}, @var{h}, @var{nu}] =} krylov (@var{A}, @var{V}, @var{k}, @var{eps1}, @var{pflg})
+## Construct an orthogonal basis @var{u} of block Krylov subspace
+##
+## @example
+## [v a*v a^2*v @dots{} a^(k+1)*v]
+## @end example
+##
+## @noindent
+## Using Householder reflections to guard against loss of orthogonality.
+##
+## If @var{V} is a vector, then @var{h} contains the Hessenberg matrix
+## such that @nospell{@xcode{a*u == u*h+rk*ek'}}, in which @code{rk =
+## a*u(:,k)-u*h(:,k)}, and @nospell{@xcode{ek'}} is the vector
+## @code{[0, 0, @dots{}, 1]} of length @code{k}. Otherwise, @var{h} is
+## meaningless.
+##
+## If @var{V} is a vector and @var{k} is greater than
+## @code{length(A)-1}, then @var{h} contains the Hessenberg matrix such
+## that @code{a*u == u*h}.
+##
+## The value of @var{nu} is the dimension of the span of the Krylov
+## subspace (based on @var{eps1}).
+##
+## If @var{b} is a vector and @var{k} is greater than @var{m-1}, then
+## @var{h} contains the Hessenberg decomposition of @var{A}.
+##
+## The optional parameter @var{eps1} is the threshold for zero. The
+## default value is 1e-12.
+##
+## If the optional parameter @var{pflg} is nonzero, row pivoting is used
+## to improve numerical behavior. The default value is 0.
+##
+## Reference: A. Hodel, P. Misra, @cite{Partial Pivoting in the Computation of
+## Krylov Subspaces of Large Sparse Systems}, Proceedings of the 42nd IEEE
+## Conference on Decision and Control, December 2003.
+## @end deftypefn
+
+## Author: A. Scottedward Hodel <a.s.hodel@eng.auburn.edu>
+
+function [Uret, H, nu] = krylov (A, V, k, eps1, pflg);
+
+ if (isa (A, "single") || isa (V, "single"))
+ defeps = 1e-6;
+ else
+ defeps = 1e-12;
+ endif
+
+ if (nargin < 3 || nargin > 5)
+ print_usage ();
+ elseif (nargin < 5)
+ ## Default permutation flag.
+ pflg = 0;
+ endif
+
+ if(nargin < 4)
+ ## Default tolerance parameter.
+ eps1 = defeps;
+ endif
+
+ if (isempty (eps1))
+ eps1 = defeps;
+ endif
+
+ if (! issquare (A) || isempty (A))
+ error ("krylov: A(%d x %d) must be a non-empty square matrix", rows (A), columns (A));
+ endif
+ na = rows (A);
+
+ [m, kb] = size (V);
+ if (m != na)
+ error ("krylov: A(%d x %d), V(%d x %d): argument dimensions do not match",
+ na, na, m, kb);
+ endif
+
+ if (! isscalar (k))
+ error ("krylov: K must be a scalar integer");
+ endif
+
+ Vnrm = norm (V, Inf);
+
+ ## check for trivial solution.
+ if (Vnrm == 0)
+ Uret = [];
+ H = [];
+ nu = 0;
+ return;
+ endif
+
+ ## Identify trivial null space.
+ abm = max (abs ([A, V]'));
+ zidx = find (abm == 0);
+
+ ## Set up vector of pivot points.
+ pivot_vec = 1:na;
+
+ iter = 0;
+ alpha = [];
+ nh = 0;
+ while (length(alpha) < na) && (columns(V) > 0) && (iter < k)
+ iter++;
+
+ ## Get orthogonal basis of V.
+ jj = 1;
+ while (jj <= columns (V) && length (alpha) < na)
+ ## Index of next Householder reflection.
+ nu = length(alpha)+1;
+
+ short_pv = pivot_vec(nu:na);
+ q = V(:,jj);
+ short_q = q(short_pv);
+
+ if (norm (short_q) < eps1)
+ ## Insignificant column; delete.
+ nv = columns (V);
+ if (jj != nv)
+ [V(:,jj), V(:,nv)] = swap (V(:,jj), V(:,nv));
+ ## FIXME -- H columns should be swapped too. Not done
+ ## since Block Hessenberg structure is lost anyway.
+ endif
+ V = V(:,1:(nv-1));
+ ## One less reflection.
+ nu--;
+ else
+ ## New householder reflection.
+ if (pflg)
+ ## Locate max magnitude element in short_q.
+ asq = abs (short_q);
+ maxv = max (asq);
+ maxidx = find (asq == maxv, 1);
+ pivot_idx = short_pv(maxidx);
+
+ ## See if need to change the pivot list.
+ if (pivot_idx != pivot_vec(nu))
+ swapidx = maxidx + (nu-1);
+ [pivot_vec(nu), pivot_vec(swapidx)] = ...
+ swap (pivot_vec(nu), pivot_vec(swapidx));
+ endif
+ endif
+
+ ## Isolate portion of vector for reflection.
+ idx = pivot_vec(nu:na);
+ jdx = pivot_vec(1:nu);
+
+ [hv, av, z] = housh (q(idx), 1, 0);
+ alpha(nu) = av;
+ U(idx,nu) = hv;
+
+ ## Reduce V per the reflection.
+ V(idx,:) = V(idx,:) - av*hv*(hv' * V(idx,:));
+ if(iter > 1)
+ ## FIXME -- not done correctly for block case.
+ H(nu,nu-1) = V(pivot_vec(nu),jj);
+ endif
+
+ ## Advance to next column of V.
+ jj++;
+ endif
+ endwhile
+
+ ## Check for oversize V (due to full rank).
+ if ((columns (V) > na) && (length (alpha) == na))
+ ## Trim to size.
+ V = V(:,1:na);
+ elseif (columns(V) > na)
+ krylov_V = V;
+ krylov_na = na;
+ krylov_length_alpha = length (alpha);
+ error ("krylov: this case should never happen; submit a bug report");
+ endif
+
+ if (columns (V) > 0)
+ ## Construct next Q and multiply.
+ Q = zeros (size (V));
+ for kk = 1:columns (Q)
+ Q(pivot_vec(nu-columns(Q)+kk),kk) = 1;
+ endfor
+
+ ## Apply Householder reflections.
+ for ii = nu:-1:1
+ idx = pivot_vec(ii:na);
+ hv = U(idx,ii);
+ av = alpha(ii);
+ Q(idx,:) = Q(idx,:) - av*hv*(hv'*Q(idx,:));
+ endfor
+ endif
+
+ ## Multiply to get new vector.
+ V = A*Q;
+ ## Project off of previous vectors.
+ nu = length (alpha);
+ for i = 1:nu
+ hv = U(:,i);
+ av = alpha(i);
+ V = V - av*hv*(hv'*V);
+ H(i,nu-columns(V)+(1:columns(V))) = V(pivot_vec(i),:);
+ endfor
+
+ endwhile
+
+ ## Back out complete U matrix.
+ ## back out U matrix.
+ j1 = columns (U);
+ for i = j1:-1:1;
+ idx = pivot_vec(i:na);
+ hv = U(idx,i);
+ av = alpha(i);
+ U(:,i) = zeros (na, 1);
+ U(idx(1),i) = 1;
+ U(idx,i:j1) = U(idx,i:j1)-av*hv*(hv'*U(idx,i:j1));
+ endfor
+
+ nu = length (alpha);
+ Uret = U;
+ if (max (max (abs (Uret(zidx,:)))) > 0)
+ warning ("krylov: trivial null space corrupted; set pflg = 1 or eps1 > %e",
+ eps1);
+ endif
+
+endfunction
+
+
+function [a1, b1] = swap (a, b)
+
+ a1 = b;
+ b1 = a;
+
+endfunction
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff