SHARE
TWEET

demonio

a guest Nov 19th, 2019 85 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. % [H,Q] = hessred(A)
  2. %
  3. % Compute the Hessenberg decomposition H = Q’*A*Q using
  4. % Householder transformations.
  5. %
  6. function [H,Q] = hessred(A)
  7.  
  8. n = length(A);
  9. Q = eye(n); % Orthogonal transform so far
  10. H = A; % Transformed matrix so far
  11.  
  12. for j = 1:n-2
  13.  
  14. % -- Find W = I-2vv’ to put zeros below H(j+1,j)
  15. u = H(j+1:end,j);
  16. u(1) = u(1) + sign(u(1))*norm(u);
  17. v = u/norm(u);
  18.  
  19. % -- H := WHW’, Q := QW
  20. H(j+1:end,:) = H(j+1:end,:)-2*v*(v’*H(j+1:end,:));
  21. H(:,j+1:end) = H(:,j+1:end)-(H(:,j+1:end)*(2*v))*v’;
  22. Q(:,j+1:end) = Q(:,j+1:end)-(Q(:,j+1:end)*(2*v))*v’;
  23.  
  24. end
  25.  
  26. end
  27. % [H] = hessqr_basic(H)
  28. %
  29. % Compute one basic (unshifted) implicit Hessenberg QR step via
  30. % Householder transformations.
  31. %
  32. function H = hessqr_basic(H)
  33.  
  34. n = length(H);
  35. V = zeros(2,n-1);
  36.  
  37. % Compute the QR factorization
  38. for j = 1:n-1
  39.  
  40. % -- Find W_j = I-2vv’ to put zero into H(j+1,j)
  41. u = H(j:j+1,j);
  42. u(1) = u(1) + sign(u(1))*norm(u);
  43. v = u/norm(u);
  44. V(:,j) = v;
  45.  
  46. % -- H := W_j H
  47. H(j:j+1,:) = H(j:j+1,:)-2*v*(v’*H(j:j+1,:));
  48.  
  49. end
  50.  
  51. % Compute RQ
  52. for j = 1:n-1
  53.  
  54. % -- H := WHW’, Q := QW
  55. v = V(:,j);
  56. H(:,j:j+1) = H(:,j:j+1)-(H(:,j:j+1)*(2*v))*v’;
  57.  
  58. end
  59.  
  60. end
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top