home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts / linear-algebra / krylovb.m < prev    next >
Encoding:
Text File  |  1998-11-05  |  2.7 KB  |  94 lines
  1. # Copyright (C) 1993, 1998 A. Scottedward Hodel
  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 the
  7. # Free Software Foundation; either version 2, or (at your option) any
  8. # later version.
  9. #
  10. # Octave is distributed in the hope that it will be useful, but WITHOUT
  11. # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12. # FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13. # 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, write to the Free
  17. # Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  
  19. function [Uret,Ucols] = krylovb(A,V,k,eps1);
  20.   # function [U,Ucols] = krylovb(A,V,k[,eps1]);
  21.   # construct orthogonal basis U of block Krylov subspace;
  22.   #     [V AV A^2*V ... A^(k+1)*V];
  23.   # method used: householder reflections to guard against loss of
  24.   # orthogonality
  25.   # eps1: threshhold for 0 (default: 1e-12)
  26.   #
  27.   # outputs:
  28.   #   returned basis U is orthogonal matrix; due to "zeroed"
  29.   #   columns of product, may not satisfy A U = U H identity
  30.   # Ucols: dimension of span of krylov subspace (based on eps1)
  31.   # if k > m-1, krylov returns the Hessenberg decompostion of A.
  32. #
  33. # A. S. Hodel Feb 1993
  34. # $Revision$
  35. # $Log$
  36.  
  37.   if(nargin == 3)
  38.     eps1 = 1e-12;
  39.   endif
  40.   na = is_square(A);
  41.   if( !na )
  42.     error("krylov: first argument must be a square matrix")
  43.   endif
  44.  
  45.   [m,kb] = size(V); 
  46.   if(m != na);
  47.     error("krylov: argument dimensions do not match")
  48.   endif
  49.  
  50.   if( !is_scalar(k) )
  51.     error("krylov: third argument must be a scalar integer")
  52.   endif
  53.  
  54.   Vnrm = norm(V,Inf);
  55.  
  56.   # compute factored QR
  57.   [U,alpha,kb] = qrhouse(V,eps1*Vnrm);
  58.   Q = krygetq(U,alpha,kb);
  59.   Uret = Q;
  60.   Ucols = kb;
  61.  
  62.   iter = 0;
  63.   while (Ucols < na) & (kb > 0) & (iter < k)
  64.     iter++;
  65.     # multiply to get new vector;
  66.     V = A*Q;
  67.  
  68.     # project off of previous vectors
  69.     nzv = [];    # set of reflection indices used so far
  70.     nj = length(alpha);
  71.     for i=1:nj
  72.      hv = U(:,i);
  73.      nzidx = find(hv != 0);    # extract nonzero entries
  74.      nzv = union(nzv,nzidx(1)); # save for call to qrhouse
  75.      hv = hv(nzidx);
  76.      vr = V(nzidx,:); 
  77.      av = alpha(i);
  78.      V(nzidx,:) = vr - (av*hv)*(hv'*vr);
  79.     end
  80.     V(nzv,:) = 0;        # clear entries covered by earlier reflections
  81.  
  82.     [hvk,alphk,kb] = qrhouse(V,eps1*Vnrm);    # now get new QR factorization
  83.     if(kb > 0)
  84.       U = [U,hvk];        # append new data to Householder data structure
  85.       alpha = [alpha;alphk];
  86.       Q2 = krygetq(U,alpha,kb);
  87.       Uret = [Uret,Q2];
  88.  
  89.       Ucols = Ucols + kb;
  90.       Q = Q2;
  91.     end
  92.   end
  93. endfunction
  94.