eiglancz.m

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

function [lmin,lmax,deltaeig,k]=eiglancz(A,m,toll)
n=size(A); V=[0*[1:n]',[1,0*[1:n-1]]'];
beta(1)=0; normb=1; k=1; deltaeig(1)=1;
while  k <= m & normb >= eps & deltaeig(k) < toll
  vk = V(:,k+1);   w = A*vk-beta(k)*V(:,k);
  alpha(k)= w'*vk; w = w - alpha(k)*vk;
  normb = norm(w,2); beta(k+1)=normb;
  if normb ~= 0
    V=[V,w/normb];
    if k==1
      lmin(1)=alpha; lmax(1)=alpha;
      k=k+1; deltaeig(k)=1;
    else
      d=alpha; b=beta(2:length(beta)-1);
      [ak,bk,ck,nch,niter]=givsturm(d,b,1,toll);
      lmax(k)=(ak(niter)+bk(niter))/2;
      [ak,bk,ck,nch,niter]=givsturm(d,b,k,toll);
      lmin(k)=(ak(niter)+bk(niter))/2;
      deltaeig(k+1)=max(abs(lmin(k)-lmin(k-1)),abs(lmax(k)-lmax(k-1)));
      k=k+1;
    end
  else
    disp('Breakdown');
    d=alpha; b=beta(2:length(beta)-1);
    [ak,bk,ck,nch,niter]=givsturm(d,b,1,toll);
    lmax(k)=(ak(niter)+bk(niter))/2;
    [ak,bk,ck,nch,niter]=givsturm(d,b,k,toll);
    lmin(k)=(ak(niter)+bk(niter))/2;
    deltaeig(k+1)=max(abs(lmin(k)-lmin(k-1)),abs(lmax(k)-lmax(k-1)));
    k=k+1;
  end
end
k=k-1;
return