www.pudn.com > RobustSF.zip > update_gbound.m, change:1998-09-09,size:2599b

```function anorm = update_gbound(anorm,alpha,beta,j)
%UPDATE_GBOUND   Update Gerscgorin estimate of 2-norm
%  ANORM = UPDATE_GBOUND(ANORM,ALPHA,BETA,J) updates the Gersgorin bound
%  for the tridiagonal in the Lanczos process after the J'th step.
%  Applies Gerscgorins circles to T_K'*T_k instead of T_k itself
%  since this gives a tighter bound.

if j==1 % Apply Gerscgorin circles to T_k'*T_k to estimate || A ||_2
i=j;
% scale to avoid overflow
scale = max(abs(alpha(i)),abs(beta(i+1)));
alpha(i) = alpha(i)/scale;
beta(i) = beta(i)/scale;
anorm = 1.01*scale*sqrt(alpha(i)^2+beta(i+1)^2 + abs(alpha(i)*beta(i+1)));
elseif j==2
i=1;
% scale to avoid overflow
scale = max(max(abs(alpha(1:2)),max(abs(beta(2:3)))));
alpha(1:2) = alpha(1:2)/scale;
beta(2:3) = beta(2:3)/scale;

anorm = max(anorm, scale*sqrt(alpha(i)^2+beta(i+1)^2 + ...
abs(alpha(i)*beta(i+1) + alpha(i+1)*beta(i+1)) + ...
abs(beta(i+1)*beta(i+2))));
i=2;
anorm = max(anorm,scale*sqrt(abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
abs(alpha(i)*beta(i+1))) );
elseif j==3
% scale to avoid overflow
scale = max(max(abs(alpha(1:3)),max(abs(beta(2:4)))));
alpha(1:3) = alpha(1:3)/scale;
beta(2:4) = beta(2:4)/scale;
i=2;
anorm = max(anorm,scale*sqrt(abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
abs(alpha(i)*beta(i+1) + alpha(i+1)*beta(i+1)) + ...
abs(beta(i+1)*beta(i+2))) );
i=3;
anorm = max(anorm,scale*sqrt(abs(beta(i)*beta(i-1)) + ...
abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
abs(alpha(i)*beta(i+1))) );
else
% scale to avoid overflow
%  scale = max(max(abs(alpha(j-2:j)),max(abs(beta(j-2:j+1)))));
%  alpha(j-2:j) = alpha(j-2:j)/scale;
%  beta(j-2:j+1) = beta(j-2:j+1)/scale;

% Avoid scaling, which is slow. At j>3 the estimate is usually quite good
% so just make sure that anorm is not made infinite by overflow.
i = j-1;
anorm1 = sqrt(abs(beta(i)*beta(i-1)) + ...
abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
abs(alpha(i)*beta(i+1) + alpha(i+1)*beta(i+1)) + ...
abs(beta(i+1)*beta(i+2)));
if isfinite(anorm1)
anorm = max(anorm,anorm1);
end
i = j;
anorm1 = sqrt(abs(beta(i)*beta(i-1)) + ...
abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
abs(alpha(i)*beta(i+1)));
if isfinite(anorm1)
anorm = max(anorm,anorm1);
end
end
```