text: 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
  