function est = mynormest1 (L, U, P, Q) %MYNORMEST1 estimate norm(A,1), using LU factorization (L*U = P*A*Q). % % Example: % est = mynormest1 (L, U, P, Q) % See also: testall % CXSparse, Copyright (c) 2006-2022, Timothy A. Davis. All Rights Reserved. % SPDX-License-Identifier: LGPL-2.1+ n = size (L,1) ; est = 0 ; S = zeros (n,1) ; for k = 1:5 if k == 1 x = ones (n,1) / n ; else j = find (abs (x) == max (abs (x))) ; j = j (1) ; x = zeros (n,1) ; x (j) = 1 ; % fprintf ('eka: k %d j %d est %g\n', k, j, est) ; end % x=A\x, but use the existing P*A*Q=L*U factorization x = Q * (U \ (L \ (P*x))) ; est_old = est ; est = sum (abs (x)) ; unchanged = 1 ; for i = 1:n if (x (i) >= 0) s = 1 ; else s = -1 ; end if (s ~= S (i)) S (i) = s ; unchanged = 0 ; end end if (any (S ~= signum (x))) S' %#ok signum(x)' %#ok error ('Hey!') ; end if k > 1 & (est <= est_old | unchanged) %#ok break ; end x = S ; % x=A'\x, but use the existing P*A*Q=L*U factorization x = P' * (L' \ (U' \ (Q'*x))) ; if k > 1 jnew = find (abs (x) == max (abs (x))) ; if (jnew == j) break ; end end end for k = 1:n x (k) = power (-1, k+1) * (1 + ((k-1)/(n-1))) ; end % x=A\x, but use the existing P*A*Q=L*U factorization x = Q * (U \ (L \ (P*x))) ; est_new = 2 * sum (abs (x)) / (3 * n) ; if (est_new > est) est = est_new ; end function s = signum (x) %SIGNUM compute sign of x s = ones (length (x),1) ; s (find (x < 0)) = -1 ; %#ok