octave_packages/m/linear-algebra/krylov.m

   1 ## Copyright (C) 1993-2012 Auburn University.  All rights reserved.
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   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.
   9 ##
  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.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {[@var{u}, @var{h}, @var{nu}] =} krylov (@var{A}, @var{V}, @var{k}, @var{eps1}, @var{pflg})
  21 ## Construct an orthogonal basis @var{u} of block Krylov subspace
  22 ##
  23 ## @example
  24 ## [v a*v a^2*v @dots{} a^(k+1)*v]
  25 ## @end example
  26 ##
  27 ## @noindent
  28 ## Using Householder reflections to guard against loss of orthogonality.
  29 ##
  30 ## If @var{V} is a vector, then @var{h} contains the Hessenberg matrix
  31 ## such that @nospell{@xcode{a*u == u*h+rk*ek'}}, in which @code{rk =
  32 ## a*u(:,k)-u*h(:,k)}, and @nospell{@xcode{ek'}} is the vector
  33 ## @code{[0, 0, @dots{}, 1]} of length @code{k}.  Otherwise, @var{h} is
  34 ## meaningless.
  35 ##
  36 ## If @var{V} is a vector and @var{k} is greater than
  37 ## @code{length(A)-1}, then @var{h} contains the Hessenberg matrix such
  38 ## that @code{a*u == u*h}.
  39 ##
  40 ## The value of @var{nu} is the dimension of the span of the Krylov
  41 ## subspace (based on @var{eps1}).
  42 ##
  43 ## If @var{b} is a vector and @var{k} is greater than @var{m-1}, then
  44 ## @var{h} contains the Hessenberg decomposition of @var{A}.
  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.
  51 ##
  52 ## Reference: A. Hodel, P. Misra, @cite{Partial Pivoting in the Computation of
  53 ## Krylov Subspaces of Large Sparse Systems}, Proceedings of the 42nd IEEE
  54 ## Conference on Decision and Control, December 2003.
  55 ## @end deftypefn
  56
  57 ## Author: A. Scottedward Hodel <a.s.hodel@eng.auburn.edu>
  58
  59 function [Uret, H, nu] = krylov (A, V, k, eps1, pflg);
  60
  61   if (isa (A, "single") || isa (V, "single"))
  62     defeps = 1e-6;
  63   else
  64     defeps = 1e-12;
  65   endif
  66
  67   if (nargin < 3 || nargin > 5)
  68     print_usage ();
  69   elseif (nargin < 5)
  70     ## Default permutation flag.
  71     pflg = 0;
  72   endif
  73
  74   if(nargin < 4)
  75     ## Default tolerance parameter.
  76     eps1 = defeps;
  77   endif
  78
  79   if (isempty (eps1))
  80     eps1 = defeps;
  81   endif
  82
  83   if (! issquare (A) || isempty (A))
  84     error ("krylov: A(%d x %d) must be a non-empty square matrix", rows (A), columns (A));
  85   endif
  86   na = rows (A);
  87
  88   [m, kb] = size (V);
  89   if (m != na)
  90     error ("krylov: A(%d x %d), V(%d x %d): argument dimensions do not match",
  91           na, na, m, kb);
  92   endif
  93
  94   if (! isscalar (k))
  95     error ("krylov: K must be a scalar integer");
  96   endif
  97
  98   Vnrm = norm (V, Inf);
  99
 100   ## check for trivial solution.
 101   if (Vnrm == 0)
 102     Uret = [];
 103     H = [];
 104     nu = 0;
 105     return;
 106   endif
 107
 108   ## Identify trivial null space.
 109   abm = max (abs ([A, V]'));
 110   zidx = find (abm == 0);
 111
 112   ## Set up vector of pivot points.
 113   pivot_vec = 1:na;
 114
 115   iter = 0;
 116   alpha = [];
 117   nh = 0;
 118   while (length(alpha) < na) && (columns(V) > 0) && (iter < k)
 119     iter++;
 120
 121     ## Get orthogonal basis of V.
 122     jj = 1;
 123     while (jj <= columns (V) && length (alpha) < na)
 124       ## Index of next Householder reflection.
 125       nu = length(alpha)+1;
 126
 127       short_pv = pivot_vec(nu:na);
 128       q = V(:,jj);
 129       short_q = q(short_pv);
 130
 131       if (norm (short_q) < eps1)
 132         ## Insignificant column; delete.
 133         nv = columns (V);
 134         if (jj != nv)
 135           [V(:,jj), V(:,nv)] = swap (V(:,jj), V(:,nv));
 136           ## FIXME -- H columns should be swapped too.  Not done
 137           ## since Block Hessenberg structure is lost anyway.
 138         endif
 139         V = V(:,1:(nv-1));
 140         ## One less reflection.
 141         nu--;
 142       else
 143         ## New householder reflection.
 144         if (pflg)
 145           ## Locate max magnitude element in short_q.
 146           asq = abs (short_q);
 147           maxv = max (asq);
 148           maxidx = find (asq == maxv, 1);
 149           pivot_idx = short_pv(maxidx);
 150
 151           ## See if need to change the pivot list.
 152           if (pivot_idx != pivot_vec(nu))
 153             swapidx = maxidx + (nu-1);
 154             [pivot_vec(nu), pivot_vec(swapidx)] = ...
 155                 swap (pivot_vec(nu), pivot_vec(swapidx));
 156           endif
 157         endif
 158
 159         ## Isolate portion of vector for reflection.
 160         idx = pivot_vec(nu:na);
 161         jdx = pivot_vec(1:nu);
 162
 163         [hv, av, z] = housh (q(idx), 1, 0);
 164         alpha(nu) = av;
 165         U(idx,nu) = hv;
 166
 167         ## Reduce V per the reflection.
 168         V(idx,:) = V(idx,:) - av*hv*(hv' * V(idx,:));
 169         if(iter > 1)
 170           ## FIXME -- not done correctly for block case.
 171           H(nu,nu-1) = V(pivot_vec(nu),jj);
 172         endif
 173
 174         ## Advance to next column of V.
 175         jj++;
 176       endif
 177     endwhile
 178
 179     ## Check for oversize V (due to full rank).
 180     if ((columns (V) > na) && (length (alpha) == na))
 181       ## Trim to size.
 182       V = V(:,1:na);
 183     elseif (columns(V) > na)
 184       krylov_V = V;
 185       krylov_na = na;
 186       krylov_length_alpha = length (alpha);
 187       error ("krylov: this case should never happen; submit a bug report");
 188     endif
 189
 190     if (columns (V) > 0)
 191       ## Construct next Q and multiply.
 192       Q = zeros (size (V));
 193       for kk = 1:columns (Q)
 194         Q(pivot_vec(nu-columns(Q)+kk),kk) = 1;
 195       endfor
 196
 197       ## Apply Householder reflections.
 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
 204     endif
 205
 206     ## Multiply to get new vector.
 207     V = A*Q;
 208     ## Project off of previous vectors.
 209     nu = length (alpha);
 210     for i = 1:nu
 211       hv = U(:,i);
 212       av = alpha(i);
 213       V = V - av*hv*(hv'*V);
 214       H(i,nu-columns(V)+(1:columns(V))) = V(pivot_vec(i),:);
 215     endfor
 216
 217   endwhile
 218
 219   ## Back out complete U matrix.
 220   ## back out U matrix.
 221   j1 = columns (U);
 222   for i = j1:-1:1;
 223     idx = pivot_vec(i:na);
 224     hv = U(idx,i);
 225     av = alpha(i);
 226     U(:,i) = zeros (na, 1);
 227     U(idx(1),i) = 1;
 228     U(idx,i:j1) = U(idx,i:j1)-av*hv*(hv'*U(idx,i:j1));
 229   endfor
 230
 231   nu = length (alpha);
 232   Uret = U;
 233   if (max (max (abs (Uret(zidx,:)))) > 0)
 234     warning ("krylov: trivial null space corrupted; set pflg = 1 or eps1 > %e",
 235              eps1);
 236   endif
 237
 238 endfunction
 239
 240
 241 function [a1, b1] = swap (a, b)
 242
 243   a1 = b;
 244   b1 = a;
 245
 246 endfunction