www.gusucode.com > matlab矩阵函数工具箱 > matlab矩阵函数工具箱/matlab矩阵函数工具箱/signm_newton.m

    function [X,k] = signm_newton(A,scal)
%SIGNM_NEWTON   Matrix sign function by Newton iteration.
%   [S,k] = SIGNM_NEWTON(A,SCAL) computes the matrix sign
%   function S of A using the scaled Newton iteration, with
%   scale factors specified by SCAL:
%   SCAL = 0: no scaling.
%   SCAL = 1: determinantal scaling (default).
%   SCAL = 2: spectral scaling.
%   SCAL = 3: norm scaling.
%   k is the number of iterations.

if nargin < 2, scal = 1; end

tol = mft_tolerance(A);
accel_tol = 1e-2;  %  Precise value not important.
need_accel = 1;
n = length(A);
maxit = 16;

X = A;
reldiff = inf;

for k = 1:maxit

      Xold = X;
      Xinv = inv(X);
      if need_accel && scal > 0
         % In practice should estimate spectral radius and 2-norms;
         % here they are computed exactly.
         switch scal
         case 1
              g = abs(det(X))^(-1/n);
         case 2
              s1 = max(abs(eig(Xinv)));
              s2 = max(abs(eig(X)));
              g = sqrt(s1/s2);
         case 3
              g = sqrt( norm(Xinv) / (norm(X)) );
         end
         X = g*X; Xinv = Xinv/g;
      end
      X = 0.5*(X + Xinv);
      diff_F = norm(X-Xold,'fro');
      reldiff_old = reldiff;
      reldiff = diff_F/norm(X,'fro');

      if need_accel && (reldiff < accel_tol), need_accel = false; end
      cged = (diff_F <= sqrt( tol*norm(X)/norm(Xinv) ) || ...
              reldiff > reldiff_old/2 && ~need_accel);
      if cged, return, end

end
error('Not converged after %2.0f iterations', maxit)