home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts.fat / la / krylov.m < prev    next >
Text File  |  1999-04-29  |  5KB  |  163 lines
  1. # Copyright (C) 1996,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 [U,H,k1] = krylov(A,v,k,eps1,pflg);
  20.   # function [U,H,k1] = krylov(A,v,k{,eps1});
  21.   # construct orthogonal basis U of Krylov subspace;
  22.   #     span([v Av A^2*v ... A^(k+1)*v]);
  23.   # via householder reflections; reflections are multiplied to obtain U
  24.   # inputs:
  25.   #   A: square matrix
  26.   #   v: column vector
  27.   #   k: desired krylov subspace dimension (as above)
  28.   #   eps1: threshhold for 0 relative to norm of current column (default: 1e-12)
  29.   #   pflg: permutation flag (default 0): avoid using zero rows of
  30.   #      [A,v] as householder pivots; this avoids spurious non-zero entries
  31.   #      in returned orthogonal matrix U, but destroys the Householder matrix
  32.   #      structure of H.
  33.   # outputs:
  34.   #   U: (n x k+1) orthogonal basis of Krylov subspace. A*U(:,1:k) = U*H
  35.   #   H: (pflg=0): Hessenberg matrix satisfying A*U(:,1:k) = U*H
  36.   #      (pflg=1): Workspace; does not satisfy above equation.
  37.   # k1: dimension of span of krylov subspace (based on eps1)
  38.   # if k > m-1, krylov returns the Hessenberg decompostion of A.
  39.  
  40.   # Written by A. S. Hodel 1992
  41.   # $Revision: 1.2 $
  42.   # $Log$
  43.  
  44.   save_empty_list_elements_ok = empty_list_elements_ok;
  45.   empty_list_elements_ok = 1;
  46.  
  47.   if    (nargin > 5)   usage("[U,H,k1] = krylov(A,v,k{,eps1,pflg})");
  48.   elseif(nargin < 5)   pflg = 0;
  49.   elseif(nargin < 4)   eps1 = 1e-12; endif
  50.   na = is_sqr(A);
  51.   if(!na) error("krylov: A(%dx%d) must be square",na,na); endif
  52.   [m,n] = size(v);
  53.   if(!is_vec(v)) error("krylov: v(%dx%d) must be a vector",m,n);
  54.   elseif(length(v) != na)
  55.     error("krylov: A(%dx%d), v(%dx1); mismatch",na,na,length(v));
  56.   endif
  57.   v = vec(v);    # make sure it's a column vector
  58.   if(!is_scal(k))
  59.     error("krylov: k(%dx%d) must be a scalar",rows(k), columns(k));
  60.   elseif( k > m)
  61.     warning("krylov: k is too large; reducing to %d",m-1);
  62.     k = m-1;
  63.   endif
  64.  
  65.   # check for zero input vector
  66.   if(norm(v) == 0) U = []; H = []; k1 = 0; return endif
  67.  
  68.   # indices of legal pivot points (identifies trivial null space).
  69.   abm = max(abs([A,v]')); nzidx = find(abm != 0); zidx = find(abm == 0);
  70.  
  71.   # check permutation flag
  72.   if(pflg)
  73.     # permute zero rows of [A,v] to trailing rows
  74.     permvec = [vec(nzidx); vec(zidx)];
  75.     pmat = eye(na); pmat = pmat(permvec,:);
  76.     ipermvec = pmat'*vec(1:na);
  77.     A = A(permvec,permvec);
  78.     v = v(permvec);
  79.   endif
  80.  
  81.   k1 = k+1;        # Assume subspace basis has full rank
  82.   [m,n] = size(A);
  83.   [hv,alpha(1),z] = housh(v,1,0);
  84.  
  85.   # initial orthogonal vector
  86.   q = zeros(n,1);
  87.   q(1) = 1;
  88.   q = q - alpha*hv*hv'*q;
  89.   normq = norm(q);
  90.   normres = normq;
  91.  
  92.   U(:,1) = hv;
  93.   j = 1;
  94.   while(j <= k & normres > eps1*normq)
  95.     # multiply to get new vector;
  96.     q = A*q;
  97.     normq = norm(q);
  98.  
  99.     # multiply by householder reflections to obtain projected vector and the
  100.     # next column of H
  101.     for i=1:j
  102.      hv = U(i:n,i);
  103.      av = alpha(i);
  104.      q(i:n,1) = q(i:n,1)-av*hv*(hv'*q(i:n,1));
  105.     endfor
  106.  
  107.     i = j+1;
  108.     # compute and apply next householder vector;
  109.     if(i <= n)
  110.       [hv,av,z] = housh(q(i:n,1),1,0);
  111.       alpha(i) = av;
  112.       q(i:n,1) = q(i:n,1)-av*hv*(hv'*q(i:n,1));
  113.       U(i:n,i) = hv;
  114.       H(1:i,j) = q(1:i);
  115.     else
  116.       av = 0;
  117.       H(:,j) = q;    # complete Hessenberg decomposition
  118.     endif
  119.  
  120.     # see if done yet
  121.     normres = norm(q(i:n));
  122.     if(normres > eps1*normq)
  123.       j = j+1;
  124.     else   
  125.       k1 = min(k1,j);    # time to quit; norm of residual is small
  126.     endif
  127.   
  128.     # back out new q vector for next pass;
  129.     j1 = columns(U);
  130.     q = zeros(n,1);
  131.     q(j1) = 1;
  132.     for i=j1:-1:1;
  133.       hv = U(i:n,i);
  134.       av = alpha(i);
  135.       q(i:n,1) = q(i:n,1)-av*hv*(hv'*q(i:n,1));
  136.     endfor
  137.   endwhile
  138.  
  139.   # back out U matrix ;
  140.   j1 = columns(U);
  141.   for i=j1:-1:1;
  142.    hv = U(i:n,i);
  143.    av = alpha(i);
  144.    U(:,i) = zeros(n,1);
  145.    U(i,i) = 1;
  146.    U(i:n,i:j1) = U(i:n,i:j1)-av*hv*(hv'*U(i:n,i:j1));
  147.   endfor
  148.  
  149.   # check permutation flag
  150.   if(pflg)
  151.     # permute rows of U back to original coordinates
  152.     U = U(ipermvec,:);
  153.   endif
  154.  
  155.   # check for spurious nonzero entries
  156.   if( max(max( abs(U(zidx,:)) )) > eps1 )
  157.     warning("krylov: trivial null space corrupted; set pflg=1 or eps1>%e",eps1);
  158.   endif
  159.  
  160.   empty_list_elements_ok = save_empty_list_elements_ok;
  161.  
  162. endfunction
  163.