home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts / linear-algebra / krylov.m < prev    next >
Text File  |  1999-10-13  |  6KB  |  202 lines
  1. # Copyright (C) 1993, 1998, 1999 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 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
  18.  
  19. function [Uret,H,nu] = krylov(A,V,k,eps1,pflg);
  20.   # function [U,H,nu] = krylov(A,V,k{,eps1,pflg});
  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.   # pflg: flag to use row pivoting  (improves numerical behavior)
  27.   #   0 [default]: no pivoting; prints a warning message if trivial
  28.   #                null space is corrupted
  29.   #   1          : pivoting performed
  30.   #
  31.   # outputs:
  32.   #   Uret: orthogonal basis of block krylov subspace
  33.   #   H: Hessenberg matrix; if V is a vector then A U = U H
  34.   #      otherwise H is meaningless
  35.   # nu: dimension of span of krylov subspace (based on eps1)
  36.   # if B is a vector and k > m-1, krylov returns H = the Hessenberg
  37.   # decompostion of A.
  38.   #
  39.   # Reference: Hodel and Misra, "Partial Pivoting in the Computation of
  40.   #     Krylov Subspaces", to be submitted to Linear Algebra and its
  41.   #     Applications
  42.   # written by A. Scottedward Hodel a.s.hodel@eng.auburn.edu
  43.  
  44.   defeps = 1e-12;
  45.   if(nargin < 3 | nargin > 5)
  46.     usage("[U,nu] = krylov(A,V,k{,eps1,pflg})")
  47.   elseif(nargin < 5)
  48.     pflg = 0;        # default permutation flag
  49.   endif
  50.   if(nargin < 4)
  51.     eps1 = defeps;    # default tolerance parameter
  52.   endif
  53.   if(isempty(eps1)) eps1 = defeps; endif
  54.  
  55.   na = is_square(A);
  56.   if( !na ) error("A(%d x %d) must be square",rows(A),columns(A)); endif
  57.  
  58.   [m,kb] = size(V);
  59.   if(m != na);
  60.     error("A(%d x %d), V(%d x %d): argument dimensions do not match", ...
  61.       na,na,m,kb)
  62.   endif
  63.  
  64.   if( !is_scalar(k) )
  65.     error("krylov: third argument must be a scalar integer")
  66.   endif
  67.  
  68.   Vnrm = norm(V,Inf);
  69.  
  70.   # check for trivial solution
  71.   if(Vnrm == 0)
  72.     Uret = []; nu = 0;  return;
  73.   endif
  74.  
  75.   # identify trivial null space
  76.   abm = max(abs([A,V]'));  nzidx = find(abm != 0);  zidx = find(abm == 0);
  77.  
  78.   # set up vector of pivot points
  79.   pivot_vec = 1:na;
  80.  
  81.   iter = 0;
  82.   alpha = [];
  83.   nh = 0;
  84.   while (length(alpha) < na) & (columns(V) > 0) & (iter < k)
  85.     iter++;
  86.  
  87.     # get orthogonal basis of V
  88.     jj = 1;
  89.     while(jj <= columns(V) & length(alpha) < na)
  90.       nu = length(alpha)+1;   # index of next Householder reflection
  91.  
  92.       short_pv = pivot_vec(nu:na);
  93.       q = V(:,jj);
  94.       short_q = q(short_pv);
  95.  
  96.       if(norm(short_q) < eps1)
  97.         # insignificant column; delete
  98.         nv = columns(V);
  99.         if(jj != nv)
  100.           [V(:,jj),V(:,nv)] = swap(V(:,jj),V(:,nv));
  101.           # FIX ME: H columns should be swapped too.  Not done since
  102.           # Block Hessenberg structure is lost anyway.
  103.         endif
  104.         V = V(:,1:(nv-1));
  105.         nu = nu - 1;    # one less reflection
  106.  
  107.       else
  108.         # new householder reflection
  109.         if(pflg)
  110.           # locate max magnitude element in short_q
  111.           asq = abs(short_q);
  112.       maxv = max(asq);
  113.           maxidx = find(asq == maxv);
  114.           pivot_idx = short_pv(maxidx(1));
  115.  
  116.           # see if need to change the pivot list
  117.           if(pivot_idx != pivot_vec(nu))
  118.             swapidx = maxidx(1) + (nu-1);
  119.             [pivot_vec(nu),pivot_vec(swapidx)] = ...
  120.               swap(pivot_vec(nu),pivot_vec(swapidx));
  121.           endif
  122.         endif
  123.  
  124.         # isolate portion of vector for reflection
  125.         idx = pivot_vec(nu:na);
  126.         jdx = pivot_vec(1:nu);
  127.  
  128.         [hv,av,z] = housh(q(idx),1,0);
  129.         alpha(nu) = av;
  130.         U(idx,nu) = hv;
  131.  
  132.         # reduce V per the reflection
  133.         V(idx,:) = V(idx,:) - av*hv*(hv' * V(idx,:));
  134.         if(iter > 1)
  135.           # FIX ME: not done correctly for block case
  136.           H(nu,nu-1) = V(pivot_vec(nu),jj);
  137.         endif
  138.  
  139.         # advance to next column of V
  140.         jj=jj+1;
  141.       endif
  142.     endwhile
  143.  
  144.     # check for oversize V (due to full rank)
  145.     if( ( columns(V) > na ) & ( length(alpha) == na ) )
  146.       # trim to size
  147.       V = V(:,1:na);
  148.     elseif( columns(V) > na )
  149.       krylov_V = V
  150.       krylov_na = na
  151.       krylov_length_alpha = length(alpha)
  152.       error("This case should never happen; submit bug report.");
  153.     endif
  154.  
  155.     if(columns(V) > 0)
  156.       # construct next Q and multiply
  157.       Q = zeros(size(V));
  158.       for kk=1:columns(Q)
  159.         Q(pivot_vec(nu-columns(Q)+kk),kk) = 1;
  160.       endfor
  161.  
  162.       # apply Householder reflections
  163.       for ii = nu:-1:1
  164.         idx = pivot_vec(ii:na);
  165.         hv = U(idx,ii);
  166.         av = alpha(ii);
  167.         Q(idx,:) = Q(idx,:) - av*hv*(hv'*Q(idx,:));
  168.       endfor
  169.     endif
  170.  
  171.     # multiply to get new vector;
  172.     V = A*Q;
  173.     # project off of previous vectors
  174.     nu = length(alpha);
  175.     for i=1:nu
  176.       hv = U(:,i);  av = alpha(i);
  177.       V = V - av*hv*(hv'*V);
  178.       H(i,nu-columns(V)+(1:columns(V))) = V(pivot_vec(i),:);
  179.     end
  180.  
  181.   endwhile
  182.  
  183.   # Back out complete U matrix
  184.   # back out U matrix ;
  185.   j1 = columns(U);
  186.   for i=j1:-1:1;
  187.     idx = pivot_vec(i:na);
  188.     hv = U(idx,i);
  189.     av = alpha(i);
  190.     U(:,i) = zeros(na,1);
  191.     U(idx(1),i) = 1;
  192.     U(idx,i:j1) = U(idx,i:j1)-av*hv*(hv'*U(idx,i:j1));
  193.   endfor
  194.  
  195.   nu = length(alpha);
  196.   Uret = U;
  197.   if( max(max( abs(Uret(zidx,:)) )) > 0)
  198.     warning("krylov: trivial null space corrupted; set pflg=1 or eps1>%e",eps1);
  199.   endif
  200.  
  201. endfunction
  202.