# HG changeset patch # User Antonio Pino Robles # Date 1432885716 -7200 # Node ID a587712dcf5fec3ff12915a43d29344020829f54 # Parent e0817ccb28345e66c03430b949898afbeea8a12a funm_atom.m: rename fun_atom to funm_atom * funm_atom.m: rename fun_atom to funm_atom. diff -r e0817ccb2834 -r a587712dcf5f funm_files/fun_atom.m --- a/funm_files/fun_atom.m Tue May 26 18:14:54 2015 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,78 +0,0 @@ -function [F,n_terms] = fun_atom(T,fun,tol,prnt) -%FUNM_ATOM Function of triangular matrix with nearly constant diagonal. -% [F, N_TERMS] = FUNM_ATOM(T, FUN, TOL, PRNT) evaluates function -% FUN at the upper triangular matrix T, where T has nearly constant -% diagonal. A Taylor series is used. -% FUN(X,K) must return the K'th derivative of -% the function represented by FUN evaluated at the vector X. -% TOL is a convergence tolerance for the Taylor series, -% defaulting to EPS. -% If PRNT ~= 0 trace information is printed. -% N_TERMS is the number of terms taken in the Taylor series. -% N_TERMS = -1 signals lack of convergence. - -if nargin < 3 | isempty(tol), tol = eps; end -if nargin < 4, prnt = 0; end - -if isequal(fun,@fun_log) % LOG is special case. - [F,n_terms] = logm_isst(T,prnt); - return -end - -itmax = 500; - -n = length(T); -if n == 1, F = feval(fun,T,0); n_terms = 1; return, end - -lambda = sum(diag(T))/n; -F = eye(n)*feval(fun,lambda,0); -f_deriv_max = zeros(itmax+n-1,1); -N = T - lambda*eye(n); -mu = norm( (eye(n)-abs(triu(T,1)))\ones(n,1),inf ); - -P = N; -max_d = 1; - -for k = 1:itmax - - f = feval(fun,lambda,k); - F_old = F; - F = F + P*f; - rel_diff = norm(F - F_old,inf)/(tol+norm(F_old,inf)); - if prnt - fprintf('%3.0f: coef = %5.0e', k, abs(f)/factorial(k)); - fprintf(' N^k/k! = %7.1e', norm(P,inf)); - fprintf(' rel_d = %5.0e',rel_diff); - fprintf(' abs_d = %5.0e',norm(F - F_old,inf)); - end - P = P*N/(k+1); - - if rel_diff <= tol - - % Approximate the maximum of derivatives in convex set containing - % eigenvalues by maximum of derivatives at eigenvalues. - for j = max_d:k+n-1 - f_deriv_max(j) = norm(feval(fun,diag(T),j),inf); - end - max_d = k+n; - omega = 0; - for j = 0:n-1 - omega = max(omega,f_deriv_max(k+j)/factorial(j)); - end - - trunc = norm(P,inf)*mu*omega; % norm(F) moved to RHS to avoid / 0. - if prnt - fprintf(' [trunc,test] = [%5.0e %5.0e]', ... - trunc, tol*norm(F,inf)) - end - if prnt == 5, trunc = 0; end % Force simple stopping test. - if trunc <= tol*norm(F,inf) - n_terms = k; - if prnt, fprintf('\n'), end, return - end - end - - if prnt, fprintf('\n'), end - -end -n_terms = -1; diff -r e0817ccb2834 -r a587712dcf5f funm_files/funm_atom.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/funm_files/funm_atom.m Fri May 29 09:48:36 2015 +0200 @@ -0,0 +1,78 @@ +function [F,n_terms] = fun_atom(T,fun,tol,prnt) +%FUNM_ATOM Function of triangular matrix with nearly constant diagonal. +% [F, N_TERMS] = FUNM_ATOM(T, FUN, TOL, PRNT) evaluates function +% FUN at the upper triangular matrix T, where T has nearly constant +% diagonal. A Taylor series is used. +% FUN(X,K) must return the K'th derivative of +% the function represented by FUN evaluated at the vector X. +% TOL is a convergence tolerance for the Taylor series, +% defaulting to EPS. +% If PRNT ~= 0 trace information is printed. +% N_TERMS is the number of terms taken in the Taylor series. +% N_TERMS = -1 signals lack of convergence. + +if nargin < 3 | isempty(tol), tol = eps; end +if nargin < 4, prnt = 0; end + +if isequal(fun,@fun_log) % LOG is special case. + [F,n_terms] = logm_isst(T,prnt); + return +end + +itmax = 500; + +n = length(T); +if n == 1, F = feval(fun,T,0); n_terms = 1; return, end + +lambda = sum(diag(T))/n; +F = eye(n)*feval(fun,lambda,0); +f_deriv_max = zeros(itmax+n-1,1); +N = T - lambda*eye(n); +mu = norm( (eye(n)-abs(triu(T,1)))\ones(n,1),inf ); + +P = N; +max_d = 1; + +for k = 1:itmax + + f = feval(fun,lambda,k); + F_old = F; + F = F + P*f; + rel_diff = norm(F - F_old,inf)/(tol+norm(F_old,inf)); + if prnt + fprintf('%3.0f: coef = %5.0e', k, abs(f)/factorial(k)); + fprintf(' N^k/k! = %7.1e', norm(P,inf)); + fprintf(' rel_d = %5.0e',rel_diff); + fprintf(' abs_d = %5.0e',norm(F - F_old,inf)); + end + P = P*N/(k+1); + + if rel_diff <= tol + + % Approximate the maximum of derivatives in convex set containing + % eigenvalues by maximum of derivatives at eigenvalues. + for j = max_d:k+n-1 + f_deriv_max(j) = norm(feval(fun,diag(T),j),inf); + end + max_d = k+n; + omega = 0; + for j = 0:n-1 + omega = max(omega,f_deriv_max(k+j)/factorial(j)); + end + + trunc = norm(P,inf)*mu*omega; % norm(F) moved to RHS to avoid / 0. + if prnt + fprintf(' [trunc,test] = [%5.0e %5.0e]', ... + trunc, tol*norm(F,inf)) + end + if prnt == 5, trunc = 0; end % Force simple stopping test. + if trunc <= tol*norm(F,inf) + n_terms = k; + if prnt, fprintf('\n'), end, return + end + end + + if prnt, fprintf('\n'), end + +end +n_terms = -1;