home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts.fat / la / qrhouse.m < prev    next >
Text File  |  1999-04-29  |  3KB  |  88 lines
  1. # Copyright (C) 1992, 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 [hv,alph,kb] = qrhouse(VV,eps1)
  20. # function [hv,alph,kb] = qrhouse(VV{,eps1})
  21. # construct orthogonal basis of span(VV) with Householder vectors
  22. # Q R = VV; Q may be obtained via routine krygetq; R is upper triangular
  23. #   if all rows of VV are nonzero; otherwise it's a permuted uppert
  24. #   triangular matrix with zero rows matching those of VV.
  25. # Inputs: 
  26. #   VV: matrix
  27. #   eps1: zero threshhold; default: 0.  Used to check if a column
  28. #     of reduced R is already upper triangular; entries with magnitude
  29. #     <= eps1 are considered to be 0.
  30. # Outputs:
  31. #   hv: matrix of householder reflection vectors as returned by housh
  32. #   alpha: vector of householder reflection values as returned by housh
  33. #   kb: computed rank of matrix
  34. # qrhouse is used in krylovb for block Arnoldi iteration
  35. #
  36. # Reference: Golub and Van Loan, MATRIX COMPUTATIONS, 2nd ed.
  37.  
  38. # Written by A. S. Hodel, 1992
  39.  
  40. if(nargin < 1 | nargin > 2)
  41.   usage("[hv,alph,kb] = qrhouse(VV{,eps1})");
  42. elseif(nargin == 1)     # default value for eps set to 0
  43.   eps1 = 0;
  44. endif
  45.  
  46.  
  47. # extract only those rows of VV that are nonzero
  48. if(is_vec(VV))    nzidx = find(abs(VV') > 0);
  49. else            nzidx = find(max(abs(VV')) > 0);    endif
  50. VVlen = rows(VV);    # remember original size
  51.  
  52. if(is_vec(VV))    VV = VV(nzidx);
  53. else            VV = VV(nzidx,:);                   endif
  54.  
  55. [Vr,Vc] = size(VV);    nits   = min(Vr,Vc);
  56. for ii = 1:nits;
  57.   # permute maximum row entry to (ii,ii) position
  58.   Vrowi = VV(ii,1:Vc);      # pivot maximum entry in this row to lead position
  59.   Vrm = max(abs(Vrowi));
  60.   Vmidx = min(find(abs(Vrowi) == Vrm));
  61.   if(Vmidx > eps1)
  62.     kb = ii;        # update computed rank
  63.     idx = kb-1;
  64.     if(Vmidx != ii)
  65.       [VV(:,kb),VV(:,Vmidx)] = swap(VV(:,kb),VV(:,Vmidx));
  66.     endif
  67.     hh = VV(:,ii);    # extract next column of VV; ignore items 1:(ii-1).
  68.     [hv(kb:Vr,kb),alph(kb),z] = housh(hh(kb:Vr),1,0);
  69.     if(kb>1)
  70.       hv(1:idx,kb) = 0;                 # zero top of hv for safety
  71.     endif
  72.     # project off of current Householder vector
  73.     VV = VV - alph(kb)*hv(:,kb)*(hv(:,kb)'*VV);
  74.   else
  75.     break;
  76.   endif
  77. endfor
  78. if(kb <=0)
  79.   hv = [];
  80.   alph = [];
  81. else
  82.   hvs = hv(:,1:kb);    # remove extraneous vectors, expand to original size
  83.   hv = zeros(VVlen,kb);
  84.   hv(nzidx,:) = hvs;
  85.   alph = alph(1:kb);
  86. end
  87. endfunction
  88.