invpower.m

来自「various matlab programs to slove various」· M 代码 · 共 33 行

function [sigma,x,niter,err]=invpower(A,z0,mu,toll,nmax)
n=max(size(A)); 
M=A-mu*eye(n); 
[L,U,P]=lu(M);
q=z0/norm(z0); 
q2=q'; 
err=[]; 
sigma=[]; 
res=toll+1; 
niter=0; 
while (res >= toll & niter <= nmax)
 niter=niter+1; 
 b=P*q; 
 y=L \ b; z = U \ y;
 q=z/norm(z); 
 z=A*q; 
 lam=q'*z;
 b=q2'; 
 y=U' \ b; 
 w=L' \ y;
 q2=(P'*w)'; 
 q2=q2/norm(q2); 
 costheta=abs(q2*q);
 if (costheta >= 5e-2),
    res=norm(z-lam*q)/costheta; 
    err=[err; res]; 
    sigma=[sigma; lam]; 
 else, 
    disp(' Multiple eigenvalue '); break; 
 end
x=q;
end
return