X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Fm%2Flinear-algebra%2Fkrylov.m;fp=octave_packages%2Fm%2Flinear-algebra%2Fkrylov.m;h=b322f7d8c255251bd27c842d2848fe0f82381ef8;hp=0000000000000000000000000000000000000000;hb=1c0469ada9531828709108a4882a751d2816994a;hpb=63de9f36673d49121015e3695f2c336ea92bc278

diff --git a/octave_packages/m/linear-algebra/krylov.m b/octave_packages/m/linear-algebra/krylov.m
new file mode 100644
index 0000000..b322f7d
--- /dev/null
+++ b/octave_packages/m/linear-algebra/krylov.m
@@ -0,0 +1,246 @@
+## 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