Re: [討論] matlab內建指令的.m檔原始碼

看板MATLAB作者snow3804 (snow3804)時間11年前 (2014/08/15 03:51)推噓0(0推 0噓 0→)

留言0則, 0人參與討論串2/2 (看更多)

※ 引述《snow3804 (snow3804)》之銘言： : matlab有hess(A)可以將矩陣A轉換成Hessenberg矩陣 : 但我想知道這中間是怎麼計算的 : 搜尋電腦裡有hess.m檔只是打開卻只看到註解而已 : http://tinyurl.com/m6v7p96 : 有沒有辦法可以知道是怎麼實作的推 iHakka: matlab是用lapack去實作的，可以追回去看看 08/13 19:27 推 infernodimon: lapack裡有印象中是用householder 08/14 19:50 我講一下原由,我在這裡查到Hessenberg的matlab程式碼 http://www.auburn.edu/~tamtiny/lecture%2010.pdf function [H,Q]=houshess(A) n=max(size(A)); Q=eye(n); H=A; for k=1:(n-2) [v,beta]=vhouse(H(k+1:n,k)); I=eye(k); N=zeros(k,n-k); m=length(v); R=eye(m)-beta*v*v'; HH=H(k+1:n,k:n); H(k+1:n,k:n)=R*H(k+1:n,k:n); HH=H(1:n,k+1:n); H(1:n,k+1:n)=H(1:n,k+1:n)*R; P=[I, N; N', R]; Q=Q*P; end return function [v,beta]=vhouse(x) n=length(x); x=x/norm(x); s=x(2:n)'*x(2:n); v=[1; x(2:n)]; if (s==0), beta=0; else mu=sqrt(x(1)^2+s); if (x(1) <= 0) v(1)=x(1)-mu; else v(1)=-s/(x(1)+mu); end beta=2*v(1)^2/(s+v(1)^2); v=v/v(1); end return 但執行的結果和matla內建的hess()結果不一樣 http://www.mathworks.com/help/matlab/ref/hess.html A=[-149 -50 -154;537 180 546;-27 -9 -25]; H=houshess(A) H = -149.0000 -42.2037 156.3165 537.6783 152.5511 -554.9272 0 0.0728 2.4489 H=hess(A) H = -149.0000 42.2037 -156.3165 -537.6783 152.5511 -554.9272 0 0.0728 2.4489 看得出來其實只是有些地方差個負號而已當然這兩個都是Hessenberg矩陣但我還是想知道程式碼是差在哪裡造成正負號不同 lapack我努力去找到程式碼,只是lapack是用fortran寫的,但我看不懂程式碼 http://www.netlib.org/lapack/explore-3.1.1-html/sgehrd.f.html -- ※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 114.42.124.128 ※ 文章網址: http://www.ptt.cc/bbs/MATLAB/M.1408045892.A.54A.html ※ 編輯: snow3804 (114.42.124.128), 08/15/2014 04:00:43