time   calls  line
                  1 function [est, v, w, iter] = normest1(A, t, X, varargin)
                  2 %NORMEST1 Estimate of 1-norm of matrix by block 1-norm power method.
                  3 %   C = NORMEST1(A) returns an estimate C of norm(A,1), where A is N-by-N.
                  4 %   A can be an explicit matrix or a function AFUN such that
                  5 %   FEVAL(@AFUN,FLAG,X) for the following values of
                  6 %     FLAG       returns
                  7 %     'dim'      N
                  8 %     'real'     1 if A is real, 0 otherwise
                  9 %     'notransp' A*X
                 10 %     'transp'   A'*X
                 11 %
                 12 %   C = NORMEST1(A,T) changes the number of columns in the iteration matrix
                 13 %   from the default 2.  Choosing T <= N/4 is recommended, otherwise it should
                 14 %   be cheaper to form the norm exactly from the elements of A, as is done
                 15 %   when N <= 4 or T == N.  If T < 0 then ABS(T) columns are used and trace
                 16 %   information is printed.  If T is given as the empty matrix ([]) then the
                 17 %   default T is used.
                 18 %
                 19 %   C = NORMEST1(A,T,X0) specifies a starting matrix X0 with columns of unit
                 20 %   1-norm and by default is random for T > 1.  If X0 is given as the empty
                 21 %   matrix ([]) then the default X0 is used.
                 22 %
                 23 %   C = NORMEST1(AFUN,T,X0,P1,P2,...) passes extra inputs P1,P2,... to
                 24 %   FEVAL(@AFUN,FLAG,X,P1,P2,...).
                 25 %
                 26 %   [C,V] = NORMEST1(A,...) and [C,V,W] = NORMEST1(A,...) also return vectors
                 27 %   V and W such that W = A*V and NORM(W,1) = C*NORM(V,1).
                 28 %
                 29 %   [C,V,W,IT] = NORMEST1(A,...) also returns a vector IT such that
                 30 %   IT(1) is the number of iterations,
                 31 %   IT(2) is the number of products of N-by-N by N-by-T matrices.
                 32 %   On average, IT(2) = 4.
                 33 %
                 34 %   Note: NORMEST1 invokes RAND.  If repeatable results are required,  then
                 35 %   see RAND for details on how to set the default stream state.
                 36 %
                 37 %   See also CONDEST, COND, NORM, RAND.
                 38 
                 39 %   Subfunctions: MYSIGN, UNDUPLI, NORMAPP.
                 40 
                 41 %   Reference: 
                 42 %   [1] Nicholas J. Higham and Fran\c{c}oise Tisseur, 
                 43 %       A Block Algorithm for Matrix 1-Norm Estimation 
                 44 %       with an Application to 1-Norm Pseudospectra, 
                 45 %       SIAM J. Matrix Anal. App. 21, 1185-1201, 2000. 
                 46 
                 47 %   Nicholas J. Higham
                 48 %   Copyright 1984-2010 The MathWorks, Inc.
                 49 %   $Revision: 1.8.4.9 $  $Date: 2010/11/17 11:29:17 $
                 50 
                 51 % Determine whether A is a matrix or a function.
           240   52 if isa(A,'double') 
                 53    n = max(size(A));
                 54    A_is_real = isreal(A);
  0.01     240   55 elseif isa(A,'function_handle') 
  0.01     240   56    n = normapp(A,'dim',[],varargin{:}); 
  0.03     240   57    A_is_real = normapp(A,'real',[],varargin{:}); 
                 58 else
                 59    error(message('MATLAB:normest1:ANotMatrixOrFunction'));
                 60 end
                 61 
           240   62 if nargin < 2 || isempty(t), t = 2; end 
           240   63 prnt = (t < 0); 
           240   64 t = abs(t); 
           240   65 if t ~= round(t) || t < 1 || t > max(n,2) 
                 66    error(message('MATLAB:normest1:TOutOfRange'))
                 67 end
           240   68 rpt_S = 0; rpt_e = 0; 
                 69 
           240   70 if t == n || n <= 4 
           240   71    if isa(A,'double') 
                 72       Y = A;
           240   73    else 
           240   74       X = eye(n); 
           240   75       Y = normapp(A,'notransp',X,varargin{:}); 
           240   76    end 
           240   77    [temp,m] = sort( sum(abs(Y)) ); 
           240   78    est = full(temp(n)); v = zeros(n,1); v(m(n)) = 1; w = Y(:,m(n)); 
           240   79    iter = [0 1]; 
           240   80    if prnt, fprintf('No iteration - norm computed exactly.\n'), end 
           240   81    return 
                 82 end
                 83 
                 84 if nargin < 3 || isempty(X)
                 85    X = ones(n,t);
                 86    X(:,2:t) = mysign(2*rand(n,t-1) - ones(n,t-1));
                 87    X = undupli(X, [], prnt);
                 88    X = X/n;
                 89 end
                 90 
                 91 if size(X,2) ~= t 
                 92   error(message('MATLAB:normest1:WrongColNum', int2str( t )));
                 93 end
                 94 
                 95 itmax = 5;  % Maximum number of iterations.
                 96 it = 0; nmv = 0;
                 97 
                 98 ind = zeros(t,1); S = zeros(n,t);
                 99 est_old = 0;
                100 
                101 while 1
                102     it = it + 1;
                103 
                104     if isa(A,'double')
                105        Y = A*X;
                106     else
                107        Y = normapp(A,'notransp',X,varargin{:});
                108     end
                109     nmv = nmv + 1;
                110 
                111     vals = sum(abs(Y));
                112     [~, m] = sort(vals); m = m(t:-1:1);
                113     vals = vals(m); vals_ind = ind(m);
                114     est = vals(1);
                115 
                116     if est > est_old || it == 2, est_j = vals_ind(1); w = Y(:,m(1)); end
                117 
                118     if prnt
                119        fprintf('%g: ', it)
                120        for i = 1:t, fprintf(' (%g, %6.2e)', vals_ind(i), vals(i)), end
                121        fprintf('\n')
                122     end
                123 
                124     if it >= 2 && est <= est_old
                125        est = est_old;
                126        info = 2; break
                127     end
                128     est_old = est;
                129 
                130     if it > itmax, it = itmax; info = 1; break, end
                131 
                132     S_old = S;
                133     S = mysign(Y);
                134     if A_is_real
                135         SS = S_old'*S; np = sum(max(abs(SS)) == n);
                136         if np == t, info = 3; break, end
                137         % Now check/fix cols of S parallel to cols of S or S_old.
                138         [S, r] = undupli(S, S_old, prnt); rpt_S = rpt_S + r;
                139     end
                140 
                141     if isa(A,'double')
                142        Z = A'*S;
                143     else
                144        Z = normapp(A,'transp',S,varargin{:});
                145     end
                146     nmv = nmv + 1;
                147 
                148     % Faster version of `for i=1:n, Zvals(i) = norm(Z(i,:), inf); end':
                149     Zvals = max(abs(Z),[],2);
                150 
                151     if it >= 2
                152        if max(Zvals) == Zvals(est_j), info = 4; break, end
                153     end
                154 
                155     [~, m] = sort(Zvals); m = m(n:-1:1);
                156     imax = t; % Number of new unit vectors; may be reduced below (if it > 1).
                157     if it == 1 
                158        ind = m(1:t);
                159        ind_hist = ind;
                160     else
                161        rep = sum(ismember(m(1:t), ind_hist));
                162        rpt_e = rpt_e + rep;
                163        if rep && prnt, fprintf('     rep e_j = %g\n',rep), end
                164        if rep == t, info = 5; break, end
                165        j = 1;
                166        for i = 1:t
                167            if j > n, imax = i-1; break, end
                168            while any( ind_hist == m(j) )
                169               j = j+1;
                170               if j > n, imax = i-1; break, end
                171            end
                172            if j > n, break, end
                173            ind(i) = m(j);
                174            j = j+1;
                175        end
                176        ind_hist = [ind_hist; ind(1:imax)];
                177     end
                178 
                179     X = zeros(n,t);
                180     for j=1:imax, X(ind(j),j) = 1; end
                181 end
                182 
                183 if prnt
                184    switch info
                185         case 1, fprintf('*Terminate: iteration limit reached.\n')
                186         case 2, fprintf('*Terminate: estimate not increased.\n')
                187         case 3, fprintf('*Terminate: repeated sign matrix.\n')
                188         case 4, fprintf('*Terminate: power method convergence test.\n')
                189         case 5, fprintf('*Terminate: repeated unit vectors.\n')
                190    end
                191 end
                192 
                193 iter = [it; nmv]; red = [rpt_S; rpt_e];
                194 
                195 v = zeros(n,1); v(est_j) = 1;
                196 if ~prnt, return, end
                197 
                198 if A_is_real, fprintf('Parallel cols: %g\n', red(1)), end
                199 fprintf('Repeated unit vectors: %g\n', red(2))

Other subfunctions in this file are not included in this listing.