home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPARFUN.DI$ / CONDEST.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  2.1 KB  |  85 lines
  1. function [c, v] = condest(A)
  2. %CONDEST Estimate the 1-norm condition number.
  3. %    Higham's modification of Hager's method.
  4. %    [C, V] = CONDEST(A) computes a lower bound C for the 1-norm condition
  5. %    number of A, and a vector V such that NORM(A*V)=NORM(A)*NORM(V)/C.
  6. %    V is an approximate null vector if C is large.
  7.  
  8. %    Reference: N.J. Higham, FORTRAN codes for estimating the one-norm of
  9. %    a real or complex matrix, with applications to condition estimation,
  10. %    ACM Trans. Math. Soft., 14 (1988), pp. 381-396.
  11.  
  12. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  13.  
  14. n = max(size(A));
  15. one = ones(n,1);
  16. maxits = 5;                    % Arbitrary limit on the number of iterations.
  17.  
  18. [L,U] = lu(A);
  19. k = find(abs(diag(U))==0);
  20. if any(k)
  21.    c = Inf;
  22.    v = zeros(n,1);
  23.    k = min(k);
  24.    v(k) = 1;
  25.    if k > 1
  26.       v(1:k-1) = U(1:k-1,1:k-1)\(-U(1:k-1,k));
  27.    end
  28.    v = v/norm(v,1);
  29.    return
  30. end
  31.  
  32. if n == 1
  33.    % A nonzero scalar is perfectly conditioned.
  34.    c = 1;
  35.    v = 1;
  36.    return
  37. end 
  38.  
  39. x = ones(n,1)/n;
  40.  
  41. v = U\(L\x);      % Av=x
  42. gamma = norm(v, 1);
  43. xi = sign(v); f = find(xi==0); xi(f) = one(f);   % set sign(0) := 1
  44. lastsign = xi;
  45. x = L'\(U'\xi);   % A'x=xi
  46. x = real(x);      % Needed only for the complex case.
  47.  
  48. for k=2:maxits
  49.  
  50.     % Get smallest index for which max(abs(z)) is attained.
  51.     [big, j] = max(abs(x)); f = find(abs(x)==big); j = f(1);
  52.     x = zeros(n,1); x(j)=1;
  53.  
  54.     v = U\(L\x);      % Av=x
  55.     xi = sign(v); f = find(xi==0); xi(f) = one(f);   % set sign(0) := 1
  56.  
  57.     vold = gamma;
  58.     gamma = norm(v, 1);
  59.  
  60.     if xi == lastsign | gamma <= vold
  61.        break
  62.     end
  63.     lastsign = xi;
  64.  
  65.     x = L'\(U'\xi);   % A'x=xi
  66.     x = real(x);      % Needed only for the complex case.
  67.  
  68.     if norm(x, inf) == x(j), break, end
  69.     if k == maxits
  70.        disp(['Warning: Not converged after ' int2str(maxits) ' iterations.'])
  71.     end
  72. end
  73.  
  74. x = (0:n-1)';
  75. x = (1-2*rem(x,2)) .* (1 + x/(n-1));
  76. x = U\(L\x);
  77. temp = 2*norm(x,1)/(3*n);
  78. if temp > gamma
  79.    v = x; gamma = temp;
  80. end
  81.  
  82. anorm = norm(abs(A),1);  % abs for complex case so true 1-norm used.
  83. c = gamma*anorm;
  84. v = v/norm(v,1);
  85.