home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts.fat / la / krygetq.m < prev    next >
Text File  |  1999-04-29  |  2KB  |  61 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 QQ = krygetq(HH,alpha,kb)
  20. # function QQ = krygetq(HH,alpha,kb)
  21. # used internally by krylovb
  22. # construct last kb columns of orthogonal transform returned by qrhouse
  23. # Inputs:
  24. #   HH,alpha: set of nh Householder reflections set returned by qrhouse
  25. #   kb: desired number of columns
  26. # Outputs: 
  27. #   QQ: n x kb orthogonal matrix; basis of orthogonal subspace provided
  28. #      by final kb reflections in HH.
  29. #   Note: due to details in qrhouse and krylovb, QQ is *not* necessarily 
  30. #      the last kb columns defined by HH.  span(QQ) is orthogonal to the
  31. #      subspace spanned by the first (nh-kb) reflections in HH.
  32.  
  33. # Written by A. S. Hodel 1993--1995
  34. # $Revision$
  35. # $Log$
  36.  
  37. [n,m] = size(HH);
  38. kb1 = m+1-kb;
  39.  
  40. # construct permuted identity on the right rows:
  41. # since qrhouse removes zero rows, we've got to check the rows to which
  42. # the current householder reflections actually reflected their columns
  43. Hs = HH(:,kb1:m);
  44. if(is_vec(Hs))
  45.   nzidx = find(Hs' != 0);
  46. else
  47.   nzidx = find(max(abs(Hs')) != 0);
  48. endif
  49. nzidx = nzidx(1:kb);
  50. QQ = zeros(n,kb);
  51. QQ(nzidx,1:kb) = eye(kb);
  52.  
  53. # back out QQ matrix 
  54. for i=m:-1:1
  55.  hv = HH(:,i);
  56.  av = alpha(i);
  57.  QQ = QQ-av*hv*(hv'*QQ);
  58. end;
  59.  
  60. endfunction
  61.